Changes to soil water and coppicing/shooting

src/phenology/shooting.jl

@@ -8,21 +8,48 @@
 
     shoot_max => 1e4 ~ preserve(parameter, u"kg/ha")
 
-    shooting(F, Rf, shoot_max, shoot, WS) => begin
-        (F >= Rf) && (shoot_max >= shoot) && (WS <= shoot_max) 
+    shooting(F, Rf, shoot_max, shoot, WS, coppice_days,shooting_interval) => begin
+        (F >= Rf) && (shoot_max >= shoot) && (WS <= shoot_max) && (coppice_days>0u"d" && coppice_days<shooting_interval)
     end ~ flag
 
+    # Days after coppicing where reshooting can occur
+    shooting_interval => begin
+	365
+    end ~ preserve(parameter,u"d")
+
+    root_shoot_ratio(WR,WS,WF) => begin
+        (WR)/(WS+WF)
+    end ~ track()
+
+    min_RSR:min_root_shoot_ratio => begin
+         .1
+    end ~ preserve(parameter)
+
+    percent_NSC => begin
+	.3
+    end ~ preserve(parameter)
+
+    non_structural_carbon(WR,percent_NSC) => begin
+        percent_NSC*WR
+    end ~ track(u"kg/ha")
+
     ShD(T_air, T_shoot, T_shoot_opt): shooting_degrees => begin
         min(T_air, T_shoot_opt) - T_shoot 
     end ~ track(when=shooting, min=0, u"K")
 
     shoot_rate => 2 ~ preserve(parameter, u"kg/ha/hr/K")
 
     # Maximum shoot growth available for coppicing based on available root drymass.
-    dShoot_max(WR, step) => WR / step ~ track(u"kg/ha/hr")
+    dShoot_max(WR, step,non_structural_carbon) => non_structural_carbon / step ~ track(u"kg/ha/hr")
 
     # Hourly shoot growth rate.
-    dShoot(shoot_rate, ShD) => shoot_rate * ShD ~ track(max=dShoot_max, u"kg/ha/hr")
+    dShoot(shoot_rate, ShD, root_shoot_ratio,non_structural_carbon,shoot,min_RSR) => begin 
+        if(root_shoot_ratio<min_RSR || shoot>=non_structural_carbon) 
+	   0
+        else
+           shoot_rate * ShD 
+        end
+    end ~ track(max=dShoot_max, u"kg/ha/hr")
 
     # Accumulated shoot growth for the season, resets every year.
     shoot(dShoot) ~ accumulate(reset=senescent, u"kg/ha")

src/physiology/gasexchange/gasexchange.jl

@@ -6,9 +6,11 @@ include("irradiance.jl")
 include("stomata.jl")
 include("../../atmosphere/atmosphere.jl")
 
+
 @system GasExchange(Atmosphere, BoundaryLayer, Calendar, Stomata, IntercellularSpace, Irradiance, EnergyBalance, C3) begin
     PPFD: photosynthetic_photon_flux_density ~ track(u"μmol/m^2/s" #= Quanta =#, override)
     LAI: leaf_area_index ~ track(override)
+    drought_factor ~ track(override)
 
     A_net_total(A_net, LAI): net_photosynthesis_total => A_net * LAI ~ track(u"μmol/m^2/s" #= CO2 =#)
     A_gross_total(A_gross, LAI): gross_photosynthesis_total => A_gross * LAI ~ track(u"μmol/m^2/s" #= CO2 =#)

src/physiology/gasexchange/stomata.jl

@@ -14,6 +14,7 @@
 end
 
 @system Stomata(StomataBase) begin
+    drought_factor ~ hold()
     g0 => 0.01 ~ preserve(u"mol/m^2/s/bar" #= H2O =#, parameter)
     g1 => 9.670307198008624 ~ preserve(parameter)
 
@@ -27,9 +28,9 @@ end
         g0 + (g1 * (1-s) * m * (A_net * hs / Cs))
     end ~ track(u"mol/m^2/s/bar" #= H2O =#, min=g0)
 
-    m: transpiration_reduction_factor => begin
+    m(drought_factor): transpiration_reduction_factor => begin
         #TODO: implement soil water module
-        1.0
+        drought_factor
     end ~ track
 
     s: senescence_reduction_factor ~ track(override)

src/physiology/nitrogen.jl

@@ -75,18 +75,18 @@
 
     # Relative drought factor from CROPGRO. Used for N_uptake_conversion_factor.
     "Relative drought factor"
-    drought_factor(ASW, minASW, field_capacity, maxASW) => begin
-        if ASW > field_capacity
-            1 - (ASW - field_capacity) / (maxASW - field_capacity)
+    drought_factor(SW, minSW, field_capacity, soil_saturation,WP) => begin
+        if SW > field_capacity
+            1.0 - (SW - field_capacity) / (soil_saturation - field_capacity)
         else
-            ((ASW - minASW) / (field_capacity - minASW))
+            1 * ((SW - WP) / (field_capacity - WP))
         end
     end ~ track(min=0.1, max=1) 
 
     # Nitrogen uptake conversion factor.
     # How much kg/ha of nitrogen for mg/cm of nitrogen (root)?
     N_uptake_conversion_factor(RLD, drought_factor, soil_depth) => begin
-        RLD * sqrt(drought_factor) * soil_depth
+        RLD * sqrt(drought_factor) * soil_depth 
     end ~ track(u"kg*cm/mg/ha")
 
     "Amount of NO3 that stays in soil"

src/physiology/photosynthesis.jl

@@ -9,12 +9,20 @@ Photosynthesis
     Gas-exchange Model
     =================#
 
+#     "Gas exchange model for sunlit leaves"
+
+#     sunlit_gasexchange(context, PPFD=Q_sun, LAI=LAI_sunlit, w=leaf_width,drought_factor =drought_factor) ~ ::GasExchange
+
+#     "Gas exchange model for shaded leaves"
+#     shaded_gasexchange(context, PPFD=Q_sh, LAI=LAI_shaded, w=leaf_width, drought_factor = drought_factor) ~ ::GasExchange
+
     "Gas exchange model for sunlit leaves"
     sunlit_gasexchange(context, PPFD=Q_sun, LAI=LAI_sunlit, w=leaf_width, s=s) ~ ::GasExchange
 
     "Gas exchange model for shaded leaves"
     shaded_gasexchange(context, PPFD=Q_sh, LAI=LAI_shaded, w=leaf_width, s=s) ~ ::GasExchange
 
+
     #=================
     =================#
 

src/rhizosphere/waterbalance.jl

@@ -8,19 +8,36 @@ Transpiration
     Parameters
     =========#
 
-    "Initial available soil water"
-    iASW => 200 ~ preserve(parameter, u"mm")
+    "Initial soil water"
+    iSW => 200 ~ preserve(parameter, u"mm")
 
-    "Maximum available soil water"
-    maxASW => 200 ~ preserve(parameter, u"mm") # (Field capacity - Wilting point) * depth
-    #maxASW => 200 ~ preserve(parameter, u"mm") # previous
+
+    "Maximum soil water/saturation"
+    soil_saturation => 500 ~ preserve(parameter, u"mm")
+
+    "Minimum soil water"
+    minSW => 175 ~ preserve(parameter, u"mm")
+
+    "Wilting point"
+    WP(saturation,soil_table,soil_class):wilting_point => begin
+        175
+        #soil_table[soil_class].wilting_point
+    end ~ preserve(parameter, u"mm")
+
+    "Irrigation"
+    irrigation => 0 ~ preserve(parameter, u"mm/hr")
+
+#     "Maximum available soil water"
+#     maxASW => 200 ~ preserve(parameter, u"mm") # (Field capacity - Wilting point) * depth
+#     #maxASW => 200 ~ preserve(parameter, u"mm") # previous
 
-    "Minimum available soil water"
-    minASW => 0 ~ preserve(parameter, u"mm") # need to be updated based on VWC by soil types
+#     "Minimum available soil water"
+#     minASW => 0 ~ preserve(parameter, u"mm") # need to be updated based on VWC by soil types
 
     #"Irrigation"
     #irrigation => 0 ~ preserve(parameter, u"mm/hr")
 
+
     "Fraction of excess water pooled"
     pool_fraction => 0 ~ preserve(parameter)
 
@@ -44,8 +61,12 @@ Transpiration
     #     ((Int(soil_class) > 0) ? (11 - 2 * Int(c)) : (SWpower0))
     # end ~ preserve
 
-    fc => 0.5 ~ preserve(parameter)
-    field_capacity(fc, maxASW, minASW) => fc * (maxASW + minASW) ~ preserve(u"mm")
+###
+    field_capacity(soil_saturation) => 327 ~ preserve(u"mm")
+
+#     fc => 0.5 ~ preserve(parameter)
+#     field_capacity(fc, maxASW, minASW) => fc * (maxASW + minASW) ~ preserve(u"mm")
+
 
     "Proportion of rain intercepted"
     interception(LAI, maxInterception, LAImaxInterception) => begin
@@ -58,16 +79,16 @@ Transpiration
     "Canopy transpiration" # Looks like missing soil surface evaporation 
     evapotranspiration(transpiration, rainInterception) => begin
         transpiration + rainInterception
-    end ~ track(u"mm/hr", max=ASWhour)
+    end ~ track(u"mm/hr", max=SWhour)
 
     "Hourly excess soil water"
-    excessSW(ASWhour, maxASWhour, evapotranspiration, irrigation, rain, poolHour) => begin
-        ASWhour + rain + poolHour - evapotranspiration + irrigation - maxASWhour
+    excessSW(SWhour, maxSWhour, evapotranspiration, irrigation, rain, poolHour) => begin
+        SWhour + rain + poolHour - evapotranspiration + irrigation - maxSWhour
     end ~ track(u"mm/hr", min=0u"mm/hr")
 
     "Hourly loss in water pool"
-    lossPool(ASWhour, maxASWhour, evapotranspiration, irrigation, rain, poolHour) => begin
-        maxASWhour - (ASWhour - evapotranspiration + irrigation + rain + poolHour)
+    lossPool(SWhour, maxSWhour, evapotranspiration, irrigation, rain, poolHour) => begin
+        maxSWhour - (SWhour - evapotranspiration + irrigation + rain + poolHour)
     end ~ track(u"mm/hr", min=0u"mm/hr", max=poolHour)
 
     "Hourly gain in water pool"
@@ -83,6 +104,10 @@ Transpiration
         (1 - pool_fraction) * excessSW
     end ~ track(u"mm/hr")
 
+###
+#     dSW(dPool, evapotranspiration#=, irrigation=#, rain) => begin
+#          -dPool - evapotranspiration#= + irrigation=# + rain
+
     "Hourly change in avilable soil water"
     dASW(dPool, evapotranspiration, irrigation, rain) => begin
          -dPool - evapotranspiration + irrigation + rain
@@ -95,13 +120,13 @@ Transpiration
         evapotranspiration / (transpiration + rainInterception)
     end ~ track(when=flag_transpiration, init=1)
 
-    ASW(dASW) ~ accumulate(u"mm", init=iASW, min=minASW, max=maxASW)
+    SW(dSW) ~ accumulate(u"mm", init=iSW, min=minSW, max=soil_saturation)
     pool(dPool) ~ accumulate(u"mm")
     runoff(dRunoff) ~ accumulate(u"mm")
 
     # ASW and MaxASW and pool in mm/hr
-    ASWhour(ASW) => ASW / u"d" ~ track(u"mm/hr")
-    maxASWhour(maxASW) => maxASW / u"d" ~ track(u"mm/hr")
+    SWhour(SW) => SW / u"d" ~ track(u"mm/hr")
+    maxSWhour(soil_saturation) => soil_saturation / u"d" ~ track(u"mm/hr")
     poolHour(pool) => pool / u"d" ~ track(u"mm/hr")
 
     "Irrigation based on profiling VWC for Slit Loam"

src/silviculture/coppicing.jl

@@ -8,13 +8,22 @@
         (time in coppicing_date) && (dormant)
     end ~ flag
 
+    # True after coppicing has occurred for the first time
+    first_coppice(coppice) ~ flag(once)
+
     # Don't have to worry about foliage during dormancy
     coppicing(step, WS, growthStem, deathStem, thinning_WS, dBud) => begin
         (WS / step) - (growthStem - deathStem - thinning_WS - dBud)
     end ~ track(when=coppice, u"kg/ha/hr")
 
+   # Days since last coppiced 
+   coppice_days(coppice,first_coppice) => begin
+        1
+   end ~ accumulate(when=first_coppice, reset=coppice,init=0, u"d")
+
+
     # Coppiced when stem biomass is zero.
-    coppiced(WS, W) => begin
-        WS == 0u"kg/ha" && W != 0u"kg/ha"
+    coppiced(WS, W, first_coppice) => begin
+        WS == 0u"kg/ha" && W != 0u"kg/ha" && first_coppice
     end ~ flag
 end