From 156f7bdef52818f7e3f06b7a3e485c5076215e58 Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 16 Feb 2022 18:21:59 -0700 Subject: [PATCH 1/6] some changes to allow tracking mortality carbon flux vars --- biogeochem/EDCohortDynamicsMod.F90 | 4 +-- biogeochem/EDPatchDynamicsMod.F90 | 17 +++++---- main/EDInitMod.F90 | 12 +++---- main/EDTypesMod.F90 | 16 ++++----- main/FatesHistoryInterfaceMod.F90 | 57 +++++++++++++++++++++--------- 5 files changed, 67 insertions(+), 39 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 5446c94dda..cb11005dc4 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -825,13 +825,13 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ currentSite%term_nindivs_canopy(currentCohort%size_class,currentCohort%pft) = & currentSite%term_nindivs_canopy(currentCohort%size_class,currentCohort%pft) + currentCohort%n - currentSite%term_carbonflux_canopy = currentSite%term_carbonflux_canopy + & + currentSite%term_carbonflux_canopy(currentCohort%pft) = currentSite%term_carbonflux_canopy(currentCohort%pft) + & currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c) else currentSite%term_nindivs_ustory(currentCohort%size_class,currentCohort%pft) = & currentSite%term_nindivs_ustory(currentCohort%size_class,currentCohort%pft) + currentCohort%n - currentSite%term_carbonflux_ustory = currentSite%term_carbonflux_ustory + & + currentSite%term_carbonflux_ustory(currentCohort%pft) = currentSite%term_carbonflux_ustory(currentCohort%pft) + & currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c) end if diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 5f89ba0b07..230b8ede73 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -194,7 +194,7 @@ subroutine disturbance_rates( site_in, bc_in) ! first calculate the fractino of the site that is primary land call get_frac_site_primary(site_in, frac_site_primary) - site_in%harvest_carbon_flux = 0._r8 + site_in%harvest_carbon_flux(:) = 0._r8 currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -233,7 +233,8 @@ subroutine disturbance_rates( site_in, bc_in) ! estimate the wood product (trunk_product_site) if (currentCohort%canopy_layer>=1) then - site_in%harvest_carbon_flux = site_in%harvest_carbon_flux + & + site_in%harvest_carbon_flux(currentCohort%pft) = & + site_in%harvest_carbon_flux(currentCohort%pft) + & currentCohort%lmort_direct * currentCohort%n * & ( currentCohort%prt%GetState(sapw_organ, all_carbon_elements) + & currentCohort%prt%GetState(struct_organ, all_carbon_elements)) * & @@ -780,7 +781,8 @@ subroutine spawn_patches( currentSite, bc_in) nc%n * ED_val_understorey_death / hlm_freq_day - currentSite%imort_carbonflux = currentSite%imort_carbonflux + & + currentSite%imort_carbonflux(currentCohort%pft) = & + currentSite%imort_carbonflux(currentCohort%pft) + & (nc%n * ED_val_understorey_death / hlm_freq_day ) * & total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 @@ -856,7 +858,8 @@ subroutine spawn_patches( currentSite, bc_in) currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) + & nc%n * currentCohort%fire_mort / hlm_freq_day - currentSite%fmort_carbonflux_canopy = currentSite%fmort_carbonflux_canopy + & + currentSite%fmort_carbonflux_canopy(currentCohort%pft) = & + currentSite%fmort_carbonflux_canopy(currentCohort%pft) + & (nc%n * currentCohort%fire_mort) * & total_c * g_per_kg * days_per_sec * ha_per_m2 @@ -865,7 +868,8 @@ subroutine spawn_patches( currentSite, bc_in) currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) + & nc%n * currentCohort%fire_mort / hlm_freq_day - currentSite%fmort_carbonflux_ustory = currentSite%fmort_carbonflux_ustory + & + currentSite%fmort_carbonflux_ustory(currentCohort%pft) = & + currentSite%fmort_carbonflux_ustory(currentCohort%pft) + & (nc%n * currentCohort%fire_mort) * & total_c * g_per_kg * days_per_sec * ha_per_m2 end if @@ -997,7 +1001,8 @@ subroutine spawn_patches( currentSite, bc_in) nc%n * currentPatch%fract_ldist_not_harvested * & logging_coll_under_frac / hlm_freq_day - currentSite%imort_carbonflux = currentSite%imort_carbonflux + & + currentSite%imort_carbonflux(currentCohort%pft) = & + currentSite%imort_carbonflux(currentCohort%pft) + & (nc%n * currentPatch%fract_ldist_not_harvested * & logging_coll_under_frac/ hlm_freq_day ) * & total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index bb83ffb668..1c43542894 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -203,7 +203,7 @@ subroutine zero_site( site_in ) ! Disturbance rates tracking site_in%primary_land_patchfusion_error = 0.0_r8 - site_in%harvest_carbon_flux = 0.0_r8 + site_in%harvest_carbon_flux(:) = 0.0_r8 site_in%potential_disturbance_rates(:) = 0.0_r8 site_in%disturbance_rates_secondary_to_secondary(:) = 0.0_r8 site_in%disturbance_rates_primary_to_secondary(:) = 0.0_r8 @@ -225,15 +225,15 @@ subroutine zero_site( site_in ) ! termination and recruitment info site_in%term_nindivs_canopy(:,:) = 0._r8 site_in%term_nindivs_ustory(:,:) = 0._r8 - site_in%term_carbonflux_canopy = 0._r8 - site_in%term_carbonflux_ustory = 0._r8 + site_in%term_carbonflux_canopy(:) = 0._r8 + site_in%term_carbonflux_ustory(:) = 0._r8 site_in%recruitment_rate(:) = 0._r8 site_in%imort_rate(:,:) = 0._r8 - site_in%imort_carbonflux = 0._r8 + site_in%imort_carbonflux(:) = 0._r8 site_in%fmort_rate_canopy(:,:) = 0._r8 site_in%fmort_rate_ustory(:,:) = 0._r8 - site_in%fmort_carbonflux_canopy = 0._r8 - site_in%fmort_carbonflux_ustory = 0._r8 + site_in%fmort_carbonflux_canopy(:) = 0._r8 + site_in%fmort_carbonflux_ustory(:) = 0._r8 site_in%fmort_rate_cambial(:,:) = 0._r8 site_in%fmort_rate_crown(:,:) = 0._r8 diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 9ce4d5fe44..8275f38cb7 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -797,15 +797,16 @@ module EDTypesMod real(r8), allocatable :: term_nindivs_ustory(:,:) ! number of understory individuals that were in cohorts which ! were terminated this timestep, on size x pft - real(r8) :: term_carbonflux_canopy ! carbon flux from live to dead pools associated + real(r8) :: term_carbonflux_canopy(1:maxpft) ! carbon flux from live to dead pools associated ! with termination mortality, per canopy level - real(r8) :: term_carbonflux_ustory ! carbon flux from live to dead pools associated + real(r8) :: term_carbonflux_ustory(1:maxpft) ! carbon flux from live to dead pools associated ! with termination mortality, per canopy level - real(r8) :: demotion_carbonflux ! biomass of demoted individuals from canopy to understory [kgC/ha/day] - real(r8) :: promotion_carbonflux ! biomass of promoted individuals from understory to canopy [kgC/ha/day] - real(r8) :: imort_carbonflux ! biomass of individuals killed due to impact mortality per year. [kgC/ha/day] - real(r8) :: fmort_carbonflux_canopy ! biomass of canopy indivs killed due to fire per year. [gC/m2/sec] - real(r8) :: fmort_carbonflux_ustory ! biomass of understory indivs killed due to fire per year [gC/m2/sec] + real(r8) :: demotion_carbonflux ! biomass of demoted individuals from canopy to understory [kgC/ha/day] + real(r8) :: promotion_carbonflux ! biomass of promoted individuals from understory to canopy [kgC/ha/day] + real(r8) :: imort_carbonflux(1:maxpft) ! biomass of individuals killed due to impact mortality per year. [kgC/ha/day] + real(r8) :: fmort_carbonflux_canopy(1:maxpft) ! biomass of canopy indivs killed due to fire per year. [gC/m2/sec] + real(r8) :: fmort_carbonflux_ustory(1:maxpft) ! biomass of understory indivs killed due to fire per year [gC/m2/sec] + real(r8) :: harvest_carbon_flux(1:maxpft) ! diagnostic site level flux of carbon as harvested plants [kg C / m2 / day] real(r8) :: recruitment_rate(1:maxpft) ! number of individuals that were recruited into new cohorts real(r8), allocatable :: demotion_rate(:) ! rate of individuals demoted from canopy to understory per FATES timestep @@ -839,7 +840,6 @@ module EDTypesMod real(r8) :: disturbance_rates_secondary_to_secondary(N_DIST_TYPES) ! actual disturbance rates from secondary patches to secondary patches [m2/m2/day] real(r8) :: potential_disturbance_rates(N_DIST_TYPES) ! "potential" disturb rates (i.e. prior to the "which is most" logic) [m2/m2/day] real(r8) :: primary_land_patchfusion_error ! error term in total area of primary patches associated with patch fusion [m2/m2/day] - real(r8) :: harvest_carbon_flux ! diagnostic site level flux of carbon as harvested plants [kg C / m2 / day] end type ed_site_type diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 20918251d2..358eb0e932 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -509,6 +509,7 @@ module FatesHistoryInterfaceMod integer :: ih_nindivs_si_pft integer :: ih_recruitment_si_pft integer :: ih_mortality_si_pft + integer :: ih_mortality_carbonflux_si_pft integer :: ih_crownarea_si_pft integer :: ih_canopycrownarea_si_pft integer :: ih_gpp_si_pft @@ -1841,6 +1842,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_nindivs_si_pft => this%hvars(ih_nindivs_si_pft)%r82d, & hio_recruitment_si_pft => this%hvars(ih_recruitment_si_pft)%r82d, & hio_mortality_si_pft => this%hvars(ih_mortality_si_pft)%r82d, & + hio_mortality_carbonflux_si_pft => this%hvars(ih_mortality_carbonflux_si_pft)%r82d, & hio_crownarea_si_pft => this%hvars(ih_crownarea_si_pft)%r82d, & hio_canopycrownarea_si_pft => this%hvars(ih_canopycrownarea_si_pft)%r82d, & hio_gpp_si_pft => this%hvars(ih_gpp_si_pft)%r82d, & @@ -2179,7 +2181,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_potential_disturbance_rate_si(io_si) = sum(sites(s)%potential_disturbance_rates(1:N_DIST_TYPES)) * days_per_year ! harvest carbon flux in [kgC/m2/d] -> [kgC/m2/yr] - hio_harvest_carbonflux_si(io_si) = sites(s)%harvest_carbon_flux * & + hio_harvest_carbonflux_si(io_si) = sum(sites(s)%harvest_carbon_flux(:)) * & days_per_year ! Loop through patches to sum up diagonistics @@ -2590,6 +2592,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) ccohort%frmort*ccohort%n / m2_per_ha hio_m9_si_scls(io_si,scls) = hio_m9_si_scls(io_si,scls) + ccohort%smort*ccohort%n / m2_per_ha + !C13 discrimination if(gpp_cached + ccohort%gpp_acc_hold > 0.0_r8)then hio_c13disc_si_scpf(io_si,scpf) = ((hio_c13disc_si_scpf(io_si,scpf) * gpp_cached) + & @@ -2617,6 +2620,12 @@ subroutine update_history_dyn(this,nc,nsites,sites) alive_m = leaf_m + fnrt_m + sapw_m total_m = alive_m + store_m + struct_m + hio_mortality_carbonflux_si_pft(io_si,ccohort%pft) = hio_mortality_carbonflux_si_pft(io_si,ccohort%pft) + & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + & + ccohort%frmort + ccohort%smort + ccohort%asmort) * & + total_m * ccohort%n * days_per_sec * years_per_day * ha_per_m2 + & + (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * total_m * & + ccohort%n * ha_per_m2 ! number density by size and biomass hio_agb_si_scls(io_si,scls) = hio_agb_si_scls(io_si,scls) + & @@ -3028,14 +3037,6 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_mortality_understory_si_scls(io_si,i_scls) = hio_mortality_understory_si_scls(io_si,i_scls) + & sites(s)%fmort_rate_ustory(i_scls, i_pft) / m2_per_ha - ! - ! carbon flux associated with mortality of trees dying by fire - hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & - sites(s)%fmort_carbonflux_canopy / g_per_kg - - hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - sites(s)%fmort_carbonflux_ustory / g_per_kg - ! ! for scag variables, also treat as happening in the newly-disurbed patch @@ -3051,18 +3052,34 @@ subroutine update_history_dyn(this,nc,nsites,sites) end do end do + ! + ! carbon flux associated with mortality of trees dying by fire + hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & + sum(sites(s)%fmort_carbonflux_canopy(:)) / g_per_kg + + hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & + sum(sites(s)%fmort_carbonflux_ustory(:)) / g_per_kg + ! treat carbon flux from imort the same way hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - sites(s)%imort_carbonflux / g_per_kg + sum(sites(s)%imort_carbonflux(:)) / g_per_kg ! + do i_pft = 1, numpft + hio_mortality_carbonflux_si_pft(io_si,ccohort%pft) = hio_mortality_carbonflux_si_pft(io_si,ccohort%pft) + & + (sites(s)%fmort_carbonflux_canopy(i_pft) + & + sites(s)%fmort_carbonflux_ustory(i_pft) + & + sites(s)%imort_carbonflux(i_pft) + & + sites(s)%harvest_carbon_flux(i_pft)) / g_per_kg !cdk + end do + sites(s)%term_nindivs_canopy(:,:) = 0._r8 sites(s)%term_nindivs_ustory(:,:) = 0._r8 - sites(s)%imort_carbonflux = 0._r8 + sites(s)%imort_carbonflux(:) = 0._r8 sites(s)%imort_rate(:,:) = 0._r8 sites(s)%fmort_rate_canopy(:,:) = 0._r8 sites(s)%fmort_rate_ustory(:,:) = 0._r8 - sites(s)%fmort_carbonflux_canopy = 0._r8 - sites(s)%fmort_carbonflux_ustory = 0._r8 + sites(s)%fmort_carbonflux_canopy(:) = 0._r8 + sites(s)%fmort_carbonflux_ustory(:) = 0._r8 sites(s)%fmort_rate_cambial(:,:) = 0._r8 sites(s)%fmort_rate_crown(:,:) = 0._r8 sites(s)%growthflux_fusion(:,:) = 0._r8 @@ -3462,14 +3479,14 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! mortality-associated carbon fluxes hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & - sites(s)%term_carbonflux_canopy * days_per_sec * ha_per_m2 + sum(sites(s)%term_carbonflux_canopy(:)) * days_per_sec * ha_per_m2 hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - sites(s)%term_carbonflux_ustory * days_per_sec * ha_per_m2 + sum(sites(s)%term_carbonflux_ustory(:)) * days_per_sec * ha_per_m2 ! and zero the site-level termination carbon flux variable - sites(s)%term_carbonflux_canopy = 0._r8 - sites(s)%term_carbonflux_ustory = 0._r8 + sites(s)%term_carbonflux_canopy(:) = 0._r8 + sites(s)%term_carbonflux_ustory(:) = 0._r8 ! ! add the site-level disturbance-associated cwd and litter input fluxes to thir respective flux fields @@ -5438,6 +5455,12 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, & index = ih_understory_mortality_carbonflux_si) + call this%set_history_var(vname='FATES_MORTALITY_CFLUX_PF', units='kg m-2 s-1', & + long='PFT-level flux of biomass carbon from live to dead pool from mortality', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_mortality_carbonflux_si_pft) + ! size class by age dimensioned variables call this%set_history_var(vname='FATES_NPLANT_SZAP', units = 'm-2', & From 9ccbcb6fcbaed08d1f733513caafe5098abb3102 Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 17 Feb 2022 13:22:50 -0700 Subject: [PATCH 2/6] undid changes to harvest_carbon_flux and added restart variables --- biogeochem/EDPatchDynamicsMod.F90 | 5 ++-- main/EDInitMod.F90 | 2 +- main/EDTypesMod.F90 | 2 +- main/FatesHistoryInterfaceMod.F90 | 5 ++-- main/FatesRestartInterfaceMod.F90 | 47 +++++++++++++++++++------------ 5 files changed, 35 insertions(+), 26 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 230b8ede73..e8e0190e5e 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -194,7 +194,7 @@ subroutine disturbance_rates( site_in, bc_in) ! first calculate the fractino of the site that is primary land call get_frac_site_primary(site_in, frac_site_primary) - site_in%harvest_carbon_flux(:) = 0._r8 + site_in%harvest_carbon_flux = 0._r8 currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -233,8 +233,7 @@ subroutine disturbance_rates( site_in, bc_in) ! estimate the wood product (trunk_product_site) if (currentCohort%canopy_layer>=1) then - site_in%harvest_carbon_flux(currentCohort%pft) = & - site_in%harvest_carbon_flux(currentCohort%pft) + & + site_in%harvest_carbon_flux = site_in%harvest_carbon_flux + & currentCohort%lmort_direct * currentCohort%n * & ( currentCohort%prt%GetState(sapw_organ, all_carbon_elements) + & currentCohort%prt%GetState(struct_organ, all_carbon_elements)) * & diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 1c43542894..63ec3fc280 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -203,7 +203,7 @@ subroutine zero_site( site_in ) ! Disturbance rates tracking site_in%primary_land_patchfusion_error = 0.0_r8 - site_in%harvest_carbon_flux(:) = 0.0_r8 + site_in%harvest_carbon_flux = 0.0_r8 site_in%potential_disturbance_rates(:) = 0.0_r8 site_in%disturbance_rates_secondary_to_secondary(:) = 0.0_r8 site_in%disturbance_rates_primary_to_secondary(:) = 0.0_r8 diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 8275f38cb7..3c9e44b5b8 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -806,7 +806,7 @@ module EDTypesMod real(r8) :: imort_carbonflux(1:maxpft) ! biomass of individuals killed due to impact mortality per year. [kgC/ha/day] real(r8) :: fmort_carbonflux_canopy(1:maxpft) ! biomass of canopy indivs killed due to fire per year. [gC/m2/sec] real(r8) :: fmort_carbonflux_ustory(1:maxpft) ! biomass of understory indivs killed due to fire per year [gC/m2/sec] - real(r8) :: harvest_carbon_flux(1:maxpft) ! diagnostic site level flux of carbon as harvested plants [kg C / m2 / day] + real(r8) :: harvest_carbon_flux ! diagnostic site level flux of carbon as harvested plants [kg C / m2 / day] real(r8) :: recruitment_rate(1:maxpft) ! number of individuals that were recruited into new cohorts real(r8), allocatable :: demotion_rate(:) ! rate of individuals demoted from canopy to understory per FATES timestep diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 358eb0e932..5c784a6677 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -2181,7 +2181,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_potential_disturbance_rate_si(io_si) = sum(sites(s)%potential_disturbance_rates(1:N_DIST_TYPES)) * days_per_year ! harvest carbon flux in [kgC/m2/d] -> [kgC/m2/yr] - hio_harvest_carbonflux_si(io_si) = sum(sites(s)%harvest_carbon_flux(:)) * & + hio_harvest_carbonflux_si(io_si) = sites(s)%harvest_carbon_flux * & days_per_year ! Loop through patches to sum up diagonistics @@ -3068,8 +3068,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_mortality_carbonflux_si_pft(io_si,ccohort%pft) = hio_mortality_carbonflux_si_pft(io_si,ccohort%pft) + & (sites(s)%fmort_carbonflux_canopy(i_pft) + & sites(s)%fmort_carbonflux_ustory(i_pft) + & - sites(s)%imort_carbonflux(i_pft) + & - sites(s)%harvest_carbon_flux(i_pft)) / g_per_kg !cdk + sites(s)%imort_carbonflux(i_pft) ) / g_per_kg !cdk end do sites(s)%term_nindivs_canopy(:,:) = 0._r8 diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 0605767cd6..e2824a716b 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -208,9 +208,9 @@ module FatesRestartInterfaceMod integer :: ir_termcflux_usto_si integer :: ir_democflux_si integer :: ir_promcflux_si - integer :: ir_imortcflux_si - integer :: ir_fmortcflux_cano_si - integer :: ir_fmortcflux_usto_si + integer :: ir_imortcflux_sipft + integer :: ir_fmortcflux_cano_sipft + integer :: ir_fmortcflux_usto_sipft integer :: ir_cwdagin_flxdg integer :: ir_cwdbgin_flxdg integer :: ir_leaflittin_flxdg @@ -1230,20 +1230,20 @@ subroutine define_restart_vars(this, initialize_variables) units='indiv/ha/da', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_promrate_sisc) - call this%set_restart_var(vname='fates_imortcflux', vtype=site_r8, & + call this%set_restart_var(vname='fates_imortcflux', vtype=cohort_r8, & long_name='biomass of indivs killed due to impact mort', & units='kgC/ha/day', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_imortcflux_si) + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_imortcflux_sipft) - call this%set_restart_var(vname='fates_fmortcflux_canopy', vtype=site_r8, & + call this%set_restart_var(vname='fates_fmortcflux_canopy', vtype=cohort_r8, & long_name='fates diagnostic biomass of canopy fire', & units='gC/m2/sec', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmortcflux_cano_si) + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmortcflux_cano_sipft) - call this%set_restart_var(vname='fates_fmortcflux_ustory', vtype=site_r8, & + call this%set_restart_var(vname='fates_fmortcflux_ustory', vtype=cohort_r8, & long_name='fates diagnostic biomass of understory fire', & units='gC/m2/sec', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmortcflux_usto_si) + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmortcflux_usto_sipft) call this%set_restart_var(vname='fates_termcflux_canopy', vtype=site_r8, & long_name='fates diagnostic term carbon flux canopy', & @@ -2224,9 +2224,16 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_termcflux_usto_si(io_idx_si) = sites(s)%term_carbonflux_ustory rio_democflux_si(io_idx_si) = sites(s)%demotion_carbonflux rio_promcflux_si(io_idx_si) = sites(s)%promotion_carbonflux - rio_imortcflux_si(io_idx_si) = sites(s)%imort_carbonflux - rio_fmortcflux_cano_si(io_idx_si) = sites(s)%fmort_carbonflux_canopy - rio_fmortcflux_usto_si(io_idx_si) = sites(s)%fmort_carbonflux_ustory + + io_idx_si_pft = io_idx_co_1st + do i_pft = 1, numpft + rio_fmortcflux_cano_sipft(io_idx_si, io_idx_si_pft) = sites(s)%fmort_carbonflux_canopy(i_pft) + rio_fmortcflux_usto_sipft(io_idx_si, io_idx_si_pft) = sites(s)%fmort_carbonflux_ustory(i_pft) + + rio_imortcflux_sipft(io_idx_si, io_idx_si_pft) = sites(s)%imort_carbonflux(i_pft) + + io_idx_si_pft = io_idx_si_pft + 1 + end do rio_cd_status_si(io_idx_si) = sites(s)%cstatus rio_dd_status_si(io_idx_si) = sites(s)%dstatus @@ -2663,9 +2670,9 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_termcflux_usto_si => this%rvars(ir_termcflux_usto_si)%r81d, & rio_democflux_si => this%rvars(ir_democflux_si)%r81d, & rio_promcflux_si => this%rvars(ir_promcflux_si)%r81d, & - rio_imortcflux_si => this%rvars(ir_imortcflux_si)%r81d, & - rio_fmortcflux_cano_si => this%rvars(ir_fmortcflux_cano_si)%r81d, & - rio_fmortcflux_usto_si => this%rvars(ir_fmortcflux_usto_si)%r81d) + rio_imortcflux_sipft => this%rvars(ir_imortcflux_sipft)%r81d, & + rio_fmortcflux_cano_sipft => this%rvars(ir_fmortcflux_cano_sipft)%r81d, & + rio_fmortcflux_usto_sipft => this%rvars(ir_fmortcflux_usto_sipft)%r81d) totalcohorts = 0 @@ -3082,10 +3089,14 @@ subroutine get_restart_vectors(this, nc, nsites, sites) sites(s)%term_carbonflux_ustory = rio_termcflux_usto_si(io_idx_si) sites(s)%demotion_carbonflux = rio_democflux_si(io_idx_si) sites(s)%promotion_carbonflux = rio_promcflux_si(io_idx_si) - sites(s)%imort_carbonflux = rio_imortcflux_si(io_idx_si) - sites(s)%fmort_carbonflux_canopy = rio_fmortcflux_cano_si(io_idx_si) - sites(s)%fmort_carbonflux_ustory = rio_fmortcflux_usto_si(io_idx_si) + io_idx_si_pft = io_idx_co_1st + do i_pft = 1, numpft + sites(s)%fmort_carbonflux_canopy(i_pft) = rio_fmortcflux_cano_sipft(io_idx_si, io_idx_si_pft) + sites(s)%fmort_carbonflux_ustory(i_pft) = rio_fmortcflux_usto_sipft(io_idx_si, io_idx_si_pft) + sites(s)%imort_carbonflux(i_pft) = rio_imortcflux_sipft(io_idx_si, io_idx_si_pft) + io_idx_si_pft = io_idx_si_pft + 1 + end do ! Site level phenology status flags From 96aac267371c43bbf27f361aa2fa4898109e6498 Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 17 Feb 2022 18:18:14 -0700 Subject: [PATCH 3/6] fixes to restart file --- main/FatesRestartInterfaceMod.F90 | 46 +++++++++++++++---------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index e2824a716b..6a3297714f 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -204,8 +204,8 @@ module FatesRestartInterfaceMod integer :: ir_growflx_fusion_siscpf integer :: ir_demorate_sisc integer :: ir_promrate_sisc - integer :: ir_termcflux_cano_si - integer :: ir_termcflux_usto_si + integer :: ir_termcflux_cano_sipft + integer :: ir_termcflux_usto_sipft integer :: ir_democflux_si integer :: ir_promcflux_si integer :: ir_imortcflux_sipft @@ -1245,15 +1245,15 @@ subroutine define_restart_vars(this, initialize_variables) units='gC/m2/sec', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fmortcflux_usto_sipft) - call this%set_restart_var(vname='fates_termcflux_canopy', vtype=site_r8, & + call this%set_restart_var(vname='fates_termcflux_canopy', vtype=cohort_r8, & long_name='fates diagnostic term carbon flux canopy', & units='', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_termcflux_cano_si ) + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_termcflux_cano_sipft ) - call this%set_restart_var(vname='fates_termcflux_ustory', vtype=site_r8, & + call this%set_restart_var(vname='fates_termcflux_ustory', vtype=cohort_r8, & long_name='fates diagnostic term carbon flux understory', & units='', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_termcflux_usto_si ) + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_termcflux_usto_sipft ) call this%set_restart_var(vname='fates_democflux', vtype=site_r8, & long_name='fates diagnostic demotion carbon flux', & @@ -1828,13 +1828,13 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_growflx_fusion_siscpf => this%rvars(ir_growflx_fusion_siscpf)%r81d, & rio_demorate_sisc => this%rvars(ir_demorate_sisc)%r81d, & rio_promrate_sisc => this%rvars(ir_promrate_sisc)%r81d, & - rio_termcflux_cano_si => this%rvars(ir_termcflux_cano_si)%r81d, & - rio_termcflux_usto_si => this%rvars(ir_termcflux_usto_si)%r81d, & + rio_termcflux_cano_sipft => this%rvars(ir_termcflux_cano_sipft)%r81d, & + rio_termcflux_usto_sipft => this%rvars(ir_termcflux_usto_sipft)%r81d, & rio_democflux_si => this%rvars(ir_democflux_si)%r81d, & rio_promcflux_si => this%rvars(ir_promcflux_si)%r81d, & - rio_imortcflux_si => this%rvars(ir_imortcflux_si)%r81d, & - rio_fmortcflux_cano_si => this%rvars(ir_fmortcflux_cano_si)%r81d, & - rio_fmortcflux_usto_si => this%rvars(ir_fmortcflux_usto_si)%r81d) + rio_imortcflux_sipft => this%rvars(ir_imortcflux_sipft)%r81d, & + rio_fmortcflux_cano_sipft => this%rvars(ir_fmortcflux_cano_sipft)%r81d, & + rio_fmortcflux_usto_sipft => this%rvars(ir_fmortcflux_usto_sipft)%r81d) totalCohorts = 0 @@ -2220,17 +2220,17 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_si_sc = io_idx_si_sc + 1 end do - rio_termcflux_cano_si(io_idx_si) = sites(s)%term_carbonflux_canopy - rio_termcflux_usto_si(io_idx_si) = sites(s)%term_carbonflux_ustory rio_democflux_si(io_idx_si) = sites(s)%demotion_carbonflux rio_promcflux_si(io_idx_si) = sites(s)%promotion_carbonflux io_idx_si_pft = io_idx_co_1st do i_pft = 1, numpft - rio_fmortcflux_cano_sipft(io_idx_si, io_idx_si_pft) = sites(s)%fmort_carbonflux_canopy(i_pft) - rio_fmortcflux_usto_sipft(io_idx_si, io_idx_si_pft) = sites(s)%fmort_carbonflux_ustory(i_pft) + rio_termcflux_cano_sipft(io_idx_si_pft) = sites(s)%term_carbonflux_canopy(i_pft) + rio_termcflux_usto_sipft(io_idx_si_pft) = sites(s)%term_carbonflux_ustory(i_pft) + rio_fmortcflux_cano_sipft(io_idx_si_pft) = sites(s)%fmort_carbonflux_canopy(i_pft) + rio_fmortcflux_usto_sipft(io_idx_si_pft) = sites(s)%fmort_carbonflux_ustory(i_pft) - rio_imortcflux_sipft(io_idx_si, io_idx_si_pft) = sites(s)%imort_carbonflux(i_pft) + rio_imortcflux_sipft(io_idx_si_pft) = sites(s)%imort_carbonflux(i_pft) io_idx_si_pft = io_idx_si_pft + 1 end do @@ -2666,8 +2666,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_growflx_fusion_siscpf => this%rvars(ir_growflx_fusion_siscpf)%r81d, & rio_demorate_sisc => this%rvars(ir_demorate_sisc)%r81d, & rio_promrate_sisc => this%rvars(ir_promrate_sisc)%r81d, & - rio_termcflux_cano_si => this%rvars(ir_termcflux_cano_si)%r81d, & - rio_termcflux_usto_si => this%rvars(ir_termcflux_usto_si)%r81d, & + rio_termcflux_cano_sipft => this%rvars(ir_termcflux_cano_sipft)%r81d, & + rio_termcflux_usto_sipft => this%rvars(ir_termcflux_usto_sipft)%r81d, & rio_democflux_si => this%rvars(ir_democflux_si)%r81d, & rio_promcflux_si => this%rvars(ir_promcflux_si)%r81d, & rio_imortcflux_sipft => this%rvars(ir_imortcflux_sipft)%r81d, & @@ -3085,16 +3085,16 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_si_sc = io_idx_si_sc + 1 end do - sites(s)%term_carbonflux_canopy = rio_termcflux_cano_si(io_idx_si) - sites(s)%term_carbonflux_ustory = rio_termcflux_usto_si(io_idx_si) sites(s)%demotion_carbonflux = rio_democflux_si(io_idx_si) sites(s)%promotion_carbonflux = rio_promcflux_si(io_idx_si) io_idx_si_pft = io_idx_co_1st do i_pft = 1, numpft - sites(s)%fmort_carbonflux_canopy(i_pft) = rio_fmortcflux_cano_sipft(io_idx_si, io_idx_si_pft) - sites(s)%fmort_carbonflux_ustory(i_pft) = rio_fmortcflux_usto_sipft(io_idx_si, io_idx_si_pft) - sites(s)%imort_carbonflux(i_pft) = rio_imortcflux_sipft(io_idx_si, io_idx_si_pft) + sites(s)%term_carbonflux_canopy(i_pft) = rio_termcflux_cano_sipft(io_idx_si_pft) + sites(s)%term_carbonflux_ustory(i_pft) = rio_termcflux_usto_sipft(io_idx_si_pft) + sites(s)%fmort_carbonflux_canopy(i_pft) = rio_fmortcflux_cano_sipft(io_idx_si_pft) + sites(s)%fmort_carbonflux_ustory(i_pft) = rio_fmortcflux_usto_sipft(io_idx_si_pft) + sites(s)%imort_carbonflux(i_pft) = rio_imortcflux_sipft(io_idx_si_pft) io_idx_si_pft = io_idx_si_pft + 1 end do From 73c6c53b42352f24b05fb4290b9711999400ed6e Mon Sep 17 00:00:00 2001 From: ckoven Date: Fri, 4 Mar 2022 10:55:40 -0700 Subject: [PATCH 4/6] added drought and c starvation mortality metrics --- main/FatesHistoryInterfaceMod.F90 | 36 +++++++++++++++++++++++++++++-- 1 file changed, 34 insertions(+), 2 deletions(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 5c784a6677..2a0975d51d 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -510,6 +510,9 @@ module FatesHistoryInterfaceMod integer :: ih_recruitment_si_pft integer :: ih_mortality_si_pft integer :: ih_mortality_carbonflux_si_pft + integer :: ih_hydraulicmortality_carbonflux_si_pft + integer :: ih_cstarvmortality_carbonflux_si_pft + integer :: ih_firemortality_carbonflux_si_pft integer :: ih_crownarea_si_pft integer :: ih_canopycrownarea_si_pft integer :: ih_gpp_si_pft @@ -1843,6 +1846,9 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_recruitment_si_pft => this%hvars(ih_recruitment_si_pft)%r82d, & hio_mortality_si_pft => this%hvars(ih_mortality_si_pft)%r82d, & hio_mortality_carbonflux_si_pft => this%hvars(ih_mortality_carbonflux_si_pft)%r82d, & + hio_cstarvmortality_carbonflux_si_pft => this%hvars(ih_cstarvmortality_carbonflux_si_pft)%r82d, & + hio_hydraulicmortality_carbonflux_si_pft => this%hvars(ih_hydraulicmortality_carbonflux_si_pft)%r82d, & + hio_firemortality_carbonflux_si_pft => this%hvars(ih_firemortality_carbonflux_si_pft)%r82d, & hio_crownarea_si_pft => this%hvars(ih_crownarea_si_pft)%r82d, & hio_canopycrownarea_si_pft => this%hvars(ih_canopycrownarea_si_pft)%r82d, & hio_gpp_si_pft => this%hvars(ih_gpp_si_pft)%r82d, & @@ -2627,6 +2633,12 @@ subroutine update_history_dyn(this,nc,nsites,sites) (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * total_m * & ccohort%n * ha_per_m2 + hio_hydraulicmortality_carbonflux_si_pft(io_si,ccohort%pft) = hio_hydraulicmortality_carbonflux_si_pft(io_si,ccohort%pft) + & + ccohort%hmort * total_m * ccohort%n * days_per_sec * years_per_day * ha_per_m2 + + hio_cstarvmortality_carbonflux_si_pft(io_si,ccohort%pft) = hio_cstarvmortality_carbonflux_si_pft(io_si,ccohort%pft) + & + ccohort%cmort * total_m * ccohort%n * days_per_sec * years_per_day * ha_per_m2 + ! number density by size and biomass hio_agb_si_scls(io_si,scls) = hio_agb_si_scls(io_si,scls) + & total_m * ccohort%n * prt_params%allom_agb_frac(ccohort%pft) * AREA_INV @@ -3063,12 +3075,14 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! treat carbon flux from imort the same way hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & sum(sites(s)%imort_carbonflux(:)) / g_per_kg - ! + do i_pft = 1, numpft - hio_mortality_carbonflux_si_pft(io_si,ccohort%pft) = hio_mortality_carbonflux_si_pft(io_si,ccohort%pft) + & + hio_mortality_carbonflux_si_pft(io_si,i_pft) = hio_mortality_carbonflux_si_pft(io_si,i_pft) + & (sites(s)%fmort_carbonflux_canopy(i_pft) + & sites(s)%fmort_carbonflux_ustory(i_pft) + & sites(s)%imort_carbonflux(i_pft) ) / g_per_kg !cdk + + hio_firemortality_carbonflux_si_pft(io_si,i_pft) = sites(s)%fmort_carbonflux_canopy(i_pft) / g_per_kg end do sites(s)%term_nindivs_canopy(:,:) = 0._r8 @@ -5460,6 +5474,24 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_mortality_carbonflux_si_pft) + call this%set_history_var(vname='FATES_MORTALITY_FIRE_CFLUX_PF', units='kg m-2 s-1', & + long='PFT-level flux of biomass carbon from live to dead pool from fire mortality', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_firemortality_carbonflux_si_pft) + + call this%set_history_var(vname='FATES_MORTALITY_HYDRAULIC_CFLUX_PF', units='kg m-2 s-1', & + long='PFT-level flux of biomass carbon from live to dead pool from hydraulic failure mortality', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_hydraulicmortality_carbonflux_si_pft) + + call this%set_history_var(vname='FATES_MORTALITY_CSTARV_CFLUX_PF', units='kg m-2 s-1', & + long='PFT-level flux of biomass carbon from live to dead pool from carbon starvation mortality', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_cstarvmortality_carbonflux_si_pft) + ! size class by age dimensioned variables call this%set_history_var(vname='FATES_NPLANT_SZAP', units = 'm-2', & From c4d96b58038c51a4f52cdd2e4f0d5050c0b3ab36 Mon Sep 17 00:00:00 2001 From: Marcos Longo Date: Tue, 3 May 2022 11:54:38 -0700 Subject: [PATCH 5/6] Revised the carbon allocation routine to ensure allocation to fine roots and leaves are truly proportional to demand. Description: 1. Bug fix in DailyPRTAllometricCarbon (parteh/PRTAllometricCarbonMod.F90). When allocating to different tissues, the code was subtracting allocation of tissues before calculating the amount for the next tissue, potentially under-allocating carbon to fine roots. 2. Minor code updates to simplify allocation calculations. It now uses a single allocation factor based on availability and total need, which is applied to all tissues. 3. Replaced a few if statements with select case, which simplifies adding other hypotheses in the future (and it's safer for selecting cases). This pull request addresses the bug discussed in #784 and supersedes pull request #800. Additional minor changes in former PR #800 will be added as subsequent pull requests. --- biogeochem/FatesSoilBGCFluxMod.F90 | 5 +- biogeophys/FatesPlantRespPhotosynthMod.F90 | 1 + main/EDPftvarcon.F90 | 11 +- parteh/PRTAllometricCarbonMod.F90 | 150 +++++++++++---------- parteh/PRTParamsFATESMod.F90 | 43 +++--- 5 files changed, 113 insertions(+), 97 deletions(-) diff --git a/biogeochem/FatesSoilBGCFluxMod.F90 b/biogeochem/FatesSoilBGCFluxMod.F90 index d14ce7b005..cebab25728 100644 --- a/biogeochem/FatesSoilBGCFluxMod.F90 +++ b/biogeochem/FatesSoilBGCFluxMod.F90 @@ -219,7 +219,8 @@ subroutine UnPackNutrientAquisitionBCs(sites, bc_in) end do ! We can exit if this is a c-only simulation - if(hlm_parteh_mode.eq.prt_carbon_allom_hyp) then + select case (hlm_parteh_mode) + case (prt_carbon_allom_hyp) ! These can now be zero'd do s = 1, nsites bc_in(s)%plant_nh4_uptake_flux(:,:) = 0._r8 @@ -227,7 +228,7 @@ subroutine UnPackNutrientAquisitionBCs(sites, bc_in) bc_in(s)%plant_p_uptake_flux(:,:) = 0._r8 end do return - end if + end select do s = 1, nsites diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index adffe67d85..08607ec9f2 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -500,6 +500,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) end select + ! MLO - Shouldn't these numbers be parameters too? lmr25top = 2.525e-6_r8 * (1.5_r8 ** ((25._r8 - 20._r8)/10._r8)) lmr25top = lmr25top * lnc_top / (umolC_to_kgC * g_per_kg) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index a149132a8d..b65cad5cc1 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -1502,7 +1502,8 @@ subroutine FatesCheckParams(is_master) if(.not.is_master) return - if (hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) then + select case (hlm_parteh_mode) + case (prt_cnp_flex_allom_hyp) ! Check to see if either RD/ECA/MIC is turned on @@ -1538,14 +1539,18 @@ subroutine FatesCheckParams(is_master) end if end if - elseif (hlm_parteh_mode .ne. prt_carbon_allom_hyp) then + case (prt_carbon_allom_hyp) + ! No additional checks needed for now. + continue + + case default write(fates_log(),*) 'FATES Plant Allocation and Reactive Transport has' write(fates_log(),*) 'only 2 modules supported, allometric carbon and CNP.' write(fates_log(),*) 'fates_parteh_mode must be set to 1 or 2 in the namelist' write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + end select ! If any PFTs are specified as either prescribed N or P uptake ! then they all must be ! diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index 5bdf624502..814d3daa36 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -51,6 +51,9 @@ module PRTAllometricCarbonMod use PRTParametersMod , only : prt_params + use EDTypesMod , only : leaves_on + use EDTypesMod , only : leaves_off + implicit none private @@ -114,7 +117,7 @@ module PRTAllometricCarbonMod ! ------------------------------------------------------------------------------------- - type, public, extends(prt_vartypes) :: callom_prt_vartypes + type, public, extends(prt_vartypes) :: callom_prt_vartypes contains @@ -144,7 +147,7 @@ module PRTAllometricCarbonMod public :: InitPRTGlobalAllometricCarbon -contains + contains subroutine InitPRTGlobalAllometricCarbon() @@ -313,6 +316,9 @@ subroutine DailyPRTAllometricCarbon(this) real(r8) :: struct_below_target ! dead (structural) biomass below target amount [kgC] real(r8) :: total_below_target ! total biomass below the allometric target [kgC] + real(r8) :: allocation_factor ! allocation factor (relative to demand) to + ! reconstruct tissues + real(r8) :: flux_adj ! adjustment made to growth flux term to minimize error [kgC] real(r8) :: store_target_fraction ! ratio between storage and leaf biomass when on allometry [kgC] @@ -451,12 +457,13 @@ subroutine DailyPRTAllometricCarbon(this) call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) ! Target leaf biomass according to allometry and trimming - if(leaf_status==2) then + select case (leaf_status) + case (leaves_on) call bleaf(dbh,ipft,canopy_trim,target_leaf_c) - else + case (leaves_off) target_leaf_c = 0._r8 - end if - + end select + ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) @@ -466,43 +473,42 @@ subroutine DailyPRTAllometricCarbon(this) ! ----------------------------------------------------------------------------------- ! III. Prioritize some amount of carbon to replace leaf/root turnover - ! Make sure it isnt a negative payment, and either pay what is available + ! Make sure it isn't a negative payment, and either pay what is available ! or forcefully pay from storage. ! ----------------------------------------------------------------------------------- - if( prt_params%evergreen(ipft) ==1 ) then + if( prt_params%evergreen(ipft) == itrue ) then leaf_c_demand = max(0.0_r8, & prt_params%leaf_stor_priority(ipft)*sum(this%variables(leaf_c_id)%turnover(:))) else leaf_c_demand = 0.0_r8 end if - fnrt_c_demand = max(0.0_r8, & + fnrt_c_demand = max(0.0_r8, & prt_params%leaf_stor_priority(ipft)*this%variables(fnrt_c_id)%turnover(icd)) total_c_demand = leaf_c_demand + fnrt_c_demand - - if (total_c_demand> nearzero ) then + + if (total_c_demand > nearzero) then ! We pay this even if we don't have the carbon ! Just don't pay so much carbon that storage+carbon_balance can't pay for it + allocation_factor = max(0.0_r8,min(1.0_r8,(store_c+carbon_balance)/total_c_demand)) - leaf_c_flux = min(leaf_c_demand, & - max(0.0_r8,(store_c+carbon_balance)* & - (leaf_c_demand/total_c_demand))) - - ! Add carbon to the youngest age pool (i.e iexp_leaf = index 1) - carbon_balance = carbon_balance - leaf_c_flux - leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux - - ! If we are testing b4b, then we pay this even if we don't have the carbon - fnrt_c_flux = min(fnrt_c_demand, & - max(0.0_r8, (store_c+carbon_balance)* & - (fnrt_c_demand/total_c_demand))) + ! MLO. Edited the code to switch the order of operations. The previous code would + ! subtract leaf flux from carbon balance before estimating the fine root flux, + ! potentially allowing less fluxes to fine roots than possible. + leaf_c_flux = leaf_c_demand * allocation_factor + fnrt_c_flux = fnrt_c_demand * allocation_factor - carbon_balance = carbon_balance - fnrt_c_flux - fnrt_c = fnrt_c + fnrt_c_flux + ! Add carbon to the youngest age pool (i.e iexp_leaf = index 1) and fine roots + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux + fnrt_c = fnrt_c + fnrt_c_flux + ! Remove fluxes from carbon balance. In case we may have drawn carbon from storage, + ! carbon_balance will become negative, in which case we will deplete carbon from + ! storage in the next step. + carbon_balance = carbon_balance - ( leaf_c_flux + fnrt_c_flux ) end if ! ----------------------------------------------------------------------------------- @@ -511,19 +517,24 @@ subroutine DailyPRTAllometricCarbon(this) ! ----------------------------------------------------------------------------------- if( carbon_balance < 0.0_r8 ) then - + + ! Store_c_flux will be negative, so store_c will be depleted store_c_flux = carbon_balance carbon_balance = carbon_balance - store_c_flux store_c = store_c + store_c_flux else - store_below_target = max(target_store_c - store_c,0.0_r8) + ! Accumulate some carbon in storage. If storage is completely depleted, aim to + ! increase storage, but not to replenish completely so we can still use some + ! carbon for growth. + store_below_target = max(0.0_r8,target_store_c - store_c) store_target_fraction = max(0.0_r8, store_c/target_store_c ) store_c_flux = min(store_below_target,carbon_balance * & max(exp(-1.*store_target_fraction**4._r8) - exp( -1.0_r8 ),0.0_r8)) + ! Move carbon from carbon balance to storage carbon_balance = carbon_balance - store_c_flux store_c = store_c + store_c_flux @@ -533,24 +544,29 @@ subroutine DailyPRTAllometricCarbon(this) ! V. If carbon is still available, prioritize some allocation to replace ! the rest of the leaf/fineroot deficit ! carbon balance is guaranteed to be >=0 beyond this point + ! MLO. Renamed demand with below target to make it consistent with the + ! definitions at the variable declaration part. ! ----------------------------------------------------------------------------------- - - leaf_c_demand = max(0.0_r8,(target_leaf_c - sum(leaf_c(1:nleafage)))) - fnrt_c_demand = max(0.0_r8,(target_fnrt_c - fnrt_c)) - total_c_demand = leaf_c_demand + fnrt_c_demand - - if( (carbon_balance > nearzero ) .and. (total_c_demand>nearzero)) then + leaf_below_target = max(0.0_r8,target_leaf_c - sum(leaf_c(1:nleafage))) + fnrt_below_target = max(0.0_r8,target_fnrt_c - fnrt_c) + + total_below_target = leaf_below_target + fnrt_below_target + + if ( (carbon_balance > nearzero) .and. (total_below_target > nearzero) ) then + ! Find fraction of carbon to be allocated to leaves and fine roots + allocation_factor = min(1.0_r8, carbon_balance / total_below_target) + + ! MLO. Edited the code to switch the order of operations. The previous code would + ! subtract leaf flux from carbon balance before estimating the fine root flux, + ! potentially allowing less fluxes to fine roots than possible. + leaf_c_flux = leaf_below_target * allocation_factor + fnrt_c_flux = fnrt_below_target * allocation_factor - leaf_c_flux = min(leaf_c_demand, & - carbon_balance*(leaf_c_demand/total_c_demand)) - carbon_balance = carbon_balance - leaf_c_flux leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux - - fnrt_c_flux = min(fnrt_c_demand, & - carbon_balance*(fnrt_c_demand/total_c_demand)) - carbon_balance = carbon_balance - fnrt_c_flux - fnrt_c = fnrt_c + fnrt_c_flux + fnrt_c = fnrt_c + fnrt_c_flux + + carbon_balance = carbon_balance - ( leaf_c_flux + fnrt_c_flux ) end if @@ -571,31 +587,22 @@ subroutine DailyPRTAllometricCarbon(this) sapw_below_target + store_below_target if ( total_below_target > nearzero ) then - - if( total_below_target > carbon_balance) then - leaf_c_flux = carbon_balance * leaf_below_target/total_below_target - fnrt_c_flux = carbon_balance * fnrt_below_target/total_below_target - sapw_c_flux = carbon_balance * sapw_below_target/total_below_target - store_c_flux = carbon_balance * store_below_target/total_below_target - else - leaf_c_flux = leaf_below_target - fnrt_c_flux = fnrt_below_target - sapw_c_flux = sapw_below_target - store_c_flux = store_below_target - end if - - carbon_balance = carbon_balance - leaf_c_flux - leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux - - carbon_balance = carbon_balance - fnrt_c_flux - fnrt_c = fnrt_c + fnrt_c_flux - - carbon_balance = carbon_balance - sapw_c_flux - sapw_c = sapw_c + sapw_c_flux - - carbon_balance = carbon_balance - store_c_flux - store_c = store_c + store_c_flux - + ! Find allocation factor based on available carbon and total demand to meet target. + allocation_factor = min(1.0_r8, carbon_balance / total_below_target) + + ! Find fluxes to individual pools + leaf_c_flux = leaf_below_target * allocation_factor + fnrt_c_flux = fnrt_below_target * allocation_factor + sapw_c_flux = sapw_below_target * allocation_factor + store_c_flux = store_below_target * allocation_factor + + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux + fnrt_c = fnrt_c + fnrt_c_flux + sapw_c = sapw_c + sapw_c_flux + store_c = store_c + store_c_flux + + carbon_balance = carbon_balance - & + ( leaf_c_flux + fnrt_c_flux + sapw_c_flux + store_c_flux ) end if end if @@ -690,11 +697,12 @@ subroutine DailyPRTAllometricCarbon(this) c_pool(dbh_id) = dbh ! Only grow leaves if we are in a "leaf-on" status - if(leaf_status==2) then - c_mask(leaf_c_id) = grow_leaf - else - c_mask(leaf_c_id) = .false. - end if + select case (leaf_status) + case (leaves_on) + c_mask(leaf_c_id) = grow_leaf + case default + c_mask(leaf_c_id) = .false. + end select c_mask(fnrt_c_id) = grow_fnrt c_mask(sapw_c_id) = grow_sapw c_mask(store_c_id) = grow_store diff --git a/parteh/PRTParamsFATESMod.F90 b/parteh/PRTParamsFATESMod.F90 index 208ff848fb..c017b4e138 100644 --- a/parteh/PRTParamsFATESMod.F90 +++ b/parteh/PRTParamsFATESMod.F90 @@ -1014,12 +1014,12 @@ subroutine PRTCheckParams(is_master) ! Check to make sure the organ ids are valid if this is the ! cnp_flex_allom_hypothesis - if ((hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) .or. & - (hlm_parteh_mode .eq. prt_carbon_allom_hyp) ) then + select case (hlm_parteh_mode) + case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp) do io = 1,norgans if(prt_params%organ_id(io) == repro_organ) then - write(fates_log(),*) 'with flexible cnp or c-only alloc hypothesese' + write(fates_log(),*) 'with flexible cnp or c-only alloc hypotheses' write(fates_log(),*) 'reproductive tissues are a special case' write(fates_log(),*) 'and therefore should not be included in' write(fates_log(),*) 'the parameter file organ list' @@ -1027,16 +1027,16 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) 'Aborting' end if if(prt_params%organ_id(io) == store_organ) then - write(fates_log(),*) 'with flexible cnp or c-only alloc hypothesese' + write(fates_log(),*) 'with flexible cnp or c-only alloc hypotheses' write(fates_log(),*) 'storage is a special case' write(fates_log(),*) 'and therefore should not be included in' write(fates_log(),*) 'the parameter file organ list' write(fates_log(),*) 'fates_prt_organ_id: ',prt_params%organ_id(:) write(fates_log(),*) 'Aborting' end if - + end do - end if + end select pftloop: do ipft = 1,npft @@ -1128,9 +1128,8 @@ subroutine PRTCheckParams(is_master) ! should not be re-translocating mass upon turnover. ! Note to advanced users. Feel free to remove these checks... ! ------------------------------------------------------------------- - - if ((hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) .or. & - (hlm_parteh_mode .eq. prt_carbon_allom_hyp) ) then + select case (hlm_parteh_mode) + case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp) do i = 1,norgans io = prt_params%organ_id(i) @@ -1165,10 +1164,11 @@ subroutine PRTCheckParams(is_master) end if end do - end if - - if (hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) then - + end select + + select case (hlm_parteh_mode) + case (prt_cnp_flex_allom_hyp) + ! Make sure nutrient storage fractions are positive if( prt_params%nitr_store_ratio(ipft) < 0._r8 ) then write(fates_log(),*) 'With parteh allometric CNP hypothesis' @@ -1235,9 +1235,9 @@ subroutine PRTCheckParams(is_master) end if end do - - end if - + + end select + ! Growth respiration ! if (parteh_mode .eq. prt_carbon_allom_hyp) then @@ -1258,8 +1258,8 @@ subroutine PRTCheckParams(is_master) ! end if ! end if - if ((hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) .or. & - (hlm_parteh_mode .eq. prt_carbon_allom_hyp) ) then + select case (hlm_parteh_mode) + case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp) ! The first nitrogen stoichiometry is used in all cases if ( (any(prt_params%nitr_stoich_p1(ipft,:) < 0.0_r8)) .or. & (any(prt_params%nitr_stoich_p1(ipft,:) >= 1.0_r8))) then @@ -1269,9 +1269,10 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) ' Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end select - if(hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) then + select case (hlm_parteh_mode) + case (prt_cnp_flex_allom_hyp) do i = 1,norgans if ( (prt_params%nitr_stoich_p1(ipft,i) < 0._r8) .or. & @@ -1309,7 +1310,7 @@ subroutine PRTCheckParams(is_master) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end select ! Check turnover time-scales From 0bcfd4fc2eb6e5a3b44aff3408fe25e2ee6995b0 Mon Sep 17 00:00:00 2001 From: Gregory Lemieux Date: Wed, 15 Jun 2022 10:08:22 -0700 Subject: [PATCH 6/6] change flux terms to static allocation by maxpft to dynamic allocation by numpft --- main/EDInitMod.F90 | 8 +++++- main/EDTypesMod.F90 | 63 ++++++++++++++++++++++----------------------- 2 files changed, 38 insertions(+), 33 deletions(-) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index b9f551b5be..34f1f49d20 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -127,7 +127,13 @@ subroutine init_site_vars( site_in, bc_in, bc_out ) allocate(site_in%growthflux_fusion(1:nlevsclass,1:numpft)) allocate(site_in%mass_balance(1:num_elements)) allocate(site_in%flux_diags(1:num_elements)) - + + allocate(site_in%term_carbonflux_canopy(1:numpft)) + allocate(site_in%term_carbonflux_ustory(1:numpft)) + allocate(site_in%imort_carbonflux(1:numpft)) + allocate(site_in%fmort_carbonflux_canopy(1:numpft)) + allocate(site_in%fmort_carbonflux_ustory(1:numpft)) + site_in%nlevsoil = bc_in%nlevsoil allocate(site_in%rootfrac_scr(site_in%nlevsoil)) allocate(site_in%zi_soil(0:site_in%nlevsoil)) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 3c9e44b5b8..8475844c83 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -792,42 +792,41 @@ module EDTypesMod ! TERMINATION, RECRUITMENT, DEMOTION, and DISTURBANCE - real(r8), allocatable :: term_nindivs_canopy(:,:) ! number of canopy individuals that were in cohorts which - ! were terminated this timestep, on size x pft - real(r8), allocatable :: term_nindivs_ustory(:,:) ! number of understory individuals that were in cohorts which - ! were terminated this timestep, on size x pft - - real(r8) :: term_carbonflux_canopy(1:maxpft) ! carbon flux from live to dead pools associated - ! with termination mortality, per canopy level - real(r8) :: term_carbonflux_ustory(1:maxpft) ! carbon flux from live to dead pools associated - ! with termination mortality, per canopy level - real(r8) :: demotion_carbonflux ! biomass of demoted individuals from canopy to understory [kgC/ha/day] - real(r8) :: promotion_carbonflux ! biomass of promoted individuals from understory to canopy [kgC/ha/day] - real(r8) :: imort_carbonflux(1:maxpft) ! biomass of individuals killed due to impact mortality per year. [kgC/ha/day] - real(r8) :: fmort_carbonflux_canopy(1:maxpft) ! biomass of canopy indivs killed due to fire per year. [gC/m2/sec] - real(r8) :: fmort_carbonflux_ustory(1:maxpft) ! biomass of understory indivs killed due to fire per year [gC/m2/sec] - real(r8) :: harvest_carbon_flux ! diagnostic site level flux of carbon as harvested plants [kg C / m2 / day] - - real(r8) :: recruitment_rate(1:maxpft) ! number of individuals that were recruited into new cohorts - real(r8), allocatable :: demotion_rate(:) ! rate of individuals demoted from canopy to understory per FATES timestep - - real(r8), allocatable :: promotion_rate(:) ! rate of individuals promoted from understory to canopy per FATES timestep - - real(r8), allocatable :: imort_rate(:,:) ! rate of individuals killed due to impact mortality per year. on size x pft array + real(r8), allocatable :: term_nindivs_canopy(:,:) ! number of canopy individuals that were in cohorts which + ! were terminated this timestep, on size x pft + real(r8), allocatable :: term_nindivs_ustory(:,:) ! number of understory individuals that were in cohorts which + ! were terminated this timestep, on size x pft + + real(r8), allocatable :: term_carbonflux_canopy(:) ! carbon flux from live to dead pools associated + ! with termination mortality, per canopy level + real(r8), allocatable :: term_carbonflux_ustory(:) ! carbon flux from live to dead pools associated + ! with termination mortality, per canopy level + real(r8), allocatable :: imort_carbonflux(:) ! biomass of individuals killed due to impact mortality per year. [kgC/ha/day] + real(r8), allocatable :: fmort_carbonflux_canopy(:) ! biomass of canopy indivs killed due to fire per year. [gC/m2/sec] + real(r8), allocatable :: fmort_carbonflux_ustory(:) ! biomass of understory indivs killed due to fire per year [gC/m2/sec] + + real(r8) :: demotion_carbonflux ! biomass of demoted individuals from canopy to understory [kgC/ha/day] + real(r8) :: promotion_carbonflux ! biomass of promoted individuals from understory to canopy [kgC/ha/day] + real(r8) :: harvest_carbon_flux ! diagnostic site level flux of carbon as harvested plants [kg C / m2 / day] + + real(r8) :: recruitment_rate(1:maxpft) ! number of individuals that were recruited into new cohorts + real(r8), allocatable :: demotion_rate(:) ! rate of individuals demoted from canopy to understory per FATES timestep + real(r8), allocatable :: promotion_rate(:) ! rate of individuals promoted from understory to canopy per FATES timestep + real(r8), allocatable :: imort_rate(:,:) ! rate of individuals killed due to impact mortality per year. on size x pft array - real(r8), allocatable :: fmort_rate_canopy(:,:) ! rate of canopy individuals killed due to fire mortality per year. - ! on size x pft array (1:nlevsclass,1:numpft) - real(r8), allocatable :: fmort_rate_ustory(:,:) ! rate of understory individuals killed due to fire mortality per year. - ! on size x pft array (1:nlevsclass,1:numpft) + real(r8), allocatable :: fmort_rate_canopy(:,:) ! rate of canopy individuals killed due to fire mortality per year. + ! on size x pft array (1:nlevsclass,1:numpft) + real(r8), allocatable :: fmort_rate_ustory(:,:) ! rate of understory individuals killed due to fire mortality per year. + ! on size x pft array (1:nlevsclass,1:numpft) - real(r8), allocatable :: fmort_rate_cambial(:,:) ! rate of individuals killed due to fire mortality - ! from cambial damage per year. on size x pft array - real(r8), allocatable :: fmort_rate_crown(:,:) ! rate of individuals killed due to fire mortality - ! from crown damage per year. on size x pft array + real(r8), allocatable :: fmort_rate_cambial(:,:) ! rate of individuals killed due to fire mortality + ! from cambial damage per year. on size x pft array + real(r8), allocatable :: fmort_rate_crown(:,:) ! rate of individuals killed due to fire mortality + ! from crown damage per year. on size x pft array - real(r8), allocatable :: growthflux_fusion(:,:) ! rate of individuals moving into a given size class bin - ! due to fusion in a given day. on size x pft array + real(r8), allocatable :: growthflux_fusion(:,:) ! rate of individuals moving into a given size class bin + ! due to fusion in a given day. on size x pft array