minor

.gitignore

@@ -2,3 +2,4 @@ data/inpara
 data
 Manifest.toml
 *debug*
+docs

src/Param/get_pftpar.jl

@@ -25,7 +25,7 @@ end
   + 1  `beta`: the coefficient for calculation of interceptions
   + 2  `D50` : the depth above which 50% of the root mas is located, mm
   + 3  `c`   : the shape parameters of logistic dose-response root distribution model
-  + 4  `Zr`  : root depth (m)
+  + 4  `Hc`  : root depth (m)
 
 Standard numbering rule for different PFTs based on MCD12, IGBP
 - `0` : 'water'
@@ -49,7 +49,7 @@ Standard numbering rule for different PFTs based on MCD12, IGBP
 function get_pftpar(LC::Int)
   # The PFTpar look up table
   PLUT = [
-    # :β, :D50, :D95, :Zr
+    # :β, :D50, :D95, :Hc
     0.06 1771 3998 20;   # 1  Evergreen Needleleaf Forest   ------ ENF
     0.02 2187 4316 18;   # 2  Evergreen Broadleaf  Forest   ------ EBF
     0.06 1668 3936 18;   # 3  Deciduous Needleleaf Forest   ------ DNF
@@ -76,5 +76,5 @@ function get_pftpar(LC::Int)
   else
     pftpar = PLUT[Typenum, :]
   end
-  @LArray pftpar (:β, :D50, :D95, :Zr)
+  @LArray pftpar (:β, :D50, :D95, :Hc)
 end

src/Soil/swb_case1.jl

@@ -41,7 +41,7 @@ function swb_case1(wa, IWS, pEc, pEs, s_tem, s_vod, soilpar, pftpar, wet, zm, zg
   Tr_p1_g = Tr_p1 * ((zm[1] - d1) * θ_sat) / (d1 * wa1_unsat + (zm[1] - d1) * θ_sat)
 
   # Moisture constraints
-  f_sm1, f_sm_s1 = swc_stress(wa1, soilpar, pEc, pftpar)
+  f_sm1, f_sm_s1 = swc_stress(wa1, pEc, soilpar, pftpar)
 
   # Actual transpiration
   Tr1_u = clamp( f_sm1 * s_vod * s_tem * Tr_p1_u, 0, d1 * (wa1_unsat - θ_wp))

src/Soil/swb_case2.jl

@@ -55,8 +55,8 @@ function swb_case2(wa, IWS, pEc, pEs, s_tem, s_vod, soilpar, pftpar, wet, zm, zg
   Tr_p2_g = Tr_p2 * ((zm[2] - d2) * θ_sat) / (d2 * wa2_unsat + (zm[2] - d2) * θ_sat)
 
   # Moisture constraints
-  f_sm1, f_sm_s1 = swc_stress(wa1, soilpar, pEc, pftpar)
-  f_sm2, _ = swc_stress(wa2_unsat, soilpar, pEc, pftpar)
+  f_sm1, f_sm_s1 = swc_stress(wa1, pEc, soilpar, pftpar)
+  f_sm2, _ = swc_stress(wa2_unsat, pEc, soilpar, pftpar)
 
   # Actual transpiration
   Tr1 = f_sm1 * s_vod * s_tem * Tr_p1

src/Soil/swb_case3.jl

@@ -64,9 +64,9 @@ function swb_case3(wa, IWS, pEc, pEs, s_tem, s_vod, soilpar, pftpar, wet, zm, zg
   Tr_p3_u = Tr_p3 * (d3 * wa3_unsat) / (d3 * wa3_unsat + (zm[3] - d3) * θ_sat)
   Tr_p3_g = Tr_p3 * ((zm[3] - d3) * θ_sat) / (d3 * wa3_unsat + (zm[3] - d3) * θ_sat)
 
-  f_sm1, f_sm_s1 = swc_stress(wa1, soilpar, pEc, pftpar)
-  f_sm2, _ = swc_stress(wa2, soilpar, pEc, pftpar)
-  f_sm3, _ = swc_stress(wa3_unsat, soilpar, pEc, pftpar)
+  f_sm1, f_sm_s1 = swc_stress(wa1, pEc, soilpar, pftpar)
+  f_sm2, _ = swc_stress(wa2, pEc, soilpar, pftpar)
+  f_sm3, _ = swc_stress(wa3_unsat, pEc, soilpar, pftpar)
 
   Tr1 = f_sm1 * s_vod * s_tem * Tr_p1
   Tr2 = f_sm2 * s_vod * s_tem * Tr_p2

src/Soil/swb_case4.jl

@@ -74,9 +74,9 @@ function swb_case4(wa, IWS, pEc, pEs, s_tem, s_vod, soilpar, pftpar, wet, zm, zg
   Tr_p1, Tr_p2, Tr_p3 = pTr_partition(pEc, wa1, wa2, wa3, soilpar, pftpar, wet, zm)
 
   # Calculate the moisture constrains for plant and soil in unsaturated zone
-  f_sm1, f_sm_s1 = swc_stress(wa1, soilpar, pEc, pftpar)
-  f_sm2, _ = swc_stress(wa2, soilpar, pEc, pftpar)
-  f_sm3, _ = swc_stress(wa3, soilpar, pEc, pftpar)
+  f_sm1, f_sm_s1 = swc_stress(wa1, pEc, soilpar, pftpar)
+  f_sm2, _ = swc_stress(wa2, pEc, soilpar, pftpar)
+  f_sm3, _ = swc_stress(wa3, pEc, soilpar, pftpar)
 
   # Actual transpiration
   Tr1 = f_sm1 * s_vod * s_tem * Tr_p1

src/Soil/swc_stress.jl

@@ -1,62 +1,63 @@
 # Soil moisture Constrains
-function swc_stress(wa, soilpar, pET, pftpar)
-  # INPUT:
-  # wa      : The antecedent soil water content expressed
-  #           as a function of the WHC in that layer
-  # soilpar : Soil parameters according to Soil type
-  # OUTPUT:
-  # S_plant : Soil moisture stress to plant transpiration
-  # S_soil  : Soil moisture stress to soil evaporation
-  # --------
-  # Stress function for plant transpiration and soil evaporation:
-  #                  (θ_c-θ_wp)
-  # wc =  ------------------------------    (about 0.4; Choudhury & Digirolamo, 1998)
-  #                  (θ_fc-θ_wp)
-  # where θ_c : the critical soil water content at which plant stress start
-  # Stree Function (Martens et al., 2017)
-  #                      θ_c-θ
-  # S_plant   =   1- (-------------------)^2   =   1-(1-w/wc)^2
-  #                     θ_c-θ_wp
-  #
-  #                    θ_c-θ
-  # S_soil    =   1- -----------------    =   1-(1-w/wc)=w/wc
-  #                    θ_c-θ_r
-  # -------------------------------------------------------------------------
-  (; θ_fc, θ_wp) = soilpar
-  (; Zr) = pftpar # root depth
+"""
+    swc_stress(wa::T, pET::T, soilpar, pftpar) where {T<:Real}
+
+# INPUT
+- `wa`      : The antecedent soil water content expressed as a function of the WHC in that layer
+- `soilpar` : Soil parameters according to Soil type
 
-  k = Zr^0.5
-  # scale [1, 25] to [1, 5]
-  k = 4 * ((k - 0.7) / 4.3) + 1
+# OUTPUT
+- `S_plant` : Soil moisture stress to plant transpiration
+- `S_soil`  : Soil moisture stress to soil evaporation
+ 
+--------
+Stress function for plant transpiration and soil evaporation:
+                (θ_c-θ_wp)
+wc =  ------------------------------    (about 0.4; Choudhury & Digirolamo, 1998)
+                (θ_fc-θ_wp)
+where θ_c : the critical soil water content at which plant stress start
+Stree Function (Martens et al., 2017)
+                    θ_c-θ
+S_plant   =   1- (-------------------)^2   =   1-(1-w/wc)^2
+                    θ_c-θ_wp
 
-  a = 0.1
-  p = 1 / (1 + pET) - a * (1 / (1 + Zr))
+                  θ_c-θ
+S_soil    =   1- -----------------    =   1-(1-w/wc)=w/wc
+                  θ_c-θ_r
+"""
+function swc_stress(wa::T, pET::T, soilpar, pftpar) where {T<:Real}
+  (; θ_fc, θ_wp) = soilpar
+  (; Hc) = pftpar # canopy height, Zhang 2022
 
+  k = Hc^0.5
+  k = 4 * ((k - 0.7) / 4.3) + 1 # scale [1, 25] to [1, 5]
+
+  b = 0.1
+  p = 1 / (1 + pET) - b / (1 + Hc) # Zhang 2022, Eq. 9
   θ_wpCH = θ_wp / k
 
   # critical soil moisture for different PFTs
   θ_c = (1 - p) * (θ_fc - θ_wpCH) + θ_wpCH
   θ_c = clamp(θ_c, θ_wpCH, θ_fc)
 
   if wa <= θ_wpCH
-    f_sm = 0
+    f_sm = 0.0
   elseif wa >= θ_c
-    f_sm = 1
+    f_sm = 1.0
   else
-    f_sm = 1 - ((θ_c - wa) / (θ_c - θ_wpCH))^k
+    f_sm = 1.0 - ((θ_c - wa) / (θ_c - θ_wpCH))^k
   end
 
   # constraint for soil evaporation
   θ_wp_soil = 0
   if wa <= θ_wp_soil
-    f_sm_s = 0
+    f_sm_s = 0.0
   elseif wa >= θ_fc
-    f_sm_s = 1
+    f_sm_s = 1.0
   else
     # f_sm_s = ((wa - θ_wp) / (θ_fc - θ_wp))^1
     f_sm_s = (wa - θ_wp_soil) / (θ_fc - θ_wp_soil)  # for soil evaporation only
   end
-
   return f_sm, f_sm_s
 end