Merge pull request #1149 from JessicaNeedham/jfn-leafmr-vert-scaling

Alternative vertical scaling of leaf maintenance respiration

biogeochem/EDPhysiologyMod.F90

@@ -45,7 +45,7 @@ module EDPhysiologyMod
   use FatesAllometryMod   , only : tree_lai
   use FatesAllometryMod   , only : tree_sai
   use FatesAllometryMod   , only : leafc_from_treelai
-  use FatesAllometryMod   , only : decay_coeff_kn
+  use FatesAllometryMod   , only : decay_coeff_vcmax
   use FatesLitterMod      , only : litter_type
   use EDTypesMod          , only : site_massbal_type
   use EDTypesMod          , only : numlevsoil_max
@@ -765,7 +765,10 @@ subroutine trim_canopy( currentSite )
 
                 ! Calculate sla_levleaf following the sla profile with overlying leaf area
                 ! Scale for leaf nitrogen profile
-                kn = decay_coeff_kn(ipft,currentCohort%vcmax25top)
+                kn = decay_coeff_vcmax(currentCohort%vcmax25top, &
+                     prt_params%leafn_vert_scaler_coeff1(ipft), &
+                     prt_params%leafn_vert_scaler_coeff2(ipft))
+
                 ! Nscaler value at leaf level z
                 nscaler_levleaf = exp(-kn * cumulative_lai)
                 ! Sla value at leaf level z after nitrogen profile scaling (m2/gC)

biogeochem/FatesAllometryMod.F90

@@ -118,7 +118,7 @@ module FatesAllometryMod
   public :: bdead_allom   ! Generic bdead wrapper
   public :: carea_allom   ! Generic crown area wrapper
   public :: bstore_allom  ! Generic maximum storage carbon wrapper
-  public :: decay_coeff_kn
+  public :: decay_coeff_vcmax ! vertical canopy decay rate, scaled on vcmax
   public :: ForceDBH      ! Method to set DBH to sync with structure
                           ! or fineroot biomass
   public :: CheckIntegratedAllometries
@@ -710,8 +710,10 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25
        end if
 
        ! Coefficient for exponential decay of 1/sla with canopy depth:
-       kn = decay_coeff_kn(pft,vcmax25top)
-
+       kn = decay_coeff_vcmax(vcmax25top, &
+                              prt_params%leafn_vert_scaler_coeff1(pft), &
+                              prt_params%leafn_vert_scaler_coeff2(pft))
+
        ! take PFT-level maximum SLA value, even if under a thick canopy (which has units of m2/gC),
        ! and put into units of m2/kgC
        sla_max = g_per_kg*prt_params%slamax(pft)
@@ -904,8 +906,11 @@ real(r8) function leafc_from_treelai( treelai, pft, c_area, nplant, cl, vcmax25t
     ! convert PFT-level canopy top and maximum SLA values and convert from m2/gC to m2/kgC
     slat = g_per_kg * prt_params%slatop(pft)
     sla_max = g_per_kg * prt_params%slamax(pft)
+
     ! Coefficient for exponential decay of 1/sla with canopy depth:
-     kn = decay_coeff_kn(pft,vcmax25top)
+    kn = decay_coeff_vcmax(vcmax25top, &
+                           prt_params%leafn_vert_scaler_coeff1(pft), &
+                           prt_params%leafn_vert_scaler_coeff2(pft))
 
     if(treelai > 0.0_r8)then 
        ! Leafc_per_unitarea at which sla_max is reached due to exponential sla profile in canopy:
@@ -2743,21 +2748,24 @@ end subroutine jackson_beta_root_profile
   ! =====================================================================================
 
 
-  real(r8) function decay_coeff_kn(pft,vcmax25top)
+  real(r8) function decay_coeff_vcmax(vcmax25top,slope_param,intercept_param)
 
     ! ---------------------------------------------------------------------------------
     ! This function estimates the decay coefficient used to estimate vertical
     ! attenuation of properties in the canopy.
     !
     ! Decay coefficient (kn) is a function of vcmax25top for each pft.
     !
-    ! Currently, this decay is applied to vcmax attenuation, and SLA (optionally)
+    ! Currently, this decay is applied to vcmax attenuation, SLA (optionally)
+    ! and leaf respiration (optionally w/ Atkin)
     !
     ! ---------------------------------------------------------------------------------
 
     !ARGUMENTS
-    integer, intent(in) :: pft
+
     real(r8),intent(in) :: vcmax25top
+    real(r8),intent(in) :: slope_param      ! multiplies vcmax25top
+    real(r8),intent(in) :: intercept_param  ! adds to vcmax25top
 
 
     !LOCAL VARIABLES
@@ -2766,12 +2774,14 @@ real(r8) function decay_coeff_kn(pft,vcmax25top)
     ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used
     ! kn = 0.11. Here, we derive kn from vcmax25 as in Lloyd et al 
     ! (2010) Biogeosciences, 7, 1833-1859
+    ! This function is also used to vertically scale leaf maintenance
+    ! respiration.
 
-    decay_coeff_kn = exp(0.00963_r8 * vcmax25top - 2.43_r8)
+    decay_coeff_vcmax = exp(slope_param * vcmax25top - intercept_param)
 
     return
-  end function decay_coeff_kn
-
+  end function decay_coeff_vcmax
+  
   ! =====================================================================================
 
   subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d, h, bdead, bl )

biogeophys/FatesPlantRespPhotosynthMod.F90

@@ -72,7 +72,8 @@ module FATESPlantRespPhotosynthMod
   use FatesRadiationMemMod, only : ipar
   use FatesTwoStreamUtilsMod, only : FatesGetCohortAbsRad
   use FatesAllometryMod     , only : VegAreaLayer
-
+  use FatesAllometryMod, only : decay_coeff_vcmax
+
   ! CIME Globals
   use shr_log_mod , only      : errMsg => shr_log_errMsg
 
@@ -146,7 +147,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
     use FatesAllometryMod, only : bleaf, bstore_allom
     use FatesAllometryMod, only : storage_fraction_of_target
     use FatesAllometryMod, only : set_root_fraction
-    use FatesAllometryMod, only : decay_coeff_kn
+
 
     use DamageMainMod, only : GetCrownReduction
 
@@ -544,7 +545,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                                  ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al 
                                  ! (2010) Biogeosciences, 7, 1833-1859
 
-                                 kn = decay_coeff_kn(ft,currentCohort%vcmax25top)
+                                 kn = decay_coeff_vcmax(currentCohort%vcmax25top, &
+                                                        prt_params%leafn_vert_scaler_coeff1(ft), &
+                                                        prt_params%leafn_vert_scaler_coeff2(ft))
 
                                  ! Scale for leaf nitrogen profile
                                  nscaler = exp(-kn * cumulative_lai)
@@ -593,7 +596,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                                  case (lmrmodel_atkin_etal_2017)
 
                                     call LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, &  ! in
-                                         nscaler,                            &  ! in
+                                         cumulative_lai,                     &  ! in
+                                         currentCohort%vcmax25top,           &  ! in
                                          ft,                                 &  ! in
                                          bc_in(s)%t_veg_pa(ifp),             &  ! in
                                          currentPatch%tveg_lpa%GetMean(),    &  ! in
@@ -2184,10 +2188,11 @@ end subroutine LeafLayerMaintenanceRespiration_Ryan_1991
   ! ====================================================================================   
 
   subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, &
-       nscaler,   &
-       ft,        &
-       veg_tempk, &
-       tgrowth,   &
+       cumulative_lai, &
+       vcmax25top,     &
+       ft,             &
+       veg_tempk,      &
+       tgrowth,        &
        lmr)
 
     use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm
@@ -2201,33 +2206,49 @@ subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, &
     use EDPftvarcon      , only : EDPftvarcon_inst
 
     ! Arguments
-    real(r8), intent(in)  :: lnc_top      ! Leaf nitrogen content per unit area at canopy top [gN/m2]
-    integer,  intent(in)  :: ft           ! (plant) Functional Type Index
-    real(r8), intent(in)  :: nscaler      ! Scale for leaf nitrogen profile
-    real(r8), intent(in)  :: veg_tempk    ! vegetation temperature  (degrees K)
-    real(r8), intent(in)  :: tgrowth      ! lagged vegetation temperature averaged over acclimation timescale (degrees K)
-    real(r8), intent(out) :: lmr          ! Leaf Maintenance Respiration  (umol CO2/m**2/s)
-
+    real(r8), intent(in)  :: lnc_top          ! Leaf nitrogen content per unit area at canopy top [gN/m2]
+    integer,  intent(in)  :: ft               ! (plant) Functional Type Index
+    real(r8), intent(in)  :: vcmax25top       ! top of canopy vcmax
+    real(r8), intent(in)  :: cumulative_lai   ! cumulative lai above the current leaf layer
+    real(r8), intent(in)  :: veg_tempk        ! vegetation temperature  (degrees K)
+    real(r8), intent(in)  :: tgrowth          ! lagged vegetation temperature averaged over acclimation timescale (degrees K)
+    real(r8), intent(out) :: lmr              ! Leaf Maintenance Respiration  (umol CO2/m**2/s)
+
     ! Locals
     real(r8) :: lmr25   ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s)
     real(r8) :: r_0     ! base respiration rate, PFT-dependent (umol CO2/m**2/s)
     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)
 
+    real(r8) :: rdark_scaler ! negative exponential scaling of rdark
+    real(r8) :: kn           ! decay coefficient
+
     ! 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
     ! Shrubs              2.0749
     ! C3 herbs/grasses    2.1956
     ! In the absence of better information, we use the same value for C4 grasses as C3 grasses.
 
-    ! note that this code uses the relationship between leaf N and respiration from Atkin et al 
-    ! for the top of the canopy, but then assumes proportionality with N through the canopy.
-
     ! 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 = max( 0._r8, nscaler * (r_0 + lmr_r_1 * lnc_top + lmr_r_2 * max(0._r8, (tgrowth - tfrz) )) )
+
+    ! This code uses the relationship between leaf N and respiration from Atkin et al 
+    ! for the top of the canopy, but then scales through the canopy based on a rdark_scaler.
+    ! To assume proportionality with N through the canopy following Lloyd et al. 2010, use the
+    ! default parameter value of 2.43, which results in the scaling of photosynthesis and respiration
+    ! being proportional through the canopy. To have a steeper decrease in respiration than photosynthesis
+    ! this number can be smaller. There is some observational evidence for this being the case
+    ! in Lamour et al. 2023. 
+
+    kn = decay_coeff_vcmax(vcmax25top, &
+                           EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff1(ft), &
+                           EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff2(ft))
+
+    rdark_scaler = exp(-kn * cumulative_lai)
+
+    r_t_ref = max(0._r8, rdark_scaler * (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))

main/EDPftvarcon.F90

@@ -92,8 +92,7 @@ module EDPftvarcon
 
      real(r8), allocatable :: maintresp_leaf_ryan1991_baserate(:)  ! leaf maintenance respiration per Ryan et al 1991
 
-     real(r8), allocatable :: leafn_vert_scaler_coeff1(:)  ! Coefficient one for decrease of leaf N through the canopy
-     real(r8), allocatable :: leafn_vert_scaler_coeff2(:)  ! Coefficient two for decrease of leaf N through the canopy 
+
 
      real(r8), allocatable :: maintresp_leaf_vert_scaler_coeff1(:) ! leaf maintenance respiration decrease through the canopy param 1
                                                                    ! only with Atkin et al. 2017 respiration model
@@ -478,13 +477,7 @@ subroutine Register_PFT(this, fates_params)
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
         dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
-    name = 'fates_leafn_vert_scaler_coeff1'
-    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
-        dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
-    name = 'fates_leafn_vert_scaler_coeff2'
-    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
-         dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
     name = 'fates_maintresp_leaf_vert_scaler_coeff1'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
@@ -935,14 +928,6 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=this%maintresp_leaf_ryan1991_baserate)
 
-    name = 'fates_leafn_vert_scaler_coeff1'
-    call fates_params%RetrieveParameterAllocate(name=name, &
-         data=this%leafn_vert_scaler_coeff1)
-
-    name = 'fates_leafn_vert_scaler_coeff2'
-    call fates_params%RetrieveParameterAllocate(name=name, &
-         data=this%leafn_vert_scaler_coeff2)
-
     name = 'fates_maintresp_leaf_vert_scaler_coeff1'
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=this%maintresp_leaf_vert_scaler_coeff1)
@@ -1776,8 +1761,6 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmt0) 'hydro_vg_alpha_node  = ',EDPftvarcon_inst%hydr_vg_alpha_node
         write(fates_log(),fmt0) 'hydro_vg_m_node  = ',EDPftvarcon_inst%hydr_vg_m_node
         write(fates_log(),fmt0) 'hydro_vg_n_node  = ',EDPftvarcon_inst%hydr_vg_n_node
-        write(fates_log(),fmt0) 'leafn_vert_scaler_coeff1 = ',EDPftvarcon_inst%leafn_vert_scaler_coeff1
-        write(fates_log(),fmt0) 'leafn_vert_scaler_coeff2 = ',EDPftvarcon_inst%leafn_vert_scaler_coeff2
         write(fates_log(),fmt0) 'maintresp_leaf_vert_scaler_coeff1 = ',EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff1
         write(fates_log(),fmt0) 'maintresp_leaf_vert_scaler_coeff2 = ',EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff2
         write(fates_log(),*) '-------------------------------------------------'

parteh/PRTParametersMod.F90

@@ -60,6 +60,11 @@ module PRTParametersMod
      real(r8), allocatable :: turnover_nitr_retrans(:,:) ! nitrogen re-translocation fraction (pft x organ)
      real(r8), allocatable :: turnover_phos_retrans(:,:) ! phosphorus re-translocation fraction (pft x organ)
                                                          ! Parameters dimensioned by PFT and leaf age
+
+     ! These vertical n profile scalers affect crown allometry and sla, thus they
+     ! are here in the PRT module
+     real(r8), allocatable :: leafn_vert_scaler_coeff1(:)  ! Coefficient one for decrease of leaf N through the canopy
+     real(r8), allocatable :: leafn_vert_scaler_coeff2(:)  ! Coefficient two for decrease of leaf N through the canopy 
 
      real(r8), allocatable :: grperc(:)                  ! Growth respiration per unit Carbon gained
                                                          ! One value for whole plant

parteh/PRTParamsFATESMod.F90

@@ -237,6 +237,14 @@ subroutine PRTRegisterPFT(fates_params)
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
+    name = 'fates_leafn_vert_scaler_coeff1'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
+    name = 'fates_leafn_vert_scaler_coeff2'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
     name = 'fates_recruit_seed_alloc_mature'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
@@ -540,6 +548,14 @@ subroutine PRTReceivePFT(fates_params)
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=prt_params%root_long)
 
+    name = 'fates_leafn_vert_scaler_coeff1'
+    call fates_params%RetrieveParameterAllocate(name=name, &
+         data=prt_params%leafn_vert_scaler_coeff1)
+
+    name = 'fates_leafn_vert_scaler_coeff2'
+    call fates_params%RetrieveParameterAllocate(name=name, &
+         data=prt_params%leafn_vert_scaler_coeff2)
+
     name = 'fates_recruit_seed_alloc_mature'
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=prt_params%seed_alloc_mature)
@@ -1040,6 +1056,8 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmti) 'organ_id = ',prt_params%organ_id
         write(fates_log(),fmt0) 'nitr_store_ratio = ',prt_params%nitr_store_ratio
         write(fates_log(),fmt0) 'phos_store_ratio = ',prt_params%phos_store_ratio
+        write(fates_log(),fmt0) 'leafn_vert_scaler_coeff1 = ',prt_params%leafn_vert_scaler_coeff1
+        write(fates_log(),fmt0) 'leafn_vert_scaler_coeff2 = ',prt_params%leafn_vert_scaler_coeff2
         write(fates_log(),*) '-------------------------------------------------'
 
      end if