minor

src/clang/module.jl

@@ -35,10 +35,11 @@ function photosynthesis_c(temp_leaf_p::Cdouble, rad_leaf::Cdouble, e_air::Cdoubl
 end
 
 function netRadiation_c(shortRad_global, CosZs,
-  temp_o, temp_u, temp_g,
+  T::Layer3{FT},
   lai_o, lai_u, lai_os, lai_us, lai::Leaf, Ω, temp_air, rh,
   α_snow_v, α_snow_n, α_v::Layer3{FT}, α_n::Layer3{FT},
-  percArea_snow_o, percArea_snow_u, perc_snow_g, Rn_Leaf::Leaf, Rns_Leaf::Leaf, Rnl_Leaf::Leaf, Ra::Radiation)
+  percArea_snow_o, percArea_snow_u, perc_snow_g, 
+  Rn_Leaf::Leaf, Rns_Leaf::Leaf, Rnl_Leaf::Leaf, Ra::Radiation)
 
   netRad_o = init_dbl()
   netRad_u = init_dbl()
@@ -47,7 +48,8 @@ function netRadiation_c(shortRad_global, CosZs,
   ccall((:netRadiation, libbeps), Cvoid,
     (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Leaf, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble,
       Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Leaf}, Ptr{Leaf}),
-    shortRad_global, CosZs, temp_o, temp_u, temp_g, lai_o, lai_u, lai_os, lai_us, lai, Ω, temp_air, rh, α_snow_v, α_snow_n, percArea_snow_o, percArea_snow_u, perc_snow_g,
+    shortRad_global, CosZs, T.o, T.u, T.g, 
+    lai_o, lai_u, lai_os, lai_us, lai, Ω, temp_air, rh, α_snow_v, α_snow_n, percArea_snow_o, percArea_snow_u, perc_snow_g,
     α_v.o, α_n.o, α_v.u, α_n.u, α_v.g, α_n.g,
     netRad_o, netRad_u, netRad_g, Ref(Rn_Leaf), Ref(Rns_Leaf))
 

src/evaporation_soil.jl

@@ -1,5 +1,4 @@
 """
-
 # Return
 - `Ewater`: evaporation from water
 - `Esoil`: evaporation from soil

src/inter_prg.jl

@@ -45,9 +45,6 @@ function inter_prg_jl(
   radiation_u = 0.0
   radiation_g = 0.0
 
-  Tco = 0.0
-  Tcu = 0.0
-
   # /*****  Vcmax-Nitrogen calculations，by G.Mo，Apr. 2011  *****/
   VCmax25, N_leaf, slope = param[37], param[47], param[48]
   Vcmax_sunlit, Vcmax_shaded = VCmax(lai, Ω, CosZs, VCmax25, N_leaf, slope)
@@ -100,7 +97,7 @@ function inter_prg_jl(
 
   m_snow_pre = state.m_snow# mass_snow
   m_water_pre = state.m_water
-  
+
   z_snow = soil.z_snow
   z_water = soil.z_water
   (z_water < 0.001) && (z_water = 0.0)
@@ -110,7 +107,7 @@ function inter_prg_jl(
   f_snow = Layer3(0.0)    # perc_snow
   m_snow = Layer3(0.0) # mass_snow
   A_snow = Layer2()
-  
+
   # the mass of intercepted liquid water and snow, overstory 
   f_water = Layer2() # perc_water
   m_water = Layer2() # mass_water
@@ -126,8 +123,9 @@ function inter_prg_jl(
   ρ_snow = init_dbl()
   α_v_sw = init_dbl()
   α_n_sw = init_dbl()
+  Tc = Layer3() # temperatures of o, u, g
 
-  # /*****  Ten time intervals in a hourly time step.6min or 360s per loop  ******/
+  # ten time intervals in a hourly time step.6min or 360s per loop
   @inbounds for k = 2:kloop+1
     ρ_snow[] = 0.0 # TODO: debug C, might error
     α_v_sw[], α_n_sw[] = 0.0, 0.0
@@ -163,7 +161,7 @@ function inter_prg_jl(
     perc_snow_o = A_snow.o / lai_o / 2 # area
     perc_snow_u = A_snow.u / lai_u / 2 # area
 
-    temp_grd = Ta   # ground temperature substituted by air temperature
+    Tc.g = Ta   # ground temperature substituted by air temperature
 
     num = 0
     while true # iteration for BWB equation until results converge
@@ -181,17 +179,18 @@ function inter_prg_jl(
         1.0 / (1.0 / Ga_u + 1.0 / Gb_u + 100))
 
       # temperatures of overstorey and understorey canopies
-      Tco = (Tc_old.o_sunlit * PAI.o_sunlit + Tc_old.o_shaded * PAI.o_shaded) /
-            (PAI.o_sunlit + PAI.o_shaded)
-      Tcu = (Tc_old.u_sunlit * PAI.u_sunlit + Tc_old.u_shaded * PAI.u_shaded) /
-            (PAI.u_sunlit + PAI.u_shaded)
+      Tc.o = (Tc_old.o_sunlit * PAI.o_sunlit + Tc_old.o_shaded * PAI.o_shaded) /
+             (PAI.o_sunlit + PAI.o_shaded)
+      Tc.u = (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(Rs, CosZs, Tco, Tcu, temp_grd,
+      radiation_o, radiation_u, radiation_g = netRadiation_jl(
+        Rs, CosZs, Tc, 
         lai_o, lai_u, lai_o + stem_o, lai_u + stem_u, PAI,
         Ω, Ta, RH,
         α_v_sw[], α_n_sw[], α_v, α_n,
-        perc_snow_o, perc_snow_o, f_snow.g, 
+        perc_snow_o, perc_snow_u, f_snow.g,
         Rn, Rns, Rnl, Ra)
 
       # /*****  Photosynthesis module by B. Chen  *****/
@@ -261,7 +260,7 @@ function inter_prg_jl(
     mass_water_g = ρ_w * z_water
 
     var.Evap_soil[k], var.Evap_SW[k], var.Evap_SS[k], z_water, z_snow =
-      evaporation_soil_jl(temp_grd, var.Ts0[k-1], RH, radiation_g, Gheat_g,
+      evaporation_soil_jl(Tc.g, var.Ts0[k-1], RH, radiation_g, Gheat_g,
         f_snow,
         z_water, z_snow, mass_water_g, m_snow, # Ref
         ρ_snow[], soil.θ_prev[1], soil.θ_sat[1])
@@ -277,7 +276,7 @@ function inter_prg_jl(
 
     var.Cs[1, k] = soil.Cs[1]
     var.Cs[2, k] = soil.Cs[1]
-    var.Tc_u[k] = Tcu
+    var.Tc_u[k] = Tc.u
     lambda[2] = soil.lambda[1]
     dz[2] = soil.dz[1]
 
@@ -336,7 +335,7 @@ function inter_prg_jl(
   state.Qhc_o = var.Qhc_o[k]
   set!(state.m_water, m_water)
   set!(state.m_snow, m_snow)
-  
+
   mid_res.Net_Rad = radiation_o + radiation_u + radiation_g
 
   OutputET!(mid_ET,
@@ -346,7 +345,7 @@ function inter_prg_jl(
     var.Evap_soil, var.Evap_SW, var.Evap_SS, var.Qhc_o, var.Qhc_u, var.Qhg, k)
   update_ET!(mid_ET, mid_res, Ta)
 
-  mid_res.gpp_o_sunlit = GPP.o_sunlit   # umol C/m2/s
+  mid_res.gpp_o_sunlit = GPP.o_sunlit   # [umol C/m2/s], !note
   mid_res.gpp_u_sunlit = GPP.u_sunlit
   mid_res.gpp_o_shaded = GPP.o_shaded
   mid_res.gpp_u_shaded = GPP.u_shaded
@@ -355,7 +354,7 @@ function inter_prg_jl(
   mid_res.z_snow = z_snow
   mid_res.ρ_snow = ρ_snow[]
 
-  # total GPP . gC/m2/step
-  mid_res.GPP = (GPP.o_sunlit + GPP.o_shaded + GPP.u_sunlit + GPP.u_shaded) * 12 * step * 0.000001
+  GPP = GPP.o_sunlit + GPP.o_shaded + GPP.u_sunlit + GPP.u_shaded
+  mid_res.GPP = GPP * 12 * step * 1e-6 # [umol m-2 s-1] -> [gC m-2]
   nothing
 end

src/netRadiation.jl

@@ -1,13 +1,19 @@
+"""
+## Arguments
+- `T`                 : temperature of o, u, g
+- `α_v, α_n`          : albedo of visible, near infrared
+- `percentArea_snow_o`: percentage of snow on overstorey (by area)
+- `percentArea_snow_u`: percentage of snow on understorey (by area)
+- `percent_snow_g`    : percentage of snow on ground (by mass)
+"""
 function netRadiation_jl(Rs_global::FT, CosZs::FT,
-  temp_o::FT, temp_u::FT, temp_g::FT,
+  T::Layer3{FT}, 
   lai_o::FT, lai_u::FT, lai_os::FT, lai_us::FT,
   lai::Leaf,
   Ω::FT, Tair::FT, RH::FT,
   α_snow_v::FT, α_snow_n::FT, α_v::Layer3{FT}, α_n::Layer3{FT},
   percArea_snow_o::FT, percArea_snow_u::FT, perc_snow_g::FT,
-  Rn_Leaf::Leaf,
-  Rns_Leaf::Leaf,
-  Rnl_Leaf::Leaf, Ra::Radiation)
+  Rn_Leaf::Leaf, Rns_Leaf::Leaf, Rnl_Leaf::Leaf, Ra::Radiation)
 
   # calculate α of canopy in this step
   α_v_os::FT = α_v.o * (1.0 - percArea_snow_o) + α_snow_v * percArea_snow_o  # visible, overstory
@@ -80,7 +86,7 @@ function netRadiation_jl(Rs_global::FT, CosZs::FT,
   Rns_u = Ra.Rns_u_dir + Ra.Rns_u_df
   Rns_g = Ra.Rns_g_dir + Ra.Rns_g_df
 
-  Rnl_o, Rnl_u, Rnl_g = cal_Rln_Longwave(Tair, RH, temp_o, temp_u, temp_g, lai_o, lai_u, Ω)
+  Rnl_o, Rnl_u, Rnl_g = cal_Rln_Longwave(Tair, RH, T, lai_o, lai_u, Ω)
 
   # 计算植被和地面的总净辐射
   Rn_o = Rns_o + Rnl_o
@@ -128,8 +134,8 @@ function netRadiation_jl(Rs_global::FT, CosZs::FT,
 end
 
 
-function cal_Rln_Longwave(temp_air::FT, rh::FT, temp_o::FT, 
-  temp_u::FT, temp_g::FT, lai_o::FT, lai_u::FT, Ω::FT) where {FT<:Real}
+function cal_Rln_Longwave(temp_air::FT, rh::FT, T::Layer3{FT}, 
+  lai_o::FT, lai_u::FT, Ω::FT) where {FT<:Real}
 
   # indicators to describe leaf distribution angles in canopy. slightly related with LAI
   cosQ_o::FT = 0.537 + 0.025 * lai_o  # Luo2018, A10, a representative zenith angle for diffuse radiation transmission
@@ -149,9 +155,9 @@ function cal_Rln_Longwave(temp_air::FT, rh::FT, temp_o::FT,
 
   # 计算植被和地面的净长波辐射
   Rl_air = cal_Rln(ϵ_air, temp_air)
-  Rl_o = cal_Rln(ϵ_o, temp_o)
-  Rl_u = cal_Rln(ϵ_u, temp_u)
-  Rl_g = cal_Rln(ϵ_g, temp_g)
+  Rl_o = cal_Rln(ϵ_o, T.o)
+  Rl_u = cal_Rln(ϵ_u, T.u)
+  Rl_g = cal_Rln(ϵ_g, T.g)
 
   Rnl_o = (ϵ_o * (Rl_air + Rl_u * (1.0 - τ_u_df) + Rl_g * τ_u_df) - 2 * Rl_o) *
           (1.0 - τ_o_df) +

test/modules/test-radiation.jl

@@ -3,16 +3,12 @@ using Test, BEPS
 @testset "netRadiation_jl" begin
   Rs_global = 500.0  # Global solar radiation
   CosZs = 0.5  # Cosine of the solar zenith angle
-  temp_o = 300.0  # Overstory temperature in Kelvin
-  temp_u = 295.0  # Understory temperature in Kelvin
-  temp_g = 290.0  # Ground temperature in Kelvin
   lai_o = 2.0  # Leaf area index of overstory
   lai_u = 1.0  # Leaf area index of understory
   lai_os = 1.5  # Leaf area index of overstory snow
   lai_us = 0.5  # Leaf area index of understory snow
   lai = Leaf(2.0)  # Leaf area index
   Ω = 0.5  # Clumping index
-  Tair = 298.0  # Air temperature in Kelvin
   RH = 0.6  # Relative humidity
   α_snow_v = 0.8  # Albedo of snow in visible spectrum
   α_snow_n = 0.7  # Albedo of snow in near-infrared spectrum
@@ -21,23 +17,24 @@ using Test, BEPS
   perc_snow_g = 0.1  # Percentage area of snow on ground
   α_v = Layer3(0.2, 0.25, 0.15)
   α_n = Layer3(0.3, 0.35, 0.25)
+
+  Tair = 298.0  # Air temperature in Kelvin
+  T = Layer3(300.0, 295.0, 290.0)  # Temperature
 
   Rn_Leaf = Leaf(0.0)  # Net radiation for leaf
   Rns_Leaf = Leaf(0.0)  # Net shortwave radiation for leaf
   Rnl_Leaf = Leaf(0.0)  # Net longwave radiation for leaf
   Ra = Radiation()  # Radiation object
 
   # Call the function with the test inputs
-  r_jl = netRadiation_jl(Rs_global, CosZs, 
-    temp_o, temp_u, temp_g, 
+  r_jl = netRadiation_jl(Rs_global, CosZs, T, 
     lai_o, lai_u, lai_os, lai_us, lai, Ω, 
     Tair, RH, 
     α_snow_v, α_snow_n, α_v, α_n,
     perc_snow_o, perc_snow_u, perc_snow_g, 
 
     Rn_Leaf, Rns_Leaf, Rnl_Leaf, Ra)
-  r_c = clang.netRadiation_c(Rs_global, CosZs,
-    temp_o, temp_u, temp_g,
+  r_c = clang.netRadiation_c(Rs_global, CosZs, T,
     lai_o, lai_u, lai_os, lai_us, lai, Ω,
     Tair, RH,
     α_snow_v, α_snow_n, α_v, α_n,