Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 15 commits into from
Aug 14, 2024
Merged
Show file tree
Hide file tree
Changes from 11 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
178 changes: 176 additions & 2 deletions biogeochem/FatesAllometryMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -1036,6 +1042,35 @@ 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

rgknox marked this conversation as resolved.
Show resolved Hide resolved
! replicate the crown damage code
! Do we really need this for grass? I would think this can be helpful for
! grazing in the future. --XLG
if(crowndamage > 1)then
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

tagging @JessicaNeedham

Copy link
Contributor

@JessicaNeedham JessicaNeedham Aug 12, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agree @XiulinGao that we should think ahead about possible grazing functionality. But as written the crown damage module doesn't change plant height, only crown area and the biomass of aboveground tissues, so I'm not sure how applicable it is to grasses. Maybe we could skip this for now, but think about making the damage code more general so it could be used in the future?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Per discussion in the fates software meeting today @rgknox will make a PR to this branch to remove this and include a graceful failure.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see. In that case then I think this is not that helpful for grass grazing. I agree that we can remove this for grasses for now.

call GetCrownReduction(crowndamage, crown_reduction)
bsap = elongf_stem * (bsap - (bsap * agb_frac * branch_frac * crown_reduction))
if(present(dbsapdd))then
dbsapdd = elongf_stem * &
(dbagwdd + dbbgwdd - ((dbagwdd + dbbgwdd) * agb_frac * branch_frac * crown_reduction))
end if
else
bsap = elongf_stem * bsap
if (present(dbsapdd))then
dbsapdd = elongf_stem * (dbagwdd + dbbgwdd)
end if
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)
Expand Down Expand Up @@ -1228,7 +1263,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. &
Expand Down Expand Up @@ -1631,6 +1666,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
! =========================================================================
Expand Down Expand Up @@ -2015,6 +2126,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)

! =========================================================================
Expand Down
56 changes: 28 additions & 28 deletions parameter_files/fates_params_default.cdl
Original file line number Diff line number Diff line change
Expand Up @@ -983,55 +983,55 @@ data:
0.8, 0.8, 0.8, 0.8 ;

fates_allom_agb1 = 0.0673, 0.1364012, 0.0393057, 0.2653695, 0.0673,
0.0728698, 0.06896, 0.06896, 0.06896, 0.01, 0.01, 0.01 ;
0.0728698, 0.06896, 0.06896, 0.06896, 0.001, 0.001, 0.003 ;

fates_allom_agb2 = 0.976, 0.9449041, 1.087335, 0.8321321, 0.976, 1.0373211,
0.572, 0.572, 0.572, 0.572, 0.572, 0.572 ;
0.572, 0.572, 0.572, 1.6592, 1.6592, 1.3456 ;

fates_allom_agb3 = 1.94, 1.94, 1.94, 1.94, 1.94, 1.94, 1.94, 1.94, 1.94,
1.94, 1.94, 1.94 ;
1.248, 1.248, 1.869 ;

fates_allom_agb4 = 0.931, 0.931, 0.931, 0.931, 0.931, 0.931, 0.931, 0.931,
0.931, 0.931, 0.931, 0.931 ;
0.931, -999.9, -999.9, -999.9 ;

fates_allom_agb_frac = 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6,
0.6, 0.6 ;
fates_allom_agb_frac = 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 1,
1, 1 ;

fates_allom_amode = 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1 ;
fates_allom_amode = 3, 3, 3, 3, 3, 3, 1, 1, 1, 5, 5, 5 ;

fates_allom_blca_expnt_diff = -0.12, -0.34, -0.32, -0.22, -0.12, -0.35, 0,
0, 0, 0, 0, 0 ;
0, 0, -0.487, -0.487, -0.259 ;

fates_allom_cmode = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;

fates_allom_d2bl1 = 0.04, 0.07, 0.07, 0.01, 0.04, 0.07, 0.07, 0.07, 0.07,
0.07, 0.07, 0.07 ;
0.0004, 0.0004, 0.0012 ;

fates_allom_d2bl2 = 1.6019679, 1.5234373, 1.3051237, 1.9621397, 1.6019679,
1.3998939, 1.3, 1.3, 1.3, 1.3, 1.3, 1.3 ;
1.3998939, 1.3, 1.3, 1.3, 1.7092, 1.7092, 1.5879 ;

fates_allom_d2bl3 = 0.55, 0.55, 0.55, 0.55, 0.55, 0.55, 0.55, 0.55, 0.55,
0.55, 0.55, 0.55 ;
0.3417, 0.3417, 0.9948 ;

fates_allom_d2ca_coefficient_max = 0.2715891, 0.3693718, 1.0787259,
0.0579297, 0.2715891, 1.1553612, 0.6568464, 0.6568464, 0.6568464,
0.6568464, 0.6568464, 0.6568464 ;
0.0408, 0.0408, 0.0862 ;

fates_allom_d2ca_coefficient_min = 0.2715891, 0.3693718, 1.0787259,
0.0579297, 0.2715891, 1.1553612, 0.6568464, 0.6568464, 0.6568464,
0.6568464, 0.6568464, 0.6568464 ;
0.6568464, 0.0408, 0.0862 ;
rgknox marked this conversation as resolved.
Show resolved Hide resolved

fates_allom_d2h1 = 78.4087704, 306.842667, 106.8745821, 104.3586841,
78.4087704, 31.4557047, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64 ;
78.4087704, 31.4557047, 0.64, 0.64, 0.64, 0.1812, 0.1812, 0.3353 ;

fates_allom_d2h2 = 0.8124383, 0.752377, 0.9471302, 1.1146973, 0.8124383,
0.9734088, 0.37, 0.37, 0.37, 0.37, 0.37, 0.37 ;
0.9734088, 0.37, 0.37, 0.37, 0.6384, 0.6384, 0.4235 ;

fates_allom_d2h3 = 47.6666164, 196.6865691, 93.9790461, 160.6835089,
47.6666164, 16.5928174, -999.9, -999.9, -999.9, -999.9, -999.9, -999.9 ;

fates_allom_dbh_maxheight = 1000, 1000, 1000, 1000, 1000, 1000, 3, 3, 2,
0.35, 0.35, 0.35 ;
20, 20, 30 ;

fates_allom_dmode = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;

Expand All @@ -1049,21 +1049,21 @@ data:

fates_allom_h2cd2 = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;

fates_allom_hmode = 5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 1, 1 ;
fates_allom_hmode = 5, 5, 5, 5, 5, 5, 1, 1, 1, 3, 3, 3 ;

fates_allom_l2fr = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;
fates_allom_l2fr = 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.67, 0.67, 1.41 ;

fates_allom_la_per_sa_int = 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8,
0.8, 0.8, 0.8 ;

fates_allom_la_per_sa_slp = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;

fates_allom_lmode = 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1 ;
fates_allom_lmode = 2, 2, 2, 2, 2, 2, 1, 1, 1, 5, 5, 5 ;

fates_allom_sai_scaler = 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1, 0.1 ;

fates_allom_smode = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;
fates_allom_smode = 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2 ;

fates_allom_stmode = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;

Expand Down Expand Up @@ -1306,10 +1306,10 @@ data:
495 ;

fates_leaf_slamax = 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.012,
0.03, 0.03, 0.03, 0.03, 0.03 ;
0.03, 0.03, 0.05, 0.05, 0.05 ;

fates_leaf_slatop = 0.012, 0.005, 0.024, 0.009, 0.03, 0.03, 0.012, 0.03,
0.03, 0.03, 0.03, 0.03 ;
0.03, 0.05, 0.05, 0.05 ;

fates_leaf_stomatal_intercept = 10000, 10000, 10000, 10000, 10000, 10000,
10000, 10000, 10000, 10000, 10000, 40000 ;
Expand Down Expand Up @@ -1463,22 +1463,22 @@ data:
0.001, 0.001, 0.12, 0.12, 0.12 ;

fates_recruit_height_min = 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 0.2, 0.2, 0.2,
0.125, 0.125, 0.125 ;
0.2, 0.2, 0.2 ;

fates_recruit_init_density = 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2, 0.2, 0.2 ;

fates_recruit_prescribed_rate = 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02,
0.02, 0.02, 0.02, 0.02, 0.02 ;

fates_recruit_seed_alloc = 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1, 0.1 ;
fates_recruit_seed_alloc = 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0,
0, 0 ;

fates_recruit_seed_alloc_mature = 0, 0, 0, 0, 0, 0, 0.9, 0.9, 0.9, 0.9, 0.9,
0.9 ;
fates_recruit_seed_alloc_mature = 0, 0, 0, 0, 0, 0, 0.9, 0.9, 0.9, 0.25, 0.25,
0.2 ;

fates_recruit_seed_dbh_repro_threshold = 90, 80, 80, 80, 90, 80, 3, 3, 2,
0.35, 0.35, 0.35 ;
3, 3, 3 ;

fates_recruit_seed_germination_rate = 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
0.5, 0.5, 0.5, 0.5, 0.5 ;
Expand Down