Skip to content

Commit

Permalink
Merge pull request #1 from rgknox/new-arctic-shrubPFT
Browse files Browse the repository at this point in the history 
fixing pft map 10/9 removing alternative cdl
  • Loading branch information
@jenniferholm
jenniferholm authored Aug 27, 2024
2 parents 4f25d3c + f388426 commit 0a68ba0
Show file tree
Hide file tree
Showing 30 changed files with 2,973 additions and 1,228 deletions.
256 changes: 137 additions & 119 deletions biogeochem/EDCanopyStructureMod.F90
View file

Large diffs are not rendered by default.

209 changes: 143 additions & 66 deletions biogeochem/EDLoggingMortalityMod.F90
View file

Large diffs are not rendered by default.

7 changes: 4 additions & 3 deletions biogeochem/EDMortalityFunctionsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -281,7 +281,7 @@ end subroutine mortality_rates

subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, btran_ft, &
mean_temp, land_use_label, age_since_anthro_disturbance, &
frac_site_primary, harvestable_forest_c, harvest_tag)
frac_site_primary, frac_site_secondary, harvestable_forest_c, harvest_tag)

!
! !DESCRIPTION:
Expand All @@ -301,6 +301,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, btran_ft, &
integer, intent(in) :: land_use_label
real(r8), intent(in) :: age_since_anthro_disturbance
real(r8), intent(in) :: frac_site_primary
real(r8), intent(in) :: frac_site_secondary

real(r8), intent(in) :: harvestable_forest_c(:) ! total carbon available for logging, kgC site-1
integer, intent(out) :: harvest_tag(:) ! tag to record the harvest status
Expand Down Expand Up @@ -329,7 +330,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, btran_ft, &
!if trees are in the canopy, then their death is 'disturbance'. This probably needs a different terminology
call mortality_rates(currentCohort,bc_in,btran_ft, mean_temp, &
cmort,hmort,bmort,frmort, smort, asmort, dgmort)
call LoggingMortality_frac(ipft, currentCohort%dbh, currentCohort%canopy_layer, &
call LoggingMortality_frac(currentSite, bc_in, ipft, currentCohort%dbh, currentCohort%canopy_layer, &
currentCohort%lmort_direct, &
currentCohort%lmort_collateral, &
currentCohort%lmort_infra, &
Expand All @@ -339,7 +340,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, btran_ft, &
bc_in%hlm_harvest_units, &
land_use_label, &
age_since_anthro_disturbance, &
frac_site_primary, harvestable_forest_c, harvest_tag)
frac_site_primary, frac_site_secondary, harvestable_forest_c, harvest_tag)

if (currentCohort%canopy_layer > 1)then
! Include understory logging mortality rates not associated with disturbance
Expand Down
933 changes: 775 additions & 158 deletions biogeochem/EDPatchDynamicsMod.F90
View file

Large diffs are not rendered by default.

43 changes: 31 additions & 12 deletions biogeochem/EDPhysiologyMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ module EDPhysiologyMod
use FatesInterfaceTypesMod, only : hlm_parteh_mode
use FatesInterfaceTypesMod, only : hlm_use_fixed_biogeog
use FatesInterfaceTypesMod, only : hlm_use_nocomp
use EDParamsMod , only : crop_lu_pft_vector
use FatesInterfaceTypesMod, only : hlm_nitrogen_spec
use FatesInterfaceTypesMod, only : hlm_phosphorus_spec
use FatesInterfaceTypesMod, only : hlm_use_tree_damage
Expand All @@ -34,6 +35,8 @@ module EDPhysiologyMod
use FatesConstantsMod, only : g_per_kg
use FatesConstantsMod, only : ndays_per_year
use FatesConstantsMod, only : nocomp_bareground
use FatesConstantsMod, only : nocomp_bareground_land
use FatesConstantsMod, only : is_crop
use FatesConstantsMod, only : area_error_2
use EDPftvarcon , only : EDPftvarcon_inst
use PRTParametersMod , only : prt_params
Expand Down Expand Up @@ -140,7 +143,8 @@ module EDPhysiologyMod
use FatesParameterDerivedMod, only : param_derived
use FatesPlantHydraulicsMod, only : InitHydrCohort
use PRTInitParamsFatesMod, only : NewRecruitTotalStoichiometry

use FatesInterfaceTypesMod , only : hlm_use_luh

implicit none
private

Expand Down Expand Up @@ -2492,20 +2496,35 @@ subroutine recruitment(currentSite, currentPatch, bc_in)
real(r8) :: seedling_layer_smp ! soil matric potential at seedling rooting depth [mm H2O suction]
integer, parameter :: recruitstatus = 1 ! whether the newly created cohorts are recruited or initialized
integer :: ilayer_seedling_root ! the soil layer at seedling rooting depth

logical :: use_this_pft ! logical flag for whether or not to allow a given PFT to recruit
!---------------------------------------------------------------------------

do ft = 1, numpft

! The following if block is for the prescribed biogeography and/or nocomp modes.
! Since currentSite%use_this_pft is a site-level quantity and thus only limits whether a given PFT
! is permitted on a given gridcell or not, it applies to the prescribed biogeography case only.
! If nocomp is enabled, then we must determine whether a given PFT is allowed on a given patch or not.
! The following if block is for the prescribed biogeography and/or nocomp modes and/or crop land use types
! Since currentSite%use_this_pft is a site-level quantity and thus only limits whether a given PFT
! is permitted on a given gridcell or not, it applies to the prescribed biogeography case only.
! If nocomp is enabled, then we must determine whether a given PFT is allowed on a given patch or not.
! Whether or not nocomp or prescribed biogeography is enabled, if land use change is enabled, then we only want to
! allow crop PFTs on patches with crop land use types

use_this_pft = .false.
if(currentSite%use_this_pft(ft).eq.itrue &
.and. ((hlm_use_nocomp .eq. ifalse) .or. (ft .eq. currentPatch%nocomp_pft_label)))then
use_this_pft = .true.
end if

if (currentSite%use_this_pft(ft) .eq. itrue .and. &
((hlm_use_nocomp .eq. ifalse) .or. &
(ft .eq. currentPatch%nocomp_pft_label))) then
if ( currentPatch%land_use_label .ne. nocomp_bareground_land ) then ! cdk
if ((hlm_use_luh .eq. itrue) .and. (is_crop(currentPatch%land_use_label))) then
if ( crop_lu_pft_vector(currentPatch%land_use_label) .eq. ft ) then
use_this_pft = .true.
else
use_this_pft = .false.
end if
end if
endif

use_this_pft_if: if(use_this_pft) then
height = EDPftvarcon_inst%hgt_min(ft)
stem_drop_fraction = prt_params%phen_stem_drop_fraction(ft)
fnrt_drop_fraction = prt_params%phen_fnrt_drop_fraction(ft)
Expand Down Expand Up @@ -2749,11 +2768,11 @@ subroutine recruitment(currentSite, currentPatch, bc_in)
currentSite%recruitment_rate(ft) = currentSite%recruitment_rate(ft) + cohort_n

endif any_recruits
endif !use_this_pft
endif use_this_pft_if
enddo !pft loop
end subroutine recruitment

! ======================================================================================
! ======================================================================================

subroutine CWDInput( currentSite, currentPatch, litt, bc_in)

Expand Down Expand Up @@ -3022,7 +3041,7 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in)
SF_val_CWD_frac_adj(c) * dead_n_dlogging * &
prt_params%allom_agb_frac(pft)

site_mass%wood_product = site_mass%wood_product + &
site_mass%wood_product_harvest(pft) = site_mass%wood_product_harvest(pft) + &
trunk_wood * currentPatch%area * logging_export_frac

! Add AG wood to litter from the non-exported fraction of wood
Expand Down
173 changes: 168 additions & 5 deletions biogeochem/FatesAllometryMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,6 @@ module FatesAllometryMod
use FatesGlobals , only : endrun => fates_endrun
use FatesGlobals , only : FatesWarn,N2S,A2S,I2S
use EDParamsMod , only : nlevleaf,dinc_vai,dlower_vai
use EDParamsMod , only : nclmax
use DamageMainMod , only : GetCrownReduction

implicit none
Expand Down Expand Up @@ -412,6 +411,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 +473,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 +535,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 @@ -670,7 +675,7 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25
real(r8), intent(in) :: c_area ! areal extent of canopy (m2)
real(r8), intent(in) :: nplant ! number of individuals in cohort per ha
integer, intent(in) :: cl ! canopy layer index
real(r8), intent(in) :: canopy_lai(nclmax) ! total leaf area index of
real(r8), intent(in) :: canopy_lai(:) ! total leaf area index of
! each canopy layer
real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at canopy
! top, ref 25C
Expand Down Expand Up @@ -803,7 +808,7 @@ real(r8) function tree_sai(pft, dbh, crowndamage, canopy_trim, elongf_stem, c_ar
real(r8), intent(in) :: c_area ! crown area (m2)
real(r8), intent(in) :: nplant ! number of plants
integer, intent(in) :: cl ! canopy layer index
real(r8), intent(in) :: canopy_lai(nclmax) ! total leaf area index of
real(r8), intent(in) :: canopy_lai(:) ! total leaf area index of
! each canopy layer
real(r8), intent(in) :: treelai ! tree LAI for checking purposes only
real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at top of crown
Expand Down Expand Up @@ -1036,6 +1041,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)
Expand Down Expand Up @@ -1228,7 +1252,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 +1655,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 +2115,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
Loading

0 comments on commit 0a68ba0

Please sign in to comment.