add InterTempLeafs

examples/BEPS_compare_with_C.qmd

@@ -7,12 +7,18 @@
 
 ```{julia}
 using BEPS
+using BenchmarkTools
+
 d = deserialize("data/p1_meteo")
 lai = readdlm("examples/input/p1_lai.txt")[:]
 
 par = (lon=120.5, lat=30.5, landcover=25, clumping=0.85,
   soil_type=8, Tsoil=2.2,
   soilwater=0.4115, snowdepth=0.0)
+
+@time df_jl, df_ET_jl = besp_main(d, lai, par; version="julia");
+@btime df_jl, df_ET_jl = besp_main(d, lai, par; version="julia");
+
 ```
 
 ```{julia}

src/DataType/InterTempVars.jl

@@ -1,3 +1,69 @@
+export InterTempLeafs
+
+@with_kw mutable struct InterTempLeafs
+  x0::Float64 = 0.0
+  Cc_new::Leaf = Leaf(x0)
+  Cs_old::Leaf = Leaf(x0)
+  Cs_new::Leaf = Leaf(x0)
+  Ci_old::Leaf = Leaf(x0)
+  Tc_old::Leaf = Leaf(x0)
+  Tc_new::Leaf = Leaf(x0)
+  Gs_old::Leaf = Leaf(x0)
+
+  # to the reference height above the canopy
+  Gc::Leaf = Leaf(x0)  # total conductance for CO2 from the intercellular space of the leaves
+  Gh::Leaf = Leaf(x0)  # total conductance for heat transfer from the leaf surface 
+  Gw::Leaf = Leaf(x0)  # total conductance for water from the intercellular space of the leaves
+  Gww::Leaf = Leaf(x0) # total conductance for water from the surface of the leaves
+
+  Gs_new::Leaf = Leaf(x0)
+  Ac::Leaf = Leaf(x0)
+  Ci_new::Leaf = Leaf(x0)
+
+  Rn::Leaf = Leaf(x0)
+  Rns::Leaf = Leaf(x0)
+  Rnl::Leaf = Leaf(x0)
+
+  leleaf::Leaf = Leaf(x0)
+  GPP::Leaf = Leaf(x0)
+  LAI::Leaf = Leaf(x0)
+  PAI::Leaf = Leaf(x0)
+end
+
+InterTempLeafs(x0) = InterTempLeafs(; x0)
+
+function reset!(l::InterTempLeafs)
+  # reset!(l.Cc_new)
+  # reset!(l.Cs_old)
+  # reset!(l.Cs_new)
+  # reset!(l.Ci_old)
+  # reset!(l.Tc_old)
+  # reset!(l.Tc_new)
+  # reset!(l.Gs_old)
+  # reset!(l.Gc)
+  # reset!(l.Gh)
+  # reset!(l.Gw)
+  # reset!(l.Gww)
+  # reset!(l.Gs_new)
+  # reset!(l.Ac)
+  # reset!(l.Ci_new)
+  # reset!(l.Rn)
+  # reset!(l.Rns)
+  # reset!(l.Rnl)
+  # reset!(l.leleaf)
+  # reset!(l.GPP)
+  # reset!(l.LAI)
+  # reset!(l.PAI)
+
+  # names = fieldnames(InterTempLeafs)[2:end]
+  # for name in names
+  #   x = getfield(l, name)
+  #   reset(x)
+  # end
+end
+
+
+
 @with_kw mutable struct InterTempVars
   Tc_u::Vector{FT} = zeros(MAX_Loop)
   Ts0::Vector{FT} = zeros(MAX_Loop)
@@ -36,10 +102,12 @@
   lambda_snow::Vector{FT} = zeros(MAX_Loop)
 
   # 记录土壤温度
-
   Cs::Matrix{FT} = zeros(layer + 2, MAX_Loop)
   Tm::Matrix{FT} = zeros(layer + 2, MAX_Loop)
   G::Matrix{FT} = zeros(layer + 2, MAX_Loop)
+
+  # Leafs
+  TempLeafs::InterTempLeafs = InterTempLeafs(0.0)
 end
 
 function init_vars!(x::InterTempVars)
@@ -66,7 +134,6 @@ function init_vars!(x::InterTempVars)
 
   # x.Ac_snow_o .= 0.
   # x.Ac_snow_u .= 0.  
-
   # x.alpha_v_sw .= 0.
   # x.alpha_n_sw .= 0.
   # x.r_rain_g .= 0.

src/DataType/Leaf.jl

@@ -1,3 +1,10 @@
+export @reset, reset!
+macro reset(x)
+  quote
+    reset!($(esc(x)))
+  end
+end
+
 @with_kw mutable struct Leaf
   o_sunlit::Cdouble = 0.0
   o_shaded::Cdouble = 0.0
@@ -26,6 +33,8 @@ for (m, f) in Base_ops
 end
 
 
+reset!(x::Leaf) = init_leaf_dbl(x, 0.0)
+
 function init_leaf_struct(x::Leaf, replacement::Leaf)
   x.o_sunlit = replacement.o_sunlit
   x.o_shaded = replacement.o_shaded
@@ -39,6 +48,7 @@ function init_leaf_dbl(x::Leaf, replacement::Float64)
   x.o_shaded = replacement
   x.u_sunlit = replacement
   x.u_shaded = replacement
+  nothing
   # ccall((:init_leaf_dbl, libbeps), Cvoid, (Ptr{Leaf}, Cdouble), Ref(x), replacement)
 end
 

src/beps_inter_prg.jl

@@ -25,35 +25,40 @@ function inter_prg_jl(
   # var = var2
   # var = InterTempVars()
   init_vars!(var)
+  reset!(var.TempLeafs)
+  @unpack Cc_new, Cs_old, Cs_new, Ci_old, 
+    Tc_old, Tc_new, Gs_old, Gc, Gh, Gw, Gww, 
+    Gs_new, Ac, Ci_new, Rn, Rns, Rnl, 
+    leleaf, GPP, LAI, PAI = var.TempLeafs
+
+  # Cc_new = Leaf()
+  # Cs_old = Leaf()
+  # Cs_new = Leaf()
+  # Ci_old = Leaf()
+  # Tc_old = Leaf()
+  # Tc_new = Leaf()
+  # Gs_old = Leaf()
+  # # to the reference height above the canopy
+  # Gc = Leaf()  # total conductance for CO2 from the intercellular space of the leaves
+  # Gh = Leaf()  # total conductance for heat transfer from the leaf surface 
+  # Gw = Leaf()  # total conductance for water from the intercellular space of the leaves
+  # Gww = Leaf() # total conductance for water from the surface of the leaves
+
+  # Gs_new = Leaf()
+  # Ac = Leaf()
+  # Ci_new = Leaf()
+
+  # Rn = Leaf()
+  # Rns = Leaf()
+  # Rnl = Leaf()
+  # leleaf = Leaf()
+  # GPP = Leaf()
+  # LAI = Leaf()
+  # PAI = Leaf()
 
   d_soil = zeros(layer + 1)
   lambda = zeros(layer + 2)
 
-  Cc_new = Leaf()
-  Cs_old = Leaf()
-  Cs_new = Leaf()
-  Ci_old = Leaf()
-  Tc_old = Leaf()
-  Tc_new = Leaf()
-  Gs_old = Leaf()
-
-  # to the reference height above the canopy
-  Gc = Leaf()  # total conductance for CO2 from the intercellular space of the leaves
-  Gh = Leaf()  # total conductance for heat transfer from the leaf surface 
-  Gw = Leaf()  # total conductance for water from the intercellular space of the leaves
-  Gww = Leaf() # total conductance for water from the surface of the leaves
-
-  Gs_new = Leaf()
-  Ac = Leaf()
-  Ci_new = Leaf()
-
-  Rn = Leaf()
-  Rns = Leaf()
-  Rnl = Leaf()
-
-  leleaf = Leaf()
-  GPP = Leaf()
-
   ra_o = 0.0
   ra_u = 0.0
   ra_g = 0.0
@@ -69,11 +74,6 @@ function inter_prg_jl(
   Tco = 0.0
   Tcu = 0.0
 
-  # # parameters for Vcmax-Nitrogen calculations
-  Kn = 0.3  #0.713/2.4
-  G_theta = 0.5
-  gamma = 0.066
-
   alpha_sat = param[24+1]       # albedo of saturated/dry soil for module rainfall 1
   alpha_dry = param[25+1]       # the albedo of dry soil
   canopyh_o = param[29+1]       # to be used for module aerodynamic_conductance
@@ -84,8 +84,12 @@ function inter_prg_jl(
 
   # /*****  Vcmax-Nitrogen calculations，by G.Mo，Apr. 2011  *****/
   if (CosZs > 0) # day time
-    K = G_theta * clumping / CosZs  # G_theta = 0.5 assuming a spherical leaf angle distribution
+    # parameters for Vcmax-Nitrogen calculations
+    G_theta = 0.5 # assuming a spherical leaf angle distribution
+    Kn = 0.3      # 0.713/2.4
+    K = G_theta * clumping / CosZs
     Vcmax0 = param[36+1]
+
     expr1 = 1 - exp(-K * lai)
     expr2 = 1 - exp(-lai * (Kn + K))
     expr3 = 1 - exp(-Kn * lai)
@@ -113,19 +117,16 @@ function inter_prg_jl(
   else
     lai_u = 1.18 * exp(-0.99 * lai_o)
   end
-
   (lai_u > lai_o) && (lai_u = 0.01)
 
   stem_o = param[8+1] * 0.2    # parameter[8].LAI max overstory
   stem_u = param[9+1] * 0.2    # parameter[9].LAI max understory
 
   # lai2: separate lai into sunlit and shaded portions
-  LAI = Leaf()
-  PAI = Leaf()
   lai2(clumping, CosZs, stem_o, stem_u, lai_o, lai_u, LAI, PAI)
 
   # /*****  Initialization of this time step  *****/
-  Ks = meteo.Srad
+  Rs = meteo.Srad
   rh_air = meteo.rh
   wind_sp = meteo.wind
   prcp = meteo.rain / step  # precipitation in meters
@@ -136,10 +137,11 @@ function inter_prg_jl(
   VPD_air = es - ea                        # water vapor deficit at the reference height
 
   q_ca = 0.622 * ea / (101.35 - 0.378 * ea)  # in g/g, unitless
-  Cp_ca = Cpd * (1 + 0.84 * q_ca)
+  cp = Cpd * (1 + 0.84 * q_ca)
   slope = cal_slope(Ta)
+  gamma = 0.066
 
-  if (Ks <= 0)
+  if (Rs <= 0)
     alpha_v_o = 0.0
     alpha_n_o = 0.0
     alpha_v_u = 0.0
@@ -210,8 +212,8 @@ function inter_prg_jl(
     init_leaf_dbl(Ci_old, 0.7 * CO2_air)
     init_leaf_dbl2(Gs_old, 1.0 / 200.0, 1.0 / 300.0)
 
-    percentArea_snow_o = var.Ac_snow_o[kkk] / lai_o / 2
-    percentArea_snow_u = var.Ac_snow_u[kkk] / lai_u / 2
+    percArea_snow_o = var.Ac_snow_o[kkk] / lai_o / 2
+    percArea_snow_u = var.Ac_snow_u[kkk] / lai_u / 2
 
     temp_grd = Ta   # ground temperature substituted by air temperature
 
@@ -236,18 +238,17 @@ function inter_prg_jl(
       Tcu = (Tc_old.u_sunlit * PAI.u_sunlit + Tc_old.u_shaded * PAI.u_shaded) / (PAI.u_sunlit + PAI.u_shaded)
 
       # /*****  Net radiation at canopy and leaf level module by X.Luo  *****/
-      radiation_o, radiation_u, radiation_g = netRadiation_jl(Ks, CosZs, Tco, Tcu, temp_grd,
+      radiation_o, radiation_u, radiation_g = netRadiation_jl(Rs, CosZs, Tco, Tcu, temp_grd,
         lai_o, lai_u, lai_o + stem_o, lai_u + stem_u, PAI,
         clumping, Ta, rh_air, var.alpha_v_sw[kkk], var.alpha_n_sw[kkk],
-        percentArea_snow_o, percentArea_snow_u,
+        percArea_snow_o, percArea_snow_u,
         var.Xg_snow[kkk],
         alpha_v_o, alpha_n_o, alpha_v_u, alpha_n_u,
         alpha_v_g, alpha_n_g, Rn, Rns, Rnl, Ra)
 
       # /*****  Photosynthesis module by B. Chen  *****/
-      # conductance for water
-      update_Gw!(Gw, Gs_old, Ga_o, Ga_u, Gb_o, Gb_u)
-      latent_heat!(leleaf, Gw, VPD_air, slope, Tc_old, Ta, rho_a, Cp_ca, gamma)
+      update_Gw!(Gw, Gs_old, Ga_o, Ga_u, Gb_o, Gb_u) # conductance for water
+      latent_heat!(leleaf, Gw, VPD_air, slope, Tc_old, Ta, rho_a, cp, gamma)
 
       if (CosZs > 0)
         photosynthesis(Tc_old, Rns, Ci_old, leleaf,
@@ -271,13 +272,13 @@ function inter_prg_jl(
       update_Gc!(Gc, Gs_new, Ga_o, Ga_u, Gb_o, Gb_u)
 
       # /***** Leaf temperatures module by L. He  *****/
-      Leaf_Temperatures_jl(Ta, slope, gamma, VPD_air, Cp_ca,
+      Leaf_Temperatures_jl(Ta, slope, gamma, VPD_air, cp,
         Gw, Gww, Gh,
         var.Xcs_o[kkk], var.Xcl_o[kkk], var.Xcs_u[kkk], var.Xcl_u[kkk],
         Rn, Tc_new)
 
-      H_o_sunlit = (Tc_new.o_sunlit - Ta) * rho_a * Cp_ca * Gh.o_sunlit
-      H_o_shaded = (Tc_new.o_shaded - Ta) * rho_a * Cp_ca * Gh.o_shaded
+      H_o_sunlit = (Tc_new.o_sunlit - Ta) * rho_a * cp * Gh.o_sunlit
+      H_o_shaded = (Tc_new.o_shaded - Ta) * rho_a * cp * Gh.o_shaded
       GH_o = H_o_sunlit * PAI.o_sunlit + H_o_shaded * PAI.o_shaded  # for next num aerodynamic conductance calculation
 
       if (abs(Tc_new.o_sunlit - Tc_old.o_sunlit) < 0.02 &&

src/beps_main.jl

@@ -50,11 +50,10 @@ function besp_main(d::DataFrame, lai::Vector, par::NamedTuple;
       end
 
       CosZs = s_coszs(jday, rstep - 1, par.lat, par.lon) # cos_solar zenith angle
+
       # /***** start simulation modules *****/
-
       fun(jday, rstep - 1, _lai, par.clumping, param, meteo, CosZs, var_o, var_n, p_soil,
         Ra, mid_res, mid_ET, var; debug)
-      # inter_prg_jl(jday, rstep - 1, _lai, par.clumping, parameter, meteo, CosZs, var_o, var_n, p_soil, mid_res)
       # Store updated variables array in temp array
       v2last .= var_n