diff --git a/src/RootDepth.jl b/src/RootDepth.jl index 1b7c715..31d29fe 100644 --- a/src/RootDepth.jl +++ b/src/RootDepth.jl @@ -163,6 +163,7 @@ function rootdepth_main( qlatflux = zeros(Float64, nzg + 2) wtd_old = wtd[i, j] + dθ = zeros(Float64, nzg) # pre-initialize so line 260 can access it outside the itime loop # 时间子循环 for itime in 1:round(Int, steps) @@ -182,13 +183,13 @@ function rootdepth_main( et_c_daily[i, j] += pet_c - watdef * 1.0e3 # 土壤水流计算 - # updated_o18, transpo18step + # updated_o18, transpo18step (isotope tracking removed) et_s_step, runoff, rechstep, flux_step, qrfcorrect, updated_smoi = soilfluxes(nzg, Δt / steps, z₋ₕ, dz, soiltxt[1, i, j], θ_wtd[i, j], dθ, dsmoideep, θ[:, i, j], wtd[i, j], ppdrip_step, pet_s, fdepth[i, j], qlatstep, qrfstep, floodstep, icefac - ; freedrain + ; freedrain=Bool(freedrain) ) # θ_eq[:, i, j], o18[:, i, j], o18ratiopp[i, j], tempsfc[i, j], # qlato18step, transpo18step @@ -198,14 +199,14 @@ function rootdepth_main( qsrun[i, j] += max(runoff - floodstep, 0.0) # m rech[i, j] += rechstep * 1.0e3 et_s[i, j] += et_s_step - transpo18[i, j] += transpo18step * 1.0e3 + # transpo18[i, j] += transpo18step * 1.0e3 # isotope tracking removed et_s_daily[i, j] += et_s_step ppacum[i, j] += ppdrip_step qrf[i, j] += qrfcorrect - # 更新土壤含水量和同位素 + # 更新土壤含水量 θ[:, i, j] .= updated_smoi - o18[:, i, j] .= updated_o18 + # o18[:, i, j] .= updated_o18 # isotope tracking removed # 更新浅层地下水位 wtd[i, j], rech_additional = updatewtd_shallow(nzg, freedrain, z₋ₕ, dz, diff --git a/src/SoilFluxes.jl b/src/SoilFluxes.jl index 2e9e239..33c77f3 100644 --- a/src/SoilFluxes.jl +++ b/src/SoilFluxes.jl @@ -202,6 +202,19 @@ function soilfluxes( cc[1] = D_mid[2] / Δz₊ₕ[2] bb[1] = -(cc[1] + Δz[1] / dt) rr[1] = -θ[1] * Δz[1] / dt - K_mid[2] + K_mid[1] + transp[1] / dt + + # k=2 must be set explicitly since the main loop starts at max(jwt,3)=3 + k = 2 + aa[k] = D_mid[k] / Δz₊ₕ[k] + if k < nzg + cc[k] = D_mid[k+1] / Δz₊ₕ[k+1] + bb[k] = -(aa[k] + cc[k] + Δz[k] / dt) + rr[k] = -θ[k] * Δz[k] / dt - K_mid[k+1] + K_mid[k] + transp[k] / dt + else + cc[k] = 0.0 + bb[k] = -(aa[k] + Δz[k] / dt) + rr[k] = -θ[k] * Δz[k] / dt + K_mid[k] + transp[k] / dt + end end # 求解三对角系统 @@ -283,7 +296,7 @@ function soilfluxes( # 处理最上层的干燥限制 k = nzg soil = get_soil_params(soiltxt) - θ_cp = soil.ρb * max(min(exp((z[k] + 1.5) / fdepth), 1.0), 0.1) + θ_cp = soil.θ_cp * max(min(exp((z[k] + 1.5) / fdepth), 1.0), 0.1) et_s = pet_s if θ_old[k] < θ_cp @@ -307,7 +320,7 @@ function soilfluxes( # 继续处理下层 for k in (nzg-1):-1:1 soil = get_soil_params(soiltxt) - θ_cp = soil.ρb * max(min(exp((z[k] + 1.5) / fdepth), 1.0), 0.1) + θ_cp = soil.θ_cp * max(min(exp((z[k] + 1.5) / fdepth), 1.0), 0.1) if θ_old[k] < θ_cp dθ = max((θ_cp - θ_old[k]) * Δz[k], 0.0) diff --git a/src/extraction.jl b/src/extraction.jl index 545838e..2e45888 100644 --- a/src/extraction.jl +++ b/src/extraction.jl @@ -111,7 +111,7 @@ function extraction( if toteasy == 0.0 rootactivity = zeros(Float64, nzg) else - rootactivity = clamp((easy .* dz2) ./ toteasy, 0.0, 1.0) + rootactivity = clamp.((easy .* dz2) ./ toteasy, 0.0, 1.0) end # 更新非活跃天数 @@ -141,7 +141,7 @@ function extraction( end # 计算土壤水分胁迫因子 - fswp = clamp(θ_root / rootfc, 0.0, 1.0) + fswp = rootfc == 0.0 ? 0.0 : clamp(θ_root / rootfc, 0.0, 1.0) # 计算冠层和土壤阻力 rs_c = fswp == 0.0 ? 5000.0 : min(rs_c_factor / fswp, 5000.0) diff --git a/test/runtests.jl b/test/runtests.jl index 208e934..bdd0ea0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,7 +5,10 @@ using ASAP, Test include("test_soil_initialization.jl") include("test_evapotranspiration.jl") include("test_interception.jl") -# include("test_soil_fluxes.jl") +include("test_soil_parameters.jl") +include("test_soil_fluxes.jl") +include("test_extraction.jl") +include("test_updatewtd_shallow.jl") +include("test_rootdepth.jl") include("wtable/runtests.jl") -# include("test_soil_parameters.jl") diff --git a/test/test_extraction.jl b/test/test_extraction.jl new file mode 100644 index 0000000..a8030f7 --- /dev/null +++ b/test/test_extraction.jl @@ -0,0 +1,123 @@ +using ASAP, Test + +# 典型参数设置（Penman-Monteith 方程参数） +const Δ_test = 145.0 # 饱和水汽压斜率 (Pa/K) +const γ_test = 66.0 # 湿度计常数 (Pa/K) +const λ_test = 2.45e6 # 汽化潜热 (J/kg) + +function make_extraction_inputs(nzg; lai=2.0, θ_val=0.25, icefac_val=Int8(0)) + z₋ₕ, dz = initializesoildepth(nzg) + Δt = 3600.0 + soiltxt = 1 + wtd = z₋ₕ[1] - 0.5 # 水位在所有层之下（不饱和） + θ = fill(θ_val, nzg) + θ_wtd = 0.3 + ra_a = 100.0 + ra_c = 50.0 + rs_c_factor = 200.0 + R_a = (Δ_test + γ_test) * ra_a * 3.0 # 近似辐射驱动力 + R_s = (Δ_test + γ_test) * ra_a * 1.0 + petfactor_s = R_s + petfactor_c = R_a + inactivedays = zeros(Int, nzg + 2) # 全部活跃 + maxinactivedays = 10 + hhveg = 10.0 + fdepth = 2.0 + icefac = fill(icefac_val, nzg) + return (nzg, z₋ₕ, dz, Δt, soiltxt, wtd, θ, θ_wtd, + Δ_test, γ_test, λ_test, + lai, ra_a, ra_c, rs_c_factor, + R_a, R_s, petfactor_s, petfactor_c, + inactivedays, maxinactivedays, + hhveg, fdepth, icefac) +end + + +@testset "extraction基本功能" begin + args = make_extraction_inputs(10) + pet_s, pet_c, watdef, dθ, dθ_deep = extraction(args...) + + @test pet_s >= 0.0 + @test pet_c >= 0.0 + @test watdef >= 0.0 + @test length(dθ) == 10 + @test dθ_deep == 0.0 # 深层无提取 + @test all(dθ .>= 0.0) # 每层提取量非负 +end + +@testset "零LAI时无冠层蒸腾" begin + args = make_extraction_inputs(10; lai=0.0) + pet_s, pet_c, watdef, dθ, dθ_deep = extraction(args...) + + # LAI=0 时，冠层蒸腾强制为零 + @test pet_c == 0.0 + @test pet_s >= 0.0 +end + +@testset "全部冰冻时无提取" begin + args = make_extraction_inputs(10; icefac_val=Int8(1)) + pet_s, pet_c, watdef, dθ, dθ_deep = extraction(args...) + + # 冰冻情况下根区水分提取便利性为零 + @test all(dθ .== 0.0) + @test watdef >= 0.0 # 水分亏缺 ≥ 0 +end + +@testset "干旱土壤产生水分亏缺" begin + # 非常干的土壤 (θ << θ_wilt) + nzg = 10 + args_dry = make_extraction_inputs(nzg; θ_val=0.05) + pet_s, pet_c, watdef_dry, dθ_dry, _ = extraction(args_dry...) + + args_wet = make_extraction_inputs(nzg; θ_val=0.35) + _, _, watdef_wet, dθ_wet, _ = extraction(args_wet...) + + # 湿润土壤水分亏缺应该比干燥时少 + @test watdef_wet <= watdef_dry + # 湿润土壤提取总量应该更多 + @test sum(dθ_wet) >= sum(dθ_dry) +end + +@testset "inactivedays正确更新" begin + nzg = 10 + args = make_extraction_inputs(nzg; θ_val=0.25) + # inactivedays初始为全0（全部活跃） + pet_s, pet_c, watdef, dθ, dθ_deep = extraction(args...) + + # 函数会修改 inactivedays（若根区活跃则清零，否则+1） + @test pet_s >= 0.0 + @test pet_c >= 0.0 +end + +@testset "蒸发量与时间步长成正比" begin + nzg = 10 + + # 短时间步长（30 min） + args1 = make_extraction_inputs(nzg) + # 修改 Δt 为 1800 s + args1 = (args1[1], args1[2], args1[3], 1800.0, args1[5:end]...) + pet_s1, pet_c1, _, _, _ = extraction(args1...) + + # 长时间步长（2 h） + args2 = make_extraction_inputs(nzg) + args2 = (args2[1], args2[2], args2[3], 7200.0, args2[5:end]...) + pet_s2, pet_c2, _, _, _ = extraction(args2...) + + # 时间步长加倍，PET近似翻倍 + @test pet_s2 > pet_s1 + @test pet_c2 > pet_c1 +end + +@testset "各层提取量之和等于冠层蒸腾量（无亏缺时）" begin + nzg = 10 + # 使用高含水量，确保无水分亏缺 + args = make_extraction_inputs(nzg; θ_val=0.38) + pet_s, pet_c, watdef, dθ, dθ_deep = extraction(args...) + + if watdef == 0.0 + # 无水分亏缺：提取量之和应等于 pet_c（m） + @test sum(dθ) ≈ pet_c * 1e-3 atol=1e-10 + else + @test sum(dθ) + watdef ≈ pet_c * 1e-3 atol=1e-10 + end +end diff --git a/test/test_rootdepth.jl b/test/test_rootdepth.jl new file mode 100644 index 0000000..861a75f --- /dev/null +++ b/test/test_rootdepth.jl @@ -0,0 +1,141 @@ +using ASAP, Test + +# rootdepth_main使用nzg=40（标准层数，与icefactor[26:40]索引匹配） +const NZG_RD = 40 + +# 构建rootdepth_main所需的最小参数集（网格大小3×3，内层循环只访问(2,2)） +function make_rootdepth_inputs(; freedrain=1, veg_val=1.0, landmask_val=1, nzg=NZG_RD) + is, ie, js, je = 1, 3, 1, 3 + nx, ny = ie, je # 3×3 grid + z₋ₕ, dz = initializesoildepth(nzg) + Δt = 3600.0 + + landmask = zeros(Int, nx, ny); landmask[2, 2] = landmask_val + veg = fill(veg_val, nx, ny) + hveg = fill(1.0, nx, ny) + soiltxt = ones(Int, 2, nx, ny) + wind = fill(3.0, nx, ny) + temp = fill(293.0, nx, ny) + qair = fill(0.01, nx, ny) + press = fill(101325.0, nx, ny) + netrad = fill(100.0, nx, ny) + rshort = fill(150.0, nx, ny) + lai = fill(2.0, nx, ny) + precip = fill(5.0, nx, ny) + qsrun = zeros(Float64, nx, ny) + θ = fill(0.25, nzg, nx, ny) + θ_eq = fill(0.20, nzg, nx, ny) + θ_wtd = fill(0.35, nx, ny) + wtd = fill(-1.0, nx, ny) + waterdeficit = zeros(Float64, nx, ny) + watext = zeros(Float64, nzg, nx, ny) + watextdeep = zeros(Float64, nx, ny) + rech = zeros(Float64, nx, ny) + deeprech = zeros(Float64, nx, ny) + et_s = zeros(Float64, nx, ny) + et_i = zeros(Float64, nx, ny) + et_c = zeros(Float64, nx, ny) + intercepstore = zeros(Float64, nx, ny) + ppacum = zeros(Float64, nx, ny) + pInfiltDepth = zeros(Float64, nx, ny) + pInfiltDepthK_old = zeros(Int8, nx, ny) + qlat = zeros(Float64, nx, ny) + qlatsum = zeros(Float64, nx, ny) + qsprings = zeros(Float64, nx, ny) + inactivedays = zeros(Int, nzg+2, nx, ny) + maxinactivedays = 10 + fdepth = fill(2.0, nx, ny) + steps = 1.0 + floodheight = zeros(Float64, nx, ny) + qrf = zeros(Float64, nx, ny) + delsfcwat = zeros(Float64, nx, ny) + icefactor = zeros(Int8, nx, ny, nzg) + wtdflux = zeros(Float64, nx, ny) + et_s_daily = zeros(Float64, nx, ny) + et_c_daily = zeros(Float64, nx, ny) + transptop = zeros(Float64, nx, ny) + infilk = zeros(Int8, nx, ny) + infilflux = zeros(Float64, nzg, nx, ny) + infilfluxday = zeros(Float64, nzg, nx, ny) + infilcounter = zeros(Int16, nzg, nx, ny) + hour = 0 + o18 = zeros(Float64, nzg, nx, ny) + o18ratiopp = zeros(Float64, nx, ny) + tempsfc = fill(293.0, nx, ny) + qlato18 = zeros(Float64, nx, ny) + transpo18 = zeros(Float64, nx, ny) + upflux = zeros(Float64, nzg, nx, ny) + + return (freedrain, is, ie, js, je, nzg, + z₋ₕ, dz, Δt, + landmask, veg, hveg, + soiltxt, wind, temp, + qair, press, netrad, + rshort, lai, precip, + qsrun, θ, θ_eq, + θ_wtd, wtd, waterdeficit, + watext, + watextdeep, rech, + deeprech, + et_s, et_i, et_c, + intercepstore, ppacum, + pInfiltDepth, pInfiltDepthK_old, qlat, + qlatsum, qsprings, inactivedays, + maxinactivedays, fdepth, + steps, floodheight, qrf, + delsfcwat, icefactor, wtdflux, + et_s_daily, et_c_daily, transptop, + infilk, infilflux, infilfluxday, + infilcounter, hour, + o18, o18ratiopp, tempsfc, + qlato18, transpo18, upflux) +end + +@testset "rootdepth_main: 全零陆地掩码（无活跃格点）" begin + # landmask全零：内层循环立即跳过，函数应返回nothing且不修改任何数组 + args = make_rootdepth_inputs(landmask_val=0) + et_s_before = copy(args[33]) # et_s is the 33rd arg + result = rootdepth_main(args...) + @test result === nothing + # 所有输出应保持不变 + @test args[33] == et_s_before +end + +@testset "rootdepth_main: 水体格点（veg=1）" begin + # veg=1 (水体)：执行潜在蒸散发计算，但跳过植被相关代码 + args = make_rootdepth_inputs(veg_val=1.0, landmask_val=1) + et_s = args[33] + et_s_init = et_s[2, 2] + + result = rootdepth_main(args...) + @test result === nothing + # 水体格点应更新土壤蒸发 + @test et_s[2, 2] >= et_s_init +end + +@testset "rootdepth_main: 自由排水+植被格点（veg=7，短草）" begin + # freedrain=1（自由排水），veg=7（短草）：执行全部子模块（extraction、soilfluxes、updatewtd_shallow） + args = make_rootdepth_inputs(freedrain=1, veg_val=7.0, landmask_val=1) + et_s = args[33] + et_c = args[35] + rech = args[30] + + result = rootdepth_main(args...) + @test result === nothing + # 植被格点应更新蒸腾和补给 + @test isfinite(et_s[2, 2]) + @test isfinite(et_c[2, 2]) + @test isfinite(rech[2, 2]) +end + +@testset "rootdepth_main: 非自由排水+植被格点（veg=7）" begin + # freedrain=0（非自由排水）：测试soilfluxes的!freedrain分支在全流程中的集成 + args = make_rootdepth_inputs(freedrain=0, veg_val=7.0, landmask_val=1) + et_s = args[33] + et_c = args[35] + + result = rootdepth_main(args...) + @test result === nothing + @test isfinite(et_s[2, 2]) + @test isfinite(et_c[2, 2]) +end diff --git a/test/test_soil_fluxes.jl b/test/test_soil_fluxes.jl index 7f76bfc..b9b12d0 100644 --- a/test/test_soil_fluxes.jl +++ b/test/test_soil_fluxes.jl @@ -2,350 +2,270 @@ using ASAP, Test # 测试基本功能 @testset "soilfluxes基本功能" begin - # 设置测试参数 nzg = 10 - i, j = 1, 1 - freedrain = 0 - dtll = 3600.0 # 1小时时间步长 - soiltxt = 1 # 土壤类型1 + dt = 3600.0 # 1小时时间步长 + soiltxt = 1 # 沙土 - # 初始化土壤层 - slz, dz = initializesoildepth(nzg) + z₋ₕ, Δz = initializesoildepth(nzg) - # 初始化输入数据 - smoiwtd = 0.3 + θ_wtd = 0.3 transp = zeros(Float64, nzg) transpdeep = 0.0 - θ = fill(0.25, nzg) # 使用希腊字母θ + θ = fill(0.25, nzg) wtd = -2.0 - precip = 5.0 # mm - pet_s = 2.0 # mm + precip = 5.0 + pet_s = 2.0 fdepth = 2.0 qlat = 0.0 qrf = 0.0 flood = 0.0 icefactor = zeros(Int8, nzg) - θ_eq = fill(0.2, nzg) # 使用希腊字母θ - o18 = θ .* (-5.0) # 初始同位素值 - precipo18 = -8.0 - tempsfc = 293.15 - qlato18 = 0.0 - transpo18 = 0.0 - - # 调用soilfluxes函数 - et_s, runoff, rech, flux, qrfcorrect, transpo18_out, θ_new, o18_new = - soilfluxes(nzg, freedrain, dtll, slz, dz, soiltxt, - smoiwtd, transp, transpdeep, θ, wtd, precip, pet_s, - fdepth, qlat, qrf, flood, icefactor, θ_eq, - o18, precipo18, tempsfc, qlato18, transpo18) - - # 基本检查 - @test et_s >= 0.0 # 蒸发应该非负 - @test runoff >= 0.0 # 径流应该非负 - @test length(θ_new) == nzg # 输出土壤含水量长度正确 - @test length(o18_new) == nzg # 输出同位素长度正确 - @test length(flux) == nzg + 1 # 通量长度正确 - - # 检查土壤含水量在合理范围内 + + et_s, runoff, rech, flux, qrfcorrect, θ_new = + soilfluxes(nzg, dt, z₋ₕ, Δz, soiltxt, + θ_wtd, transp, transpdeep, copy(θ), wtd, + precip, pet_s, fdepth, qlat, qrf, flood, icefactor) + + @test et_s >= 0.0 + @test runoff >= 0.0 + @test length(θ_new) == nzg + @test length(flux) == nzg + 1 soil = get_soil_params(soiltxt) for k in 1:nzg - @test θ_new[k] >= 0.0 # 含水量非负 - @test θ_new[k] <= soil.θ_sat * 1.1 # 不超过饱和含水量太多 + @test θ_new[k] >= 0.0 + @test θ_new[k] <= soil.θ_sat * 1.01 end - - # 检查同位素值非负 - @test all(o18_new .>= 0.0) end # 测试不同降水条件 @testset "不同降水条件测试" begin nzg = 5 - slz, dz = initializesoildepth(nzg) - - # 基础参数 - i, j = 1, 1 - freedrain = 0 - dtll = 3600.0 + z₋ₕ, Δz = initializesoildepth(nzg) + dt = 3600.0 soiltxt = 1 - smoiwtd = 0.3 + θ_wtd = 0.3 transp = zeros(Float64, nzg) transpdeep = 0.0 θ = fill(0.25, nzg) wtd = -2.0 pet_s = 2.0 fdepth = 2.0 - qlat = 0.0 - qrf = 0.0 - flood = 0.0 icefactor = zeros(Int8, nzg) - θ_eq = fill(0.2, nzg) - o18 = θ .* (-5.0) - precipo18 = -8.0 - tempsfc = 293.15 - qlato18 = 0.0 - transpo18 = 0.0 - - # 测试无降水 - et_s1, runoff1, _, _, _, _, _, _ = - soilfluxes(nzg, freedrain, dtll, slz, dz, soiltxt, - smoiwtd, transp, transpdeep, copy(θ), wtd, 0.0, pet_s, - fdepth, qlat, qrf, flood, icefactor, θ_eq, - copy(o18), precipo18, tempsfc, qlato18, transpo18) - - # 测试大降水 - et_s2, runoff2, _, _, _, _, _, _ = - soilfluxes(nzg, freedrain, dtll, slz, dz, soiltxt, - smoiwtd, transp, transpdeep, copy(θ), wtd, 50.0, pet_s, - fdepth, qlat, qrf, flood, icefactor, θ_eq, - copy(o18), precipo18, tempsfc, qlato18, transpo18) + + # 无降水 + _, runoff1, _, _, _, _ = + soilfluxes(nzg, dt, z₋ₕ, Δz, soiltxt, + θ_wtd, transp, transpdeep, copy(θ), wtd, + 0.0, pet_s, fdepth, 0.0, 0.0, 0.0, icefactor) + + # 大降水 + _, runoff2, _, _, _, _ = + soilfluxes(nzg, dt, z₋ₕ, Δz, soiltxt, + θ_wtd, transp, transpdeep, copy(θ), wtd, + 50.0, pet_s, fdepth, 0.0, 0.0, 0.0, icefactor) # 大降水应该产生更多径流 @test runoff2 >= runoff1 end # 测试自由排水条件 -@testset "自由排水条件测试" begin +@testset "自由排水条件" begin nzg = 5 - slz, dz = initializesoildepth(nzg) - - # 基础参数 - i, j = 1, 1 - freedrain = 1 # 自由排水 - dtll = 3600.0 + z₋ₕ, Δz = initializesoildepth(nzg) + dt = 3600.0 soiltxt = 1 - smoiwtd = 0.3 + θ_wtd = 0.3 transp = zeros(Float64, nzg) - transpdeep = 0.0 θ = fill(0.25, nzg) wtd = -2.0 - precip = 10.0 - pet_s = 2.0 - fdepth = 2.0 - qlat = 0.0 - qrf = 0.0 - flood = 0.0 icefactor = zeros(Int8, nzg) - θ_eq = fill(0.2, nzg) - o18 = θ .* (-5.0) - precipo18 = -8.0 - tempsfc = 293.15 - qlato18 = 0.0 - transpo18 = 0.0 - - et_s, runoff, rech, flux, qrfcorrect, transpo18_out, θ_new, o18_new = - soilfluxes(nzg, freedrain, dtll, slz, dz, soiltxt, - smoiwtd, transp, transpdeep, θ, wtd, precip, pet_s, - fdepth, qlat, qrf, flood, icefactor, θ_eq, - o18, precipo18, tempsfc, qlato18, transpo18) - - # 自由排水条件下应该有补给 - @test abs(rech) >= 0.0 # 补给可能为正或负 -end -# 测试有蒸腾的情况 -@testset "有蒸腾测试" begin - nzg = 5 - slz, dz = initializesoildepth(nzg) + et_s, runoff, rech, flux, qrfcorrect, θ_new = + soilfluxes(nzg, dt, z₋ₕ, Δz, soiltxt, + θ_wtd, transp, 0.0, copy(θ), wtd, + 10.0, 2.0, 2.0, 0.0, 0.0, 0.0, icefactor; + freedrain=true) - # 基础参数 - i, j = 1, 1 - freedrain = 0 - dtll = 3600.0 - soiltxt = 1 - smoiwtd = 0.3 - transp = fill(0.001, nzg) # 每层有蒸腾 - transpdeep = 0.0 - θ = fill(0.3, nzg) # 较高含水量 - wtd = -2.0 - precip = 0.0 # 无降水 - pet_s = 2.0 - fdepth = 2.0 - qlat = 0.0 - qrf = 0.0 - flood = 0.0 - icefactor = zeros(Int8, nzg) - θ_eq = fill(0.2, nzg) - o18 = θ .* (-5.0) - precipo18 = -8.0 - tempsfc = 293.15 - qlato18 = 0.0 - transpo18 = 0.0 - - θ_original = copy(θ) - - et_s, runoff, rech, flux, qrfcorrect, transpo18_out, θ_new, o18_new = - soilfluxes(nzg, freedrain, dtll, slz, dz, soiltxt, - smoiwtd, transp, transpdeep, θ, wtd, precip, pet_s, - fdepth, qlat, qrf, flood, icefactor, θ_eq, - o18, precipo18, tempsfc, qlato18, transpo18) - - # 有蒸腾应该降低土壤含水量 - @test sum(θ_new) <= sum(θ_original) - @test transpo18_out >= 0.0 # 蒸腾同位素输出应该非负 + @test et_s >= 0.0 + @test runoff >= 0.0 + @test length(θ_new) == nzg end # 测试冰冻因子影响 @testset "冰冻因子测试" begin nzg = 5 - slz, dz = initializesoildepth(nzg) - - # 基础参数 - i, j = 1, 1 - freedrain = 0 - dtll = 3600.0 + z₋ₕ, Δz = initializesoildepth(nzg) + dt = 3600.0 soiltxt = 1 - smoiwtd = 0.3 + θ_wtd = 0.3 transp = zeros(Float64, nzg) - transpdeep = 0.0 θ = fill(0.25, nzg) wtd = -2.0 precip = 10.0 - pet_s = 2.0 - fdepth = 2.0 - qlat = 0.0 - qrf = 0.0 - flood = 0.0 - θ_eq = fill(0.2, nzg) - o18 = θ .* (-5.0) - precipo18 = -8.0 - tempsfc = 293.15 - qlato18 = 0.0 - transpo18 = 0.0 # 无冰冻 - icefactor1 = zeros(Int8, nzg) - et_s1, runoff1, _, _, _, _, _, _ = - soilfluxes(nzg, freedrain, dtll, slz, dz, soiltxt, - smoiwtd, transp, transpdeep, copy(θ), wtd, precip, pet_s, - fdepth, qlat, qrf, flood, icefactor1, θ_eq, - copy(o18), precipo18, tempsfc, qlato18, transpo18) - - # 部分冰冻 - icefactor2 = ones(Int8, nzg) # 全部冰冻 - et_s2, runoff2, _, _, _, _, _, _ = - soilfluxes(nzg, freedrain, dtll, slz, dz, soiltxt, - smoiwtd, transp, transpdeep, copy(θ), wtd, precip, pet_s, - fdepth, qlat, qrf, flood, icefactor2, θ_eq, - copy(o18), precipo18, tempsfc, qlato18, transpo18) - - # 冰冻会影响水分运动，通常导致更多径流 + icefactor_none = zeros(Int8, nzg) + _, runoff1, _, _, _, _ = + soilfluxes(nzg, dt, z₋ₕ, Δz, soiltxt, + θ_wtd, transp, 0.0, copy(θ), wtd, + precip, 2.0, 2.0, 0.0, 0.0, 0.0, icefactor_none) + + # 全部冰冻 + icefactor_all = ones(Int8, nzg) + _, runoff2, _, _, _, _ = + soilfluxes(nzg, dt, z₋ₕ, Δz, soiltxt, + θ_wtd, transp, 0.0, copy(θ), wtd, + precip, 2.0, 2.0, 0.0, 0.0, 0.0, icefactor_all) + + # 冰冻会阻碍水分运动，导致更多径流 @test runoff2 >= runoff1 end -# 测试不同土壤类型 -@testset "不同土壤类型测试" begin +# 测试饱和土壤产生径流 +@testset "饱和土壤径流" begin nzg = 5 - slz, dz = initializesoildepth(nzg) + z₋ₕ, Δz = initializesoildepth(nzg) + # 使用粘土 (type 11, Ksat=0.000001283 m/s ≈ 4.6mm/hr) — 50mm降水大幅超过入渗能力 + soiltxt = 11 + soil = get_soil_params(soiltxt) + icefactor = zeros(Int8, nzg) - # 基础参数 - i, j = 1, 1 - freedrain = 0 - dtll = 3600.0 - smoiwtd = 0.3 + θ_sat = fill(soil.θ_sat, nzg) + + _, runoff, _, _, _, _ = + soilfluxes(nzg, 3600.0, z₋ₕ, Δz, soiltxt, + soil.θ_sat, zeros(nzg), 0.0, θ_sat, -2.0, + 50.0, 2.0, 2.0, 0.0, 0.0, 0.0, icefactor) + + @test runoff > 0.0 +end + +# 测试非自由排水条件：水位在土柱以下（jwt <= 2） +@testset "非自由排水-水位在土柱以下" begin + nzg = 10 + z₋ₕ, Δz = initializesoildepth(nzg) + dt = 3600.0 + soiltxt = 4 # 壤土 + θ_wtd = 0.4 transp = zeros(Float64, nzg) - transpdeep = 0.0 θ = fill(0.25, nzg) - wtd = -2.0 - precip = 10.0 - pet_s = 2.0 - fdepth = 2.0 - qlat = 0.0 - qrf = 0.0 - flood = 0.0 + # wtd低于最底层边界，jwt=1 → jwt<=2分支 + wtd = z₋ₕ[1] - 0.5 icefactor = zeros(Int8, nzg) - θ_eq = fill(0.2, nzg) - o18 = θ .* (-5.0) - precipo18 = -8.0 - tempsfc = 293.15 - qlato18 = 0.0 - transpo18 = 0.0 - - results = [] - - # 测试多种土壤类型 - for soiltxt in [1, 5, 10] - et_s, runoff, rech, flux, qrfcorrect, transpo18_out, θ_new, o18_new = - soilfluxes(nzg, freedrain, dtll, slz, dz, soiltxt, - smoiwtd, transp, transpdeep, copy(θ), wtd, precip, pet_s, - fdepth, qlat, qrf, flood, icefactor, θ_eq, - copy(o18), precipo18, tempsfc, qlato18, transpo18) - - push!(results, (soiltxt=soiltxt, et_s=et_s, runoff=runoff)) - - # 基本合理性检查 - @test et_s >= 0.0 - @test runoff >= 0.0 - end - # 不同土壤类型应该产生不同的结果 - @test length(unique([r.runoff for r in results])) > 1 + et_s, runoff, rech, flux, qrfcorrect, θ_new = + soilfluxes(nzg, dt, z₋ₕ, Δz, soiltxt, + θ_wtd, transp, 0.0, copy(θ), wtd, + 5.0, 1.0, 2.0, 0.0, 0.0, 0.0, icefactor; + freedrain=false) + + @test et_s >= 0.0 + @test runoff >= 0.0 + @test length(θ_new) == nzg + @test length(flux) == nzg + 1 + @test all(isfinite.(θ_new)) + # 非自由排水且水位极低，底层无渗漏 + @test flux[1] == 0.0 end -# 测试三对角矩阵求解器 -@testset "tridag! 函数测试" begin - n = 5 - a = [0.0, 1.0, 1.0, 1.0, 1.0] # 下对角线 - b = [2.0, 2.0, 2.0, 2.0, 2.0] # 主对角线 - c = [1.0, 1.0, 1.0, 1.0, 0.0] # 上对角线 - r = [1.0, 2.0, 3.0, 4.0, 5.0] # 右端向量 - u = zeros(Float64, n) # 解向量 - - # 调用三对角求解器 - tridag!(a, b, c, r, u) - - # 检查解的合理性 - @test length(u) == n - @test all(isfinite.(u)) # 解应该是有限的 -end +# 测试非自由排水条件：水位在中间层（jwt > 3） +@testset "非自由排水-水位在中间层" begin + nzg = 10 + z₋ₕ, Δz = initializesoildepth(nzg) + dt = 3600.0 + soiltxt = 4 # 壤土 + soil = get_soil_params(soiltxt) + θ_wtd = soil.θ_sat + transp = zeros(Float64, nzg) + # 土柱接近饱和，水位在层5和层6之间（wtd=-0.55 → jwt=6） + θ = fill(soil.θ_sat * 0.9, nzg) + wtd = 0.5 * (z₋ₕ[5] + z₋ₕ[6]) # between layers 5 and 6 + icefactor = zeros(Int8, nzg) -# 边界条件测试 -@testset "边界条件测试" begin - nzg = 3 # 使用较少的层数简化测试 - slz, dz = initializesoildepth(nzg) + et_s, runoff, rech, flux, qrfcorrect, θ_new = + soilfluxes(nzg, dt, z₋ₕ, Δz, soiltxt, + θ_wtd, transp, 0.0, copy(θ), wtd, + 3.0, 1.0, 2.0, 0.0, 0.0, 0.0, icefactor; + freedrain=false) + + @test et_s >= 0.0 + @test runoff >= 0.0 + @test length(θ_new) == nzg + @test all(isfinite.(θ_new)) + # 非自由排水，底层通量为零 + @test flux[1] == 0.0 +end - # 基础参数 - i, j = 1, 1 - freedrain = 0 - dtll = 3600.0 - soiltxt = 1 - smoiwtd = 0.3 +# 对比自由排水与非自由排水的底层渗漏差异 +@testset "自由排水与非自由排水底层渗漏对比" begin + nzg = 10 + z₋ₕ, Δz = initializesoildepth(nzg) + dt = 3600.0 + soiltxt = 1 # 沙土（Ksat大，容易产生底层渗漏） + θ_wtd = 0.35 transp = zeros(Float64, nzg) - transpdeep = 0.0 - θ = fill(0.25, nzg) + θ = fill(0.3, nzg) wtd = -2.0 - pet_s = 2.0 - fdepth = 2.0 - qlat = 0.0 - qrf = 0.0 - flood = 0.0 icefactor = zeros(Int8, nzg) - θ_eq = fill(0.2, nzg) - o18 = θ .* (-5.0) - precipo18 = -8.0 - tempsfc = 293.15 - qlato18 = 0.0 - transpo18 = 0.0 - - # 测试极端干燥条件 - θ_dry = fill(0.05, nzg) # 很干的土壤 - et_s, runoff, _, _, _, _, θ_new, _ = - soilfluxes(nzg, freedrain, dtll, slz, dz, soiltxt, - smoiwtd, transp, transpdeep, θ_dry, wtd, 0.0, pet_s, - fdepth, qlat, qrf, flood, icefactor, θ_eq, - copy(o18), precipo18, tempsfc, qlato18, transpo18) - - # 干燥条件下蒸发应该受限 - @test et_s <= pet_s - - # 测试饱和条件 - soil = get_soil_params(soiltxt) - θ_sat = fill(soil.θ_sat, nzg) - et_s2, runoff2, _, _, _, _, θ_new2, _ = - soilfluxes(nzg, freedrain, dtll, slz, dz, soiltxt, - smoiwtd, transp, transpdeep, θ_sat, wtd, 50.0, pet_s, - fdepth, qlat, qrf, flood, icefactor, θ_eq, - copy(o18), precipo18, tempsfc, qlato18, transpo18) - - # 饱和条件下应该有显著径流 - @test runoff2 > 0.0 + + # 自由排水 + _, _, rech_free, flux_free, _, _ = + soilfluxes(nzg, dt, z₋ₕ, Δz, soiltxt, + θ_wtd, transp, 0.0, copy(θ), wtd, + 5.0, 1.0, 2.0, 0.0, 0.0, 0.0, icefactor; + freedrain=true) + + # 非自由排水 + _, _, rech_nonfree, flux_nonfree, _, _ = + soilfluxes(nzg, dt, z₋ₕ, Δz, soiltxt, + θ_wtd, transp, 0.0, copy(θ), wtd, + 5.0, 1.0, 2.0, 0.0, 0.0, 0.0, icefactor; + freedrain=false) + + # 非自由排水时底层通量为零（无重力排水） + @test flux_nonfree[1] == 0.0 + @test rech_nonfree == 0.0 + # 自由排水时底层有向下渗漏（Q[1]<0 表示向下排水） + @test rech_free <= 0.0 end + +# 测试三对角矩阵求解器 +@testset "tridag! 函数测试" begin + # 测试简单3x3系统: [2 1 0; 1 2 1; 0 1 2] * x = [5; 8; 5] → x = [1; 2; 1]+adjust + # 使用一个有解析解的简单系统 + n = 3 + a = [0.0, -1.0, -1.0] # 下对角线 (a[1]未使用) + b = [2.0, 2.0, 2.0] # 主对角线 + c = [-1.0, -1.0, 0.0] # 上对角线 (c[n]未使用) + r = [1.0, 2.0, 3.0] # 右端向量 + u = zeros(Float64, n) + + # 方程组: 2u1 - u2 = 1; -u1 + 2u2 - u3 = 2; -u2 + 2u3 = 3 + # 解: u1 = 2.5, u2 = 4.0, u3 = 3.5 + ASAP.tridag!(a, b, c, r, u) + + @test u[1] ≈ 2.5 atol=1e-10 + @test u[2] ≈ 4.0 atol=1e-10 + @test u[3] ≈ 3.5 atol=1e-10 + + # 测试更大的系统 + n2 = 5 + a2 = [0.0, 1.0, 1.0, 1.0, 1.0] + b2 = [3.0, 3.0, 3.0, 3.0, 3.0] + c2 = [1.0, 1.0, 1.0, 1.0, 0.0] + r2 = [1.0, 2.0, 3.0, 4.0, 5.0] + u2 = zeros(Float64, n2) + + ASAP.tridag!(a2, b2, c2, r2, u2) + + @test length(u2) == n2 + @test all(isfinite.(u2)) + # 验证 A*u ≈ r + Av = [b2[1]*u2[1] + c2[1]*u2[2], + a2[2]*u2[1] + b2[2]*u2[2] + c2[2]*u2[3], + a2[3]*u2[2] + b2[3]*u2[3] + c2[3]*u2[4], + a2[4]*u2[3] + b2[4]*u2[4] + c2[4]*u2[5], + a2[5]*u2[4] + b2[5]*u2[5]] + @test Av ≈ r2 atol=1e-10 +end + diff --git a/test/test_soil_parameters.jl b/test/test_soil_parameters.jl index bd84326..c31922c 100644 --- a/test/test_soil_parameters.jl +++ b/test/test_soil_parameters.jl @@ -1,37 +1,66 @@ +using ASAP, Test + # 测试获取土壤参数 @testset "获取土壤参数" begin soil = get_soil_params(1) - @test soil.slmsts ≈ 0.395 - @test soil.soilcp ≈ 0.050 - @test soil.slbs ≈ 4.05 - @test soil.slcons ≈ 0.000176 - @test soil.slpots ≈ -0.121 - @test soil.klatfactor ≈ 2.0 + @test soil.θ_sat ≈ 0.395 + @test soil.θ_cp ≈ 0.050 + @test soil.b ≈ 4.05 + @test soil.Ksat ≈ 0.000176 + @test soil.ψsat ≈ -0.121 + @test soil.K_latfactor ≈ 2.0 + + # 测试所有13种土壤类型 + for i in 1:13 + soil_i = get_soil_params(i) + @test soil_i.θ_sat > 0.0 + @test soil_i.θ_cp > 0.0 + @test soil_i.θ_sat > soil_i.θ_cp + @test soil_i.Ksat > 0.0 + @test soil_i.b > 0.0 + @test soil_i.K_latfactor > 0.0 + @test soil_i.ψsat < 0.0 # 基质势为负值 + end # 测试边界条件 - @test_throws BoundsError get_soil_params(0) - @test_throws BoundsError get_soil_params(14) + @test_throws ErrorException get_soil_params(0) + @test_throws ErrorException get_soil_params(14) end # 测试初始化土壤参数 @testset "初始化土壤参数" begin nzg = 10 - fieldcp, slwilt = init_soil_param(nzg) - @test size(fieldcp) == (nzg, NSTYP) - @test length(slwilt) == NSTYP - @test all(slwilt .> 0) + fieldcp, θ_wilt = init_soil_param(nzg) + @test size(fieldcp) == (nzg, ASAP.NSTYP) + @test length(θ_wilt) == ASAP.NSTYP + @test all(θ_wilt .> 0) + # 凋萎点应小于饱和含水量 + for i in 1:ASAP.NSTYP + soil = get_soil_params(i) + @test θ_wilt[i] < soil.θ_sat + end end # 测试导水率计算 @testset "导水率计算" begin - smoi = 0.3 - nsoil = 1 - k = cal_K(smoi, nsoil) - @test k > 0.0 - @test k ≈ SLCONS[nsoil] * (smoi / SLMSTS[nsoil])^(2.0 * SLBS[nsoil] + 3.0) + soil = get_soil_params(1) - # 测试边界条件 - @test_throws BoundsError cal_K(0.3, 0) - @test_throws BoundsError cal_K(0.3, 14) + # 完全饱和时，K = Ksat + k_sat = cal_K(soil.θ_sat, soil.θ_sat, soil.Ksat, soil.b) + @test k_sat ≈ soil.Ksat + + # 含水量越高，导水率越大 + k1 = cal_K(0.2, soil.θ_sat, soil.Ksat, soil.b) + k2 = cal_K(0.3, soil.θ_sat, soil.Ksat, soil.b) + @test k2 > k1 + + # 导水率应为正值 + @test k1 > 0.0 + @test k2 > 0.0 + + # 验证公式: K = Ksat * (θ/θ_sat)^(2b+3) + θ = 0.25 + k_expected = soil.Ksat * (θ / soil.θ_sat)^(2.0 * soil.b + 3.0) + @test cal_K(θ, soil.θ_sat, soil.Ksat, soil.b) ≈ k_expected end diff --git a/test/test_updatewtd_shallow.jl b/test/test_updatewtd_shallow.jl new file mode 100644 index 0000000..12956f2 --- /dev/null +++ b/test/test_updatewtd_shallow.jl @@ -0,0 +1,114 @@ +using ASAP, Test + +@testset "find_jwt基本功能" begin + # 使用 nzg=10 的标准层深度 + z₋ₕ, _ = initializesoildepth(10) + # z₋ₕ 从最深到表面 (负值到0) + + # 水位深于所有层（默认返回1） + jwt_deep = find_jwt(z₋ₕ[1] - 1.0, z₋ₕ) + @test jwt_deep == 1 + + # 水位在顶层（接近地表） + jwt_top = find_jwt(z₋ₕ[10] + 0.01, z₋ₕ) + @test jwt_top >= 1 + + # 水位正好在第一层中心 + z_center1 = 0.5 * (z₋ₕ[1] + z₋ₕ[2]) # 第1层中心 + jwt1 = find_jwt(z_center1, z₋ₕ) + @test jwt1 == 2 # 水位在第1层，上边界是 z₋ₕ[2] + + # 水位在某个中间层 + z_center5 = 0.5 * (z₋ₕ[5] + z₋ₕ[6]) # 第5层中心 + jwt5 = find_jwt(z_center5, z₋ₕ) + # 应该返回6（z₋ₕ[6] 是第5层的上边界） + @test jwt5 == 6 +end + +@testset "find_jwt边界条件" begin + z₋ₕ = [-3.0, -2.0, -1.0, -0.5, -0.1, 0.0] # nzg=5 手工构造 + + # wtd=-1.5: z₋ₕ[3]=-1.0 > -1.5? YES → jwt=3 + @test find_jwt(-1.5, z₋ₕ) == 3 + + # wtd=-0.8: z₋ₕ[4]=-0.5 > -0.8? YES → jwt=4 + @test find_jwt(-0.8, z₋ₕ) == 4 + + # wtd=-3.5 (低于所有层边界): 无 z₋ₕ[k] > -3.5 → jwt=1 (default) + @test find_jwt(-3.5, z₋ₕ) == 1 + + # wtd=-0.05 (在最浅层内部): z₋ₕ[5]=-0.1 > -0.05? No → jwt=1 (default) + # 因为 nzg=5 的 z₋ₕ 中没有比 -0.05 更大的元素（除了z₋ₕ[6]=0，但不检查） + @test find_jwt(-0.05, z₋ₕ) == 1 +end + +@testset "updatewtd_shallow水位在有效层内" begin + nzg = 20 + z₋ₕ, dz = initializesoildepth(nzg) + soiltxt = 1 + soil = get_soil_params(soiltxt) + fdepth = 10.0 + + # 水位在某一中间层内 + # 选择一个合适的 wtd + k_target = 10 # 第10层 + wtd_init = 0.5 * (z₋ₕ[k_target] + z₋ₕ[k_target+1]) # 层中心 + + # 构造平衡含水量 (略低于饱和含水量) + θ_sat_approx = soil.θ_sat * 0.9 + θ_eq = fill(θ_sat_approx * 0.6, nzg) # 平衡含水量 < 饱和量 + + # 土壤接近饱和 + θ = fill(θ_sat_approx * 0.95, nzg) + θ_wtd = soil.θ_sat + + wtd_new, rech = updatewtd_shallow(nzg, 0, z₋ₕ, dz, soiltxt, θ_eq, θ_wtd, θ, wtd_init, fdepth) + + # 基本检查：函数应当返回有限值 + @test isfinite(wtd_new) + @test isfinite(rech) +end + +@testset "updatewtd_shallow自由排水模式" begin + nzg = 20 + z₋ₕ, dz = initializesoildepth(nzg) + soiltxt = 1 + soil = get_soil_params(soiltxt) + fdepth = 10.0 + + θ_eq = fill(soil.θ_sat * 0.6, nzg) + θ = fill(soil.θ_sat * 0.3, nzg) + θ_wtd = soil.θ_sat + + # 水位在最深层中心 + wtd_init = 0.5 * (z₋ₕ[1] + z₋ₕ[2]) + + # 自由排水 (freedrain=1) + wtd_new, rech = updatewtd_shallow(nzg, 1, z₋ₕ, dz, soiltxt, θ_eq, θ_wtd, θ, wtd_init, fdepth) + + @test isfinite(wtd_new) + @test isfinite(rech) +end + +@testset "cal_factor指数衰减函数" begin + fdepth = 2.0 + # z = -1.5 (即模型中的 z + 1.5 = 0): factor = exp(0) = 1.0 + @test ASAP.cal_factor(-1.5, fdepth) ≈ 1.0 + + # z = -0.5 (z + 1.5 = 1.0): factor = exp(0.5) > 1, 但 clamp 到 1.0 + f1 = ASAP.cal_factor(-0.5, fdepth) + @test f1 ≈ 1.0 # clamped to 1.0 + + # z = -3.5 (z + 1.5 = -1.0): factor = exp(-0.5) ≈ 0.607 + f2 = ASAP.cal_factor(-3.5, fdepth) + @test f2 ≈ exp(-1.0) atol=1e-10 + + # 非常深的层 → 接近下限 0.1 + f3 = ASAP.cal_factor(-100.0, fdepth) + @test f3 ≈ 0.1 + + # fdepth 增大 → 衰减更慢 + f_shallow = ASAP.cal_factor(-5.0, 1.0) + f_deep = ASAP.cal_factor(-5.0, 10.0) + @test f_deep > f_shallow +end diff --git a/test/wtable/runtests.jl b/test/wtable/runtests.jl index 59823ed..5380794 100644 --- a/test/wtable/runtests.jl +++ b/test/wtable/runtests.jl @@ -4,3 +4,4 @@ include("test_watertable.jl") include("test_isotope_tracing.jl") include("test_river_routing.jl") include("test_rivers_dw_flood.jl") +include("test-river.jl") diff --git a/test/wtable/test-river.jl b/test/wtable/test-river.jl index b1ceb0e..99c2c75 100644 --- a/test/wtable/test-river.jl +++ b/test/wtable/test-river.jl @@ -1,3 +1,5 @@ +using ASAP, Test + @testset "River-Groundwater Interaction" begin # 测试河流-地下水交换 imax, jmax = 3, 3