test rainfall_stage

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

src/clang/BEPS_c.jl

@@ -19,29 +19,29 @@ 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}
 
   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")

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

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

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")

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

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