improve inter_prg_jl

README.md

@@ -45,14 +45,17 @@ 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");
+@time df_jl, df_ET_jl, Tg = besp_main(d, lai, par; version="julia");
 ```
 
 > Figure1: The bias of Julia version compared with C, `bias = (Julia - C)/ C`.
-![Figure1: bias of julia version compared with c](./images/Figure1_bias_of_julia-version.png)
+![](./images/Figure1_bias_of_julia-version.png)
 
 The bias of `SH` is the largest due to double number accuracy, about 1.48%, which is acceptable.
 
+> Figure2: The variation of soil temperature at different depths.
+![](./images/Figure2_variation_of_Tg.png)
+
 See [examples/example_01.qmd](examples/example_01.qmd) for details.
 
 ## TODO

examples/example_01.qmd

@@ -16,8 +16,8 @@ par = (lon=120.5, lat=30.5, landcover=25, clumping=0.85,
 ```
 
 ```{julia}
-@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");
+@time df_jl, df_ET_jl, Tg = besp_main(d, lai, par; version="julia");
+@time df_c, df_ET_c, Tg = besp_main(d, lai, par; version="c");
 
 sum(df_jl)
 sum(df_c)
@@ -53,3 +53,24 @@ ps = map(plot_var, vars)
 plot(ps..., layout=(3, 4), size=(1400, 800))
 # write_fig("images/Figure1_bias_of_julia-version.png", 15, 8; show=false)
 ```
+
+
+```{julia}
+## 土壤温度的变化
+begin
+  using Plots
+  gr(framestyle=:box)
+
+  depths = [0.05, 0.10, 0.20, 0.40, 1.25] |> cumsum |> x -> round.(x, digits=4)
+  
+  p = plot(; size=(1400, 700), title="Soil temperature")
+  for i = 1:5
+    label = "depth: $(depths[i])"
+    plot!(Tg[:, i]; label)
+  end
+  p
+end
+
+write_fig("images/Figure2_variation_of_Tg.png", 15, 8; show=false)
+
+```

src/DataType/InterTempVars.jl

@@ -32,35 +32,13 @@ 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
+# function reset!(l::InterTempLeafs)
+#   names = fieldnames(InterTempLeafs)[2:end]
+#   for name in names
+#     x = getfield(l, name)
+#     reset!(x)
+#   end
+# end
 
 
 

src/DataType/Other_structs.jl

@@ -21,6 +21,7 @@ ClimateData(Srad, LR, temp, rh, rain, wind) =
   gpp_u_sunlit::Cdouble = 0.0
   gpp_o_shaded::Cdouble = 0.0
   gpp_u_shaded::Cdouble = 0.0
+
   plant_resp::Cdouble = 0.0
   npp_o::Cdouble = 0.0
   npp_u::Cdouble = 0.0

src/beps_inter_prg.jl

@@ -25,7 +25,7 @@ function inter_prg_jl(
   # var = var2
   # var = InterTempVars()
   init_vars!(var)
-  reset!(var.TempLeafs)
+  # 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, 
@@ -89,7 +89,7 @@ function inter_prg_jl(
     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)
@@ -126,15 +126,15 @@ function inter_prg_jl(
   lai2(clumping, CosZs, stem_o, stem_u, lai_o, lai_u, LAI, PAI)
 
   # /*****  Initialization of this time step  *****/
-  Rs = meteo.Srad
-  rh_air = meteo.rh
+  Rs      = meteo.Srad
+  rh_air  = meteo.rh
   wind_sp = meteo.wind
-  prcp = meteo.rain / step  # precipitation in meters
-  Ta = meteo.temp
+  prcp    = meteo.rain / step  # precipitation in meters
+  Ta      = meteo.temp
 
   es = cal_es(Ta)  # to estimate saturated water vapor pressure in kpa
-  ea = es * rh_air / 100                   # to be used for module photosynthesis
-  VPD_air = es - ea                        # water vapor deficit at the reference height
+  ea = es * rh_air / 100                   # used in `photosynthesis`
+  VPD = es - ea                            # water vapor deficit at the reference height
 
   q_ca = 0.622 * ea / (101.35 - 0.378 * ea)  # in g/g, unitless
   cp = Cpd * (1 + 0.84 * q_ca)
@@ -154,8 +154,8 @@ function inter_prg_jl(
   end
 
   # Ground surface temperature
-  var.Ts0[1] = clamp(var_o[3+1], Ta - 2.0, Ta + 2.0)
-  var.Tsm0[1] = clamp(var_o[5+1], Ta - 2.0, Ta + 2.0)
+  var.Ts0[1] = clamp(var_o[3+1], Ta - 2.0, Ta + 2.0)  # ground0
+  var.Tsm0[1] = clamp(var_o[5+1], Ta - 2.0, Ta + 2.0) # 
   var.Tsn0[1] = clamp(var_o[4+1], Ta - 2.0, Ta + 2.0) # snow0
   var.Tsn1[1] = clamp(var_o[6+1], Ta - 2.0, Ta + 2.0) # snow1
   var.Tsn2[1] = clamp(var_o[7+1], Ta - 2.0, Ta + 2.0) # snow2
@@ -248,7 +248,7 @@ function inter_prg_jl(
 
       # /*****  Photosynthesis module by B. Chen  *****/
       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)
+      latent_heat!(leleaf, Gw, VPD, slope, Tc_old, Ta, rho_a, cp, gamma)
 
       if (CosZs > 0)
         photosynthesis(Tc_old, Rns, Ci_old, leleaf,
@@ -272,7 +272,7 @@ 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,
+      Leaf_Temperatures_jl(Ta, slope, gamma, VPD, cp,
         Gw, Gww, Gh,
         var.Xcs_o[kkk], var.Xcl_o[kkk], var.Xcs_u[kkk], var.Xcl_u[kkk],
         Rn, Tc_new)
@@ -343,9 +343,11 @@ function inter_prg_jl(
       surface_temperature_jl(Ta, rh_air, Zsp[], Zp[],
         var.Cs[2, kkk], var.Cs[1, kkk], Gheat_g, d_soil[2], var.rho_snow[kkk], var.Tc_u[kkk],
         radiation_g, var.Evap_soil[kkk], var.Evap_SW[kkk], var.Evap_SS[kkk],
-        lambda[2], var.Xg_snow[kkk],
-
-        var.G[2, kkk], var.Ts0[kkk-1], var.Tm[2, kkk-1], var.Tm[1, kkk-1], var.Tsm0[kkk-1], 
+        lambda[2], 
+        var.Xg_snow[kkk], var.G[2, kkk], 
+        var.Ts0[kkk-1], 
+        # T_soil1_last::FT, T_any0_last::FT, T_soil0_last::FT,
+        var.Tm[2, kkk-1], var.Tm[1, kkk-1], var.Tsm0[kkk-1], 
         var.Tsn0[kkk-1], var.Tsn1[kkk-1], var.Tsn2[kkk-1])
 
     Update_temp_soil_c(soilp, var.Tm[1, kkk])

src/beps_main.jl

@@ -17,6 +17,7 @@ function besp_main(d::DataFrame, lai::Vector, par::NamedTuple;
 
   df_out = DataFrame(zeros(n, length(vars)), vars)
   df_ET = DataFrame(zeros(n, length(vars_ET)), vars_ET)
+  Tg = zeros(n, layer)
 
   var_o = zeros(41)
   var_n = zeros(41)
@@ -50,16 +51,18 @@ 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)
+
+      Tg[k, :] .= p_soil.temp_soil_c[1:layer]
       # Store updated variables array in temp array
       v2last .= var_n
 
       fill_res!(df_out, mid_res, k)
       fill_res!(df_ET, mid_ET, k)
     end # End of hourly loop
   end # End of daily loop
-  df_out, df_ET
+  df_out, df_ET, Tg
 end