Merge pull request #1086 from JessicaNeedham/jfn-atkin-lmr-cap

Prevent negative values of Rdark
NGEET · Sep 28, 2023 · 040e23a · 040e23a
2 parents 73f31ab + ddff4c4
commit 040e23a
Show file tree

Hide file tree

Showing 4 changed files with 63 additions and 15 deletions.
diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90
@@ -22,7 +22,7 @@ module FATESPlantRespPhotosynthMod
 
   use FatesGlobals,      only : endrun => fates_endrun
   use FatesGlobals,      only : fates_log
-  use FatesGlobals,      only : FatesWarn,N2S,A2S
+  use FatesGlobals,      only : FatesWarn,N2S,A2S,I2S
   use FatesConstantsMod, only : r8 => fates_r8
   use FatesConstantsMod, only : itrue
   use FatesConstantsMod, only : nearzero
@@ -2140,6 +2140,11 @@ subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, &
    use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm
    use FatesConstantsMod, only : umolC_to_kgC
    use FatesConstantsMod, only : g_per_kg
+   use FatesConstantsMod, only : lmr_b
+   use FatesConstantsMod, only : lmr_c
+   use FatesConstantsMod, only : lmr_TrefC
+   use FatesConstantsMod, only : lmr_r_1
+   use FatesConstantsMod, only : lmr_r_2
    use EDPftvarcon      , only : EDPftvarcon_inst
 
    ! Arguments
@@ -2156,15 +2161,6 @@ subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, &
    real(r8) :: r_t_ref ! acclimated ref respiration rate (umol CO2/m**2/s)
    real(r8) :: lmr25top  ! canopy top leaf maint resp rate at 25C for this pft (umol CO2/m**2/s)
 
-   ! Parameters
-   ! values from Atkin et al., 2017 https://doi.org/10.1007/978-3-319-68703-2_6
-   ! and Heskel et al., 2016 https://doi.org/10.1073/pnas.1520282113
-   real(r8), parameter :: b = 0.1012_r8       ! (degrees C**-1)
-   real(r8), parameter :: c = -0.0005_r8      ! (degrees C**-2)
-   real(r8), parameter :: TrefC = 25._r8      ! (degrees C)
-   real(r8), parameter :: r_1 = 0.2061_r8     ! (umol CO2/m**2/s / (gN/(m2 leaf))) 
-   real(r8), parameter :: r_2 = -0.0402_r8    ! (umol CO2/m**2/s/degree C)
-
    ! parameter values of r_0 as listed in Atkin et al 2017: (umol CO2/m**2/s) 
    ! Broad-leaved trees  1.7560
    ! Needle-leaf trees   1.4995
@@ -2178,9 +2174,15 @@ subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, &
    ! r_0 currently put into the EDPftvarcon_inst%dev_arbitrary_pft
    ! all figs in Atkin et al 2017 stop at zero Celsius so we will assume acclimation is fixed below that
    r_0 = EDPftvarcon_inst%maintresp_leaf_atkin2017_baserate(ft)
-   r_t_ref = nscaler * (r_0 + r_1 * lnc_top + r_2 * max(0._r8, (tgrowth - tfrz) ))
+   r_t_ref = max( 0._r8, nscaler * (r_0 + lmr_r_1 * lnc_top + lmr_r_2 * max(0._r8, (tgrowth - tfrz) )) )
+
+   if (r_t_ref .eq. 0._r8) then
+      warn_msg = 'Rdark is negative at this temperature and is capped at 0. tgrowth (C): '//trim(N2S(tgrowth-tfrz))//' pft: '//trim(I2S(ft))
+      call FatesWarn(warn_msg,index=4)            
+   end if
 
-   lmr = r_t_ref * exp(b * (veg_tempk - tfrz - TrefC) + c * ((veg_tempk-tfrz)**2 - TrefC**2))
+   lmr = r_t_ref * exp(lmr_b * (veg_tempk - tfrz - lmr_TrefC) + lmr_c * &
+        ((veg_tempk-tfrz)**2 - lmr_TrefC**2))
 
 end subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017
 

diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90
@@ -1717,12 +1717,15 @@ subroutine FatesCheckParams(is_master)
      ! -----------------------------------------------------------------------------------
     use FatesConstantsMod  , only : fates_check_param_set
     use FatesConstantsMod  , only : itrue, ifalse
+    use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm
+    use FatesConstantsMod, only : lmr_r_1
+    use FatesConstantsMod, only : lmr_r_2
     use EDParamsMod        , only : logging_mechanical_frac, logging_collateral_frac
     use EDParamsMod        , only : logging_direct_frac,logging_export_frac
     use EDParamsMod        , only : radiation_model
     use FatesInterfaceTypesMod, only : hlm_use_fixed_biogeog,hlm_use_sp, hlm_name
     use FatesInterfaceTypesMod, only : hlm_use_inventory_init
-
+    
      ! Argument
      logical, intent(in) :: is_master    ! Only log if this is the master proc
 
@@ -1737,7 +1740,11 @@ subroutine FatesCheckParams(is_master)
      integer  :: fates_pft  ! used in fixed biogeog mode
 
      real(r8) :: sumarea    ! area of PFTs in nocomp mode.
-
+     real(r8) :: neg_lmr_temp ! temperature at which lmr would got negative 
+     real(r8) :: r_0 ! base respiartion rate, PFT-dependent
+     real(r8) :: lnc_top ! leaf nitrogen content at top of canopy
+
+
      npft = size(EDPftvarcon_inst%freezetol,1)
 
      if(.not.is_master) return
@@ -2031,6 +2038,32 @@ subroutine FatesCheckParams(is_master)
      end do !ipft
 
 
+     ! Check the temperature at which Rdark would become negative for each PFT -
+     ! given their parameters
+     !------------------------------------------------------------------------------------
+     do ipft = 1,npft
+
+        r_0 = EDPftvarcon_inst%maintresp_leaf_atkin2017_baserate(ipft)
+
+        lnc_top = prt_params%nitr_stoich_p1(ipft, prt_params%organ_param_id(leaf_organ))
+
+        ! From LeafLayerMaintenanceRespiration_Atkin_etal_2017
+        ! r_t_ref = nscaler * (r_0 + r_1 * lnc_top + r_2 * max(0._r8, (tgrowth - tfrz) ))
+
+        ! find temperature at which whole term is negative
+        neg_lmr_temp = ( -1._r8 * (  r_0  + lmr_r_1 * lnc_top ) ) / lmr_r_2
+
+        write(fates_log(),*)  'PFT  ',  ipft
+        write(fates_log(),*)  'will have  negative Rdark at ', neg_lmr_temp, 'degrees C' 
+        write(fates_log(),*)  'with these values of slatop, nitrogen stoichiometry and' 
+        write(fates_log(),*)  'maintresp_leaf_atkin2017_baserate.'
+        write(fates_log(),*)  'See LeafLayerMaintenanceRespiration_Atkin_etal_2017 in '
+        write(fates_log(),*)  'FatesPlantRespPhotosynthMod'  
+
+     end do ! ipft
+
+
+
 !!    ! Checks for HYDRO
 !!    if( hlm_use_planthydro == itrue ) then
 !!

diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90
@@ -316,4 +316,17 @@ module FatesConstantsMod
   ! PI
   real(fates_r8), parameter, public :: pi_const = 3.14159265359_fates_r8
 
+  ! Rdark constants from Atkin et al., 2017 https://doi.org/10.1007/978-3-319-68703-2_6
+  ! and Heskel et al., 2016 https://doi.org/10.1073/pnas.1520282113
+  real(fates_r8), parameter, public :: lmr_b = 0.1012_fates_r8       ! (degrees C**-1)
+
+  real(fates_r8), parameter, public :: lmr_c = -0.0005_fates_r8      ! (degrees C**-2)
+
+  real(fates_r8), parameter, public :: lmr_TrefC = 25._fates_r8      ! (degrees C)
+
+  real(fates_r8), parameter, public :: lmr_r_1 = 0.2061_fates_r8     ! (umol CO2/m**2/s / (gN/(m2 leaf))) 
+
+  real(fates_r8), parameter, public :: lmr_r_2 = -0.0402_fates_r8    ! (umol CO2/m**2/s/degree C)
+
+
 end module FatesConstantsMod
diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90
@@ -1940,8 +1940,8 @@ subroutine FatesReportParameters(masterproc)
 
       call FatesReportPFTParams(masterproc)
       call FatesReportParams(masterproc)
-      call FatesCheckParams(masterproc)    ! Check general fates parameters
       call PRTDerivedParams()              ! Update PARTEH derived constants
+      call FatesCheckParams(masterproc)    ! Check general fates parameters
       call PRTCheckParams(masterproc)      ! Check PARTEH parameters
       call SpitFireCheckParams(masterproc)