From 672672424ac174c9ffecb64026b34fa3afc465ce Mon Sep 17 00:00:00 2001 From: Dongdong Kong Date: Sun, 13 Oct 2024 17:22:25 +0800 Subject: [PATCH] test rainfall_stage --- src/Soil/UpdateSoilMoisture.jl | 1 - src/clang/BEPS_c.jl | 41 +++++++++++++---------- src/rainfall_stage.jl | 50 +++++++++++++--------------- src/snowpack.jl | 1 - test/modules/modules.jl | 1 + test/modules/test-rainfall_stage1.jl | 17 ++++++++++ test/test-beps_main.jl | 4 +-- 7 files changed, 68 insertions(+), 47 deletions(-) create mode 100644 test/modules/test-rainfall_stage1.jl diff --git a/src/Soil/UpdateSoilMoisture.jl b/src/Soil/UpdateSoilMoisture.jl index dfdf3ca..c5ca72d 100644 --- a/src/Soil/UpdateSoilMoisture.jl +++ b/src/Soil/UpdateSoilMoisture.jl @@ -106,7 +106,6 @@ end # Function to calculate soil water uptake from a layer function Soil_Water_Uptake(p::Soil, Trans_o::Float64, Trans_u::Float64, Evap_soil::Float64) - ρ_w = 1025.0 Trans = Trans_o + Trans_u p.Ett[1] = Trans / ρ_w * p.dt[1] + Evap_soil / ρ_w # for the top layer diff --git a/src/clang/BEPS_c.jl b/src/clang/BEPS_c.jl index 24abe85..0e579fd 100644 --- a/src/clang/BEPS_c.jl +++ b/src/clang/BEPS_c.jl @@ -19,13 +19,13 @@ include("module.jl") # include("beps_main.jl") # include("debug_Rln.jl") -# function plantresp(LC, mid_res::Results, lai_yr, lai, temp_air, temp_soil, CosZs) +# function plantresp(LC, mid_res::Results, lai_yr, lai, Tair, temp_soil, CosZs) # ccall((:plantresp, libbeps), Cvoid, (Cint, Ptr{Results}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), -# LC, Ref(mid_res), lai_yr, lai, temp_air, temp_soil, CosZs) +# LC, Ref(mid_res), lai_yr, lai, Tair, temp_soil, CosZs) # end function inter_prg_c(jday, rstep, - lai::T, clumping::T, parameter::Vector{T}, meteo::Met, CosZs::T, + lai::T, Ω::T, parameter::Vector{T}, meteo::Met, CosZs::T, var_o::Vector{T}, var_n::Vector{T}, soilp::Soil_c, Ra::Radiation, mid_res::Results, mid_ET::OutputET, var::InterTempVars; debug=false, kw...) where {T<:Real} @@ -33,15 +33,15 @@ function inter_prg_c(jday, rstep, ccall((:inter_prg_c, libbeps), Cvoid, (Cint, Cint, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Met}, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Soil_c}, Ptr{Results}, Ptr{OutputET}), - jday, rstep, lai, clumping, parameter, + jday, rstep, lai, Ω, parameter, Ref(meteo), CosZs, var_o, var_n, Ref(soilp), Ref(mid_res), Ref(mid_ET)) end -function Vcmax_Jmax(lai_o, clumping, Vcmax0, slope_Vcmax_N, leaf_N, CosZs, Vcmax_sunlit, Vcmax_shaded, Jmax_sunlit, Jmax_shaded) +function Vcmax_Jmax(lai_o, Ω, Vcmax0, slope_Vcmax_N, leaf_N, CosZs, Vcmax_sunlit, Vcmax_shaded, Jmax_sunlit, Jmax_shaded) ccall((:Vcmax_Jmax, libbeps), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), - lai_o, clumping, Vcmax0, slope_Vcmax_N, leaf_N, CosZs, + lai_o, Ω, Vcmax0, slope_Vcmax_N, leaf_N, CosZs, Vcmax_sunlit, Vcmax_shaded, Jmax_sunlit, Jmax_shaded) end @@ -82,37 +82,44 @@ end # passed test # lai2(0.8, 0.6, 0.2, 0.4, 0.2, 0.4) -function lai2(clumping, CosZs, stem_o, stem_u, lai_o, lai_u) +function lai2(Ω, CosZs, stem_o, stem_u, lai_o, lai_u) LAI = Leaf() PAI = Leaf() - lai2(clumping, CosZs, stem_o, stem_u, lai_o, lai_u, LAI, PAI) + lai2(Ω, CosZs, stem_o, stem_u, lai_o, lai_u, LAI, PAI) LAI, PAI end ## add Leaf struct -function lai2(clumping, CosZs, stem_o, stem_u, lai_o, lai_u, LAI::Leaf, PAI::Leaf) +function lai2(Ω, CosZs, stem_o, stem_u, lai_o, lai_u, LAI::Leaf, PAI::Leaf) ccall((:lai2, libbeps), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Leaf}, Ptr{Leaf}), - clumping, CosZs, stem_o, stem_u, lai_o, lai_u, Ref(LAI), Ref(PAI)) + Ω, CosZs, stem_o, stem_u, lai_o, lai_u, Ref(LAI), Ref(PAI)) end function meteo_pack(temp, rh, meteo_pack_output) ccall((:meteo_pack, libbeps), Cvoid, (Cdouble, Cdouble, Ptr{Cdouble}), temp, rh, meteo_pack_output) end -function rainfall_stage1(temp_air, prcp, mass_water_o_last, mass_water_u_last, lai_o, lai_u, clumping) - mass_water_o = init_dbl() - mass_water_u = init_dbl() +function rainfall_stage1(Tair, prcp, m_water_o_last, m_water_u_last, lai_o, lai_u, Ω) + m_water_o = init_dbl() + m_water_u = init_dbl() percent_water_o = init_dbl() percent_water_u = init_dbl() prcp_g = init_dbl() - ccall((:rainfall_stage1, libbeps), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), temp_air, prcp, mass_water_o_last, mass_water_u_last, lai_o, lai_u, clumping, mass_water_o, mass_water_u, percent_water_o, percent_water_u, prcp_g) + ccall((:rainfall_stage1, libbeps), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), Tair, prcp, m_water_o_last, m_water_u_last, lai_o, lai_u, Ω, m_water_o, m_water_u, percent_water_o, percent_water_u, prcp_g) - mass_water_o[], mass_water_u[], percent_water_o[], percent_water_u[], prcp_g[] + m_water_o[], m_water_u[], percent_water_o[], percent_water_u[], prcp_g[] end -function rainfall_stage2(evapo_water_o, evapo_water_u, mass_water_o::Ref, mass_water_u::Ref) - ccall((:rainfall_stage2, libbeps), Cvoid, (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), evapo_water_o, evapo_water_u, mass_water_o, mass_water_u) +function rainfall_stage1(Tair, prcp, perc_water, m_water, m_water_pre, lai_o, lai_u, Ω) + # m_water_o_last, m_water_u_last, + # perc_water = Layer2{Float64}() + # m_water = Layer2{Float64}() + prcp_g = rainfall_stage1_jl(Tair, prcp, perc_water, m_water, m_water_pre, lai_o, lai_u, Ω) +end + +function rainfall_stage2(evapo_water_o, evapo_water_u, m_water_o::Ref, m_water_u::Ref) + ccall((:rainfall_stage2, libbeps), Cvoid, (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), evapo_water_o, evapo_water_u, m_water_o, m_water_u) end include("snowpack_stage.jl") diff --git a/src/rainfall_stage.jl b/src/rainfall_stage.jl index 5fb0691..e18dc6f 100644 --- a/src/rainfall_stage.jl +++ b/src/rainfall_stage.jl @@ -1,42 +1,40 @@ -function rainfall_stage1_jl(Tair::Float64, prcp::Float64, - perc_water::Layer2{Float64}, - mass_water::Layer2{Float64}, - mass_water_pre::Layer2{Float64}, - lai_o::Float64, lai_u::Float64, clumping::Float64) +# - m : [kg m-2] +# - prcp : [m m-2 s-1] +function water_change(m_water_pre, prcp, lai, Ω) + mMax_water = 0.1 * lai + τ = 1 - exp(-lai * Ω) + m_water_o = m_water_pre + prcp * kstep * ρ_w * τ + m_water_o = clamp(m_water_o, 0, mMax_water) + + Δm_water_o = max(m_water_o - m_water_pre, 0.0) + perc_water_o = min(m_water_o / mMax_water, 1.0) + m_water_o, perc_water_o, Δm_water_o +end + +# [kg m-2] -> [m s-1] +kg2m(p) = p / ρ_w / kstep +m2kg(m) = m * kstep * ρ_w +# - m_water: change +function rainfall_stage1_jl(Tair::Float64, prcp::Float64, + perc_water::Layer2{Float64}, m_water::Layer2{Float64}, m_water_pre::Layer2{Float64}, + lai_o::Float64, lai_u::Float64, Ω::Float64) # Ta > 0, otherwise it is snow fall Tair <= 0.0 && (prcp = 0.0) prcp_o = prcp - prcp_g = 0.0 # overstorey - mass_water.o = mass_water_pre.o + prcp_o * kstep * ρ_w * (1 - exp(-lai_o * clumping)) - massMax_water_o = 0.1 * lai_o - - mass_water.o = clamp(mass_water.o, 0, massMax_water_o) - massStep_water_o = max(mass_water.o - mass_water_pre.o, 0.0) - - perc_water.o = min(mass_water.o / massMax_water_o, 1.0) - + m_water.o, perc_water.o, Δm_water_o = water_change(m_water_pre.o, prcp_o, lai_o, Ω) # understorey - prcp_u = prcp_o - massStep_water_o / ρ_w / kstep - mass_water.u = mass_water_pre.u + prcp_u * kstep * ρ_w * (1 - exp(-lai_u * clumping)) - massMax_water_u = 0.1 * lai_u - - mass_water.u = clamp(mass_water.u, 0, massMax_water_u) - massStep_water_u = max(mass_water.u - mass_water_pre.u, 0.0) - - perc_water.u = min(mass_water.u / massMax_water_u, 1.0) - prcp_g = prcp_u - massStep_water_u / ρ_w / kstep + prcp_u = prcp_o - Δm_water_o / ρ_w / kstep + m_water.u, perc_water.u, Δm_water_u = water_change(m_water_pre.u, prcp_u, lai_u, Ω) + prcp_g = prcp_u - Δm_water_u / ρ_w / kstep return prcp_g end - function rainfall_stage2_jl(evapo_water_o::Float64, evapo_water_u::Float64, mass_water::Layer2{Float64}) - - # kstep # 6min or 360s per step mass_water.o = max(mass_water.o - evapo_water_o * kstep, 0.0) mass_water.u = max(mass_water.u - evapo_water_u * kstep, 0.0) end diff --git a/src/snowpack.jl b/src/snowpack.jl index 1c02be7..032782c 100644 --- a/src/snowpack.jl +++ b/src/snowpack.jl @@ -148,7 +148,6 @@ function snowpack_stage3_jl(Tair::Float64, Tsnow::Float64, Tsnow_last::Float64, zs_sup = depth_snow # already considered sublimation ms_sup = m_snow.g - ρ_w = 1025 # [kg m-3] Δm = cal_melt(zs_sup, ρ_snow, Tsnow) # [kg m-2] con_melt = Tsnow > 0 && Tsnow_last <= 0 && ms_sup > 0 con_frozen = Tsnow <= 0 && Tsnow_last > 0 && depth_water > 0 diff --git a/test/modules/modules.jl b/test/modules/modules.jl index 756ed3f..87f6320 100644 --- a/test/modules/modules.jl +++ b/test/modules/modules.jl @@ -1,3 +1,4 @@ +include("test-rainfall_stage1.jl") include("test-surface_temperature.jl") include("test-radiation.jl") include("test-snowpack.jl") diff --git a/test/modules/test-rainfall_stage1.jl b/test/modules/test-rainfall_stage1.jl new file mode 100644 index 0000000..bc21ab4 --- /dev/null +++ b/test/modules/test-rainfall_stage1.jl @@ -0,0 +1,17 @@ +@testset "rainfall_stage1_jl" begin + Tair = 5.0 # Air temperature in Celsius + prcp = 10.0/1000 # Precipitation in mm + perc_water = Layer2{Float64}(0.1, 0.2) # Percentage water + m_water = Layer2{Float64}(0.01, 0.02) # Current water mass + m_water_pre = Layer2{Float64}(0.05, 0.005) # Previous water mass + lai_o = 1.0 # Leaf area index of overstory + lai_u = 1.0 # Leaf area index of understory + Ω = 0.5 # Clumping index + + # Call the function with the test inputs + prcp_g = rainfall_stage1_jl(Tair, prcp, perc_water, m_water, m_water_pre, lai_o, lai_u, Ω) + + r2 = clang.rainfall_stage1(Tair, prcp, perc_water, m_water, m_water_pre, lai_o, lai_u, Ω) + + @test prcp_g == r2 +end diff --git a/test/test-beps_main.jl b/test/test-beps_main.jl index 6077f13..f3bb4d3 100644 --- a/test/test-beps_main.jl +++ b/test/test-beps_main.jl @@ -16,7 +16,7 @@ par = (lon=120.5, lat=30.5, landcover=25, clumping=0.85, @testset "besp_main julia" begin d = deserialize(path_proj("data/p1_meteo")) d.tem = d.tem .- 5.0 - + @time df_jl, df_ET_jl = besp_main(d, lai, par; version="julia") @time df_c, df_ET_c = besp_main(d, lai, par; version="c") r = sum(df_jl) @@ -35,7 +35,7 @@ par = (lon=120.5, lat=30.5, landcover=25, clumping=0.85, @show l @show nanmax(l) @test true - + # @test nanmax(l) < 1.5 # SH, 1.48%的误差, current 0.09% @test nanmax(l) < 2.5 # GPP, 2.38%的误差, current 0.09% end