Merge pull request #1132 from mpaiao/mpaiao-pr-cstarve

Exponential-based carbon starvation mortality function

biogeochem/EDMortalityFunctionsMod.F90

@@ -12,7 +12,11 @@ module EDMortalityFunctionsMod
    use FatesCohortMod        , only : fates_cohort_type
    use EDTypesMod            , only : ed_site_type
    use EDParamsMod,            only : maxpft
+   use EDParamsMod           , only : mort_cstarvation_model
    use FatesConstantsMod     , only : itrue,ifalse
+   use FatesConstantsMod     , only : cstarvation_model_lin
+   use FatesConstantsMod     , only : cstarvation_model_exp
+   use FatesConstantsMod     , only : nearzero
    use FatesConstantsMod     , only : ihard_stress_decid
    use FatesConstantsMod     , only : isemi_stress_decid
    use FatesConstantsMod     , only : leaves_off
@@ -190,24 +194,44 @@ subroutine mortality_rates( cohort_in,bc_in, btran_ft, mean_temp,             &
        ! Carbon Starvation induced mortality.
        if ( cohort_in%dbh  >  0._r8 ) then
 
-          ! We compare storage with leaf biomass if plant were fully flushed, otherwise
-          ! mortality would be underestimated for plants that lost all leaves and have no
-          ! storage to flush new ones.
-          ! MLO. Why isn't this comparing with storage allometry (i.e., accounting for
-          !      cushion)?
+          ! Find the current ratio between storage biomass and leaf biomass, which will be
+          ! used to define carbon starvation mortality.  The reference leaf biomass is 
+          ! always for when plants are fully flushed (but accounting for damage and 
+          ! trimming).
           call bleaf(cohort_in%dbh,cohort_in%pft,cohort_in%crowndamage,cohort_in%canopy_trim, &
                1.0_r8, target_leaf_c)
           store_c = cohort_in%prt%GetState(store_organ,carbon12_element)
-
           call storage_fraction_of_target(target_leaf_c, store_c, frac)
-          if( frac .lt. 1._r8) then
-             cmort = max(0.0_r8,EDPftvarcon_inst%mort_scalar_cstarvation(cohort_in%pft) * &
-                  (1.0_r8 - frac))
-          else
+
+          ! Select the carbon starvation mortality model (linear or exponential)s.
+          select case (mort_cstarvation_model)
+          case (cstarvation_model_lin)
+             ! Linear model. Carbon starvation mortality will be zero when fraction of
+             ! storage is greater than or equal to mort_upthresh_cstarvation, and will
+             ! increase to the maximum mortality (mort_scalar_cstarvation) when frac = 0.
+             cmort = EDPftvarcon_inst%mort_scalar_cstarvation(cohort_in%pft) * &
+                max(0.0_r8, (EDPftvarcon_inst%mort_upthresh_cstarvation(cohort_in%pft)-frac) / &
+                            EDPftvarcon_inst%mort_upthresh_cstarvation(cohort_in%pft) )
+
+          case (cstarvation_model_exp)
+             ! Exponential model.  Maximum carbon starvation mortality 
+             ! (mort_scalar_cstarvation) occurs when frac=0. Parameter 
+             ! mort_upthresh_cstarvation controls the the e-folding factor for frac. The
+             ! smaller the mort_upthresh_cstarvation, the faster the mortality will decay.
+             cmort = EDPftvarcon_inst%mort_scalar_cstarvation(cohort_in%pft) * &
+                     exp(- frac / EDPftvarcon_inst%mort_upthresh_cstarvation(cohort_in%pft))
+
+          case default
+              write(fates_log(),*) &
+                 'Invalid carbon starvation model (',mort_cstarvation_model,').'
+              call endrun(msg=errMsg(sourcefile, __LINE__))
+          end select
+
+          ! Make sure the mortality is set to zero when tiny.
+          if (cmort <= nearzero) then
              cmort = 0.0_r8
-          endif
+          end if
 
-
        else
           write(fates_log(),*) 'dbh problem in mortality_rates', &
                cohort_in%dbh,cohort_in%pft,cohort_in%n,cohort_in%canopy_layer

biogeochem/FatesAllometryMod.F90

@@ -616,20 +616,25 @@ end subroutine bleaf
   subroutine storage_fraction_of_target(c_store_target, c_store, frac)
 
     !--------------------------------------------------------------------------------
-    ! returns the storage pool as a fraction of its target (only if it is below its target)
-    ! used in both the carbon starvation mortlaity scheme as well as the optional
-    ! respiration throttling logic
+    !    This subroutine returns the ratio between the storage pool and the target 
+    ! storage.  This subroutine is used both the carbon starvation mortality scheme 
+    ! and the optional respiration throttling. We impose checks so it cannot go negative
+    ! due to truncation errors, but this function can return values greater than 1.
+    ! 
+    !    Fractions exceeding do not impact the default linear carbon starvation model
+    ! (mort_cstarvation_model=2), because mortality becomes zero, but they allow carbon 
+    ! starvation mortality rates to continue decaying when the exponential carbon 
+    ! starvation model is used (mort_cstarvation_model=2).
+    ! 
+    !    Fraction values above 1 do not impact lowstorage_maintresp_reduction either,
+    ! as that routine imposes no reduction once the fraction exceeds 1.
     !--------------------------------------------------------------------------------
 
     real(r8),intent(in)    :: c_store_target  ! target storage carbon [kg]
     real(r8),intent(in)    :: c_store         ! storage carbon [kg]
     real(r8),intent(out)   :: frac
 
-    if( c_store_target > 0._r8 .and. c_store <= c_store_target )then
-       frac = c_store/ c_store_target
-    else
-       frac = 1._r8
-    endif
+    frac = max(0._r8, c_store / max( c_store_target, nearzero) )
 
   end subroutine storage_fraction_of_target
 

main/EDParamsMod.F90

@@ -10,6 +10,7 @@ module EDParamsMod
    use FatesGlobals        , only : fates_log
    use FatesGlobals        , only : endrun => fates_endrun
    use FatesConstantsMod,    only : fates_unset_r8
+   use FatesConstantsMod,    only : cstarvation_model_lin
    use FatesConstantsMod,    only : n_landuse_cats
 
    ! CIME Globals
@@ -50,7 +51,11 @@ module EDParamsMod
                                                       ! 1=non-acclimating, 2=Kumarathunge et al., 2019
 
    integer,protected, public :: radiation_model       ! Switch betrween Norman (1) and Two-stream (2) radiation models
-
+
+   integer,protected, public :: mort_cstarvation_model ! Switch for carbon starvation mortality:
+                                                       ! 1 -- Linear model
+                                                       ! 2 -- Exponential model
+
    real(r8),protected, public :: fates_mortality_disturbance_fraction ! the fraction of canopy mortality that results in disturbance
    real(r8),protected, public :: ED_val_comp_excln                    ! weighting factor for canopy layer exclusion and promotion
    real(r8),protected, public :: ED_val_vai_top_bin_width             ! width in VAI units of uppermost leaf+stem layer scattering element
@@ -140,6 +145,7 @@ module EDParamsMod
    character(len=param_string_length),parameter,public :: name_radiation_model = "fates_rad_model"
    character(len=param_string_length),parameter,public :: ED_name_hydr_htftype_node = "fates_hydro_htftype_node"
    character(len=param_string_length),parameter,public :: ED_name_mort_disturb_frac = "fates_mort_disturb_frac"
+   character(len=param_string_length),parameter,public :: ED_name_mort_cstarvation_model = "fates_mort_cstarvation_model"
    character(len=param_string_length),parameter,public :: ED_name_comp_excln = "fates_comp_excln"
    character(len=param_string_length),parameter,public :: ED_name_vai_top_bin_width = "fates_vai_top_bin_width"
    character(len=param_string_length),parameter,public :: ED_name_vai_width_increase_factor = "fates_vai_width_increase_factor"
@@ -311,6 +317,7 @@ subroutine FatesParamsInit()
     maintresp_leaf_model                  = -9
     radiation_model                       = -9
     fates_mortality_disturbance_fraction  = nan
+    mort_cstarvation_model                = -9
     ED_val_comp_excln                     = nan
     ED_val_vai_top_bin_width              = nan
     ED_val_vai_width_increase_factor      = nan
@@ -423,6 +430,10 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=ED_name_mort_disturb_frac, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
+    ! Temporary until we add parameter to file
+    !call fates_params%RegisterParameter(name=ED_name_mort_cstarvation_model, dimension_shape=dimension_shape_scalar, &
+    !     dimension_names=dim_names_scalar)
+
     call fates_params%RegisterParameter(name=ED_name_comp_excln, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
@@ -637,6 +648,12 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetrieveParameter(name=ED_name_mort_disturb_frac, &
           data=fates_mortality_disturbance_fraction)
 
+    ! Temporary until we add parameter to file 
+    !call fates_params%RetrieveParameter(name=ED_name_mort_cstarvation_model, &
+    !     data=tmpreal)
+    !    mort_cstarvation_model = nint(tmpreal)
+    mort_cstarvation_model = cstarvation_model_lin
+
     call fates_params%RetrieveParameter(name=ED_name_comp_excln, &
          data=ED_val_comp_excln)
 

main/FatesConstantsMod.F90

@@ -141,6 +141,10 @@ module FatesConstantsMod
   integer, parameter, public :: lmrmodel_ryan_1991         = 1
   integer, parameter, public :: lmrmodel_atkin_etal_2017   = 2
 
+  ! integer labels for specifying carbon starvation model
+  integer, parameter, public :: cstarvation_model_lin = 1 ! Linear scaling
+  integer, parameter, public :: cstarvation_model_exp = 2 ! Exponential scaling
+
   ! Error Tolerances
 
   ! Allowable error in carbon allocations, should be applied to estimates