Implement grass allometric equations and update default allometric parameters for grass PFT

biogeochem/FatesAllometryMod.F90

@@ -412,6 +412,9 @@ subroutine bagw_allom(d,ipft,crowndamage, elongf_stem, bagw,dbagwdd)
       case (4) ! 3par_pwr
          call h_allom(d,ipft,h,dhdd)
          call dh2bagw_3pwr(d,h,dhdd,p1,p2,p3,wood_density,c2b,bagw,dbagwdd)
+      case (5) ! 3par_pwr_grass
+         call h_allom(d,ipft,h,dhdd)
+         call dh2bagw_3pwr_grass(d,h,dhdd,p1,p2,p3,c2b,bagw,dbagwdd)
       case DEFAULT
          write(fates_log(),*) 'An undefined AGB allometry was specified: ',allom_amode
          write(fates_log(),*) 'Aborting'
@@ -471,6 +474,9 @@ subroutine blmax_allom(d,ipft,blmax,dblmaxdd)
       case(4) ! dh2blmax_3pwr
          call h_allom(d,ipft,h,dhdd)
          call dh2blmax_3pwr(d,h,dhdd,p1,p2,p3,slatop,dbh_maxh,c2b,blmax,dblmaxdd)
+      case (5) ! dh2blmax_3pwr_grass
+         call h_allom(d,ipft,h,dhdd)
+         call dh2blmax_3pwr_grass(d,h,dhdd,p1,p2,p3,dbh_maxh,c2b,blmax,dblmaxdd)
       case DEFAULT
          write(fates_log(),*) 'An undefined leaf allometry was specified: ', &
               allom_lmode
@@ -530,7 +536,7 @@ subroutine carea_allom(dbh,nplant,site_spread,ipft,crowndamage,c_area,inverse)
           call carea_2pwr(dbh,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max, & 
                crowndamage, c_area, do_inverse)
           capped_allom = .false.
-       case(3)
+       case(3,5)
           dbh_eff = min(dbh,dbh_maxh)
           call carea_2pwr(dbh_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max, &
                crowndamage, c_area, do_inverse)
@@ -1036,6 +1042,25 @@ subroutine bsap_allom(d,ipft,crowndamage,canopy_trim,elongf_stem, sapw_area,bsap
           end if
        end if
 
+    case(2) ! this is a 'sapwood' function specifically for grass PFT that do not produce
+            ! dead woody biomass. So bsap = bagw. Might remove the bsap and bdead for grass
+            ! in the future as there is no need to distinguish the two for grass above- and belowground biomass
+
+       call bagw_allom(d,ipft,crowndamage,elongf_stem,bagw,dbagwdd)
+       call bbgw_allom(d,ipft, elongf_stem,bbgw,dbbgwdd)
+       bsap = bagw + bbgw
+
+       ! This is a grass-only functionnal type, no need to run crown-damage effects
+
+       bsap = elongf_stem * bsap
+       if (present(dbsapdd))then
+          dbsapdd = elongf_stem * (dbagwdd + dbbgwdd)
+       end if
+
+       if(present(dbsapdd))then
+          dbsapdd = dbagwdd + dbbgwdd
+       end if
+
     case DEFAULT
        write(fates_log(),*) 'An undefined sapwood allometry was specified: ', &
             prt_params%allom_smode(ipft)
@@ -1228,7 +1253,7 @@ subroutine bdead_allom(bagw,bbgw,bsap,ipft,bdead,dbagwdd,dbbgwdd,dbsapdd,dbdeadd
             dbdeaddd = dbagwdd/agb_fraction
          end if
 
-      case(2,3,4)
+      case(2,3,4,5)
 
          bdead = bagw + bbgw - bsap
          if(present(dbagwdd) .and. present(dbbgwdd) .and. &
@@ -1631,6 +1656,82 @@ subroutine dh2blmax_3pwr(d,h,dhdd,p1,p2,p3,slatop,dbh_maxh,c2b,blmax,dblmaxdd)
   end subroutine dh2blmax_3pwr
 
 
+
+  ! =========================================================================
+
+  subroutine dh2blmax_3pwr_grass(d,h,dhdd,p1,p2,p3,dbh_maxh,c2b,blmax,dblmaxdd)
+    !------------------------------------------------------------------------
+    !
+    ! This function calculates the maximum leaf biomass using diameter (basal
+    ! diameter for grass) and plant height based on grass leaf allometry developed
+    ! in Gao et al. 2024
+    !
+    !-------------------
+    ! References
+    !-------------------
+    ! Gao X., Koven C., and Kueppers L. 2024. Allometric relationships and trade-offs
+    ! in 11 common Mediterranean-climate grasses. Ecological Applications.
+    ! https://esajournals.onlinelibrary.wiley.com/doi/10.1002/eap.2976
+
+    !------------------
+    ! Input arguments
+    !------------------
+    ! d          -- Basal diameter for grass or any herbaceous plants             [      cm]
+    ! h          -- Plant height                                                  [       m]
+    ! dhdd       -- Height derivative with dbh                                    [    m/cm]
+    ! p1         -- Parameter 1 (log-intercept)                                   [      --]
+    ! p2         -- Parameter 2 (log-slope associated with d)                     [      --]
+    ! p3         -- Parameter 3 (log-slope associated with h)                     [      --]
+    ! dbh_maxh   -- DBH at maximum height                                         [      cm]
+    ! c2b        -- Carbon to biomass multiplier ~ 2                              [  kg/kgC]
+    !
+    !-----------------
+    ! Output arguments
+    !-----------------
+    ! blmax      -- Leaf biomass                                                  [      kgC]
+    ! dblmaxdd   -- Leaf biomass derivative                                       [   kgC/cm]
+    !
+    !-----------------
+    ! Default parameters have been updated for three FATES grass PFTs according to Gao et al. 2024
+    !---------------------------------------------------------------------------------------------
+
+    !-----Arguments
+    real(r8), intent(in)              :: d         ! plant diameter                   [    cm]
+    real(r8), intent(in)              :: h         ! plant height                     [     m]
+    real(r8), intent(in)              :: dhdd      ! height derivative                [  m/cm]
+    real(r8), intent(in)              :: p1        ! log-intercept parameter          [     -]
+    real(r8), intent(in)              :: p2        ! log-slope associated with d      [     -]
+    real(r8), intent(in)              :: p3        ! log-slope associated with h      [     -]
+    real(r8), intent(in)              :: c2b       ! carbon to biomass multiplier     [kg/kgC]
+    real(r8), intent(in)              :: dbh_maxh  ! diameter at maximum height       [    cm]
+    real(r8), intent(out)             :: blmax     ! leaf biomass                     [   kgC]
+    real(r8), intent(out), optional   :: dblmaxdd  ! leaf biomass derivative          [kgC/cm]
+    !----Local variables
+    real(r8)  :: duse
+
+
+    !----Cap diameter
+    duse = min(d, dbh_maxh)
+
+    !----Calculate leaf biomass
+    blmax = p1 * duse**p2 * h**p3 / c2b
+
+    !----Calculate leaf biomass derivative if needed
+
+    if (present(dblmaxdd))then
+       if(d .ge. dbh_maxh)then
+          dblmaxdd = 0._r8
+       else
+          dblmaxdd = blmax * (p2 / duse + p3 * dhdd / h)
+       end if
+    end if
+
+    return
+  end subroutine dh2blmax_3pwr_grass
+
+
+
+
   ! =========================================================================
   ! Diameter to height (D2H) functions
   ! =========================================================================
@@ -2015,6 +2116,69 @@ end subroutine dh2bagw_3pwr
   ! ============================================================================
 
 
+  subroutine dh2bagw_3pwr_grass(d,h,dhdd,p1,p2,p3,c2b,bagw,dbagwdd)
+    !---------------------------------------------------------------------------
+    !
+    ! This function calculates aboveground biomass (excluding leaf biomass) using
+    ! basal diamerer (cm) and plant height (m) as size references, specifically
+    ! for grass or herbaceous plants (can be used for other PFTs if supported by data)
+    !
+    !----------------
+    ! Reference
+    !----------------
+    ! Gao X., Koven C., and Kueppers L. 2024. Allometric relationships and trade-offs in 11
+    ! common Mediterranean-climate grasses. Ecological Applications.
+    ! https://esajournals.onlinelibrary.wiley.com/doi/10.1002/eap.2976
+    !
+    !----------------
+    ! Input arguments
+    !----------------
+    ! d                    -- Basal diameter                      [    cm]
+    ! h                    -- Plant height                        [     m]
+    ! dhdd                 -- Height derivative with diameter     [  m/cm]
+    ! p1                   -- Log-intercept                       [     -]
+    ! p2                   -- Log-slope associated with d         [     -]
+    ! p3                   -- Log-slope associated with h         [     -]
+    ! c2b                  -- Carbon to biomass multiplier        [kg/kgC]
+    !
+    !----------------
+    ! Output variables
+    !----------------
+    ! bagw                 -- Aboveground biomass                 [   kgC]
+    ! dbagwdd              -- Aboveground biomass derivative      [kgC/cm]
+    !
+    !---------------------------------------------------------------------------
+
+
+    !----Arguments
+    real(r8), intent(in)              :: d               ! plant diameter                      [    cm]
+    real(r8), intent(in)              :: h               ! plant height                        [     m]
+    real(r8), intent(in)              :: dhdd            ! height derivative w/ diameter       [  m/cm]
+    real(r8), intent(in)              :: p1              ! log-intercept                       [     -]
+    real(r8), intent(in)              :: p2              ! log-slope associated with d         [     -]
+    real(r8), intent(in)              :: p3              ! log-slope associated with h         [     -]
+    real(r8), intent(in)              :: c2b             ! biomass to carbon multiplier        [kg/kgC]
+    real(r8), intent(out)             :: bagw            ! aboveground biomass excluding leaf  [   kgC]
+    real(r8), intent(out),optional    :: dbagwdd         ! aboveground biomass derivative      [kgC/cm]
+
+    !----Calculate aboveground biomass
+
+    bagw = p1 * (d**p2) * (h**p3) / c2b
+
+    !----Compute the aboveground biomass derivative with basal diameter if needed
+    if (present(dbagwdd)) then
+       dbagwdd = p2 * bagw / d + p3 * bagw * dhdd / h
+    end if
+
+    return
+  end subroutine dh2bagw_3pwr_grass
+
+
+
+  ! ============================================================================
+
+
+
   subroutine d2bagw_2pwr(d,p1,p2,c2b,bagw,dbagwdd)
 
     ! =========================================================================

main/EDPftvarcon.F90

@@ -2183,6 +2183,19 @@ subroutine FatesCheckParams(is_master)
         end if
 
 
+        ! Check to make sure that if a grass sapwood allometry is used, it is not
+        ! a woody plant.
+        if ( ( prt_params%allom_smode(ipft)==2 ) .and. (prt_params%woody(ipft)==itrue) ) then
+           write(fates_log(),*) 'Allometry mode 2 is a mode that is only appropriate'
+           write(fates_log(),*) 'for a grass functional type. Sapwood allometry is set with'
+           write(fates_log(),*) 'fates_allom_smode in the parameter file. Woody versus non woody'
+           write(fates_log(),*) 'plants are set via fates_woody in the parameter file.'
+           write(fates_log(),*) 'Current settings for pft number: ',ipft
+           write(fates_log(),*) 'fates_woody: true'
+           write(fates_log(),*) 'fates_allom_smode: ',prt_params%allom_smode(ipft)
+           write(fates_log(),*) 'Please correct this discrepancy before re-running. Aborting.'
+           call endrun(msg=errMsg(sourcefile, __LINE__))
+        end if
 
         ! Check if fraction of storage to reproduction is between 0-1
         ! ----------------------------------------------------------------------------------