diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 4318ee469f..ec552f542c 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -70,6 +70,7 @@ module EDCohortDynamicsMod use FatesAllometryMod , only : tree_lai, tree_sai use FatesAllometryMod , only : set_root_fraction use PRTGenericMod, only : prt_carbon_allom_hyp + use PRTGenericMod, only : prt_csimpler_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes use PRTGenericMod, only : all_carbon_elements @@ -85,14 +86,21 @@ module EDCohortDynamicsMod use PRTGenericMod, only : SetState use PRTAllometricCarbonMod, only : callom_prt_vartypes + use PRTAllometricCarbonMod, only : csimpler_allom_prt_vartypes use PRTAllometricCarbonMod, only : ac_bc_inout_id_netdc use PRTAllometricCarbonMod, only : ac_bc_in_id_pft use PRTAllometricCarbonMod, only : ac_bc_in_id_ctrim use PRTAllometricCarbonMod, only : ac_bc_inout_id_dbh use PRTAllometricCarbonMod, only : ac_bc_in_id_lstat + use PRTAllometricCarbonMod, only : ac_bc_in_id_efleaf + use PRTAllometricCarbonMod, only : ac_bc_in_id_effnrt + use PRTAllometricCarbonMod, only : ac_bc_in_id_efstem use PRTAllometricCNPMod, only : cnp_allom_prt_vartypes use PRTAllometricCNPMod, only : acnp_bc_in_id_pft, acnp_bc_in_id_ctrim use PRTAllometricCNPMod, only : acnp_bc_in_id_lstat, acnp_bc_inout_id_dbh + use PRTAllometricCNPMod, only : acnp_bc_in_id_efleaf + use PRTAllometricCNPMod, only : acnp_bc_in_id_effnrt + use PRTAllometricCNPMod, only : acnp_bc_in_id_efstem use PRTAllometricCNPMod, only : acnp_bc_inout_id_rmaint_def, acnp_bc_in_id_netdc use PRTAllometricCNPMod, only : acnp_bc_in_id_netdnh4, acnp_bc_in_id_netdno3, acnp_bc_in_id_netdp use PRTAllometricCNPMod, only : acnp_bc_out_id_cefflux, acnp_bc_out_id_nefflux @@ -147,8 +155,9 @@ module EDCohortDynamicsMod subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & - prt, laimemory, sapwmemory, structmemory, & - status, recruitstatus,ctrim, carea, clayer, spread, bc_in) + prt, leafmemory, fnrtmemory, sapwmemory, structmemory, & + elongf_leaf, elongf_fnrt, elongf_stem, status, & + recruitstatus, ctrim, carea, clayer, spread, bc_in) ! ! !DESCRIPTION: ! create new cohort @@ -171,7 +180,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & integer, intent(in) :: clayer ! canopy status of cohort ! (1 = canopy, 2 = understorey, etc.) integer, intent(in) :: status ! growth status of plant - ! (2 = leaves on , 1 = leaves off) + ! (1 = abscissing, 2 = flushing) integer, intent(in) :: recruitstatus ! recruit status of plant ! (1 = recruitment , 0 = other) real(r8), intent(in) :: nn ! number of individuals in cohort @@ -179,9 +188,17 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & real(r8), intent(in) :: hite ! height: meters real(r8), intent(in) :: coage ! cohort age in years real(r8), intent(in) :: dbh ! dbh: cm + real(r8), intent(in) :: elongf_leaf ! leaf elongation factor (fraction) + real(r8), intent(in) :: elongf_fnrt ! fine-root "elongation factor" (fraction) + real(r8), intent(in) :: elongf_stem ! stem "elongation factor" (fraction) + ! For all elongation factors: + ! 0 means fully abscissed + ! 1 means fully flushed class(prt_vartypes),target :: prt ! The allocated PARTEH ! object - real(r8), intent(in) :: laimemory ! target leaf biomass- set from + real(r8), intent(in) :: leafmemory ! target leaf biomass- set from + ! previous year: kGC per indiv + real(r8), intent(in) :: fnrtmemory ! target fine root biomass- set from ! previous year: kGC per indiv real(r8), intent(in) :: sapwmemory ! target sapwood biomass- set from ! previous year: kGC per indiv @@ -230,6 +247,9 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & new_cohort%pft = pft new_cohort%status_coh = status + new_cohort%efleaf_coh = elongf_leaf + new_cohort%effnrt_coh = elongf_fnrt + new_cohort%efstem_coh = elongf_stem new_cohort%n = nn new_cohort%hite = hite new_cohort%dbh = dbh @@ -237,7 +257,8 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & new_cohort%canopy_trim = ctrim new_cohort%canopy_layer = clayer new_cohort%canopy_layer_yesterday = real(clayer, r8) - new_cohort%laimemory = laimemory + new_cohort%leafmemory = leafmemory + new_cohort%fnrtmemory = fnrtmemory new_cohort%sapwmemory = sapwmemory new_cohort%structmemory = structmemory @@ -400,7 +421,7 @@ subroutine InitPRTBoundaryConditions(new_cohort) select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp) ! Register boundary conditions for the Carbon Only Allometric Hypothesis @@ -409,12 +430,18 @@ subroutine InitPRTBoundaryConditions(new_cohort) call new_cohort%prt%RegisterBCIn(ac_bc_in_id_pft,bc_ival = new_cohort%pft) call new_cohort%prt%RegisterBCIn(ac_bc_in_id_ctrim,bc_rval = new_cohort%canopy_trim) call new_cohort%prt%RegisterBCIn(ac_bc_in_id_lstat,bc_ival = new_cohort%status_coh) + call new_cohort%prt%RegisterBCIn(ac_bc_in_id_efleaf,bc_rval = new_cohort%efleaf_coh) + call new_cohort%prt%RegisterBCIn(ac_bc_in_id_effnrt,bc_rval = new_cohort%effnrt_coh) + call new_cohort%prt%RegisterBCIn(ac_bc_in_id_efstem,bc_rval = new_cohort%efstem_coh) case (prt_cnp_flex_allom_hyp) call new_cohort%prt%RegisterBCIn(acnp_bc_in_id_pft,bc_ival = new_cohort%pft) call new_cohort%prt%RegisterBCIn(acnp_bc_in_id_ctrim,bc_rval = new_cohort%canopy_trim) call new_cohort%prt%RegisterBCIn(acnp_bc_in_id_lstat,bc_ival = new_cohort%status_coh) + call new_cohort%prt%RegisterBCIn(acnp_bc_in_id_efleaf,bc_rval = new_cohort%efleaf_coh) + call new_cohort%prt%RegisterBCIn(acnp_bc_in_id_effnrt,bc_rval = new_cohort%effnrt_coh) + call new_cohort%prt%RegisterBCIn(acnp_bc_in_id_efstem,bc_rval = new_cohort%efstem_coh) call new_cohort%prt%RegisterBCIn(acnp_bc_in_id_netdc, bc_rval = new_cohort%npp_acc) call new_cohort%prt%RegisterBCIn(acnp_bc_in_id_netdnh4, bc_rval = new_cohort%daily_nh4_uptake) call new_cohort%prt%RegisterBCIn(acnp_bc_in_id_netdno3, bc_rval = new_cohort%daily_no3_uptake) @@ -463,10 +490,16 @@ subroutine InitPRTObject(prt) ! Potential Extended types class(callom_prt_vartypes), pointer :: c_allom_prt + class(csimpler_allom_prt_vartypes), pointer :: csimpler_allom_prt class(cnp_allom_prt_vartypes), pointer :: cnp_allom_prt select case(hlm_parteh_mode) + case (prt_csimpler_allom_hyp) + + allocate(csimpler_allom_prt) + prt => csimpler_allom_prt + case (prt_carbon_allom_hyp) allocate(c_allom_prt) @@ -530,7 +563,10 @@ subroutine nan_cohort(cc_p) currentCohort%canopy_layer = fates_unset_int ! canopy status of cohort (1 = canopy, 2 = understorey, etc.) currentCohort%canopy_layer_yesterday = nan ! recent canopy status of cohort (1 = canopy, 2 = understorey, etc.) currentCohort%NV = fates_unset_int ! Number of leaf layers: - - currentCohort%status_coh = fates_unset_int ! growth status of plant (2 = leaves on , 1 = leaves off) + currentCohort%status_coh = fates_unset_int ! growth status of plant (2 = flushing , 1 = fully abscissed, 3 = abscissing) + currentCohort%efleaf_coh = nan ! leaf elongation factor (fraction from 0 (fully abscissed) to 1 (fully flushed) + currentCohort%effnrt_coh = nan ! fine-root "elongation factor" (fraction from 0 (fully abscissed) to 1 (fully flushed) + currentCohort%efstem_coh = nan ! stem "elongation factor" (fraction from 0 (fully abscissed) to 1 (fully flushed) currentCohort%size_class = fates_unset_int ! size class index currentCohort%size_class_lasttimestep = fates_unset_int ! size class index currentCohort%size_by_pft_class = fates_unset_int ! size by pft classification index @@ -540,8 +576,9 @@ subroutine nan_cohort(cc_p) currentCohort%n = nan ! number of individuals in cohort per 'area' (10000m2 default) currentCohort%dbh = nan ! 'diameter at breast height' in cm currentCohort%coage = nan ! age of the cohort in years - currentCohort%hite = nan ! height: meters - currentCohort%laimemory = nan ! target leaf biomass- set from previous year: kGC per indiv + currentCohort%hite = nan ! height: meters + currentCohort%leafmemory = nan ! target leaf biomass- set from previous year: kGC per indiv + currentCohort%fnrtmemory = nan ! target fine-root biomass- set from previous year: kGC per indiv currentCohort%sapwmemory = nan ! target sapwood biomass- set from previous year: kGC per indiv currentCohort%structmemory = nan ! target structural biomass- set from previous year: kGC per indiv currentCohort%lai = nan ! leaf area index of cohort m2/m2 @@ -647,6 +684,9 @@ subroutine zero_cohort(cc_p) currentCohort%NV = 0 currentCohort%status_coh = 0 + currentCohort%efleaf_coh = 0._r8 + currentCohort%effnrt_coh = 0._r8 + currentCohort%efstem_coh = 0._r8 currentCohort%rdark = 0._r8 currentCohort%resp_m = 0._r8 currentCohort%resp_m_def = 0._r8 @@ -765,7 +805,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ if (currentcohort%n < min_n_safemath .and. level == 1) then terminate = itrue if ( debug ) then - write(fates_log(),*) 'terminating cohorts 0',currentCohort%n/currentPatch%area,currentCohort%dbh,call_index + write(fates_log(),*) 'terminating cohorts 0',currentCohort%n/currentPatch%area,currentCohort%dbh,currentCohort%pft,call_index endif endif @@ -778,7 +818,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ (currentCohort%dbh < 0.00001_r8 .and. store_c < 0._r8) ) then terminate = itrue if ( debug ) then - write(fates_log(),*) 'terminating cohorts 1',currentCohort%n/currentPatch%area,currentCohort%dbh,call_index + write(fates_log(),*) 'terminating cohorts 1',currentCohort%n/currentPatch%area,currentCohort%dbh,currentCohort%pft,call_index endif endif @@ -786,7 +826,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ if (currentCohort%canopy_layer > nclmax ) then terminate = itrue if ( debug ) then - write(fates_log(),*) 'terminating cohorts 2', currentCohort%canopy_layer,call_index + write(fates_log(),*) 'terminating cohorts 2', currentCohort%canopy_layer,currentCohort%pft,call_index endif endif @@ -796,7 +836,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ terminate = itrue if ( debug ) then write(fates_log(),*) 'terminating cohorts 3', & - sapw_c,leaf_c,fnrt_c,store_c,call_index + sapw_c,leaf_c,fnrt_c,store_c,currentCohort%pft,call_index endif endif @@ -805,7 +845,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ terminate = itrue if ( debug ) then write(fates_log(),*) 'terminating cohorts 4', & - struct_c,sapw_c,leaf_c,fnrt_c,store_c,call_index + struct_c,sapw_c,leaf_c,fnrt_c,store_c,currentCohort%pft,call_index endif endif @@ -1196,7 +1236,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) write(fates_log(),*) 'Cohort I, Cohort II' write(fates_log(),*) 'n:',currentCohort%n,nextc%n write(fates_log(),*) 'isnew:',currentCohort%isnew,nextc%isnew - write(fates_log(),*) 'laimemory:',currentCohort%laimemory,nextc%laimemory + write(fates_log(),*) 'leafmemory:',currentCohort%leafmemory,nextc%leafmemory + write(fates_log(),*) 'fnrtmemory:',currentCohort%fnrtmemory,nextc%fnrtmemory write(fates_log(),*) 'hite:',currentCohort%hite,nextc%hite write(fates_log(),*) 'coage:',currentCohort%coage,nextc%coage write(fates_log(),*) 'dbh:',currentCohort%dbh,nextc%dbh @@ -1234,8 +1275,11 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! ----------------------------------------------------------------- call UpdateCohortBioPhysRates(currentCohort) - currentCohort%laimemory = (currentCohort%n*currentCohort%laimemory & - + nextc%n*nextc%laimemory)/newn + currentCohort%leafmemory = (currentCohort%n*currentCohort%leafmemory & + + nextc%n*nextc%leafmemory)/newn + + currentCohort%fnrtmemory = (currentCohort%n*currentCohort%fnrtmemory & + + nextc%n*nextc%fnrtmemory)/newn currentCohort%sapwmemory = (currentCohort%n*currentCohort%sapwmemory & + nextc%n*nextc%sapwmemory)/newn @@ -1828,7 +1872,8 @@ subroutine copy_cohort( currentCohort,copyc ) n%dbh = o%dbh n%coage = o%coage n%hite = o%hite - n%laimemory = o%laimemory + n%leafmemory = o%leafmemory + n%fnrtmemory = o%fnrtmemory n%sapwmemory = o%sapwmemory n%structmemory = o%structmemory n%lai = o%lai @@ -1839,6 +1884,9 @@ subroutine copy_cohort( currentCohort,copyc ) n%canopy_layer_yesterday = o%canopy_layer_yesterday n%nv = o%nv n%status_coh = o%status_coh + n%efleaf_coh = o%efleaf_coh + n%effnrt_coh = o%effnrt_coh + n%efstem_coh = o%efstem_coh n%canopy_trim = o%canopy_trim n%excl_weight = o%excl_weight n%prom_weight = o%prom_weight diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 537e74824d..bc898236ee 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -84,7 +84,8 @@ module EDPatchDynamicsMod use PRTGenericMod, only : struct_organ use PRTLossFluxesMod, only : PRTBurnLosses use FatesInterfaceTypesMod, only : hlm_parteh_mode - use PRTGenericMod, only : prt_carbon_allom_hyp + use PRTGenericMod, only : prt_csimpler_allom_hyp + use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use SFParamsMod, only : SF_VAL_CWD_FRAC use EDParamsMod, only : logging_event_code diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 5d66f56d39..96adb14cd1 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -51,6 +51,7 @@ module EDPhysiologyMod use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : leaves_on use EDTypesMod , only : leaves_off + use EDTypesMod , only : leaves_pshed use EDTypesMod , only : min_n_safemath use PRTGenericMod , only : num_elements use PRTGenericMod , only : element_list @@ -63,6 +64,9 @@ module EDPhysiologyMod use EDTypesMod , only : phen_dstat_moistoff use EDTypesMod , only : phen_dstat_moiston use EDTypesMod , only : phen_dstat_timeon + use EDTypesMod , only : phen_dstat_pshed + use EDTypesMod , only : drgt_phen_model_smoist + use EDTypesMod , only : drgt_phen_model_gradual use EDTypesMod , only : init_recruit_trim use shr_log_mod , only : errMsg => shr_log_errMsg use FatesGlobals , only : fates_log @@ -87,6 +91,7 @@ module EDPhysiologyMod use FatesAllometryMod , only : carea_allom use FatesAllometryMod , only : CheckIntegratedAllometries use FatesAllometryMod, only : set_root_fraction + use PRTGenericMod, only : prt_csimpler_allom_hyp use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes @@ -236,7 +241,7 @@ subroutine PreDisturbanceLitterFluxes( currentSite, currentPatch, bc_in ) ! Calculate seed germination rate, the status flags prevent ! germination from occuring when the site is in a drought ! (for drought deciduous) or too cold (for cold deciduous) - call SeedGermination(litt, currentSite%cstatus, currentSite%dstatus) + call SeedGermination(litt, currentSite%cstatus, currentSite%dstatus(1:numpft)) ! Send fluxes from newly created litter into the litter pools ! This litter flux is from non-disturbance inducing mortality, as well @@ -424,9 +429,9 @@ subroutine trim_canopy( currentSite ) real(r8) :: work(workmax) ! work array real(r8) :: initial_trim ! Initial trim - real(r8) :: optimum_trim ! Optimum trim value - real(r8) :: initial_laimem ! Initial laimemory - real(r8) :: optimum_laimem ! Optimum laimemory + real(r8) :: optimum_trim ! Optimum trim value + real(r8) :: initial_leafmem ! Initial leaf memory + real(r8) :: optimum_leafmem ! Optimum leaf memory !---------------------------------------------------------------------- @@ -446,15 +451,15 @@ subroutine trim_canopy( currentSite ) currentCohort => currentPatch%tallest do while (associated(currentCohort)) - ! Save off the incoming trim and laimemory + ! Save off the incoming trim and leafmemory initial_trim = currentCohort%canopy_trim - initial_laimem = currentCohort%laimemory + initial_leafmem = currentCohort%leafmemory ! Add debug diagnstic output to determine which cohort if (debug) then write(fates_log(),*) 'Current cohort:', icohort write(fates_log(),*) 'Starting canopy trim:', initial_trim - write(fates_log(),*) 'Starting laimemory:', currentCohort%laimemory + write(fates_log(),*) 'Starting leafmemory:', currentCohort%leafmemory endif trimmed = .false. @@ -600,10 +605,6 @@ subroutine trim_canopy( currentSite ) if (currentCohort%hite > EDPftvarcon_inst%hgt_min(ipft)) then currentCohort%canopy_trim = currentCohort%canopy_trim - & EDPftvarcon_inst%trim_inc(ipft) - if (prt_params%evergreen(ipft) /= 1)then - currentCohort%laimemory = currentCohort%laimemory * & - (1.0_r8 - EDPftvarcon_inst%trim_inc(ipft)) - endif trimmed = .true. @@ -649,15 +650,15 @@ subroutine trim_canopy( currentSite ) ! optimum_trim = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_trim - optimum_laimem = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_laimem + optimum_leafmem = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_leafmem ! Determine if the optimum trim value makes sense. The smallest cohorts tend to have unrealistic fits. if (optimum_trim > 0. .and. optimum_trim < 1.) then currentCohort%canopy_trim = optimum_trim - ! If the cohort pft is not evergreen we reduce the laimemory as well + ! If the cohort pft is not evergreen we reduce the leafmemory as well if (prt_params%evergreen(ipft) /= 1) then - currentCohort%laimemory = optimum_laimem + currentCohort%leafmemory = optimum_leafmem endif trimmed = .true. @@ -696,7 +697,11 @@ subroutine phenology( currentSite, bc_in ) ! ! !USES: use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - use EDParamsMod, only : ED_val_phen_drought_threshold, ED_val_phen_doff_time + use EDParamsMod, only : phen_drought_model + use EDParamsMod, only : ED_val_phen_doff_time + use EDParamsMod, only : ED_val_phen_drought_threshold + use EDParamsMod, only : ED_val_phen_moist_threshold + use EDParamsMod, only : ED_val_phen_doff_time use EDParamsMod, only : ED_val_phen_a, ED_val_phen_b, ED_val_phen_c, ED_val_phen_chiltemp use EDParamsMod, only : ED_val_phen_mindayson, ED_val_phen_ncolddayslim, ED_val_phen_coldtemp @@ -714,46 +719,58 @@ subroutine phenology( currentSite, bc_in ) integer :: ncolddays ! no days underneath the threshold for leaf drop integer :: i_wmem ! Loop counter for water mem days integer :: i_tmem ! Loop counter for veg temp mem days - integer :: dayssincedleafon ! Days since drought-decid leaf-on started - integer :: dayssincedleafoff ! Days since drought-decid leaf-off started - integer :: dayssincecleafon ! Days since cold-decid leaf-on started - integer :: dayssincecleafoff ! Days since cold-decid leaf-off started - real(r8) :: mean_10day_liqvol ! mean liquid volume (m3/m3) over last 10 days + integer :: ft ! plant functional type index + real(r8) :: mean_10day_liqvol ! mean soil liquid water volume over last 10 days + real(r8) :: mean_10day_smp ! mean soil matric potential over last 10 days real(r8) :: leaf_c ! leaf carbon [kg] real(r8) :: fnrt_c ! fineroot carbon [kg] real(r8) :: sapw_c ! sapwood carbon [kg] real(r8) :: store_c ! storage carbon [kg] real(r8) :: struct_c ! structure carbon [kg] real(r8) :: gdd_threshold ! GDD accumulation function, - integer :: ilayer_swater ! Layer index for soil water - ! which also depends on chilling days. + real(r8) :: rootfrac_notop ! Total rooting fraction excluding the top soil layer integer :: ncdstart ! beginning of counting period for chilling degree days. integer :: gddstart ! beginning of counting period for growing degree days. - real(r8) :: temp_in_C ! daily averaged temperature in celcius + integer :: nlevroot ! Number of rooting levels to consider + real(r8) :: temp_in_C ! daily averaged temperature in celsius + real(r8) :: elongf_prev ! Elongation factor from previous time + real(r8) :: elongf_1st ! First guess for elongation factor + + logical :: smoist_below_threshold ! Is soil moisture below threshold? + logical :: recent_flush ! Last full flushing event is still very recent. + logical :: recent_abscission ! Last abscission event is still very recent. + logical :: exceed_min_on_period ! Have leaves been flushed for a minimum period of time? + logical :: exceed_min_off_period ! Have leaves been off for a minimum period of time? + logical :: prolonged_on_period ! Has leaves been flushed for too long? + logical :: prolonged_off_period ! Have leaves been abscissed for too long? + logical :: last_flush_long_ago ! Has it been a very long time since last flushing? integer, parameter :: canopy_leaf_lifespan = 365 ! Maximum lifespan of drought decid leaves - integer, parameter :: min_daysoff_dforcedflush = 30 ! THis is the number of days that must had elapsed - ! since leaves had dropped, in order to forcably - ! flush leaves again. This does not impact flushing - ! due to real moisture constraints, and will prevent - ! drought deciduous in perennially wet environments - ! that have been forced to drop their leaves, from - ! flushing them back immediately. - - real(r8),parameter :: dphen_soil_depth = 0.1 ! Use liquid soil water that is - ! closest to this depth [m] + integer, parameter :: min_daysoff_dforcedflush = 30 ! This is the number of days that must had elapsed + ! since leaves had dropped, in order to forcably + ! flush leaves again. This does not impact flushing + ! due to real moisture constraints, and will prevent + ! drought deciduous in perennially wet environments + ! that have been forced to drop their leaves, from + ! flushing them back immediately. + integer, parameter :: dd_offon_toler = 30 ! When flushing or shedding leaves, we check that + ! the dates are near last year's dates. This controls + ! the tolerance for deviating from last year. + + real(r8), parameter :: elongf_min = 0.05_r8 ! Minimum elongation factor. If elongation factor + ! reaches or falls below elongf_min, we assume + ! complete abscission. This avoids carrying out + ! a residual amount of leaves, which may create + ! computational problems. The current threshold + ! is the same used in ED-2.2. ! This is the integer model day. The first day of the simulation is 1, and it ! continues monotonically, indefinitely model_day_int = nint(hlm_model_day) - ! Use the following layer index to calculate drought conditions - ilayer_swater = minloc(abs(bc_in%z_sisl(:)-dphen_soil_depth),dim=1) - - - ! Parameter of drought decid leaf loss in mm in top layer...FIX(RF,032414) + ! Parameter of drought decid leaf loss in mm in top layer...FIX(RF,032414) ! - this is arbitrary and poorly understood. Needs work. ED_ !Parameters: defaults from Botta et al. 2000 GCB,6 709-725 !Parameters, default from from SDGVM model of senesence @@ -845,15 +862,15 @@ subroutine phenology( currentSite, bc_in ) ! not had occured yet, so set it to last year to get things rolling if (model_day_int < currentSite%cleafoffdate) then - dayssincecleafoff = model_day_int - (currentSite%cleafoffdate - 365) + currentSite%cndaysleafoff = model_day_int - (currentSite%cleafoffdate - 365) else - dayssincecleafoff = model_day_int - currentSite%cleafoffdate + currentSite%cndaysleafoff = model_day_int - currentSite%cleafoffdate end if if (model_day_int < currentSite%cleafondate) then - dayssincecleafon = model_day_int - (currentSite%cleafondate-365) + currentSite%cndaysleafon = model_day_int - (currentSite%cleafondate-365) else - dayssincecleafon = model_day_int - currentSite%cleafondate + currentSite%cndaysleafon = model_day_int - currentSite%cleafondate end if @@ -869,12 +886,12 @@ subroutine phenology( currentSite, bc_in ) if ( (currentSite%cstatus == phen_cstat_iscold .or. & currentSite%cstatus == phen_cstat_nevercold) .and. & (currentSite%grow_deg_days > gdd_threshold) .and. & - (dayssincecleafoff > ED_val_phen_mindayson) .and. & + (currentSite%cndaysleafoff > ED_val_phen_mindayson) .and. & (currentSite%nchilldays >= 1)) then currentSite%cstatus = phen_cstat_notcold ! Set to not-cold status (leaves can come on) - currentSite%cleafondate = model_day_int - dayssincecleafon = 0 - currentSite%grow_deg_days = 0._r8 ! zero GDD for the rest of the year until counting season begins. + currentSite%cleafondate = model_day_int + currentSite%cndaysleafon = 0 + currentSite%grow_deg_days = 0._r8 ! zero GDD for the rest of the year until counting season begins. if ( debug ) write(fates_log(),*) 'leaves on' endif !GDD @@ -892,7 +909,7 @@ subroutine phenology( currentSite, bc_in ) if ( (currentSite%cstatus == phen_cstat_notcold) .and. & (model_day_int > num_vegtemp_mem) .and. & (ncolddays > ED_val_phen_ncolddayslim) .and. & - (dayssincecleafon > ED_val_phen_mindayson) )then + (currentSite%cndaysleafon > ED_val_phen_mindayson) )then currentSite%grow_deg_days = 0._r8 ! The equations for Botta et al ! are for calculations of @@ -900,7 +917,8 @@ subroutine phenology( currentSite, bc_in ) ! clear this value, it will cause ! leaves to flush later in the year currentSite%cstatus = phen_cstat_iscold ! alter status of site to 'leaves off' - currentSite%cleafoffdate = model_day_int ! record leaf off date + currentSite%cleafoffdate = model_day_int ! record leaf off date + currentSite%cndaysleafoff = 0 if ( debug ) write(fates_log(),*) 'leaves off' endif @@ -912,14 +930,15 @@ subroutine phenology( currentSite, bc_in ) ! plants from re-emerging in areas without at least some cold days if( (currentSite%cstatus == phen_cstat_notcold) .and. & - (dayssincecleafoff > 400)) then ! remove leaves after a whole year - ! when there is no 'off' period. + (currentSite%cndaysleafoff > 400)) then ! remove leaves after a whole year + ! when there is no 'off' period. currentSite%grow_deg_days = 0._r8 currentSite%cstatus = phen_cstat_nevercold ! alter status of site to imply that this - ! site is never really cold enough - ! for cold deciduous - currentSite%cleafoffdate = model_day_int ! record leaf off date + ! site is never really cold enough + ! for cold deciduous + currentSite%cleafoffdate = model_day_int ! record leaf off date + currentSite%cndaysleafoff = 0 if ( debug ) write(fates_log(),*) 'leaves off' endif @@ -937,129 +956,262 @@ subroutine phenology( currentSite, bc_in ) ! to come into equlibirium. ! E*: If the soil is always wet, the leaves come on at the beginning of the window, and then ! last for their lifespan. - ! ISSUES - ! 1. It's not clear what water content we should track. Here we are tracking the top layer, - ! but we probably should track something like BTRAN, but BTRAN is defined for each PFT, - ! and there could potentially be more than one stress-dec PFT.... ? - ! 2. In the beginning, the window is set at an arbitrary time of the year, so the leaves - ! might come on in the dry season, using up stored reserves - ! for the stress-dec plants, and potentially killing them. To get around this, - ! we need to read in the 'leaf on' date from some kind of start-up file - ! but we would need that to happen for every resolution, etc. - ! 3. Will this methodology properly kill off the stress-dec trees where there is no - ! water stress? What about where the wet period coincides with the warm period? - ! We would just get them overlapping with the cold-dec trees, even though that isn't appropriate - ! Why don't the drought deciduous trees grow in the North? - ! Is cold decidousness maybe even the same as drought deciduosness there (and so does this - ! distinction actually matter??).... - - ! Accumulate surface water memory of last 10 days. - ! Liquid volume in ground layer (m3/m3) - do i_wmem = 1,numWaterMem-1 !shift memory along one - currentSite%water_memory(numWaterMem+1-i_wmem) = currentSite%water_memory(numWaterMem-i_wmem) - enddo - currentSite%water_memory(1) = bc_in%h2o_liqvol_sl(ilayer_swater) - - ! Calculate the mean water content over the last 10 days (m3/m3) - mean_10day_liqvol = sum(currentSite%water_memory(1:numWaterMem))/real(numWaterMem,r8) - - ! In drought phenology, we often need to force the leaves to stay - ! on or off as moisture fluctuates... - - ! Calculate days since leaves have come off, but make a provision - ! for the first year of simulation, we have to assume a leaf drop - ! date to start, so if that is in the future, set it to last year - - if (model_day_int < currentSite%dleafoffdate) then - dayssincedleafoff = model_day_int - (currentSite%dleafoffdate-365) - else - dayssincedleafoff = model_day_int - currentSite%dleafoffdate - endif - - ! the leaves are on. How long have they been on? - if (model_day_int < currentSite%dleafondate) then - dayssincedleafon = model_day_int - (currentSite%dleafondate-365) - else - dayssincedleafon = model_day_int - currentSite%dleafondate - endif - ! LEAF ON: DROUGHT DECIDUOUS WETNESS - ! Here, we used a window of oppurtunity to determine if we are - ! close to the time when then leaves came on last year - - ! Has it been ... - ! a) a year, plus or minus 1 month since we last had leaf-on? - ! b) Has there also been at least a nominaly short amount of "leaf-off" - ! c) is the model day at least > 10 (let soil water spin-up) - ! Note that cold-starts begin in the "leaf-on" - ! status - if ( (currentSite%dstatus == phen_dstat_timeoff .or. & - currentSite%dstatus == phen_dstat_moistoff) .and. & - (model_day_int > numWaterMem) .and. & - (dayssincedleafon >= 365-30 .and. dayssincedleafon <= 365+30 ) .and. & - (dayssincedleafoff > ED_val_phen_doff_time) ) then - - ! If leaves are off, and have been off for at least a few days - ! and the time is consistent with the correct - ! time window... test if the moisture conditions allow for leaf-on - - if ( mean_10day_liqvol >= ED_val_phen_drought_threshold ) then - currentSite%dstatus = phen_dstat_moiston ! set status to leaf-on - currentSite%dleafondate = model_day_int ! save the model day we start flushing - dayssincedleafon = 0 - endif - endif + ! Add PFT look to account for different PFT rooting depth profiles. + pft_drgt_loop: do ft=1,numpft - ! LEAF ON: DROUGHT DECIDUOUS TIME EXCEEDANCE - ! If we still haven't done budburst by end of window, then force it - ! If the status is "phen_dstat_moistoff", it means this site currently has - ! leaves off due to actual moisture limitations. - ! So we trigger bud-burst at the end of the month since - ! last year's bud-burst. If this is imposed, then we set the new - ! status to indicate bud-burst was forced by timing + ! Update soil moisture information memory (we always track the last 10 days) + do i_wmem = numWaterMem,2,-1 !shift memory to previous day, to make room for current day + currentSite%liqvol_memory(i_wmem,ft) = currentSite%liqvol_memory(i_wmem-1,ft) + currentSite%smp_memory (i_wmem,ft) = currentSite%smp_memory (i_wmem-1,ft) + end do - if( currentSite%dstatus == phen_dstat_moistoff ) then - if ( dayssincedleafon > 365+30 ) then - currentSite%dstatus = phen_dstat_timeon ! force budburst! - currentSite%dleafondate = model_day_int ! record leaf on date - dayssincedleafon = 0 + ! Find the rooting depth distribution for PFT + call set_root_fraction( currentSite%rootfrac_scr, ft, currentSite%zi_soil, & + bc_in%max_rooting_depth_index_col ) + nlevroot = max(2,min(ubound(currentSite%zi_soil,1),bc_in%max_rooting_depth_index_col)) + + ! The top most layer is typically very thin (~ 2cm) and dries rather quickly. Despite + ! being thin, it can have a non-negligible rooting fraction (e.g., using + ! exponential_2p_root_profile with default parameters make the top layer to contain + ! about 7% of the total fine root density). To avoid overestimating dryness, we + ! ignore the top layer when calculating the memory. + rootfrac_notop = sum(currentSite%rootfrac_scr(2:nlevroot)) + if ( rootfrac_notop <= nearzero ) then + ! Unlikely, but just in case all roots are in the first layer, we use the second + ! layer the second layer (to avoid FPE issues). + currentSite%rootfrac_scr(2) = 1.0_r8 + rootfrac_notop = 1.0_r8 end if - end if - ! But if leaves are off due to time, then we enforce - ! a longer cool-down (because this is a perrenially wet system) + ! Set the memory to be the weighted average of the soil properties, using the + ! root fraction of each layer (except the topmost one) as the weighting factor. + currentSite%liqvol_memory(1,ft) = sum( bc_in%h2o_liqvol_sl (2:nlevroot) * & + currentSite%rootfrac_scr(2:nlevroot) ) / rootfrac_notop + currentSite%smp_memory (1,ft) = sum( bc_in%smp_sl (2:nlevroot) * & + currentSite%rootfrac_scr(2:nlevroot) ) / rootfrac_notop + + ! Calculate the mean soil moisture ( liquid volume (m3/m3) and matric potential (mm)) + ! over the last 10 days + mean_10day_liqvol = sum(currentSite%liqvol_memory(1:numWaterMem,ft))/real(numWaterMem,r8) + mean_10day_smp = sum(currentSite%smp_memory (1:numWaterMem,ft))/real(numWaterMem,r8) + + ! Compare the moisture with the threshold. + if ( ED_val_phen_drought_threshold >= 0. ) then + ! Liquid volume in reference layer (m3/m3) + smoist_below_threshold = mean_10day_liqvol < ED_val_phen_drought_threshold + else + ! Soil matric potential in reference layer (mm) + smoist_below_threshold = mean_10day_smp < ED_val_phen_drought_threshold + end if - if(currentSite%dstatus == phen_dstat_timeoff ) then - if (dayssincedleafoff > min_daysoff_dforcedflush) then - currentSite%dstatus = phen_dstat_timeon ! force budburst! - currentSite%dleafondate = model_day_int ! record leaf on date - dayssincedleafon = 0 + ! Calculate days since last flushing and shedding event, but make a provision + ! for the first year of simulation, we have to assume leaf drop / leaf flush + ! dates to start, so if that is in the future, set it to last year + if (model_day_int < currentSite%dleafoffdate(ft)) then + currentSite%dndaysleafoff(ft) = model_day_int - (currentSite%dleafoffdate(ft)-365) + else + currentSite%dndaysleafoff(ft) = model_day_int - currentSite%dleafoffdate(ft) + end if + if (model_day_int < currentSite%dleafondate(ft)) then + currentSite%dndaysleafon(ft) = model_day_int - (currentSite%dleafondate(ft)-365) + else + currentSite%dndaysleafon(ft) = model_day_int - currentSite%dleafondate(ft) end if - end if - ! LEAF OFF: DROUGHT DECIDUOUS LIFESPAN - if the leaf gets to - ! the end of its useful life. A*, E* - ! i.e. Are the leaves rouhgly at the end of their lives? - if ( (currentSite%dstatus == phen_dstat_moiston .or. & - currentSite%dstatus == phen_dstat_timeon ) .and. & - (dayssincedleafon > canopy_leaf_lifespan) )then - currentSite%dstatus = phen_dstat_timeoff !alter status of site to 'leaves off' - currentSite%dleafoffdate = model_day_int !record leaf on date - endif - - ! LEAF OFF: DROUGHT DECIDUOUS DRYNESS - if the soil gets too dry, - ! and the leaves have already been on a while... + ! Elongation factor from the previous step. + elongf_prev = currentSite%elong_factor(ft) + + + + + + ! Find elongation factor by comparing the moisture with the thresholds. + case_drought_phen: select case (phen_drought_model) + case (drgt_phen_model_smoist) + !------ + ! Default drought deciduous phenology. + !------ + + ! In drought phenology, we often need to force the leaves to stay + ! on or off as moisture fluctuates... + + + ! Save some conditions in logical variables to simplify code below + exceed_min_on_period = & + any( currentSite%dstatus(ft) == [phen_dstat_timeon,phen_dstat_moiston] ) .and. & + (currentSite%dndaysleafon(ft) > dleafon_drycheck) + exceed_min_off_period = & + ( currentSite%dstatus(ft) == phen_dstat_timeoff ) .and. & + ( currentSite%dndaysleafoff(ft) > min_daysoff_dforcedflush ) + prolonged_on_period = & + any( currentSite%dstatus(ft) == [phen_dstat_timeon,phen_dstat_moiston] ) .and. & + ( currentSite%dndaysleafon(ft) > canopy_leaf_lifespan ) + prolonged_off_period = & + any( currentSite%dstatus(ft) == [phen_dstat_timeoff,phen_dstat_moistoff] ) .and. & + ( currentSite%dndaysleafoff(ft) > ED_val_phen_doff_time ) .and. & + ( currentSite%dndaysleafon(ft) >= 365-dd_offon_toler ) .and. & + ( currentSite%dndaysleafon(ft) <= 365+dd_offon_toler ) + last_flush_long_ago = & + ( currentSite%dstatus(ft) == phen_dstat_moistoff ) .and. & + ( currentSite%dndaysleafon(ft) > 365+dd_offon_toler ) + + + ! Revision of the conditions to simplify nested if and added an if/elseif/else + ! structure to ensure only up to one change occurs at any given time (ML 20211120). + drought_smoist_ifelse: if (model_day_int <= numWaterMem) then + ! Too early in the simulation. Do not change phenology status as we need to + ! populate the soil moisture memory. + continue + + elseif ( prolonged_off_period .and. ( .not. smoist_below_threshold ) ) then + ! LEAF ON: DROUGHT DECIDUOUS WETNESS + ! Here, we used a window of oppurtunity to determine if we are + ! close to the time when then leaves came on last year + ! The following conditions must be met + ! a) a year, plus or minus 1 month since we last had leaf-on? + ! b) Has there also been at least a nominaly short amount of "leaf-off"? + ! c) Is the soil moisture sufficiently high? + currentSite%dstatus(ft) = phen_dstat_moiston ! set status to leaf-on + currentSite%dleafondate(ft) = model_day_int ! save the model day we start flushing + currentSite%dndaysleafon(ft) = 0 + currentSite%elong_factor(ft) = 1. + + elseif ( last_flush_long_ago ) then + ! LEAF ON: DROUGHT DECIDUOUS TIME EXCEEDANCE + ! If we still haven't done budburst by end of window, then force it + + ! If the status is "phen_dstat_moistoff", it means this site currently has + ! leaves off due to actual moisture limitations. + ! So we trigger bud-burst at the end of the month since + ! last year's bud-burst. If this is imposed, then we set the new + ! status to indicate bud-burst was forced by timing + currentSite%dstatus(ft) = phen_dstat_timeon ! force budburst! + currentSite%dleafondate(ft) = model_day_int ! record leaf on date + currentSite%dndaysleafon(ft) = 0 + currentSite%elong_factor(ft) = 1. + + elseif ( exceed_min_off_period ) then + ! LEAF ON: DROUGHT DECIDUOUS EXCEEDED MINIMUM OFF PERIOD + ! Leaves were off due to time, not really moisture, so we allow them to + ! flush again as soon as they exceed a minimum off time + ! This typically occurs in a perennially wet system. + currentSite%dstatus(ft) = phen_dstat_timeon ! force budburst! + currentSite%dleafondate(ft) = model_day_int ! record leaf on date + currentSite%dndaysleafon(ft) = 0 + currentSite%elong_factor(ft) = 1. + + elseif ( prolonged_on_period ) then + ! LEAF OFF: DROUGHT DECIDUOUS LIFESPAN + ! Are the leaves rouhgly at the end of their lives? If so, shed leaves + ! even if it is not dry. + currentSite%dstatus(ft) = phen_dstat_timeoff !alter status of site to 'leaves off' + currentSite%dleafoffdate(ft) = model_day_int !record leaf on date + currentSite%dndaysleafoff(ft) = 0 + currentSite%elong_factor(ft) = 0. + + elseif ( exceed_min_on_period .and. smoist_below_threshold ) then + ! LEAF OFF: DROUGHT DECIDUOUS DRYNESS - if the soil gets too dry, + ! and the leaves have already been on a while... + currentSite%dstatus(ft) = phen_dstat_moistoff ! alter status of site to 'leaves off' + currentSite%dleafoffdate(ft) = model_day_int ! record leaf on date + currentSite%dndaysleafoff(ft) = 0 + currentSite%elong_factor(ft) = 0. + end if drought_smoist_ifelse + + + case (drgt_phen_model_gradual) + !------ + ! ED2-like deciduous. We compare the moisture with the lower and upper + ! thresholds. If the moisture is in between the thresholds, we must also + ! check whether or not the drought is developing or regressing. + !------ + + + ! First guess elongation factor + if (ED_val_phen_drought_threshold >= 0.) then + elongf_1st = elongf_min + (1.0_r8 - elongf_min ) * & + ( mean_10day_liqvol - ED_val_phen_drought_threshold ) / & + ( ED_val_phen_moist_threshold - ED_val_phen_drought_threshold ) + else + elongf_1st = elongf_min + (1.0_r8 - elongf_min ) * & + ( mean_10day_smp - ED_val_phen_drought_threshold ) / & + ( ED_val_phen_moist_threshold - ED_val_phen_drought_threshold ) + end if + elongf_1st = max(0.0_r8,min(1.0_r8,elongf_1st)) + + + + ! Save some conditions in logical variables to simplify code below + recent_flush = elongf_prev >= elongf_min .and. & + ( currentSite%dndaysleafon(ft) <= dleafon_drycheck ) + recent_abscission = elongf_prev < elongf_min .and. & + ( currentSite%dndaysleafoff(ft) <= min_daysoff_dforcedflush ) + prolonged_on_period = all( [elongf_prev,elongf_1st] >= elongf_min ) .and. & + ( currentSite%dndaysleafon(ft) > canopy_leaf_lifespan ) + last_flush_long_ago = all( [elongf_prev,elongf_1st] < elongf_min ) .and. & + ( currentSite%dndaysleafon(ft) > 365+dd_offon_toler ) + + + ! Make sure elongation factor is bounded and check for special cases. + drought_gradual_ifelse: if ( model_day_int <= numWaterMem ) then + ! Too early in the simulation, keep the same elongation factor as the day before. + currentSite%elong_factor(ft) = elongf_prev + + elseif ( prolonged_on_period ) then + ! Leaves have been on for too long and exceeded leaf lifespan. Force abscission + currentSite%elong_factor(ft) = 0.0_r8 ! Force full budburst + currentSite%dstatus(ft) = phen_dstat_timeoff ! Flag that this has been forced + currentSite%dleafoffdate(ft) = model_day_int ! Record leaf off date + currentSite%dndaysleafoff(ft) = 0 ! Reset clock + + elseif ( last_flush_long_ago ) then + ! Plant has not flushed at all for a very long time. Force flushing + currentSite%elong_factor(ft) = elongf_min ! Force minimum budburst + currentSite%dstatus(ft) = phen_dstat_timeon ! Flag that this has been forced + currentSite%dleafondate(ft) = model_day_int ! Record leaf on date + currentSite%dndaysleafon(ft) = 0 ! Reset clock + + elseif ( recent_flush .and. elongf_1st < elongf_prev ) then + ! Leaves have only recently reached flushed status. Elongation factor cannot decrease + currentSite%elong_factor(ft) = elongf_prev ! Elongation factor cannot decrease + currentSite%dstatus(ft) = phen_dstat_timeon ! Flag that this has been forced + + elseif ( recent_abscission .and. elongf_1st > elongf_min ) then + ! Leaves have only recently abscissed. Prevent plant to flush leaves. + currentSite%elong_factor(ft) = 0.0_r8 ! Elongation factor must remain 0. + currentSite%dstatus(ft) = phen_dstat_timeoff ! Flag that this has been forced + + elseif ( elongf_1st < elongf_min ) then + ! First guess of elongation factor below minimum. Impose full abscission. + currentSite%elong_factor(ft) = 0.0_r8 + + if (elongf_prev >= elongf_min ) then + ! This is the first day moisture fell below minimum. Flag change of status. + currentSite%dstatus(ft) = phen_dstat_moistoff ! Flag that this has not been forced + currentSite%dleafoffdate(ft) = model_day_int ! Record leaf off date + currentSite%dndaysleafoff(ft) = 0 ! Reset clock + end if - if ( (currentSite%dstatus == phen_dstat_moiston .or. & - currentSite%dstatus == phen_dstat_timeon ) .and. & - (model_day_int > numWaterMem) .and. & - (mean_10day_liqvol <= ED_val_phen_drought_threshold) .and. & - (dayssincedleafon > dleafon_drycheck ) ) then - currentSite%dstatus = phen_dstat_moistoff ! alter status of site to 'leaves off' - currentSite%dleafoffdate = model_day_int ! record leaf on date - endif + else + ! First guess of elongation factor is valid, use it. + currentSite%elong_factor(ft) = elongf_1st + + + if (elongf_prev < elongf_min ) then + ! This is the first day moisture allows leaves to exist. Flag change of status. + currentSite%dstatus(ft) = phen_dstat_moiston ! Flag that this has not been forced + currentSite%dleafondate(ft) = model_day_int ! Record leaf on date + currentSite%dndaysleafon(ft) = 0 ! Reset clock + elseif (elongf_1st < elongf_prev) then + currentSite%dstatus(ft) = phen_dstat_pshed ! Flag partial shedding, + ! but do not reset the clock + end if + end if drought_gradual_ifelse + end select case_drought_phen + end do pft_drgt_loop call phenology_leafonoff(currentSite) @@ -1082,22 +1234,49 @@ subroutine phenology_leafonoff(currentSite) type(ed_cohort_type), pointer :: currentCohort real(r8) :: leaf_c ! leaf carbon [kg] + real(r8) :: fnrt_c ! fine root carbon [kg] real(r8) :: sapw_c ! sapwood carbon [kg] real(r8) :: struct_c ! structural wood carbon [kg] + + real(r8) :: a_sapw0 ! target sapwood cross section area [m2] (dummy) + real(r8) :: c_agw0 ! target Above ground biomass [kgC] + real(r8) :: c_bgw0 ! target Below ground biomass [kgC] + + real(r8) :: leaf_target ! target leaf carbon (allometry scaled by elongation factor) [kg] + real(r8) :: fnrt_target ! target fine root carbon (allometry scaled by elongation factor) [kg] + real(r8) :: sapw_target ! target sapwood carbon (allometry scaled by elongation factor) [kg] + real(r8) :: struct_target ! target structural wood carbon (allometry scaled by elongation factor) [kg] + + real(r8) :: leaf_deficit ! leaf carbon deficit (relative to target) [kg] + real(r8) :: fnrt_deficit ! fine root carbon deficit (relative to target) [kg] + real(r8) :: sapw_deficit ! sapwood carbon deficit (relative to target) [kg] + real(r8) :: struct_deficit ! structural wood carbon deficit (relative to target) [kg] + real(r8) :: total_deficit ! total carbon deficit (relative to target) [kg] + + real(r8) :: eff_leaf_drop_fraction ! Effective leaf drop fraction + real(r8) :: eff_fnrt_drop_fraction ! Effective fine-root drop fraction + real(r8) :: eff_sapw_drop_fraction ! Effective sapwood drop fraction + real(r8) :: eff_struct_drop_fraction ! Effective structural wood drop fraction + real(r8) :: store_c ! storage carbon [kg] real(r8) :: store_c_transfer_frac ! Fraction of storage carbon used to flush leaves real(r8) :: totalmemory ! total memory of carbon [kg] + logical :: is_flushing_time ! Time to flush leaves + logical :: is_shedding_time ! Time to shed leaves integer :: ipft real(r8), parameter :: leaf_drop_fraction = 1.0_r8 real(r8), parameter :: carbon_store_buffer = 0.10_r8 + real(r8) :: fnrt_drop_fraction real(r8) :: stem_drop_fraction + + logical, parameter :: debug = .true. ! Print debug info? !------------------------------------------------------------------------ currentPatch => CurrentSite%oldest_patch - do while(associated(currentPatch)) + patch_loop: do while(associated(currentPatch)) currentCohort => currentPatch%tallest - do while(associated(currentCohort)) + cohort_loop: do while(associated(currentCohort)) ipft = currentCohort%pft @@ -1105,231 +1284,202 @@ subroutine phenology_leafonoff(currentSite) if(debug) call currentCohort%prt%CheckMassConservation(ipft,0) + ! Update memory variables to ensure they are up to date with recent growth. + call bleaf(currentCohort%dbh,ipft,currentCohort%canopy_trim,currentCohort%leafmemory) + call bfineroot(currentCohort%dbh,ipft,currentCohort%canopy_trim,currentCohort%fnrtmemory) + call bsap_allom(currentCohort%dbh,ipft,currentCohort%canopy_trim,a_sapw0, currentCohort%sapwmemory) + call bagw_allom(currentCohort%dbh,ipft,c_agw0) + call bbgw_allom(currentCohort%dbh,ipft,c_bgw0) + call bdead_allom(c_agw0,c_bgw0,currentCohort%sapwmemory,ipft,currentCohort%structmemory) + + ! Get current carbon pools store_c = currentCohort%prt%GetState(store_organ, all_carbon_elements) leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + fnrt_c = currentCohort%prt%GetState(fnrt_organ, all_carbon_elements) sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) + fnrt_drop_fraction = EDPftvarcon_inst%phen_fnrt_drop_fraction(ipft) stem_drop_fraction = EDPftvarcon_inst%phen_stem_drop_fraction(ipft) - ! COLD LEAF ON - ! The site level flags signify that it is no-longer too cold - ! for leaves. Time to signal flushing - - if (prt_params%season_decid(ipft) == itrue)then - if ( currentSite%cstatus == phen_cstat_notcold )then ! we have just moved to leaves being on . - if (currentCohort%status_coh == leaves_off)then ! Are the leaves currently off? - currentCohort%status_coh = leaves_on ! Leaves are on, so change status to - ! stop flow of carbon out of bstore. - - if(store_c>nearzero) then - ! flush either the amount required from the laimemory, or -most- of the storage pool - ! RF: added a criterion to stop the entire store pool emptying and triggering termination mortality - ! n.b. this might not be necessary if we adopted a more gradual approach to leaf flushing... - store_c_transfer_frac = min((EDPftvarcon_inst%phenflush_fraction(ipft)* & - currentCohort%laimemory)/store_c,(1.0_r8-carbon_store_buffer)) - - if(prt_params%woody(ipft).ne.itrue)then - totalmemory=currentCohort%laimemory+currentCohort%sapwmemory+currentCohort%structmemory - store_c_transfer_frac = min((EDPftvarcon_inst%phenflush_fraction(ipft)* & - totalmemory)/store_c, (1.0_r8-carbon_store_buffer)) - endif - - else - store_c_transfer_frac = 0.0_r8 - end if - - ! This call will request that storage carbon will be transferred to - ! leaf tissues. It is specified as a fraction of the available storage - if(prt_params%woody(ipft) == itrue) then - - call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, store_c_transfer_frac) - currentCohort%laimemory = 0.0_r8 - - else - - ! Check that the stem drop fraction is set to non-zero amount otherwise flush all carbon store to leaves - if (stem_drop_fraction .gt. 0.0_r8) then - - call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac*currentCohort%laimemory/totalmemory) - - call PRTPhenologyFlush(currentCohort%prt, ipft, sapw_organ, & - store_c_transfer_frac*currentCohort%sapwmemory/totalmemory) - - call PRTPhenologyFlush(currentCohort%prt, ipft, struct_organ, & - store_c_transfer_frac*currentCohort%structmemory/totalmemory) - - else - - call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac) - - end if - - currentCohort%laimemory = 0.0_r8 - currentCohort%structmemory = 0.0_r8 - currentCohort%sapwmemory = 0.0_r8 - - endif - endif !pft phenology - endif ! growing season - - !COLD LEAF OFF - if (currentSite%cstatus == phen_cstat_nevercold .or. & - currentSite%cstatus == phen_cstat_iscold) then ! past leaf drop day? Leaves still on tree? - - if (currentCohort%status_coh == leaves_on) then ! leaves have not dropped - - ! leaf off occur on individuals bigger than specific size for grass - if (currentCohort%dbh > EDPftvarcon_inst%phen_cold_size_threshold(ipft) & - .or. prt_params%woody(ipft)==itrue) then - - ! This sets the cohort to the "leaves off" flag - currentCohort%status_coh = leaves_off - - ! Remember what the lai was (leaf mass actually) was for next year - ! the same amount back on in the spring... - - currentCohort%laimemory = leaf_c - - ! Drop Leaves (this routine will update the leaf state variables, - ! for carbon and any other element that are prognostic. It will - ! also track the turnover masses that will be sent to litter later on) - - call PRTDeciduousTurnover(currentCohort%prt,ipft, & - leaf_organ, leaf_drop_fraction) - if(prt_params%woody(ipft).ne.itrue)then - - currentCohort%sapwmemory = sapw_c * stem_drop_fraction - - currentCohort%structmemory = struct_c * stem_drop_fraction - - call PRTDeciduousTurnover(currentCohort%prt,ipft, & - sapw_organ, stem_drop_fraction) - - call PRTDeciduousTurnover(currentCohort%prt,ipft, & - struct_organ, stem_drop_fraction) - - endif ! woody plant check - endif ! individual dbh size check - endif !leaf status - endif !currentSite status - endif !season_decid - - ! DROUGHT LEAF ON - ! Site level flag indicates it is no longer in drought condition - ! deciduous plants can flush - - if (prt_params%stress_decid(ipft) == itrue )then - - if (currentSite%dstatus == phen_dstat_moiston .or. & - currentSite%dstatus == phen_dstat_timeon )then - - ! we have just moved to leaves being on . - if (currentCohort%status_coh == leaves_off)then - - !is it the leaf-on day? Are the leaves currently off? - - currentCohort%status_coh = leaves_on ! Leaves are on, so change status to - ! stop flow of carbon out of bstore. - - if(store_c>nearzero) then - - store_c_transfer_frac = & - min((EDPftvarcon_inst%phenflush_fraction(ipft)*currentCohort%laimemory)/store_c, & - (1.0_r8-carbon_store_buffer)) - - if(prt_params%woody(ipft).ne.itrue)then - - totalmemory=currentCohort%laimemory+currentCohort%sapwmemory+currentCohort%structmemory - store_c_transfer_frac = min(EDPftvarcon_inst%phenflush_fraction(ipft)*totalmemory/store_c, & - (1.0_r8-carbon_store_buffer)) - - endif - - else - store_c_transfer_frac = 0.0_r8 - endif - - ! This call will request that storage carbon will be transferred to - ! leaf tissues. It is specified as a fraction of the available storage - if(prt_params%woody(ipft) == itrue) then - - call PRTPhenologyFlush(currentCohort%prt, ipft, & - leaf_organ, store_c_transfer_frac) - - currentCohort%laimemory = 0.0_r8 - - else - - ! Check that the stem drop fraction is set to non-zero amount otherwise flush all carbon store to leaves - if (stem_drop_fraction .gt. 0.0_r8) then - - call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac*currentCohort%laimemory/totalmemory) - - call PRTPhenologyFlush(currentCohort%prt, ipft, sapw_organ, & - store_c_transfer_frac*currentCohort%sapwmemory/totalmemory) - - call PRTPhenologyFlush(currentCohort%prt, ipft, struct_organ, & - store_c_transfer_frac*currentCohort%structmemory/totalmemory) - - else - - call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac) - - end if - - currentCohort%laimemory = 0.0_r8 - currentCohort%structmemory = 0.0_r8 - currentCohort%sapwmemory = 0.0_r8 - - endif ! woody plant check - endif !currentCohort status again? - endif !currentSite status - - !DROUGHT LEAF OFF - if (currentSite%dstatus == phen_dstat_moistoff .or. & - currentSite%dstatus == phen_dstat_timeoff) then - - if (currentCohort%status_coh == leaves_on) then ! leaves have not dropped - - ! This sets the cohort to the "leaves off" flag - currentCohort%status_coh = leaves_off - - ! Remember what the lai (leaf mass actually) was for next year - currentCohort%laimemory = leaf_c - - call PRTDeciduousTurnover(currentCohort%prt,ipft, & - leaf_organ, leaf_drop_fraction) - - if(prt_params%woody(ipft).ne.itrue)then - - currentCohort%sapwmemory = sapw_c * stem_drop_fraction - currentCohort%structmemory = struct_c * stem_drop_fraction + ! MLO. Sanity check. If this is a deciduous PFT, then leaf status ought to + ! be on or off. If it is something else, stop the run. + if ( prt_params%season_decid(ipft) == itrue .or. & + prt_params%stress_decid(ipft) == itrue ) then + select case(currentCohort%status_coh) + case (leaves_off,leaves_on,leaves_pshed) + continue + case default + write(fates_log(),'(a)') '--------------------------------------------------------------------' + write(fates_log(),'(a)') ' Odd leaf status - Leaves are not off, on, or partially abscissing. ' + write(fates_log(),'(a)') '--------------------------------------------------------------------' + write(fates_log(),'(a,1x,i5)' ) ' PFT = ',ipft + write(fates_log(),'(a,1x,i5)' ) ' Season deciduous = ',prt_params%season_decid(ipft) + write(fates_log(),'(a,1x,i5)' ) ' Stress deciduous = ',prt_params%stress_decid(ipft) + write(fates_log(),'(a,1x,i5)' ) ' Site status = ',currentSite%cstatus + write(fates_log(),'(a,1x,i5)' ) ' Cohort status = ',currentCohort%status_coh + write(fates_log(),'(a,1x,f12.5)') ' Leaf elongation = ',currentCohort%efleaf_coh + write(fates_log(),'(a,1x,f12.5)') ' Fine-root elongation = ',currentCohort%effnrt_coh + write(fates_log(),'(a,1x,f12.5)') ' Stem elongation = ',currentCohort%efstem_coh + write(fates_log(),'(a,1x,f12.5)') ' DBH = ',currentCohort%dbh + write(fates_log(),'(a,1x,f12.5)') ' Leaf_c = ',leaf_c + write(fates_log(),'(a,1x,f12.5)') ' Store_c = ',store_c + write(fates_log(),'(a,1x,f12.5)') ' Fnrt_c = ',fnrt_c + write(fates_log(),'(a,1x,f12.5)') ' Sapw_c = ',sapw_c + write(fates_log(),'(a,1x,f12.5)') ' Struct_c = ',struct_c + write(fates_log(),'(a)') '--------------------------------------------------------------------' + write(fates_log(),'(a)') '' + write(fates_log(),'(a)') '' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + end if - call PRTDeciduousTurnover(currentCohort%prt,ipft, & - sapw_organ, stem_drop_fraction) - call PRTDeciduousTurnover(currentCohort%prt,ipft, & - struct_organ, stem_drop_fraction) - endif + ! MLO. To avoid duplicating code for drought and cold deciduous PFTs, we first + ! check whether or not it's time to flush or time to shed leaves, then + ! use a common code for flushing or shedding leaves. + if (prt_params%season_decid(ipft) == itrue) then ! Cold deciduous + ! Set elongation factor to 0 or 1 (partial shedding not defined for cold deciduous) + select case (currentSite%cstatus) + case (phen_cstat_nevercold,phen_cstat_iscold) + currentCohort%efleaf_coh = 0.0_r8 + case (phen_cstat_notcold) + currentCohort%efleaf_coh = 1.0_r8 + end select + + ! A. Is this the time for COLD LEAVES to switch to ON? + is_flushing_time = ( currentSite%cstatus == phen_cstat_notcold .and. & ! We just moved to leaves being on + currentCohort%status_coh == leaves_off ) ! Leaves are currently off + ! B. Is this the time for COLD LEAVES to switch to OFF? + is_shedding_time = any(currentSite%cstatus == [phen_cstat_nevercold,phen_cstat_iscold]) .and. & ! Past leaf drop day or too cold + currentCohort%status_coh == leaves_on .and. & ! Leaves have not dropped yet + ( currentCohort%dbh > EDPftvarcon_inst%phen_cold_size_threshold(ipft) .or. & ! Grasses are big enough or... + prt_params%woody(ipft) == itrue ) ! this is a woody PFT. + + elseif (prt_params%stress_decid(ipft) == itrue ) then ! Drought deciduous + ! Use elongation factor + currentCohort%efleaf_coh = currentSite%elong_factor(ipft) + + + ! A. Is this the time for DROUGHT LEAVES to switch to ON? + is_flushing_time = any( currentSite%dstatus(ipft) == [phen_dstat_moiston,phen_dstat_timeon] ) .and. & ! Leaf flushing time (moisture or time) + any( currentCohort%status_coh == [leaves_off,leaves_pshed] ) + ! B. Is this the time for DROUGHT LEAVES to switch to OFF? + ! This will be true when leaves are abscissing (partially or fully) due to moisture or time + is_shedding_time = any( currentSite%dstatus(ipft) == [phen_dstat_moistoff,phen_dstat_timeoff,phen_dstat_pshed] ) .and. & + any( currentCohort%status_coh == [leaves_on,leaves_pshed] ) + else + ! This PFT is not deciduous. + is_flushing_time = .false. + is_shedding_time = .false. + currentCohort%efleaf_coh = 1.0_r8 + end if ! (prt_params%season_decid(ipft) == itrue) + + ! Find the effective "elongation factor" for fine roots and stems. The effective drop fraction is + ! a combination of the elongation factor (e) and the drop fraction (x), which will ensure + ! that the remaining tissue biomass will be exactly e when x=1, and exactly the original + ! biomass when x = 0. For leaves it is always assumed that the drop fraction is one. + currentCohort%effnrt_coh = 1.0_r8 - (1.0_r8 - currentCohort%efleaf_coh ) * fnrt_drop_fraction + currentCohort%efstem_coh = 1.0_r8 - (1.0_r8 - currentCohort%efleaf_coh ) * stem_drop_fraction + + ! Find the target biomass for each tissue when accounting for elongation factors. + ! Note that the target works for both flushing and shedding leaves. + leaf_target = currentCohort%efleaf_coh * currentCohort%leafmemory + fnrt_target = currentCohort%effnrt_coh * currentCohort%fnrtmemory + sapw_target = currentCohort%efstem_coh * currentCohort%sapwmemory + struct_target = currentCohort%efstem_coh * currentCohort%structmemory + + + ! A. This is time to switch to (COLD or DROUGHT) LEAF ON + flush_block: if (is_flushing_time) then + currentCohort%status_coh = leaves_on ! Leaves are on, so change status to + ! stop flow of carbon out of bstore. + + ! Transfer carbon from storage to living tissues (only if there is any carbon in storage) + transf_block: if ( store_c > nearzero ) then + ! Find the total deficit. We no longer distinguish between woody and non-woody + ! PFTs here (as sapwmemory is the same as sapw_c if this is a woody tissue). + leaf_deficit = max(0.0_r8, leaf_target - leaf_c ) + fnrt_deficit = max(0.0_r8, fnrt_target - fnrt_c ) + sapw_deficit = max(0.0_r8, sapw_target - sapw_c ) + struct_deficit = max(0.0_r8, struct_target - struct_c) + total_deficit = leaf_deficit + fnrt_deficit + sapw_deficit + struct_deficit + + ! Flush either the amount required from the memory, or -most- of the storage pool + ! RF: added a criterion to stop the entire store pool emptying and triggering termination mortality + ! n.b. this might not be necessary if we adopted a more gradual approach to leaf flushing... + store_c_transfer_frac = min( EDPftvarcon_inst%phenflush_fraction(ipft) * total_deficit / store_c, & + 1.0_r8 - carbon_store_buffer ) + + ! This call will request that storage carbon will be transferred to + ! each tissue. It is specified as a fraction of the available storage + ! MLO - Just to be safe, skip steps in the unlikely case total_deficit is zero, to avoid FPE errors. + if (total_deficit > nearzero) then + call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & + store_c_transfer_frac*leaf_deficit/total_deficit) + call PRTPhenologyFlush(currentCohort%prt, ipft, fnrt_organ, & + store_c_transfer_frac*fnrt_deficit/total_deficit) + + if ( nint(prt_params%woody(ipft)) == ifalse ) then + call PRTPhenologyFlush(currentCohort%prt, ipft, sapw_organ, & + store_c_transfer_frac*sapw_deficit/total_deficit) + call PRTPhenologyFlush(currentCohort%prt, ipft, struct_organ, & + store_c_transfer_frac*struct_deficit/total_deficit) + end if + end if + else + ! Not enough carbon to flush any living tissue. + store_c_transfer_frac = 0.0_r8 + end if transf_block + end if flush_block + + + + ! B. This is time to switch to (COLD or DROUGHT) LEAF OFF + shed_block: if (is_shedding_time) then + if ( currentCohort%efleaf_coh > 0.0_r8 ) then + ! Partial shedding + currentCohort%status_coh = leaves_pshed + else + ! Complete abscission + currentCohort%status_coh = leaves_off + end if - endif - endif !status - endif !drought dec. + ! Find the effective fraction to drop. This fraction must be calculated every time + ! because we must account for partial abscission. The simplest approach is to simply + ! use the ratio between the target and the original biomass of each pool. The + ! max(tissue_c,nearzero) is overly cautious, because leaf_c = 0 would imply that + ! leaves are already off, and this wouldn't be considered shedding time. + eff_leaf_drop_fraction = max( 0.0_r8, min( 1.0_r8,1.0_r8 - leaf_target / max( leaf_c , nearzero ) ) ) + eff_fnrt_drop_fraction = max( 0.0_r8, min( 1.0_r8,1.0_r8 - fnrt_target / max( fnrt_c , nearzero ) ) ) + eff_sapw_drop_fraction = max( 0.0_r8, min( 1.0_r8,1.0_r8 - sapw_target / max( sapw_c , nearzero ) ) ) + eff_struct_drop_fraction = max( 0.0_r8, min( 1.0_r8,1.0_r8 - struct_target / max( struct_c, nearzero ) ) ) + + ! Drop leaves + call PRTDeciduousTurnover(currentCohort%prt,ipft, leaf_organ, eff_leaf_drop_fraction) + + ! Drop fine roots + call PRTDeciduousTurnover(currentCohort%prt,ipft, fnrt_organ, eff_fnrt_drop_fraction) + + ! If plant is not woody, shed sapwood and heartwood (they may have a minimum amount of woody tissues for + ! running plant hydraulics, and it makes sense to shed them along with leaves when they should be off). + ! MLO - stem_drop_fraction is a PFT parameter, do we really need this check for woody/non-woody PFT? + if ( nint(prt_params%woody(ipft)) == ifalse ) then + ! Shed sapwood and heartwood. + call PRTDeciduousTurnover(currentCohort%prt,ipft,sapw_organ , eff_sapw_drop_fraction ) + call PRTDeciduousTurnover(currentCohort%prt,ipft,struct_organ, eff_struct_drop_fraction) + end if + end if shed_block if(debug) call currentCohort%prt%CheckMassConservation(ipft,1) currentCohort => currentCohort%shorter - enddo !currentCohort + end do cohort_loop currentPatch => currentPatch%younger - enddo !currentPatch + end do patch_loop end subroutine phenology_leafonoff @@ -1761,7 +1911,7 @@ subroutine SeedGermination( litt, cold_stat, drought_stat ) ! !ARGUMENTS type(litter_type) :: litt integer, intent(in) :: cold_stat ! Is the site in cold leaf-off status? - integer, intent(in) :: drought_stat ! Is the site in drought leaf-off status? + integer, dimension(numpft), intent(in) :: drought_stat ! Is the site in drought leaf-off status? ! ! !LOCAL VARIABLES: integer :: pft @@ -1793,7 +1943,7 @@ subroutine SeedGermination( litt, cold_stat, drought_stat ) litt%seed_germ_in(pft) = 0.0_r8 endif if ((prt_params%stress_decid(pft) == itrue ) .and. & - (any(drought_stat == [phen_dstat_timeoff,phen_dstat_moistoff]))) then + (any(drought_stat(pft) == [phen_dstat_timeoff,phen_dstat_moistoff,phen_dstat_pshed]))) then litt%seed_germ_in(pft) = 0.0_r8 end if @@ -1837,11 +1987,14 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) real(r8) :: c_leaf ! target leaf biomass [kgC] real(r8) :: c_fnrt ! target fine root biomass [kgC] real(r8) :: c_sapw ! target sapwood biomass [kgC] - real(r8) :: a_sapw ! target sapwood cross section are [m2] (dummy) + real(r8) :: a_sapw ! target sapwood cross section area [m2] (dummy) real(r8) :: c_agw ! target Above ground biomass [kgC] real(r8) :: c_bgw ! target Below ground biomass [kgC] real(r8) :: c_struct ! target Structural biomass [kgc] real(r8) :: c_store ! target Storage biomass [kgC] + real(r8) :: elongf_leaf ! leaf elongation factor [fraction] + real(r8) :: elongf_fnrt ! fine-root "elongation factor" [fraction] + real(r8) :: elongf_stem ! stem "elongation factor" [fraction] real(r8) :: m_leaf ! leaf mass (element agnostic) [kg] real(r8) :: m_fnrt ! fine-root mass (element agnostic) [kg] real(r8) :: m_sapw ! sapwood mass (element agnostic) [kg] @@ -1853,6 +2006,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) real(r8) :: mass_avail ! The mass of each nutrient/carbon available in the seed_germination pool [kg] real(r8) :: mass_demand ! Total mass demanded by the plant to achieve the stoichiometric targets ! of all the organs in the recruits. Used for both [kg per plant] and [kg per cohort] + real(r8) :: fnrt_drop_fraction real(r8) :: stem_drop_fraction !---------------------------------------------------------------------- @@ -1862,7 +2016,6 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) 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. @@ -1875,7 +2028,8 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) temp_cohort%pft = ft temp_cohort%hite = EDPftvarcon_inst%hgt_min(ft) temp_cohort%coage = 0.0_r8 - stem_drop_fraction = EDPftvarcon_inst%phen_stem_drop_fraction(ft) + fnrt_drop_fraction = EDPftvarcon_inst%phen_fnrt_drop_fraction(ft) + stem_drop_fraction = EDPftvarcon_inst%phen_stem_drop_fraction(ft) call h2d_allom(temp_cohort%hite,ft,temp_cohort%dbh) @@ -1888,47 +2042,62 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) call bdead_allom(c_agw,c_bgw,c_sapw,ft,c_struct) call bstore_allom(temp_cohort%dbh,ft,temp_cohort%canopy_trim,c_store) - ! Default assumption is that leaves are on + ! Default assumption is that leaves are on and fully flushed cohortstatus = leaves_on - temp_cohort%laimemory = 0.0_r8 - temp_cohort%sapwmemory = 0.0_r8 - temp_cohort%structmemory = 0.0_r8 + elongf_leaf = 1.0_r8 + elongf_fnrt = 1.0_r8 + elongf_stem = 1.0_r8 + + ! MLO update 20211123. The "memory" variables are now always set to the on-allometry state + ! (accounting for canopy trimming when necessary). We no longer reset them to zero when + ! leaves are flushing. This makes partial deciduousness a bit easier to implement. + temp_cohort%leafmemory = c_leaf + temp_cohort%fnrtmemory = c_fnrt + temp_cohort%sapwmemory = c_sapw + temp_cohort%structmemory = c_struct ! But if the plant is seasonally (cold) deciduous, and the site status is flagged ! as "cold", then set the cohort's status to leaves_off, and remember the leaf biomass if ((prt_params%season_decid(ft) == itrue) .and. & (any(currentSite%cstatus == [phen_cstat_nevercold,phen_cstat_iscold]))) then - temp_cohort%laimemory = c_leaf - c_leaf = 0.0_r8 - - ! If plant is not woody then set sapwood and structural biomass as well - if (prt_params%woody(ft).ne.itrue) then - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction - c_sapw = (1.0_r8 - stem_drop_fraction) * c_sapw - c_struct = (1.0_r8 - stem_drop_fraction) * c_struct - endif - cohortstatus = leaves_off - endif + elongf_leaf = 0.0_r8 + elongf_fnrt = 1.0_r8 - fnrt_drop_fraction + elongf_stem = 1.0_r8 - stem_drop_fraction + + + c_leaf = elongf_leaf * c_leaf + c_fnrt = elongf_fnrt * c_fnrt + c_sapw = elongf_stem * c_sapw + c_struct = elongf_stem * c_struct - ! Or.. if the plant is drought deciduous, and the site status is flagged as - ! "in a drought", then likewise, set the cohort's status to leaves_off, and remember leaf - ! biomass - if ((prt_params%stress_decid(ft) == itrue) .and. & - (any(currentSite%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff]))) then - temp_cohort%laimemory = c_leaf - c_leaf = 0.0_r8 - - ! If plant is not woody then set sapwood and structural biomass as well - if(prt_params%woody(ft).ne.itrue)then - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction - c_sapw = (1.0_r8 - stem_drop_fraction) * c_sapw - c_struct = (1.0_r8 - stem_drop_fraction) * c_struct - endif cohortstatus = leaves_off - endif + end if + + ! Or.. if the plant is drought deciduous, make sure leaf status is consistent with the + ! leaf elongation factor. + ! For tissues other than leaves, the actual drop fraction is a combination of the + ! elongation factor (e) and the drop fraction (x), which will ensure that the remaining + ! tissue biomass will be exactly e when x=1, and exactly the original biomass when x = 0. + if ( prt_params%stress_decid(ft) == itrue ) then + elongf_leaf = currentSite%elong_factor(ft) + elongf_fnrt = 1.0_r8 - (1.0_r8 - elongf_leaf ) * fnrt_drop_fraction + elongf_stem = 1.0_r8 - (1.0_r8 - elongf_leaf ) * stem_drop_fraction + + c_leaf = elongf_leaf * c_leaf + c_fnrt = elongf_fnrt * c_fnrt + c_sapw = elongf_stem * c_sapw + c_struct = elongf_stem * c_struct + + ! For the initial state, we always assume that leaves are flushing (instead of partially abscissing) + ! whenever the elongation factor is non-zero. If the elongation factor is zero, then leaves are in + ! the "off" state. + if ( elongf_leaf > 0.0_r8 ) then + cohortstatus = leaves_on + else + cohortstatus = leaves_off + end if + end if ! Cycle through available carbon and nutrients, find the limiting element @@ -2044,7 +2213,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) end select select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) ! Put all of the leaf mass into the first bin call SetState(prt,leaf_organ, element_id,m_leaf,1) @@ -2103,9 +2272,9 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) ! This initializes the cohort call create_cohort(currentSite,currentPatch, temp_cohort%pft, temp_cohort%n, & temp_cohort%hite, temp_cohort%coage, temp_cohort%dbh, prt, & - temp_cohort%laimemory, temp_cohort%sapwmemory, temp_cohort%structmemory, & - cohortstatus, recruitstatus, & - temp_cohort%canopy_trim,temp_cohort%c_area, & + temp_cohort%leafmemory, temp_cohort%fnrtmemory, temp_cohort%sapwmemory, & + temp_cohort%structmemory, elongf_leaf, elongf_fnrt, elongf_stem, cohortstatus, & + recruitstatus, temp_cohort%canopy_trim,temp_cohort%c_area, & currentPatch%NCL_p, currentSite%spread, bc_in) ! Note that if hydraulics is on, the number of cohorts may had diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 42264ca776..8d0772571a 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -158,6 +158,7 @@ module FatesAllometryMod ! ============================================================================ subroutine CheckIntegratedAllometries(dbh,ipft,canopy_trim, & + elongf_leaf, elongf_fnrt, elongf_stem, & bl,bfr,bsap,bstore,bdead, & grow_leaf, grow_fr, grow_sap, grow_store, grow_dead, & max_err, l_pass) @@ -173,6 +174,9 @@ subroutine CheckIntegratedAllometries(dbh,ipft,canopy_trim, & real(r8),intent(in) :: dbh ! diameter of plant [cm] integer,intent(in) :: ipft ! plant functional type index real(r8),intent(in) :: canopy_trim ! trimming function + real(r8),intent(in) :: elongf_leaf ! Leaf elongation factor + real(r8),intent(in) :: elongf_fnrt ! Fine-root elongation factor + real(r8),intent(in) :: elongf_stem ! Stem elongation factor real(r8),intent(in) :: bl ! integrated leaf biomass [kgC] real(r8),intent(in) :: bfr ! integrated fine root biomass [kgC] real(r8),intent(in) :: bsap ! integrated sapwood biomass [kgC] @@ -203,6 +207,7 @@ subroutine CheckIntegratedAllometries(dbh,ipft,canopy_trim, & if (grow_leaf) then call bleaf(dbh,ipft,canopy_trim,bl_diag) + bl_diag = bl_diag * elongf_leaf if( abs(bl_diag-bl) > max_err ) then if(verbose_logging) then write(fates_log(),*) 'disparity in integrated/diagnosed leaf carbon' @@ -217,6 +222,7 @@ subroutine CheckIntegratedAllometries(dbh,ipft,canopy_trim, & if (grow_fr) then call bfineroot(dbh,ipft,canopy_trim,bfr_diag) + bfr_diag = bfr_diag * elongf_fnrt if( abs(bfr_diag-bfr) > max_err ) then if(verbose_logging) then write(fates_log(),*) 'disparity in integrated/diagnosed fineroot carbon' @@ -231,6 +237,7 @@ subroutine CheckIntegratedAllometries(dbh,ipft,canopy_trim, & if (grow_sap) then call bsap_allom(dbh,ipft,canopy_trim,asap_diag,bsap_diag) + bsap_diag = bsap_diag * elongf_stem if( abs(bsap_diag-bsap) > max_err ) then if(verbose_logging) then write(fates_log(),*) 'disparity in integrated/diagnosed sapwood carbon' @@ -261,7 +268,9 @@ subroutine CheckIntegratedAllometries(dbh,ipft,canopy_trim, & call bsap_allom(dbh,ipft,canopy_trim,asap_diag,bsap_diag) call bagw_allom(dbh,ipft,bagw_diag) call bbgw_allom(dbh,ipft,bbgw_diag) - call bdead_allom( bagw_diag, bbgw_diag, bsap_diag, ipft, bdead_diag ) + call bdead_allom( bagw_diag, bbgw_diag, bsap_diag, ipft, bdead_diag ) + bdead_diag = bdead_diag * elongf_stem + if( abs(bdead_diag-bdead) > max_err ) then if(verbose_logging) then write(fates_log(),*) 'disparity in integrated/diagnosed structural carbon' @@ -869,6 +878,8 @@ subroutine bsap_allom(d,ipft,canopy_trim,sapw_area,bsap,dbsapdd) real(r8) :: dhdd real(r8) :: bl real(r8) :: dbldd + real(r8) :: blmax ! maximum leaf biomss per allometry + real(r8) :: dblmaxdd real(r8) :: bbgw real(r8) :: dbbgwdd real(r8) :: bagw @@ -886,7 +897,7 @@ subroutine bsap_allom(d,ipft,canopy_trim,sapw_area,bsap,dbsapdd) ! Currently only one sapwood allometry model. the slope ! of the la:sa to diameter line is zero. ! --------------------------------------------------------------------- - case(1) ! linearly related to leaf area based on target leaf biomass + case(1) ! linearly related to leaf area based on target TRIMMED leaf biomass ! and slatop (no provisions for slamax) call h_allom(d,ipft,h,dhdd) @@ -902,6 +913,28 @@ subroutine bsap_allom(d,ipft,canopy_trim,sapw_area,bsap,dbsapdd) ! (this comes into play typically in very small plants) bsap_cap = max_frac*(bagw+bbgw) + if(bsap>bsap_cap) then + bsap = bsap_cap + if(present(dbsapdd))then + dbsapdd = max_frac*(dbagwdd+dbbgwdd) + end if + end if + case(2) ! linearly related to leaf area based on target UNTRIMMED leaf biomass + ! and slatop (no provisions for slamax) + + call h_allom(d,ipft,h,dhdd) + call blmax_allom(d,ipft,blmax,dblmaxdd) + call bsap_ltarg_slatop(d,h,dhdd,blmax,dblmaxdd,ipft,sapw_area,bsap,dbsapdd) + + ! Perform a capping/check on total woody biomass + call bagw_allom(d,ipft,bagw,dbagwdd) + call bbgw_allom(d,ipft,bbgw,dbbgwdd) + + ! Force sapwood to be less than a maximum fraction of total biomass + ! We omit the sapwood area from this calculation + ! (this comes into play typically in very small plants) + bsap_cap = max_frac*(bagw+bbgw) + if(bsap>bsap_cap) then bsap = bsap_cap if(present(dbsapdd))then @@ -1010,30 +1043,35 @@ subroutine bstore_allom(d,ipft,canopy_trim,bstore,dbstoredd) real(r8),intent(in) :: canopy_trim ! Crown trimming function [0-1] real(r8),intent(out) :: bstore ! allometric target storage [kgC] real(r8),intent(out),optional :: dbstoredd ! change storage per cm [kgC/cm] - - real(r8) :: bl ! Allometric target leaf biomass - real(r8) :: dbldd ! Allometric target change in leaf biomass per cm - - + + real(r8) :: bl ! Allometric target leaf biomass (TRIMMED) + real(r8) :: dbldd ! Allometric target change in leaf biomass per cm (TRIMMED) + real(r8) :: blmax ! Allometric target leaf biomass (UNTRIMMED) + real(r8) :: dblmaxdd ! Allometric target change in leaf biomass per cm (UNTRIMMED) + + ! TODO: allom_stmode needs to be added to the parameter file - associate( allom_stmode => prt_params%allom_stmode(ipft), & cushion => prt_params%cushion(ipft) ) select case(int(allom_stmode)) - case(1) ! Storage is constant proportionality of trimmed maximum leaf + case(1) ! Storage is constant proportionality of TRIMMED maximum leaf ! biomass (ie cushion * bleaf) - call bleaf(d,ipft,canopy_trim,bl,dbldd) call bstore_blcushion(d,bl,dbldd,cushion,ipft,bstore,dbstoredd) - + + case(2) ! Storage is constant proportionality of UNTRIMMED maximum leaf + ! biomass (ie cushion * bleaf) + call blmax_allom(d,ipft,blmax,dblmaxdd) + call bstore_blcushion(d,blmax,dblmaxdd,cushion,ipft,bstore,dbstoredd) + case DEFAULT write(fates_log(),*) 'An undefined fine storage allometry was specified: ', & allom_stmode write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end select - + end associate return end subroutine bstore_allom diff --git a/biogeochem/FatesSoilBGCFluxMod.F90 b/biogeochem/FatesSoilBGCFluxMod.F90 index d14ce7b005..c6d7089802 100644 --- a/biogeochem/FatesSoilBGCFluxMod.F90 +++ b/biogeochem/FatesSoilBGCFluxMod.F90 @@ -15,6 +15,7 @@ module FatesSoilBGCFluxMod use PRTGenericMod , only : num_elements use PRTGenericMod , only : element_list use PRTGenericMod , only : element_pos + use PRTGenericMod , only : prt_csimpler_allom_hyp use PRTGenericMod , only : prt_carbon_allom_hyp use PRTGenericMod , only : prt_cnp_flex_allom_hyp use PRTGenericMod , only : prt_vartypes @@ -219,7 +220,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_csimpler_allom_hyp,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 +229,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..5148c2275b 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -39,8 +39,9 @@ module FATESPlantRespPhotosynthMod use PRTGenericMod, only : max_nleafage use EDTypesMod, only : do_fates_salinity use EDParamsMod, only : q10_mr + use PRTGenericMod, only : prt_csimpler_allom_hyp use PRTGenericMod, only : prt_carbon_allom_hyp - use PRTGenericMod, only : prt_cnp_flex_allom_hyp + use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : all_carbon_elements use PRTGenericMod, only : nitrogen_element use PRTGenericMod, only : leaf_organ @@ -417,7 +418,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if ( .not.rate_mask_z(iv,ft,cl) .or. & (hlm_use_planthydro.eq.itrue) .or. & (nleafage > 1) .or. & - (hlm_parteh_mode .ne. prt_carbon_allom_hyp ) ) then + (hlm_parteh_mode .ne. prt_csimpler_allom_hyp .and. & + hlm_parteh_mode .ne. prt_carbon_allom_hyp) ) then if (hlm_use_planthydro.eq.itrue ) then @@ -480,7 +482,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! Then scale this value at the top of the canopy for canopy depth ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp) lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft) @@ -500,6 +502,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) @@ -639,7 +642,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) fnrt_c = currentCohort%prt%GetState(fnrt_organ, all_carbon_elements) select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp) live_stem_n = prt_params%allom_agb_frac(currentCohort%pft) * & sapw_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 06bcd1858a..93c94f53e4 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -61,6 +61,7 @@ module EDInitMod use FatesAllometryMod , only : bstore_allom use PRTGenericMod , only : StorageNutrientTarget use FatesInterfaceTypesMod, only : hlm_parteh_mode + use PRTGenericMod, only : prt_csimpler_allom_hyp use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes @@ -189,16 +190,23 @@ subroutine zero_site( site_in ) ! PHENOLOGY site_in%cstatus = fates_unset_int ! are leaves in this pixel on or off? - site_in%dstatus = fates_unset_int + site_in%dstatus(:) = fates_unset_int site_in%grow_deg_days = nan ! growing degree days site_in%snow_depth = nan site_in%nchilldays = fates_unset_int site_in%ncolddays = fates_unset_int - site_in%cleafondate = fates_unset_int ! doy of leaf on - site_in%cleafoffdate = fates_unset_int ! doy of leaf off - site_in%dleafondate = fates_unset_int ! doy of leaf on drought - site_in%dleafoffdate = fates_unset_int ! doy of leaf on drought - site_in%water_memory(:) = nan + site_in%cleafondate = fates_unset_int ! doy of leaf on (cold) + site_in%cleafoffdate = fates_unset_int ! doy of leaf off (cold) + site_in%dleafondate(:) = fates_unset_int ! doy of leaf on (drought) + site_in%dleafoffdate(:) = fates_unset_int ! doy of leaf off (drought) + site_in%cndaysleafon = fates_unset_int ! days since leaf on (cold) + site_in%cndaysleafoff = fates_unset_int ! days since leaf off (cold) + site_in%dndaysleafon(:) = fates_unset_int ! days since leaf on (drought) + site_in%dndaysleafoff(:) = fates_unset_int ! days since leaf off (drought) + site_in%elong_factor(:) = nan ! Elongation factor (0 - full abscission; 1 - fully flushed) + + site_in%liqvol_memory(:,:) = nan + site_in%smp_memory(:,:) = nan site_in%vegtemp_memory(:) = nan ! record of last 10 days temperature for senescence model. ! Disturbance rates tracking @@ -278,11 +286,17 @@ subroutine set_site_properties( nsites, sites,bc_in ) real(r8) :: GDD integer :: dstat ! drought status phenology flag real(r8) :: acc_NI - real(r8) :: watermem + real(r8) :: liqvolmem + real(r8) :: smpmem + real(r8) :: elong_factor ! Elongation factor (0 - fully off; 1 - fully on) integer :: cleafon ! DOY for cold-decid leaf-on, initial guess integer :: cleafoff ! DOY for cold-decid leaf-off, initial guess integer :: dleafoff ! DOY for drought-decid leaf-off, initial guess integer :: dleafon ! DOY for drought-decid leaf-on, initial guess + integer :: cndleafon ! days since leaf on (cold), initial guess + integer :: cndleafoff ! days since leaf off (cold), initial guess + integer :: dndleafon ! days since leaf on (drought), initial guess + integer :: dndleafoff ! days since leaf off (drought), initial guess integer :: ft ! PFT loop real(r8) :: sumarea ! area of PFTs in nocomp mode. integer :: hlm_pft ! used in fixed biogeog mode @@ -300,12 +314,18 @@ subroutine set_site_properties( nsites, sites,bc_in ) GDD = 30.0_r8 cleafon = 100 cleafoff = 300 + cndleafon = 0 + cndleafoff = 0 cstat = phen_cstat_notcold ! Leaves are on acc_NI = 0.0_r8 dstat = phen_dstat_moiston ! Leaves are on - dleafoff = 300 dleafon = 100 - watermem = 0.5_r8 + dleafoff = 300 + dndleafon = 0 + dndleafoff = 0 + liqvolmem = 0.5_r8 + smpmem = 0._r8 + elong_factor = 1._r8 do s = 1,nsites sites(s)%nchilldays = 0 @@ -316,15 +336,21 @@ subroutine set_site_properties( nsites, sites,bc_in ) sites(s)%cleafondate = cleafon sites(s)%cleafoffdate = cleafoff - sites(s)%dleafoffdate = dleafoff - sites(s)%dleafondate = dleafon + sites(s)%cndaysleafon = cndleafon + sites(s)%cndaysleafoff = cndleafoff + sites(s)%dleafoffdate(1:numpft) = dleafoff + sites(s)%dleafondate(1:numpft) = dleafon + sites(s)%dndaysleafon(1:numpft) = dndleafon + sites(s)%dndaysleafoff(1:numpft) = dndleafoff sites(s)%grow_deg_days = GDD - sites(s)%water_memory(1:numWaterMem) = watermem + sites(s)%liqvol_memory(1:numWaterMem,1:numpft) = liqvolmem + sites(s)%smp_memory(1:numWaterMem,1:numpft) = smpmem sites(s)%vegtemp_memory(1:num_vegtemp_mem) = 0._r8 sites(s)%cstatus = cstat - sites(s)%dstatus = dstat + sites(s)%dstatus(1:numpft) = dstat + sites(s)%elong_factor(1:numpft) = elong_factor sites(s)%acc_NI = acc_NI sites(s)%NF = 0.0_r8 @@ -696,7 +722,11 @@ subroutine init_cohorts( site_in, patch_in, bc_in) real(r8) :: m_sapw ! Generic mass for sapwood [kg] real(r8) :: m_store ! Generic mass for storage [kg] real(r8) :: m_repro ! Generic mass for reproductive tissues [kg] - real(r8) :: stem_drop_fraction + real(r8) :: elongf_leaf ! Leaf elongation factor + real(r8) :: elongf_fnrt ! Fine-root "elongation factor" + real(r8) :: elongf_stem ! Stem "elongation factor" + real(r8) :: fnrt_drop_fraction ! Fraction of fine roots to absciss when leaves absciss + real(r8) :: stem_drop_fraction ! Fraction of stems to absciss when leaves absciss integer, parameter :: rstatus = 0 integer init @@ -787,36 +817,61 @@ subroutine init_cohorts( site_in, patch_in, bc_in) call bstore_allom(temp_cohort%dbh, pft, temp_cohort%canopy_trim, c_store) - temp_cohort%laimemory = 0._r8 - temp_cohort%sapwmemory = 0._r8 - temp_cohort%structmemory = 0._r8 - cstatus = leaves_on - + ! MLO update 20211123. The "memory" variables are now always set to the + ! on-allometry state (accounting for canopy trimming when necessary). + ! We no longer reset them to zero when leaves are flushing. This makes + ! partial deciduousness a bit easier to implement. + temp_cohort%leafmemory = c_leaf + temp_cohort%fnrtmemory = c_fnrt + temp_cohort%sapwmemory = c_sapw + temp_cohort%structmemory = c_struct + + ! Assume leaves are fully flushed, and update if needed. + cstatus = leaves_on + elongf_leaf = 1.0_r8 + elongf_fnrt = 1.0_r8 + elongf_stem = 1.0_r8 + + fnrt_drop_fraction = EDPftvarcon_inst%phen_fnrt_drop_fraction(temp_cohort%pft) stem_drop_fraction = EDPftvarcon_inst%phen_stem_drop_fraction(temp_cohort%pft) if(hlm_use_sp.eq.ifalse)then ! do not override SP vales with phenology if( prt_params%season_decid(pft) == itrue .and. & any(site_in%cstatus == [phen_cstat_nevercold,phen_cstat_iscold])) then - temp_cohort%laimemory = c_leaf - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction - c_leaf = 0._r8 - c_sapw = (1.0_r8-stem_drop_fraction) * c_sapw - c_struct = (1.0_r8-stem_drop_fraction) * c_struct - cstatus = leaves_off - endif + elongf_leaf = 0._r8 + elongf_fnrt = 1.0_r8 - fnrt_drop_fraction + elongf_stem = 1.0_r8 - stem_drop_fraction + + c_leaf = elongf_leaf * c_leaf + c_fnrt = elongf_fnrt * c_fnrt + c_sapw = elongf_stem * c_sapw + c_struct = elongf_stem * c_struct - if ( prt_params%stress_decid(pft) == itrue .and. & - any(site_in%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff])) then - temp_cohort%laimemory = c_leaf - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction - c_leaf = 0._r8 - c_sapw = (1.0_r8-stem_drop_fraction) * c_sapw - c_struct = (1.0_r8-stem_drop_fraction) * c_struct cstatus = leaves_off - endif + elseif ( prt_params%stress_decid(pft) == itrue) then + ! If the plant is drought deciduous, make sure leaf status is + ! always consistent with the leaf elongation factor. For tissues + ! other than leaves, the actual drop fraction is a combination of the + ! elongation factor (e) and the drop fraction (x), which will ensure + ! that the remaining tissue biomass will be exactly e when x=1, and + ! exactly the original biomass when x = 0. + elongf_leaf = site_in%elong_factor(pft) + elongf_fnrt = 1.0_r8 - (1.0_r8 - elongf_leaf ) * fnrt_drop_fraction + elongf_stem = 1.0_r8 - (1.0_r8 - elongf_leaf ) * stem_drop_fraction + + + c_leaf = elongf_leaf * c_leaf + c_fnrt = elongf_fnrt * c_fnrt + c_sapw = elongf_stem * c_sapw + c_struct = elongf_stem * c_struct + + if (elongf_leaf > 0.0_r8) then + cstatus = leaves_on + else + cstatus = leaves_off + end if + end if end if ! SP mode @@ -850,26 +905,26 @@ subroutine init_cohorts( site_in, patch_in, bc_in) case(nitrogen_element) - m_struct = c_struct*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(struct_organ)) - m_leaf = c_leaf*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(leaf_organ)) - m_fnrt = c_fnrt*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(fnrt_organ)) - m_sapw = c_sapw*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(sapw_organ)) + m_struct = c_struct*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(struct_organ)) + m_leaf = c_leaf*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(leaf_organ)) + m_fnrt = c_fnrt*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(fnrt_organ)) + m_sapw = c_sapw*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(sapw_organ)) m_repro = 0._r8 - m_store = StorageNutrientTarget(pft,element_id,m_leaf,m_fnrt,m_sapw,m_struct) + m_store = StorageNutrientTarget(pft,element_id,m_leaf,m_fnrt,m_sapw,m_struct) case(phosphorus_element) - m_struct = c_struct*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(struct_organ)) - m_leaf = c_leaf*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(leaf_organ)) - m_fnrt = c_fnrt*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(fnrt_organ)) - m_sapw = c_sapw*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(sapw_organ)) + m_struct = c_struct*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(struct_organ)) + m_leaf = c_leaf*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(leaf_organ)) + m_fnrt = c_fnrt*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(fnrt_organ)) + m_sapw = c_sapw*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(sapw_organ)) m_repro = 0._r8 - m_store = StorageNutrientTarget(pft,element_id,m_leaf,m_fnrt,m_sapw,m_struct) + m_store = StorageNutrientTarget(pft,element_id,m_leaf,m_fnrt,m_sapw,m_struct) end select select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) ! Put all of the leaf mass into the first bin call SetState(prt_obj,leaf_organ, element_id,m_leaf,1) @@ -893,9 +948,10 @@ subroutine init_cohorts( site_in, patch_in, bc_in) call prt_obj%CheckInitialConditions() call create_cohort(site_in, patch_in, pft, temp_cohort%n, temp_cohort%hite, & - temp_cohort%coage, temp_cohort%dbh, prt_obj, temp_cohort%laimemory, & - temp_cohort%sapwmemory, temp_cohort%structmemory, cstatus, rstatus, & - temp_cohort%canopy_trim, temp_cohort%c_area,1, site_in%spread, bc_in) + temp_cohort%coage, temp_cohort%dbh, prt_obj, temp_cohort%leafmemory, & + temp_cohort%fnrtmemory, temp_cohort%sapwmemory, temp_cohort%structmemory, & + elongf_leaf, elongf_fnrt, elongf_stem, cstatus, rstatus, & + temp_cohort%canopy_trim, temp_cohort%c_area, 1, site_in%spread, bc_in) deallocate(temp_cohort) ! get rid of temporary cohort diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 152911b59c..1335462755 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -27,6 +27,7 @@ module EDMainMod use FatesInterfaceTypesMod , only : numpft use FatesInterfaceTypesMod , only : hlm_use_nocomp use PRTGenericMod , only : prt_carbon_allom_hyp + use PRTGenericMod , only : prt_csimpler_allom_hyp use PRTGenericMod , only : prt_cnp_flex_allom_hyp use PRTGenericMod , only : nitrogen_element use PRTGenericMod , only : phosphorus_element @@ -415,7 +416,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) ! Conduct Maintenance Turnover (parteh) if(debug) call currentCohort%prt%CheckMassConservation(ft,3) - if(any(currentSite%dstatus == [phen_dstat_moiston,phen_dstat_timeon])) then + if(any(currentSite%dstatus(ft) == [phen_dstat_moiston,phen_dstat_timeon])) then is_drought = .false. else is_drought = .true. diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 51926d275a..87363f3071 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -36,7 +36,12 @@ module EDParamsMod integer,protected, public :: photo_tempsens_model ! switch for choosing the model that defines the temperature ! sensitivity of photosynthetic parameters (vcmax, jmax). ! 1=non-acclimating (NOT YET IMPLEMENTED) - + + integer,protected, public :: phen_drought_model ! Switch for chooshing the drought-deciduous phenology model + ! 0 = Default FATES ( leaves on/off based on 10-day average soil moisture threshold in rooting zone ) + ! 1 = ED2-like (partial deciduousness based on 10-day average soil moisture upper and lower threshold in rooting zone) + ! Other options coming at some point. + real(r8),protected, public :: fates_mortality_disturbance_fraction ! the fraction of canopy mortality that results in disturbance real(r8),protected, public :: ED_val_comp_excln real(r8),protected, public :: ED_val_vai_top_bin_width @@ -48,6 +53,7 @@ module EDParamsMod real(r8),protected, public :: ED_val_cwd_flig real(r8),protected, public :: ED_val_base_mr_20 real(r8),protected, public :: ED_val_phen_drought_threshold + real(r8),protected, public :: ED_val_phen_moist_threshold real(r8),protected, public :: ED_val_phen_doff_time real(r8),protected, public :: ED_val_phen_a real(r8),protected, public :: ED_val_phen_b @@ -106,6 +112,7 @@ module EDParamsMod character(len=param_string_length),parameter,public :: ED_name_photo_temp_acclim_timescale = "fates_photo_temp_acclim_timescale" character(len=param_string_length),parameter,public :: name_photo_tempsens_model = "fates_photo_tempsens_model" character(len=param_string_length),parameter,public :: name_maintresp_model = "fates_maintresp_model" + character(len=param_string_length),parameter,public :: name_phen_drought_model = "fates_phen_drought_model" character(len=param_string_length),parameter,public :: ED_name_hydr_htftype_node = "fates_hydr_htftype_node" character(len=param_string_length),parameter,public :: ED_name_mort_disturb_frac = "fates_mort_disturb_frac" character(len=param_string_length),parameter,public :: ED_name_comp_excln = "fates_comp_excln" @@ -118,6 +125,7 @@ module EDParamsMod character(len=param_string_length),parameter,public :: ED_name_cwd_flig= "fates_cwd_flig" character(len=param_string_length),parameter,public :: ED_name_base_mr_20= "fates_base_mr_20" character(len=param_string_length),parameter,public :: ED_name_phen_drought_threshold= "fates_phen_drought_threshold" + character(len=param_string_length),parameter,public :: ED_name_phen_moist_threshold= "fates_phen_moist_threshold" character(len=param_string_length),parameter,public :: ED_name_phen_doff_time= "fates_phen_doff_time" character(len=param_string_length),parameter,public :: ED_name_phen_a= "fates_phen_a" character(len=param_string_length),parameter,public :: ED_name_phen_b= "fates_phen_b" @@ -226,6 +234,7 @@ subroutine FatesParamsInit() photo_temp_acclim_timescale = nan photo_tempsens_model = -9 maintresp_model = -9 + phen_drought_model = -9 fates_mortality_disturbance_fraction = nan ED_val_comp_excln = nan ED_val_vai_top_bin_width = nan @@ -237,6 +246,7 @@ subroutine FatesParamsInit() ED_val_cwd_flig = nan ED_val_base_mr_20 = nan ED_val_phen_drought_threshold = nan + ED_val_phen_moist_threshold = nan ED_val_phen_doff_time = nan ED_val_phen_a = nan ED_val_phen_b = nan @@ -305,7 +315,10 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=name_maintresp_model,dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) - + + call fates_params%RegisterParameter(name=name_phen_drought_model,dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + call fates_params%RegisterParameter(name=name_theta_cj_c3, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -345,6 +358,9 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=ED_name_phen_drought_threshold, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) + call fates_params%RegisterParameter(name=ED_name_phen_moist_threshold, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + call fates_params%RegisterParameter(name=ED_name_phen_doff_time, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -487,7 +503,11 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=name_maintresp_model, & data=tmpreal) maintresp_model = nint(tmpreal) - + + call fates_params%RetreiveParameter(name=name_phen_drought_model, & + data=tmpreal) + phen_drought_model = nint(tmpreal) + call fates_params%RetreiveParameter(name=ED_name_mort_disturb_frac, & data=fates_mortality_disturbance_fraction) @@ -521,6 +541,9 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=ED_name_phen_drought_threshold, & data=ED_val_phen_drought_threshold) + call fates_params%RetreiveParameter(name=ED_name_phen_moist_threshold, & + data=ED_val_phen_moist_threshold) + call fates_params%RetreiveParameter(name=ED_name_phen_doff_time, & data=ED_val_phen_doff_time) @@ -647,6 +670,7 @@ subroutine FatesReceiveParams(fates_params) hydr_htftype_node(:) = nint(hydr_htftype_real(:)) deallocate(hydr_htftype_real) + end subroutine FatesReceiveParams ! ===================================================================================== @@ -667,6 +691,7 @@ subroutine FatesReportParams(is_master) write(fates_log(),fmt0) 'photo_temp_acclim_timescale = ',photo_temp_acclim_timescale write(fates_log(),fmti) 'hydr_htftype_node = ',hydr_htftype_node write(fates_log(),fmt0) 'fates_mortality_disturbance_fraction = ',fates_mortality_disturbance_fraction + write(fates_log(),fmti) 'phen_drought_model = ',phen_drought_model write(fates_log(),fmt0) 'ED_val_comp_excln = ',ED_val_comp_excln write(fates_log(),fmt0) 'ED_val_vai_top_bin_width = ',ED_val_vai_top_bin_width write(fates_log(),fmt0) 'ED_val_vai_width_increase_factor = ',ED_val_vai_width_increase_factor @@ -677,6 +702,7 @@ subroutine FatesReportParams(is_master) write(fates_log(),fmt0) 'ED_val_cwd_flig = ',ED_val_cwd_flig write(fates_log(),fmt0) 'ED_val_base_mr_20 = ', ED_val_base_mr_20 write(fates_log(),fmt0) 'ED_val_phen_drought_threshold = ',ED_val_phen_drought_threshold + write(fates_log(),fmt0) 'ED_val_phen_moist_threshold = ',ED_val_phen_moist_threshold write(fates_log(),fmt0) 'ED_val_phen_doff_time = ',ED_val_phen_doff_time write(fates_log(),fmt0) 'ED_val_phen_a = ',ED_val_phen_a write(fates_log(),fmt0) 'ED_val_phen_b = ',ED_val_phen_b diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index a149132a8d..bc72edf57d 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -8,6 +8,7 @@ module EDPftvarcon ! !USES: use EDTypesMod , only : maxSWb, ivis, inir use EDTypesMod , only : n_uptake_mode, p_uptake_mode + use EDTypesMod , only : drgt_phen_model_gradual use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : itrue, ifalse @@ -17,7 +18,9 @@ module EDPftvarcon use FatesLitterMod, only : ilabile,icellulose,ilignin use PRTGenericMod, only : leaf_organ, fnrt_organ, store_organ use PRTGenericMod, only : sapw_organ, struct_organ, repro_organ - use PRTGenericMod, only : prt_cnp_flex_allom_hyp,prt_carbon_allom_hyp + use PRTGenericMod, only : prt_csimpler_allom_hyp + use PRTGenericMod, only : prt_carbon_allom_hyp + use PRTGenericMod, only : prt_cnp_flex_allom_hyp use FatesInterfaceTypesMod, only : hlm_nitrogen_spec, hlm_phosphorus_spec use FatesInterfaceTypesMod, only : hlm_parteh_mode use FatesInterfaceTypesMod, only : hlm_nu_com @@ -167,7 +170,8 @@ module EDPftvarcon ! on bud-burst [kgC/kgC] real(r8), allocatable :: phen_cold_size_threshold(:) ! stem/leaf drop occurs on DBH size of decidious non-woody ! (coastal grass) plants larger than the threshold value - real(r8), allocatable :: phen_stem_drop_fraction(:) ! Fraction of stem dropped/senescened for decidious + real(r8), allocatable :: phen_fnrt_drop_fraction(:) ! Fraction of fine roots abscissed when leaves absciss + real(r8), allocatable :: phen_stem_drop_fraction(:) ! Fraction of stem abscissed when leaves absciss, for deciduous ! non-woody (grass) plants ! Nutrient Aquisition parameters @@ -570,6 +574,10 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_phen_fnrt_drop_fraction' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + ! Nutrient competition parameters @@ -917,6 +925,10 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%phen_stem_drop_fraction) + name = 'fates_phen_fnrt_drop_fraction' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%phen_fnrt_drop_fraction) + name = 'fates_prescribed_nuptake' call fates_params%RetreiveParameterAllocate(name=name, & data=this%prescribed_nuptake) @@ -1433,6 +1445,7 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'taus = ',EDPftvarcon_inst%taus write(fates_log(),fmt0) 'phenflush_fraction',EDpftvarcon_inst%phenflush_fraction write(fates_log(),fmt0) 'phen_cold_size_threshold = ',EDPftvarcon_inst%phen_cold_size_threshold + write(fates_log(),fmt0) 'phen_fnrt_drop_fraction',EDpftvarcon_inst%phen_fnrt_drop_fraction write(fates_log(),fmt0) 'phen_stem_drop_fraction',EDpftvarcon_inst%phen_stem_drop_fraction write(fates_log(),fmt0) 'fire_alpha_SH = ',EDPftvarcon_inst%fire_alpha_SH write(fates_log(),fmt0) 'allom_frbstor_repro = ',EDPftvarcon_inst%allom_frbstor_repro @@ -1479,8 +1492,14 @@ subroutine FatesCheckParams(is_master) ! ----------------------------------------------------------------------------------- use FatesConstantsMod , only : fates_check_param_set use FatesConstantsMod , only : itrue, ifalse - use EDParamsMod , only : logging_mechanical_frac, logging_collateral_frac, logging_direct_frac - use FatesInterfaceTypesMod , only : hlm_use_fixed_biogeog,hlm_use_sp + use EDParamsMod , only : logging_mechanical_frac + use EDParamsMod , only : logging_collateral_frac + use EDParamsMod , only : logging_direct_frac + use EDParamsMod , only : phen_drought_model + use EDParamsMod , only : ED_val_phen_drought_threshold + use EDParamsMod , only : ED_val_phen_moist_threshold + use FatesInterfaceTypesMod, only : hlm_use_fixed_biogeog + use FatesInterfaceTypesMod, only : hlm_use_sp ! Argument logical, intent(in) :: is_master ! Only log if this is the master proc @@ -1502,7 +1521,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 +1558,19 @@ subroutine FatesCheckParams(is_master) end if end if - elseif (hlm_parteh_mode .ne. prt_carbon_allom_hyp) then + case (prt_csimpler_allom_hyp,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(),*) 'only 3 modules supported, simpler allometric carbon,' + write(fates_log(),*) 'default allometric carbon, and CNP.' + write(fates_log(),*) 'fates_parteh_mode must be set to 0, 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 ! @@ -1677,8 +1702,18 @@ subroutine FatesCheckParams(is_master) write(fates_log(),*) ' Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if + if ( ( EDPftvarcon_inst%phen_fnrt_drop_fraction(ipft) < 0.0_r8 ) .or. & + ( EDPFtvarcon_inst%phen_fnrt_drop_fraction(ipft) > 1.0_r8 ) ) then + write(fates_log(),*) ' Abscission rate for fine roots must be between 0 and 1 for ' + write(fates_log(),*) ' deciduous PFTs.' + write(fates_log(),*) ' PFT#: ',ipft + write(fates_log(),*) ' evergreen flag: (shold be 0):',int(prt_params%evergreen(ipft)) + write(fates_log(),*) ' phen_fnrt_drop_fraction: ', EDPFtvarcon_inst%phen_fnrt_drop_fraction(ipft) + write(fates_log(),*) ' Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if if ( ( EDPftvarcon_inst%phen_stem_drop_fraction(ipft) < 0.0_r8 ) .or. & - ( EDPFtvarcon_inst%phen_stem_drop_fraction(ipft) > 1 ) ) then + ( EDPFtvarcon_inst%phen_stem_drop_fraction(ipft) > 1.0_r8 ) ) then write(fates_log(),*) ' Deciduous non-wood plants must keep 0-100% of their stems' write(fates_log(),*) ' during the deciduous period.' write(fates_log(),*) ' PFT#: ',ipft @@ -1784,6 +1819,40 @@ subroutine FatesCheckParams(is_master) !! end do + ! When using the the gradual (ED2-like) phenology, we must ensure that the lower + ! and upper thresholds are consistent (i.e., that both are based on either soil + ! water content or soil matric potential). + select case (phen_drought_model) + case(drgt_phen_model_gradual) + + if (ED_val_phen_drought_threshold*ED_val_phen_moist_threshold < 0._r8) then + ! In case the product of the lower and upper thresholds is negative, the + ! thresholds are inconsistent as both should be defined using the same + ! quantity. + write(fates_log(),*) ' When using gradual (ED2-like) drought deciduous phenology,' + write(fates_log(),*) ' the moist threshold should be have the same sign as' + write(fates_log(),*) ' the dry threshold. Positive = soil water content [m3/m3],' + write(fates_log(),*) ' Negative = soil matric potential [mm].' + write(fates_log(),*) ' fates_phen_drought_threshold (dry threshold) = ',ED_val_phen_drought_threshold + write(fates_log(),*) ' fates_phen_moist_threshold (moist threshold) = ',ED_val_phen_moist_threshold + write(fates_log(),*) ' Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + elseif ( ED_val_phen_drought_threshold >= ED_val_phen_moist_threshold) then + write(fates_log(),*) ' When using gradual (ED2-like) drought deciduous phenology,' + write(fates_log(),*) ' the moist threshold should be greater than the dry threshold.' + write(fates_log(),*) ' By greater we mean more positive or less negative, and' + write(fates_log(),*) ' they cannot be the identical.' + write(fates_log(),*) ' fates_phen_drought_threshold (dry threshold) = ',ED_val_phen_drought_threshold + write(fates_log(),*) ' fates_phen_moist_threshold (moist threshold) = ',ED_val_phen_moist_threshold + write(fates_log(),*) ' Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + end select + + + return end subroutine FatesCheckParams diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 9ce4d5fe44..b7fd6b6083 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -33,7 +33,7 @@ module EDTypesMod (/ 10, 4 /) !!! MUST SUM TO maxPatchesPerSite !!! integer, public :: maxCohortsPerPatch = 100 ! maximum number of cohorts per patch - integer, parameter, public :: nclmax = 2 ! Maximum number of canopy layers + integer, parameter, public :: nclmax = 3 ! Maximum number of canopy layers integer, parameter, public :: ican_upper = 1 ! Nominal index for the upper canopy integer, parameter, public :: ican_ustory = 2 ! Nominal index for diagnostics that refer ! to understory layers (all layers that @@ -90,12 +90,23 @@ module EDTypesMod ! can be approximated to be equal to the visible band + integer, parameter, public :: leaves_pshed = 3 ! Flag specifying that a deciduous plant has leaves + ! but is shedding them (partial shedding). This plant + ! should not allocate carbon towards growth or + ! reproduction. integer, parameter, public :: leaves_on = 2 ! Flag specifying that a deciduous plant has leaves ! and should be allocating to them as well integer, parameter, public :: leaves_off = 1 ! Flag specifying that a deciduous plant has dropped ! its leaves and should not be trying to allocate ! towards any growth. + integer, parameter, public :: drgt_phen_model_smoist = 0 ! Switch for the default FATES drought + ! deciduous phenology. + integer, parameter, public :: drgt_phen_model_gradual = 1 ! Switch for the ED2-like drought deciduous + ! phenology, which allows for gradual + ! abscission and flushing by using two + ! thresholds. + ! Flag to turn on/off salinity effects on the effective "btran" ! btran stress function. @@ -141,10 +152,11 @@ module EDTypesMod integer, parameter, public :: phen_cstat_iscold = 1 ! This (location/plant) is in a cold-state where leaves should have fallen integer, parameter, public :: phen_cstat_notcold = 2 ! This site is in a warm-state where leaves are allowed to flush - integer, parameter, public :: phen_dstat_timeoff = 0 ! Leaves off due to time exceedance (drought phenology) - integer, parameter, public :: phen_dstat_moistoff = 1 ! Leaves off due to moisture avail (drought phenology) - integer, parameter, public :: phen_dstat_moiston = 2 ! Leaves on due to moisture avail (drought phenology) - integer, parameter, public :: phen_dstat_timeon = 3 ! Leaves on due to time exceedance (drought phenology) + integer, parameter, public :: phen_dstat_timeoff = 0 ! Leaves off due to time exceedance (drought phenology) + integer, parameter, public :: phen_dstat_moistoff = 1 ! Leaves off due to moisture avail (drought phenology) + integer, parameter, public :: phen_dstat_moiston = 2 ! Leaves on due to moisture avail (drought phenology) + integer, parameter, public :: phen_dstat_timeon = 3 ! Leaves on due to time exceedance (drought phenology) + integer, parameter, public :: phen_dstat_pshed = 4 ! Leaves partially abscissing (drought phenology) ! SPITFIRE @@ -219,7 +231,8 @@ module EDTypesMod real(r8) :: coage ! cohort age in years real(r8) :: hite ! height: meters integer :: indexnumber ! unique number for each cohort. (within clump?) - real(r8) :: laimemory ! target leaf biomass- set from previous year: kGC per indiv + real(r8) :: leafmemory ! target leaf biomass- set from previous year: kGC per indiv + real(r8) :: fnrtmemory ! target fine-root biomass- set from previous year: kGC per indiv real(r8) :: sapwmemory ! target sapwood biomass- set from previous year: kGC per indiv real(r8) :: structmemory ! target structural biomass- set from previous year: kGC per indiv integer :: canopy_layer ! canopy status of cohort (1 = canopy, 2 = understorey, etc.) @@ -236,6 +249,11 @@ module EDTypesMod real(r8) :: prom_weight ! How much of this cohort is promoted each year, as a proportion of all cohorts:- integer :: nv ! Number of leaf layers: - integer :: status_coh ! growth status of plant (2 = leaves on , 1 = leaves off) + real(r8) :: efleaf_coh ! Elongation factor for leaves (fraction) + real(r8) :: effnrt_coh ! Elongation factor for fine roots (fraction) + real(r8) :: efstem_coh ! Elongation factor for stem (fraction) + ! For all the elongation factors, zero means fully abscissed, and + ! one means fully flushed. real(r8) :: c_area ! areal extent of canopy (m2) real(r8) :: treelai ! lai of an individual within cohort leaf area (m2) / crown area (m2) real(r8) :: treesai ! stem area index of an indiv. within cohort: stem area (m2) / crown area (m2) @@ -737,20 +755,28 @@ module EDTypesMod ! 400 days, leaves are dropped and flagged as non-cold region ! 1 = this site is in a cold-state where leaves should have fallen ! 2 = this site is in a warm-state where leaves are allowed to flush - integer :: dstatus ! are leaves in this pixel on or off for drought decid + integer :: dstatus(maxpft) ! are leaves in this pixel on or off for drought decid ! 0 = leaves off due to time exceedance ! 1 = leaves off due to moisture avail ! 2 = leaves on due to moisture avail ! 3 = leaves on due to time exceedance + ! 4 = leaves partially on (ED2-like phenology) integer :: nchilldays ! num chilling days: (for botta gdd trheshold calculation) integer :: ncolddays ! num cold days: (must exceed threshold to drop leaves) real(r8) :: vegtemp_memory(num_vegtemp_mem) ! record of last 10 days temperature for senescence model. deg C integer :: cleafondate ! model date (day integer) of leaf on (cold):- integer :: cleafoffdate ! model date (day integer) of leaf off (cold):- - integer :: dleafondate ! model date (day integer) of leaf on drought:- - integer :: dleafoffdate ! model date (day integer) of leaf off drought:- + integer :: cndaysleafon ! number of days since leaf on period started (cold) + integer :: cndaysleafoff ! number of days since leaf off period started (cold) + integer :: dleafondate(maxpft) ! model date (day integer) of leaf on drought:- + integer :: dleafoffdate(maxpft) ! model date (day integer) of leaf off drought:- + integer :: dndaysleafon(maxpft) ! number of days since leaf on period started (drought) + integer :: dndaysleafoff(maxpft) ! number of days since leaf off period started (drought) + real(r8) :: elong_factor(maxpft) ! Elongation factor (ED2-like phenology). This is zero when leaves are + ! completely off, and one when they are completely flushed. - real(r8) :: water_memory(numWaterMem) ! last 10 days of soil moisture memory... + real(r8) :: liqvol_memory(numWaterMem,maxpft) ! last 10 days of soil liquid water volume (drought phenology) + real(r8) :: smp_memory(numWaterMem,maxpft) ! last 10 days of soil matric potential (drought phenology) ! FIRE @@ -1049,7 +1075,8 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%dbh = ', ccohort%dbh write(fates_log(),*) 'co%hite = ', ccohort%hite write(fates_log(),*) 'co%coage = ', ccohort%coage - write(fates_log(),*) 'co%laimemory = ', ccohort%laimemory + write(fates_log(),*) 'co%leafmemory = ', ccohort%leafmemory + write(fates_log(),*) 'co%fnrtmemory = ', ccohort%fnrtmemory write(fates_log(),*) 'co%sapwmemory = ', ccohort%sapwmemory write(fates_log(),*) 'co%structmemory = ', ccohort%structmemory @@ -1068,6 +1095,9 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%canopy_layer_yesterday = ', ccohort%canopy_layer_yesterday write(fates_log(),*) 'co%nv = ', ccohort%nv write(fates_log(),*) 'co%status_coh = ', ccohort%status_coh + write(fates_log(),*) 'co%efleaf_coh = ', ccohort%efleaf_coh + write(fates_log(),*) 'co%effnrt_coh = ', ccohort%effnrt_coh + write(fates_log(),*) 'co%efstem_coh = ', ccohort%efstem_coh write(fates_log(),*) 'co%canopy_trim = ', ccohort%canopy_trim write(fates_log(),*) 'co%excl_weight = ', ccohort%excl_weight write(fates_log(),*) 'co%prom_weight = ', ccohort%prom_weight diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 726100a37b..693d892b06 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -210,7 +210,7 @@ module FatesConstantsMod ! Approximate molar mass of water vapor to dry air (-) real(fates_r8), parameter, public :: molar_mass_ratio_vapdry= 0.622_fates_r8 - ! Gravity constant on earth [m/s] + ! Gravity constant on earth [m/s^2] real(fates_r8), parameter, public :: grav_earth = 9.8_fates_r8 ! Megapascals to pascals diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index f43289da15..bdacf20236 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -67,9 +67,12 @@ module FatesHistoryInterfaceMod use FatesConstantsMod , only : J_per_kJ use FatesConstantsMod , only : m2_per_ha use FatesConstantsMod , only : m_per_cm + use FatesConstantsMod , only : mm_per_m use FatesConstantsMod , only : sec_per_min use FatesConstantsMod , only : umol_per_mol use FatesConstantsMod , only : pa_per_mpa + use FatesConstantsMod , only : dens_fresh_liquid_water + use FatesConstantsMod , only : grav_earth use FatesLitterMod , only : litter_type use FatesConstantsMod , only : secondaryforest @@ -78,6 +81,7 @@ module FatesHistoryInterfaceMod use PRTGenericMod , only : all_carbon_elements use PRTGenericMod , only : carbon12_element use PRTGenericMod , only : nitrogen_element, phosphorus_element + use PRTGenericMod , only : prt_csimpler_allom_hyp use PRTGenericMod , only : prt_carbon_allom_hyp implicit none @@ -221,6 +225,9 @@ module FatesHistoryInterfaceMod integer :: ih_bleaf_canopy_si_scpf integer :: ih_bleaf_understory_si_scpf + ! Size-class x PFT LAI states + integer :: ih_lai_canopy_si_scpf + integer :: ih_lai_understory_si_scpf integer :: ih_totvegn_scpf @@ -341,15 +348,11 @@ module FatesHistoryInterfaceMod integer :: ih_h2oveg_hydro_err_si integer :: ih_site_cstatus_si - integer :: ih_site_dstatus_si integer :: ih_gdd_si integer :: ih_site_nchilldays_si integer :: ih_site_ncolddays_si integer :: ih_cleafoff_si integer :: ih_cleafon_si - integer :: ih_dleafoff_si - integer :: ih_dleafon_si - integer :: ih_meanliqvol_si integer :: ih_nesterov_fire_danger_si integer :: ih_fire_nignitions_si @@ -514,6 +517,12 @@ module FatesHistoryInterfaceMod integer :: ih_canopycrownarea_si_pft integer :: ih_gpp_si_pft integer :: ih_npp_si_pft + integer :: ih_site_dstatus_si_pft + integer :: ih_dleafoff_si_pft + integer :: ih_dleafon_si_pft + integer :: ih_meanliqvol_si_pft + integer :: ih_meansmp_si_pft + integer :: ih_elong_factor_si_pft integer :: ih_nocomp_pftpatchfraction_si_pft integer :: ih_nocomp_pftnpatches_si_pft integer :: ih_nocomp_pftburnedarea_si_pft @@ -1861,6 +1870,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_bstor_understory_si_scpf => this%hvars(ih_bstor_understory_si_scpf)%r82d, & hio_bleaf_canopy_si_scpf => this%hvars(ih_bleaf_canopy_si_scpf)%r82d, & hio_bleaf_understory_si_scpf => this%hvars(ih_bleaf_understory_si_scpf)%r82d, & + hio_lai_canopy_si_scpf => this%hvars(ih_lai_canopy_si_scpf)%r82d, & + hio_lai_understory_si_scpf => this%hvars(ih_lai_understory_si_scpf)%r82d, & hio_mortality_canopy_si_scpf => this%hvars(ih_mortality_canopy_si_scpf)%r82d, & hio_mortality_understory_si_scpf => this%hvars(ih_mortality_understory_si_scpf)%r82d, & hio_nplant_canopy_si_scpf => this%hvars(ih_nplant_canopy_si_scpf)%r82d, & @@ -2011,16 +2022,18 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_mortality_canopy_si_scag => this%hvars(ih_mortality_canopy_si_scag)%r82d, & hio_mortality_understory_si_scag => this%hvars(ih_mortality_understory_si_scag)%r82d, & hio_site_cstatus_si => this%hvars(ih_site_cstatus_si)%r81d, & - hio_site_dstatus_si => this%hvars(ih_site_dstatus_si)%r81d, & hio_gdd_si => this%hvars(ih_gdd_si)%r81d, & hio_site_ncolddays_si => this%hvars(ih_site_ncolddays_si)%r81d, & hio_site_nchilldays_si => this%hvars(ih_site_nchilldays_si)%r81d, & hio_cleafoff_si => this%hvars(ih_cleafoff_si)%r81d, & hio_cleafon_si => this%hvars(ih_cleafon_si)%r81d, & - hio_dleafoff_si => this%hvars(ih_dleafoff_si)%r81d, & - hio_dleafon_si => this%hvars(ih_dleafoff_si)%r81d, & + hio_site_dstatus_si_pft => this%hvars(ih_site_dstatus_si_pft)%r82d, & + hio_dleafoff_si_pft => this%hvars(ih_dleafoff_si_pft)%r82d, & + hio_dleafon_si_pft => this%hvars(ih_dleafon_si_pft)%r82d, & + hio_meanliqvol_si_pft => this%hvars(ih_meanliqvol_si_pft)%r82d, & + hio_meansmp_si_pft => this%hvars(ih_meansmp_si_pft)%r82d, & + hio_elong_factor_si_pft => this%hvars(ih_elong_factor_si_pft)%r82d, & hio_tveg24 => this%hvars(ih_tveg24_si)%r81d, & - hio_meanliqvol_si => this%hvars(ih_meanliqvol_si)%r81d, & hio_cbal_err_fates_si => this%hvars(ih_cbal_err_fates_si)%r81d, & hio_err_fates_si => this%hvars(ih_err_fates_si)%r82d ) @@ -2061,10 +2074,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! Canopy spread index (0-1) hio_canopy_spread_si(io_si) = sites(s)%spread - ! Site statuses (stati?) for cold deciduous and drought - ! deciduous - hio_site_cstatus_si(io_si) = real(sites(s)%cstatus,r8) - hio_site_dstatus_si(io_si) = real(sites(s)%dstatus,r8) + ! Update the site status for cold deciduous (drought deciduous is now PFT dependent) + hio_site_cstatus_si(io_si) = real(sites(s)%cstatus,r8) ! Number of chill days and cold days hio_site_nchilldays_si(io_si) = real(sites(s)%nchilldays,r8) @@ -2073,17 +2084,34 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! Growing degree-days hio_gdd_si(io_si) = sites(s)%grow_deg_days - ! Model days elapsed since leaf on/off for cold- and drought-deciduous + ! Model days elapsed since leaf off/on for cold deciduous hio_cleafoff_si(io_si) = real(model_day_int - sites(s)%cleafoffdate,r8) hio_cleafon_si(io_si) = real(model_day_int - sites(s)%cleafondate,r8) - hio_dleafoff_si(io_si) = real(model_day_int - sites(s)%dleafoffdate,r8) - hio_dleafon_si(io_si) = real(model_day_int - sites(s)%dleafondate,r8) - ! Mean liquid water content (m3/m3) used for drought phenology - if(model_day_int>numWaterMem)then - hio_meanliqvol_si(io_si) = & - sum(sites(s)%water_memory(1:numWaterMem))/real(numWaterMem,r8) - end if + + ! Update drought deciduous information (now separated by PFT). + do i_pft = 1,numpft + ! Update the site-PFT status for drought deciduous + hio_site_dstatus_si_pft(io_si,i_pft) = real(sites(s)%dstatus(i_pft),r8) + + ! Model days elapsed since leaf off/on for drought deciduous + hio_dleafoff_si_pft(io_si,i_pft) = real(sites(s)%dndaysleafon (i_pft),r8) + hio_dleafon_si_pft(io_si,i_pft) = real(sites(s)%dndaysleafoff(i_pft),r8) + + ! Leaf elongation factor (0 means fully abscissed, 1 means fully flushed). + hio_elong_factor_si_pft(io_si,i_pft) = sites(s)%elong_factor(i_pft) + + if(model_day_int>numWaterMem)then + ! Mean liquid water content (m3/m3) used for drought phenology + hio_meanliqvol_si_pft(io_si,i_pft) = & + sum(sites(s)%liqvol_memory(1:numWaterMem,i_pft))/real(numWaterMem,r8) + + ! Mean soil matric potential (Pa) used for drought phenology + hio_meansmp_si_pft(io_si,i_pft) = & + sum(sites(s)%smp_memory(1:numWaterMem,i_pft))/real(numWaterMem,r8) & + * dens_fresh_liquid_water * grav_earth * mm_per_m + end if + end do ! track total wood product accumulation at the site level hio_woodproduct_si(io_si) = sites(s)%resources_management%trunk_product_site & @@ -4044,7 +4072,7 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,bc_in,dt_tstep) do j_bc = j_t,j_b vwc = bc_in(s)%h2o_liqvol_sl(j_bc) - psi = site_hydr%wrf_soil(j)%p%psi_from_th(vwc) + psi = site_hydr%wrf_soil(j)%p%psi_from_th(vwc) ! MLO: Any reason for not using smp_sl? ! cap capillary pressure ! psi = max(-1e5_r8,psi) Removing cap as that is inconstistent ! with model internals and physics. Should @@ -4397,13 +4425,6 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_site_cstatus_si) - call this%set_history_var(vname='FATES_DROUGHT_STATUS', & - units='', & - long='site-level drought status, <2 too dry for leaves, >=2 not too dry', & - use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & - upfreq=1, ivar=ivar, initialize=initialize_variables, & - index=ih_site_dstatus_si) - call this%set_history_var(vname='FATES_GDD', units='degree_Celsius', & long='site-level growing degree days', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', & @@ -4433,26 +4454,6 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_cleafon_si) - call this%set_history_var(vname='FATES_DAYSINCE_DROUGHTLEAFOFF', & - units='days', long='site level days elapsed since drought leaf drop', & - use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & - upfreq=1, ivar=ivar, initialize=initialize_variables, & - index=ih_dleafoff_si) - - call this%set_history_var(vname='FATES_DAYSINCE_DROUGHTLEAFON', & - units='days', & - long='site-level days elapsed since drought leaf flush', & - use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & - upfreq=1, ivar=ivar, initialize=initialize_variables, & - index=ih_dleafon_si) - - call this%set_history_var(vname='FATES_MEANLIQVOL_DROUGHTPHEN', & - units='m3 m-3', & - long='site-level mean liquid water volume for drought phenolgy', & - use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & - upfreq=1, ivar=ivar, initialize=initialize_variables, & - index=ih_meanliqvol_si) - call this%set_history_var(vname='FATES_CANOPY_SPREAD', units='', & long='scaling factor (0-1) between tree basal area and canopy area', & use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & @@ -4520,6 +4521,49 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_mortality_si_pft) + !MLO - Drought-deciduous phenology variables are now defined for each PFT. + call this%set_history_var(vname='FATES_DROUGHT_STATUS_PF', & + units='', & + long='PFT-level drought status, <2 too dry for leaves, >=2 not too dry', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_site_dstatus_si_pft) + + call this%set_history_var(vname='FATES_DAYSINCE_DROUGHTLEAFOFF_PF', & + units='days', long='PFT-level days elapsed since drought leaf drop', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_dleafoff_si_pft) + + call this%set_history_var(vname='FATES_DAYSINCE_DROUGHTLEAFON_PF', & + units='days', & + long='PFT-level days elapsed since drought leaf flush', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_dleafon_si_pft) + + call this%set_history_var(vname='FATES_MEANLIQVOL_DROUGHTPHEN_PF', & + units='m3 m-3', & + long='PFT-level mean liquid water volume for drought phenolgy', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_meanliqvol_si_pft) + + call this%set_history_var(vname='FATES_MEANSMP_DROUGHTPHEN_PF', & + units='Pa', & + long='PFT-level mean soil matric potential for drought phenology', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_meansmp_si_pft) + + call this%set_history_var(vname='FATES_ELONG_FACTOR_PF', & + units='1', & + long='PFT-level mean elongation factor (partial flushing/abscission)', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_elong_factor_si_pft) + + nocomp_if: if (hlm_use_nocomp .eq. itrue) then call this%set_history_var(vname='FATES_NOCOMP_NPATCHES_PF', units='', & long='number of patches per PFT (nocomp-mode-only)', & diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 3ce7589f48..d9be9c43f6 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -65,6 +65,7 @@ module FatesInterfaceMod use PRTGenericMod , only : element_list use PRTGenericMod , only : element_pos use EDParamsMod , only : eca_plant_escalar + use PRTGenericMod , only : prt_csimpler_allom_hyp use PRTGenericMod , only : prt_carbon_allom_hyp use PRTGenericMod , only : prt_cnp_flex_allom_hyp use PRTGenericMod , only : carbon12_element @@ -306,7 +307,7 @@ subroutine zero_bcs(fates,s) ! Fates -> BGC fragmentation mass fluxes select case(hlm_parteh_mode) - case(prt_carbon_allom_hyp) + case(prt_csimpler_allom_hyp,prt_carbon_allom_hyp) fates%bc_out(s)%litt_flux_cel_c_si(:) = 0._r8 fates%bc_out(s)%litt_flux_lig_c_si(:) = 0._r8 fates%bc_out(s)%litt_flux_lab_c_si(:) = 0._r8 @@ -627,7 +628,7 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in) ! Fates -> BGC fragmentation mass fluxes select case(hlm_parteh_mode) - case(prt_carbon_allom_hyp) + case(prt_csimpler_allom_hyp,prt_carbon_allom_hyp) allocate(bc_out%litt_flux_cel_c_si(nlevdecomp_in)) allocate(bc_out%litt_flux_lig_c_si(nlevdecomp_in)) allocate(bc_out%litt_flux_lab_c_si(nlevdecomp_in)) @@ -939,7 +940,7 @@ subroutine InitPARTEHGlobals() ! automatically. select case(hlm_parteh_mode) - case(prt_carbon_allom_hyp) + case(prt_csimpler_allom_hyp,prt_carbon_allom_hyp) num_elements = 1 allocate(element_list(num_elements)) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 507f01dbee..42a57949e5 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -48,6 +48,7 @@ module FatesInventoryInitMod use EDPftvarcon , only : EDPftvarcon_inst use FatesInterfaceTypesMod, only : hlm_parteh_mode use EDCohortDynamicsMod, only : InitPRTObject + use PRTGenericMod, only : prt_csimpler_allom_hyp use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes @@ -915,7 +916,11 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & real(r8) :: m_sapw ! Generic mass for sapwood [kg] real(r8) :: m_store ! Generic mass for storage [kg] real(r8) :: m_repro ! Generic mass for reproductive tissues [kg] - real(r8) :: stem_drop_fraction + real(r8) :: elongf_leaf ! Leaf elongation factor + real(r8) :: elongf_fnrt ! Fine-root "elongation factor" + real(r8) :: elongf_stem ! Stem "elongation factor" + real(r8) :: fnrt_drop_fraction ! Fine-root abscission fraction + real(r8) :: stem_drop_fraction ! Stem abscission fraction integer :: i_pft, ncohorts_to_create character(len=128),parameter :: wr_fmt = & @@ -1039,34 +1044,63 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & call bstore_allom(temp_cohort%dbh, temp_cohort%pft, temp_cohort%canopy_trim, c_store) - temp_cohort%laimemory = 0._r8 - temp_cohort%sapwmemory = 0._r8 - temp_cohort%structmemory = 0._r8 + + ! MLO update 20211123. The "memory" variables are now always set to the on-allometry state + ! (accounting for canopy trimming when necessary). We no longer reset them to zero when + ! leaves are flushing. This makes partial deciduousness a bit easier to implement. + temp_cohort%leafmemory = c_leaf + temp_cohort%fnrtmemory = c_fnrt + temp_cohort%sapwmemory = c_sapw + temp_cohort%structmemory = c_struct + cstatus = leaves_on + elongf_leaf = 1.0_r8 + elongf_fnrt = 1.0_r8 + elongf_stem = 1.0_r8 + fnrt_drop_fraction = EDPftvarcon_inst%phen_fnrt_drop_fraction(temp_cohort%pft) stem_drop_fraction = EDPftvarcon_inst%phen_stem_drop_fraction(temp_cohort%pft) if( prt_params%season_decid(temp_cohort%pft) == itrue .and. & any(csite%cstatus == [phen_cstat_nevercold,phen_cstat_iscold])) then - temp_cohort%laimemory = c_leaf - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction - c_leaf = 0._r8 - c_sapw = (1._r8 - stem_drop_fraction) * c_sapw - c_struct = (1._r8 - stem_drop_fraction) * c_struct - cstatus = leaves_off - endif - - if ( prt_params%stress_decid(temp_cohort%pft) == itrue .and. & - any(csite%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff])) then - temp_cohort%laimemory = c_leaf - temp_cohort%sapwmemory = c_sapw * stem_drop_fraction - temp_cohort%structmemory = c_struct * stem_drop_fraction - c_leaf = 0._r8 - c_sapw = (1._r8 - stem_drop_fraction) * c_sapw - c_struct = (1._r8 - stem_drop_fraction) * c_struct - cstatus = leaves_off - endif + elongf_leaf = 0.0_r8 + elongf_fnrt = 1._r8 - fnrt_drop_fraction + elongf_stem = 1._r8 - stem_drop_fraction + + c_leaf = elongf_leaf * c_leaf + c_fnrt = elongf_fnrt * c_fnrt + c_sapw = elongf_stem * c_sapw + c_struct = elongf_stem * c_struct + + cstatus = leaves_off + elseif ( prt_params%stress_decid(temp_cohort%pft) == itrue ) then + ! Drought deciduous. For the default approach, elongation factor is either + ! zero (full abscission) or one (fully flushed), but this can also be a + ! fraction in other approaches. Here we assume that leaves are "on" (i.e. + ! either fully flushed or growing) if elongation factor is not 0 for the + ! initial conditions. + ! + ! For tissues other than leaves, the actual drop fraction is a combination + ! of the elongation factor (e) and the drop fraction (x), which will ensure + ! that the remaining tissue biomass will be exactly e when x=1, and exactly + ! the original biomass when x = 0. + elongf_leaf = csite%elong_factor(temp_cohort%pft) + elongf_fnrt = 1.0_r8 - (1.0_r8 - elongf_leaf ) * fnrt_drop_fraction + elongf_stem = 1.0_r8 - (1.0_r8 - elongf_leaf ) * stem_drop_fraction + + + c_leaf = elongf_leaf * c_leaf + c_fnrt = elongf_fnrt * c_fnrt + c_sapw = elongf_stem * c_sapw + c_struct = elongf_stem * c_struct + if (elongf_leaf > 0.0_r8) then + ! Assume leaves are growing even if they are not fully flushed. + cstatus = leaves_on + else + ! Leaves are off (abscissing). + cstatus = leaves_off + end if + end if prt_obj => null() call InitPRTObject(prt_obj) @@ -1140,7 +1174,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & end select select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) ! Equally distribute leaf mass into available age-bins do iage = 1,nleafage @@ -1167,8 +1201,9 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & call create_cohort(csite, cpatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, & temp_cohort%coage, temp_cohort%dbh, & - prt_obj, temp_cohort%laimemory,temp_cohort%sapwmemory, temp_cohort%structmemory, & - cstatus, rstatus, temp_cohort%canopy_trim,temp_cohort%c_area, & + prt_obj, temp_cohort%leafmemory, temp_cohort%fnrtmemory, & + temp_cohort%sapwmemory, temp_cohort%structmemory, elongf_leaf, elongf_fnrt, & + elongf_stem, cstatus, rstatus, temp_cohort%canopy_trim,temp_cohort%c_area, & 1, csite%spread, bc_in) deallocate(temp_cohort) ! get rid of temporary cohort diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 0605767cd6..dd676a4ba9 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -80,13 +80,12 @@ module FatesRestartInterfaceMod integer :: ir_npatch_si integer :: ir_cd_status_si - integer :: ir_dd_status_si integer :: ir_nchill_days_si integer :: ir_ncold_days_si integer :: ir_leafondate_si integer :: ir_leafoffdate_si - integer :: ir_dleafondate_si - integer :: ir_dleafoffdate_si + integer :: ir_cndaysleafon_si + integer :: ir_cndaysleafoff_si integer :: ir_acc_ni_si integer :: ir_gdd_si integer :: ir_snow_depth_si @@ -100,7 +99,8 @@ module FatesRestartInterfaceMod integer :: ir_coage_co integer :: ir_g_sb_laweight_co integer :: ir_height_co - integer :: ir_laimemory_co + integer :: ir_leafmemory_co + integer :: ir_fnrtmemory_co integer :: ir_sapwmemory_co integer :: ir_structmemory_co integer :: ir_nplant_co @@ -156,6 +156,9 @@ module FatesRestartInterfaceMod integer :: ir_resp_tstep_co integer :: ir_pft_co integer :: ir_status_co + integer :: ir_efleaf_co + integer :: ir_effnrt_co + integer :: ir_efstem_co integer :: ir_isnew_co ! Litter @@ -187,7 +190,14 @@ module FatesRestartInterfaceMod integer :: ir_litter_moisture_pa_nfsc ! Site level - integer :: ir_watermem_siwm + integer :: ir_dd_status_sift + integer :: ir_dleafondate_sift + integer :: ir_dleafoffdate_sift + integer :: ir_dndaysleafon_sift + integer :: ir_dndaysleafoff_sift + integer :: ir_elong_factor_sift + integer :: ir_liqvolmem_siwmft + integer :: ir_smpmem_siwmft integer :: ir_vegtempmem_sitm integer :: ir_seed_bank_sift integer :: ir_spread_si @@ -606,10 +616,6 @@ subroutine define_restart_vars(this, initialize_variables) long_name='status flag for cold deciduous plants', units='unitless', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_cd_status_si ) - call this%set_restart_var(vname='fates_drought_dec_status', vtype=site_int, & - long_name='status flag for drought deciduous plants', units='unitless', flushval = flushinvalid, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_dd_status_si ) - call this%set_restart_var(vname='fates_chilling_days', vtype=site_int, & long_name='chilling day counter', units='unitless', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_nchill_days_si ) @@ -619,20 +625,20 @@ subroutine define_restart_vars(this, initialize_variables) hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_ncold_days_si ) call this%set_restart_var(vname='fates_leafondate', vtype=site_int, & - long_name='the day of year for leaf on', units='day of year', flushval = flushinvalid, & + long_name='the day of year for leaf on (cold deciduous)', units='day of year', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_leafondate_si ) call this%set_restart_var(vname='fates_leafoffdate', vtype=site_int, & - long_name='the day of year for leaf off', units='day of year', flushval = flushinvalid, & + long_name='the day of year for leaf off (cold deciduous)', units='day of year', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_leafoffdate_si ) - call this%set_restart_var(vname='fates_drought_leafondate', vtype=site_int, & - long_name='the day of year for drought based leaf-on', units='day of year', flushval = flushinvalid, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_dleafondate_si ) + call this%set_restart_var(vname='fates_ndaysleafon', vtype=site_int, & + long_name='number of days since leaf on (cold deciduous)', units='days', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_cndaysleafon_si ) - call this%set_restart_var(vname='fates_drought_leafoffdate', vtype=site_int, & - long_name='the day of year for drought based leaf-off', units='day of year', flushval = flushinvalid, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_dleafoffdate_si ) + call this%set_restart_var(vname='fates_ndaysleafoff', vtype=site_int, & + long_name='number of days since leaf off (cold deciduous)', units='days', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_cndaysleafoff_si ) call this%set_restart_var(vname='fates_acc_nesterov_id', vtype=site_r8, & long_name='a nesterov index accumulator', units='unitless', flushval = flushzero, & @@ -714,10 +720,15 @@ subroutine define_restart_vars(this, initialize_variables) long_name='ed cohort - plant height', units='m', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_height_co ) - call this%set_restart_var(vname='fates_laimemory', vtype=cohort_r8, & + call this%set_restart_var(vname='fates_leafmemory', vtype=cohort_r8, & long_name='ed cohort - target leaf biomass set from prev year', & units='kgC/indiv', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_laimemory_co ) + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_leafmemory_co ) + + call this%set_restart_var(vname='fates_fnrtmemory', vtype=cohort_r8, & + long_name='ed cohort - target fine-root biomass set from prev year', & + units='kgC/indiv', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fnrtmemory_co ) call this%set_restart_var(vname='fates_sapwmemory', vtype=cohort_r8, & long_name='ed cohort - target sapwood biomass set from prev year', & @@ -882,6 +893,18 @@ subroutine define_restart_vars(this, initialize_variables) long_name='ed cohort - plant phenology status', units='unitless', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_status_co ) + call this%set_restart_var(vname='fates_efleaf_coh', vtype=cohort_r8, & + long_name='ed cohort - leaf elongation factor', units='unitless', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_efleaf_co ) + + call this%set_restart_var(vname='fates_effnrt_coh', vtype=cohort_r8, & + long_name='ed cohort - fine-root elongation factor', units='unitless', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_effnrt_co ) + + call this%set_restart_var(vname='fates_efstem_coh', vtype=cohort_r8, & + long_name='ed cohort - stem elongation factor', units='unitless', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_efstem_co ) + call this%set_restart_var(vname='fates_isnew', vtype=cohort_int, & long_name='ed cohort - binary flag specifying if a plant has experienced a full day cycle', & units='0/1', flushval = flushone, & @@ -1154,10 +1177,39 @@ subroutine define_restart_vars(this, initialize_variables) ! site x time level vars ! - call this%set_restart_var(vname='fates_water_memory', vtype=cohort_r8, & + call this%set_restart_var(vname='fates_drought_dec_status', vtype=cohort_int, & + long_name='status flag for drought deciduous plants', units='unitless', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_dd_status_sift ) + + call this%set_restart_var(vname='fates_drought_leafondate', vtype=cohort_int, & + long_name='the day of year for drought based leaf-on', units='day of year', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_dleafondate_sift ) + + call this%set_restart_var(vname='fates_drought_leafoffdate', vtype=cohort_int, & + long_name='the day of year for drought based leaf-off', units='day of year', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_dleafoffdate_sift ) + + call this%set_restart_var(vname='fates_drought_ndaysleafon', vtype=cohort_int, & + long_name='number of days since leaf on (drought deciduous)', units='days', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_dndaysleafon_sift ) + + call this%set_restart_var(vname='fates_drought_ndaysleafoff', vtype=cohort_int, & + long_name='number of days since leaf off (drought deciduous)', units='days', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_dndaysleafoff_sift ) + + call this%set_restart_var(vname='fates_elong_factor', vtype=cohort_r8, & + long_name='leaf elongation factor (0 - completely abscissed; 1 - completely flushed)', units='unitless', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_elong_factor_sift ) + + call this%set_restart_var(vname='fates_liqvol_memory', vtype=cohort_r8, & long_name='last 10 days of volumetric soil water, by site x day-index', & units='m3/m3', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_watermem_siwm ) + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_liqvolmem_siwmft ) + + call this%set_restart_var(vname='fates_smp_memory', vtype=cohort_r8, & + long_name='last 10 days of soil matric potential, by site x day-index', & + units='mm', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_smpmem_siwmft ) call this%set_restart_var(vname='fates_vegtemp_memory', vtype=cohort_r8, & long_name='last 10 days of 24-hour vegetation temperature, by site x day-index', & @@ -1745,13 +1797,12 @@ subroutine set_restart_vectors(this,nc,nsites,sites) associate( rio_npatch_si => this%rvars(ir_npatch_si)%int1d, & rio_cd_status_si => this%rvars(ir_cd_status_si)%int1d, & - rio_dd_status_si => this%rvars(ir_dd_status_si)%int1d, & rio_nchill_days_si => this%rvars(ir_nchill_days_si)%int1d, & rio_ncold_days_si => this%rvars(ir_ncold_days_si)%int1d, & rio_leafondate_si => this%rvars(ir_leafondate_si)%int1d, & rio_leafoffdate_si => this%rvars(ir_leafoffdate_si)%int1d, & - rio_dleafondate_si => this%rvars(ir_dleafondate_si)%int1d, & - rio_dleafoffdate_si => this%rvars(ir_dleafoffdate_si)%int1d, & + rio_cndaysleafon_si => this%rvars(ir_cndaysleafon_si)%int1d, & + rio_cndaysleafoff_si => this%rvars(ir_cndaysleafoff_si)%int1d, & rio_acc_ni_si => this%rvars(ir_acc_ni_si)%r81d, & rio_gdd_si => this%rvars(ir_gdd_si)%r81d, & rio_snow_depth_si => this%rvars(ir_snow_depth_si)%r81d, & @@ -1769,7 +1820,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_coage_co => this%rvars(ir_coage_co)%r81d, & rio_g_sb_laweight_co => this%rvars(ir_g_sb_laweight_co)%r81d, & rio_height_co => this%rvars(ir_height_co)%r81d, & - rio_laimemory_co => this%rvars(ir_laimemory_co)%r81d, & + rio_leafmemory_co => this%rvars(ir_leafmemory_co)%r81d, & + rio_fnrtmemory_co => this%rvars(ir_fnrtmemory_co)%r81d, & rio_sapwmemory_co => this%rvars(ir_sapwmemory_co)%r81d, & rio_structmemory_co => this%rvars(ir_structmemory_co)%r81d, & rio_nplant_co => this%rvars(ir_nplant_co)%r81d, & @@ -1803,6 +1855,9 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_resp_tstep_co => this%rvars(ir_resp_tstep_co)%r81d, & rio_pft_co => this%rvars(ir_pft_co)%int1d, & rio_status_co => this%rvars(ir_status_co)%int1d, & + rio_efleaf_co => this%rvars(ir_efleaf_co)%r81d, & + rio_effnrt_co => this%rvars(ir_effnrt_co)%r81d, & + rio_efstem_co => this%rvars(ir_efstem_co)%r81d, & rio_isnew_co => this%rvars(ir_isnew_co)%int1d, & rio_gnd_alb_dif_pasb => this%rvars(ir_gnd_alb_dif_pasb)%r81d, & rio_gnd_alb_dir_pasb => this%rvars(ir_gnd_alb_dir_pasb)%r81d, & @@ -1813,7 +1868,14 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_agesinceanthrodist_pa => this%rvars(ir_agesinceanthrodist_pa)%r81d, & rio_nocomp_pft_label_pa => this%rvars(ir_nocomp_pft_label_pa)%int1d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & - rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d, & + rio_dd_status_sift => this%rvars(ir_dd_status_sift)%int1d, & + rio_dleafondate_sift => this%rvars(ir_dleafondate_sift)%int1d, & + rio_dleafoffdate_sift => this%rvars(ir_dleafoffdate_sift)%int1d, & + rio_dndaysleafon_sift => this%rvars(ir_dndaysleafon_sift)%int1d, & + rio_dndaysleafoff_sift => this%rvars(ir_dndaysleafoff_sift)%int1d, & + rio_elong_factor_sift => this%rvars(ir_elong_factor_sift)%r81d, & + rio_liqvolmem_siwmft => this%rvars(ir_liqvolmem_siwmft)%r81d, & + rio_smpmem_siwmft => this%rvars(ir_smpmem_siwmft)%r81d, & rio_vegtempmem_sitm => this%rvars(ir_vegtempmem_sitm)%r81d, & rio_recrate_sift => this%rvars(ir_recrate_sift)%r81d, & rio_use_this_pft_sift => this%rvars(ir_use_this_pft_sift)%int1d, & @@ -1997,7 +2059,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_dbh_co(io_idx_co) = ccohort%dbh rio_coage_co(io_idx_co) = ccohort%coage rio_height_co(io_idx_co) = ccohort%hite - rio_laimemory_co(io_idx_co) = ccohort%laimemory + rio_leafmemory_co(io_idx_co) = ccohort%leafmemory + rio_fnrtmemory_co(io_idx_co) = ccohort%fnrtmemory rio_sapwmemory_co(io_idx_co) = ccohort%sapwmemory rio_structmemory_co(io_idx_co) = ccohort%structmemory rio_g_sb_laweight_co(io_idx_co)= ccohort%g_sb_laweight @@ -2042,6 +2105,9 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_resp_tstep_co(io_idx_co) = ccohort%resp_tstep rio_pft_co(io_idx_co) = ccohort%pft rio_status_co(io_idx_co) = ccohort%status_coh + rio_efleaf_co(io_idx_co) = ccohort%efleaf_coh + rio_effnrt_co(io_idx_co) = ccohort%effnrt_coh + rio_efstem_co(io_idx_co) = ccohort%efstem_coh if ( ccohort%isnew ) then rio_isnew_co(io_idx_co) = new_cohort else @@ -2229,14 +2295,26 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_fmortcflux_usto_si(io_idx_si) = sites(s)%fmort_carbonflux_ustory rio_cd_status_si(io_idx_si) = sites(s)%cstatus - rio_dd_status_si(io_idx_si) = sites(s)%dstatus rio_nchill_days_si(io_idx_si) = sites(s)%nchilldays rio_ncold_days_si(io_idx_si) = sites(s)%ncolddays rio_leafondate_si(io_idx_si) = sites(s)%cleafondate rio_leafoffdate_si(io_idx_si) = sites(s)%cleafoffdate + rio_cndaysleafon_si(io_idx_si) = sites(s)%cndaysleafon + rio_cndaysleafoff_si(io_idx_si) = sites(s)%cndaysleafoff + + ! Drought-deciduous phenology are now PFT dependent + io_idx_si_pft = io_idx_co_1st + do i_pft = 1,numpft + rio_dd_status_sift(io_idx_si_pft) = sites(s)%dstatus(i_pft) + rio_dleafondate_sift(io_idx_si_pft) = sites(s)%dleafondate(i_pft) + rio_dleafoffdate_sift(io_idx_si_pft) = sites(s)%dleafoffdate(i_pft) + rio_dndaysleafon_sift(io_idx_si_pft) = sites(s)%dndaysleafon(i_pft) + rio_dndaysleafoff_sift(io_idx_si_pft) = sites(s)%dndaysleafoff(i_pft) + rio_elong_factor_sift(io_idx_si_pft) = sites(s)%elong_factor(i_pft) + + io_idx_si_pft = io_idx_si_pft + 1 + end do - rio_dleafondate_si(io_idx_si) = sites(s)%dleafondate - rio_dleafoffdate_si(io_idx_si) = sites(s)%dleafoffdate rio_acc_ni_si(io_idx_si) = sites(s)%acc_NI rio_gdd_si(io_idx_si) = sites(s)%grow_deg_days rio_snow_depth_si(io_idx_si) = sites(s)%snow_depth @@ -2248,8 +2326,11 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_npatch_si(io_idx_si) = patchespersite do i = 1,numWaterMem ! numWaterMem currently 10 - rio_watermem_siwm( io_idx_si_wmem ) = sites(s)%water_memory(i) - io_idx_si_wmem = io_idx_si_wmem + 1 + do i_pft=1,numpft + rio_liqvolmem_siwmft( io_idx_si_wmem ) = sites(s)%liqvol_memory(i,i_pft) + rio_smpmem_siwmft( io_idx_si_wmem ) = sites(s)%smp_memory(i,i_pft) + io_idx_si_wmem = io_idx_si_wmem + 1 + end do end do do i = 1, num_vegtemp_mem @@ -2576,13 +2657,12 @@ subroutine get_restart_vectors(this, nc, nsites, sites) associate( rio_npatch_si => this%rvars(ir_npatch_si)%int1d, & rio_cd_status_si => this%rvars(ir_cd_status_si)%int1d, & - rio_dd_status_si => this%rvars(ir_dd_status_si)%int1d, & rio_nchill_days_si => this%rvars(ir_nchill_days_si)%int1d, & rio_ncold_days_si => this%rvars(ir_ncold_days_si)%int1d, & rio_leafondate_si => this%rvars(ir_leafondate_si)%int1d, & rio_leafoffdate_si => this%rvars(ir_leafoffdate_si)%int1d, & - rio_dleafondate_si => this%rvars(ir_dleafondate_si)%int1d, & - rio_dleafoffdate_si => this%rvars(ir_dleafoffdate_si)%int1d, & + rio_cndaysleafon_si => this%rvars(ir_cndaysleafon_si)%int1d, & + rio_cndaysleafoff_si => this%rvars(ir_cndaysleafoff_si)%int1d, & rio_acc_ni_si => this%rvars(ir_acc_ni_si)%r81d, & rio_gdd_si => this%rvars(ir_gdd_si)%r81d, & rio_snow_depth_si => this%rvars(ir_snow_depth_si)%r81d, & @@ -2600,7 +2680,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_coage_co => this%rvars(ir_coage_co)%r81d, & rio_g_sb_laweight_co => this%rvars(ir_g_sb_laweight_co)%r81d, & rio_height_co => this%rvars(ir_height_co)%r81d, & - rio_laimemory_co => this%rvars(ir_laimemory_co)%r81d, & + rio_leafmemory_co => this%rvars(ir_leafmemory_co)%r81d, & + rio_fnrtmemory_co => this%rvars(ir_fnrtmemory_co)%r81d, & rio_sapwmemory_co => this%rvars(ir_sapwmemory_co)%r81d, & rio_structmemory_co => this%rvars(ir_structmemory_co)%r81d, & rio_nplant_co => this%rvars(ir_nplant_co)%r81d, & @@ -2634,6 +2715,9 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_resp_tstep_co => this%rvars(ir_resp_tstep_co)%r81d, & rio_pft_co => this%rvars(ir_pft_co)%int1d, & rio_status_co => this%rvars(ir_status_co)%int1d, & + rio_efleaf_co => this%rvars(ir_efleaf_co)%r81d, & + rio_effnrt_co => this%rvars(ir_effnrt_co)%r81d, & + rio_efstem_co => this%rvars(ir_efstem_co)%r81d, & rio_isnew_co => this%rvars(ir_isnew_co)%int1d, & rio_gnd_alb_dif_pasb => this%rvars(ir_gnd_alb_dif_pasb)%r81d, & rio_gnd_alb_dir_pasb => this%rvars(ir_gnd_alb_dir_pasb)%r81d, & @@ -2644,7 +2728,14 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_agesinceanthrodist_pa => this%rvars(ir_agesinceanthrodist_pa)%r81d, & rio_nocomp_pft_label_pa => this%rvars(ir_nocomp_pft_label_pa)%int1d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & - rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d, & + rio_dd_status_sift => this%rvars(ir_dd_status_sift)%int1d, & + rio_dleafondate_sift => this%rvars(ir_dleafondate_sift)%int1d, & + rio_dleafoffdate_sift => this%rvars(ir_dleafoffdate_sift)%int1d, & + rio_dndaysleafon_sift => this%rvars(ir_dndaysleafon_sift)%int1d, & + rio_dndaysleafoff_sift => this%rvars(ir_dndaysleafoff_sift)%int1d, & + rio_elong_factor_sift => this%rvars(ir_elong_factor_sift)%r81d, & + rio_liqvolmem_siwmft => this%rvars(ir_liqvolmem_siwmft)%r81d, & + rio_smpmem_siwmft => this%rvars(ir_smpmem_siwmft)%r81d, & rio_vegtempmem_sitm => this%rvars(ir_vegtempmem_sitm)%r81d, & rio_recrate_sift => this%rvars(ir_recrate_sift)%r81d, & rio_use_this_pft_sift => this%rvars(ir_use_this_pft_sift)%int1d, & @@ -2802,7 +2893,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%coage = rio_coage_co(io_idx_co) ccohort%g_sb_laweight= rio_g_sb_laweight_co(io_idx_co) ccohort%hite = rio_height_co(io_idx_co) - ccohort%laimemory = rio_laimemory_co(io_idx_co) + ccohort%leafmemory = rio_leafmemory_co(io_idx_co) + ccohort%fnrtmemory = rio_fnrtmemory_co(io_idx_co) ccohort%sapwmemory = rio_sapwmemory_co(io_idx_co) ccohort%structmemory= rio_structmemory_co(io_idx_co) ccohort%n = rio_nplant_co(io_idx_co) @@ -2843,6 +2935,9 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%resp_tstep = rio_resp_tstep_co(io_idx_co) ccohort%pft = rio_pft_co(io_idx_co) ccohort%status_coh = rio_status_co(io_idx_co) + ccohort%efleaf_coh = rio_efleaf_co(io_idx_co) + ccohort%effnrt_coh = rio_effnrt_co(io_idx_co) + ccohort%efstem_coh = rio_efstem_co(io_idx_co) ccohort%isnew = ( rio_isnew_co(io_idx_co) .eq. new_cohort ) call UpdateCohortBioPhysRates(ccohort) @@ -3017,8 +3112,11 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end if do i = 1,numWaterMem - sites(s)%water_memory(i) = rio_watermem_siwm( io_idx_si_wmem ) - io_idx_si_wmem = io_idx_si_wmem + 1 + do i_pft=1,numpft + sites(s)%liqvol_memory(i,i_pft) = rio_liqvolmem_siwmft( io_idx_si_wmem ) + sites(s)%smp_memory(i,i_pft) = rio_smpmem_siwmft( io_idx_si_wmem ) + io_idx_si_wmem = io_idx_si_wmem + 1 + end do end do do i = 1, num_vegtemp_mem @@ -3090,13 +3188,25 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! Site level phenology status flags sites(s)%cstatus = rio_cd_status_si(io_idx_si) - sites(s)%dstatus = rio_dd_status_si(io_idx_si) sites(s)%nchilldays = rio_nchill_days_si(io_idx_si) sites(s)%ncolddays = rio_ncold_days_si(io_idx_si) sites(s)%cleafondate = rio_leafondate_si(io_idx_si) sites(s)%cleafoffdate = rio_leafoffdate_si(io_idx_si) - sites(s)%dleafondate = rio_dleafondate_si(io_idx_si) - sites(s)%dleafoffdate = rio_dleafoffdate_si(io_idx_si) + sites(s)%cndaysleafon = rio_cndaysleafon_si(io_idx_si) + sites(s)%cndaysleafoff = rio_cndaysleafoff_si(io_idx_si) + + ! Fill drought-deciduous variables, which are now PFT variables. + io_idx_si_pft = io_idx_co_1st + do i_pft = 1,numpft + sites(s)%dstatus(i_pft) = rio_dd_status_sift(io_idx_si_pft) + sites(s)%dleafondate(i_pft) = rio_dleafondate_sift(io_idx_si_pft) + sites(s)%dleafoffdate(i_pft) = rio_dleafoffdate_sift(io_idx_si_pft) + sites(s)%dndaysleafon(i_pft) = rio_dndaysleafon_sift(io_idx_si_pft) + sites(s)%dndaysleafoff(i_pft) = rio_dndaysleafoff_sift(io_idx_si_pft) + sites(s)%elong_factor(i_pft) = rio_elong_factor_sift(io_idx_si_pft) + io_idx_si_pft = io_idx_si_pft + 1 + end do + sites(s)%acc_NI = rio_acc_ni_si(io_idx_si) sites(s)%grow_deg_days = rio_gdd_si(io_idx_si) sites(s)%snow_depth = rio_snow_depth_si(io_idx_si) diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 7a914fb573..5ce1481217 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -403,6 +403,9 @@ variables: double fates_phen_evergreen(fates_pft) ; fates_phen_evergreen:units = "logical flag" ; fates_phen_evergreen:long_name = "Binary flag for evergreen leaf habit" ; + double fates_phen_fnrt_drop_fraction(fates_pft) ; + fates_phen_fnrt_drop_fraction:units = "fraction" ; + fates_phen_fnrt_drop_fraction:long_name = "fraction of fine roots to drop during drought/cold" ; double fates_phen_season_decid(fates_pft) ; fates_phen_season_decid:units = "logical flag" ; fates_phen_season_decid:long_name = "Binary flag for seasonal-deciduous leaf habit" ; @@ -719,12 +722,18 @@ variables: double fates_phen_doff_time ; fates_phen_doff_time:units = "days" ; fates_phen_doff_time:long_name = "day threshold compared against days since leaves became off-allometry" ; + double fates_phen_drought_model ; + fates_phen_drought_model:units = "none" ; + fates_phen_drought_model:long_name = "which method to use for drought phenology: 0 - FATES default; 1 - Semi-deciduous (ED2-like)" ; double fates_phen_drought_threshold ; - fates_phen_drought_threshold:units = "m3/m3" ; - fates_phen_drought_threshold:long_name = "liquid volume in soil layer, threashold for drought phenology" ; + fates_phen_drought_threshold:units = "m3/m3 or mm" ; + fates_phen_drought_threshold:long_name = "threshold for drought phenology (or lower threshold when fates_phen_drought_model = 2); the quantity depends on the sign: if positive, the threshold is volumetric soil moisture (m3/m3). If negative, the threshold is soil matric potentical (mm)." ; double fates_phen_mindayson ; fates_phen_mindayson:units = "days" ; fates_phen_mindayson:long_name = "day threshold compared against days since leaves became on-allometry" ; + double fates_phen_moist_threshold ; + fates_phen_moist_threshold:units = "m3/m3 or mm" ; + fates_phen_moist_threshold:long_name = "upper threshold for drought phenology (only for fates_phen_drought_model=2); the quantity depends on the sign: if positive, the threshold is volumetric soil moisture (m3/m3). If negative, the threshold is soil matric potentical (mm)." ; double fates_phen_ncolddayslim ; fates_phen_ncolddayslim:units = "days" ; fates_phen_ncolddayslim:long_name = "day threshold exceedance for temperature leaf-drop" ; @@ -1170,6 +1179,8 @@ data: fates_phen_evergreen = 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0 ; + fates_phen_fnrt_drop_fraction = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; + fates_phen_season_decid = 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0 ; fates_phen_stem_drop_fraction = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; @@ -1463,10 +1474,14 @@ data: fates_phen_doff_time = 100 ; - fates_phen_drought_threshold = 0.15 ; + fates_phen_drought_model = 0 ; + + fates_phen_drought_threshold = -203943.2 ; fates_phen_mindayson = 90 ; + fates_phen_moist_threshold = -122365.9 ; + fates_phen_ncolddayslim = 5 ; fates_photo_temp_acclim_timescale = 30 ; diff --git a/parteh/PRTAllometricCNPMod.F90 b/parteh/PRTAllometricCNPMod.F90 index 5617d71e5d..b566a6fb9f 100644 --- a/parteh/PRTAllometricCNPMod.F90 +++ b/parteh/PRTAllometricCNPMod.F90 @@ -55,7 +55,7 @@ module PRTAllometricCNPMod use FatesConstantsMod , only : fates_unset_int use FatesConstantsMod , only : sec_per_day use PRTParametersMod , only : prt_params - use EDTypesMod , only : leaves_on,leaves_off + use EDTypesMod , only : leaves_on,leaves_off,leaves_pshed implicit none private @@ -143,16 +143,17 @@ module PRTAllometricCNPMod ! Input only Boundary Indices (These are public) ! ------------------------------------------------------------------------------------- - integer, public, parameter :: acnp_bc_in_id_pft = 1 ! Index for the PFT input BC - integer, public, parameter :: acnp_bc_in_id_ctrim = 2 ! Index for the canopy trim function - integer, public, parameter :: acnp_bc_in_id_lstat = 3 ! phenology status logical - integer, public, parameter :: acnp_bc_in_id_netdc = 4 ! Index for the net daily C input BC - integer, public, parameter :: acnp_bc_in_id_netdnh4 = 5 ! Index for the net daily NH4 input BC - integer, public, parameter :: acnp_bc_in_id_netdno3 = 6 ! Index for the net daily NO3 input BC - integer, public, parameter :: acnp_bc_in_id_netdp = 7 ! Index for the net daily P input BC - - ! 0=leaf off, 1=leaf on - integer, parameter :: num_bc_in = 7 + integer, public, parameter :: acnp_bc_in_id_pft = 1 ! Index for the PFT input BC + integer, public, parameter :: acnp_bc_in_id_ctrim = 2 ! Index for the canopy trim function + integer, public, parameter :: acnp_bc_in_id_lstat = 3 ! phenology status logical + integer, public, parameter :: acnp_bc_in_id_netdc = 4 ! Index for the net daily C input BC + integer, public, parameter :: acnp_bc_in_id_netdnh4 = 5 ! Index for the net daily NH4 input BC + integer, public, parameter :: acnp_bc_in_id_netdno3 = 6 ! Index for the net daily NO3 input BC + integer, public, parameter :: acnp_bc_in_id_netdp = 7 ! Index for the net daily P input BC + integer, public, parameter :: acnp_bc_in_id_efleaf = 8 ! Leaf elongation factor + integer, public, parameter :: acnp_bc_in_id_effnrt = 9 ! Fine-root "elongation factor" + integer, public, parameter :: acnp_bc_in_id_efstem = 10 ! Stem "elongation factor" + integer, parameter :: num_bc_in = 10 ! ------------------------------------------------------------------------------------- ! Output Boundary Indices (These are public) @@ -334,6 +335,9 @@ subroutine DailyPRTAllometricCNP(this) real(r8) :: n_gain ! Daily nitrogen uptake through fine-roots [kgN] real(r8) :: p_gain ! Daily phosphorus uptake through fine-roots [kgN] real(r8) :: canopy_trim ! The canopy trimming function [0-1] + real(r8) :: elongf_leaf ! Leaf elongation factor [0-1] + real(r8) :: elongf_fnrt ! Fine-root "elongation factor" [0-1] + real(r8) :: elongf_stem ! Stem "elongation factor" [0-1] ! Pointers to output bcs real(r8),pointer :: c_efflux ! Total plant efflux of carbon (kgC) @@ -396,6 +400,9 @@ subroutine DailyPRTAllometricCNP(this) p_gain = this%bc_in(acnp_bc_in_id_netdp)%rval; p_gain0 = p_gain canopy_trim = this%bc_in(acnp_bc_in_id_ctrim)%rval ipft = this%bc_in(acnp_bc_in_id_pft)%ival + elongf_leaf = this%bc_in(acnp_bc_in_id_efleaf)%rval + elongf_fnrt = this%bc_in(acnp_bc_in_id_effnrt)%rval + elongf_stem = this%bc_in(acnp_bc_in_id_efstem)%rval ! Output only boundary conditions c_efflux => this%bc_out(acnp_bc_out_id_cefflux)%rval; c_efflux = 0._r8 @@ -441,6 +448,17 @@ subroutine DailyPRTAllometricCNP(this) target_c(repro_id) = 0._r8 target_dcdd(repro_id) = 0._r8 + ! Correct tissue targets based on the elongation factor + target_c(leaf_id) = elongf_leaf * target_c(leaf_id) + target_c(fnrt_id) = elongf_fnrt * target_c(fnrt_id) + target_c(sapw_id) = elongf_stem * target_c(sapw_id) + target_c(struct_id) = elongf_stem * target_c(struct_id) + ! MLO - Need to check whether or not the multiplication for the growth is correct or not. + target_dcdd(leaf_id) = elongf_leaf * target_dcdd(leaf_id) + target_dcdd(fnrt_id) = elongf_fnrt * target_dcdd(fnrt_id) + target_dcdd(sapw_id) = elongf_stem * target_dcdd(sapw_id) + target_dcdd(struct_id) = elongf_stem * target_dcdd(struct_id) + ! Initialize the the state, and keep a record of this state ! as we may actuall run the allocation process twice, and ! will need this state to both reset, and measure total @@ -647,6 +665,9 @@ subroutine CNPPrioritizedReplacement(this, & real(r8) :: target_n ! Target mass of N for a given organ [kg] real(r8) :: target_p ! Target mass of P for a given organ [kg] real(r8) :: c_gain0 + real(r8) :: elongf_leaf ! Leaf elongation factor + real(r8) :: elongf_fnrt ! Fine-root "elongation factor" + real(r8) :: elongf_stem ! Stem "elongation factor" integer :: priority_code ! Index for priority level of each organ real(r8) :: sum_c_demand ! Carbon demanded to bring tissues up to allometry (kg) real(r8) :: sum_n_deficit ! The nitrogen deficit of all pools for given priority level (kg) @@ -674,6 +695,9 @@ subroutine CNPPrioritizedReplacement(this, & c_gain0 = c_gain leaf_status = this%bc_in(acnp_bc_in_id_lstat)%ival + elongf_leaf = this%bc_in(acnp_bc_in_id_efleaf)%rval + elongf_fnrt = this%bc_in(acnp_bc_in_id_effnrt)%rval + elongf_stem = this%bc_in(acnp_bc_in_id_efstem)%rval ipft = this%bc_in(acnp_bc_in_id_pft)%ival canopy_trim = this%bc_in(acnp_bc_in_id_ctrim)%rval @@ -700,7 +724,8 @@ subroutine CNPPrioritizedReplacement(this, & ! Don't allow allocation to leaves if they are in an "off" status. ! Also, dont allocate to replace turnover if this is not evergreen ! (this prevents accidental re-flushing on the day they drop) - if( ((leaf_status.eq.leaves_off) .or. (prt_params%evergreen(ipft) .ne. itrue)) & + if( ( any(leaf_status == [leaves_off,leaves_pshed]) .or. & + (prt_params%evergreen(ipft) .ne. itrue) ) & .and. (organ_list(ii).eq.leaf_organ)) cycle ! 1 is the highest priority code possible @@ -865,7 +890,7 @@ subroutine CNPPrioritizedReplacement(this, & ! Don't allow allocation to leaves if they are in an "off" status. ! (this prevents accidental re-flushing on the day they drop) - if((leaf_status.eq.leaves_off) .and. (organ_list(ii).eq.leaf_organ)) cycle + if(any(leaf_status == [leaves_off,leaves_pshed]) .and. (organ_list(ii).eq.leaf_organ)) cycle ! 1 is the highest priority code possible if( priority_code == i_pri ) then @@ -1009,7 +1034,10 @@ subroutine CNPStatureGrowth(this,c_gain, n_gain, p_gain, & real(r8), pointer :: dbh integer :: ipft real(r8) :: canopy_trim - real(r8) :: leaf_status + integer :: leaf_status + real(r8) :: elongf_leaf + real(r8) :: elongf_fnrt + real(r8) :: elongf_stem integer :: i, ii ! organ index loops (masked and unmasked) integer :: istep ! outer step iteration loop @@ -1088,6 +1116,9 @@ subroutine CNPStatureGrowth(this,c_gain, n_gain, p_gain, & integer, parameter :: p_limited = 3 leaf_status = this%bc_in(acnp_bc_in_id_lstat)%ival + elongf_leaf = this%bc_in(acnp_bc_in_id_efleaf)%rval + elongf_fnrt = this%bc_in(acnp_bc_in_id_effnrt)%rval + elongf_stem = this%bc_in(acnp_bc_in_id_efstem)%rval dbh => this%bc_inout(acnp_bc_inout_id_dbh)%rval ipft = this%bc_in(acnp_bc_in_id_pft)%ival canopy_trim = this%bc_in(acnp_bc_in_id_ctrim)%rval @@ -1105,16 +1136,21 @@ subroutine CNPStatureGrowth(this,c_gain, n_gain, p_gain, & if( c_gain <= calloc_abs_error .or. & n_gain <= 0.1_r8*calloc_abs_error .or. & p_gain <= 0.02_r8*calloc_abs_error .or. & - leaf_status.eq.leaves_off ) then + any(leaf_status == [leaves_off,leaves_pshed]) ) then return end if - intgr_params(:) = fates_unset_r8 - intgr_params(acnp_bc_in_id_ctrim) = this%bc_in(acnp_bc_in_id_ctrim)%rval - intgr_params(acnp_bc_in_id_pft) = real(this%bc_in(acnp_bc_in_id_pft)%ival) - - + intgr_params(:) = fates_unset_r8 + intgr_params(acnp_bc_in_id_ctrim) = this%bc_in(acnp_bc_in_id_ctrim)%rval + intgr_params(acnp_bc_in_id_pft) = real(this%bc_in(acnp_bc_in_id_pft)%ival,r8) + intgr_params(acnp_bc_in_id_lstat) = real(this%bc_in(acnp_bc_in_id_lstat)%ival,r8) + intgr_params(acnp_bc_in_id_efleaf) = this%bc_in(acnp_bc_in_id_efleaf)%rval + intgr_params(acnp_bc_in_id_effnrt) = this%bc_in(acnp_bc_in_id_effnrt)%rval + intgr_params(acnp_bc_in_id_efstem) = this%bc_in(acnp_bc_in_id_efstem)%rval + + + state_mask(:) = .false. mask_organs(:) = fates_unset_int @@ -1358,6 +1394,7 @@ subroutine CNPStatureGrowth(this,c_gain, n_gain, p_gain, & end do call CheckIntegratedAllometries(state_array_out(dbh_id),ipft,canopy_trim, & + elongf_leaf, elongf_fnrt, elongf_stem, & leafc_tp1, state_array_out(fnrt_id), state_array_out(sapw_id), & state_array_out(store_id), state_array_out(struct_id), & state_mask(leaf_id), state_mask(fnrt_id), state_mask(sapw_id), & @@ -1434,6 +1471,9 @@ subroutine CNPStatureGrowth(this,c_gain, n_gain, p_gain, & write(fates_log(),*) 'totalC',totalC write(fates_log(),*) 'pft: ',ipft write(fates_log(),*) 'dbh: ',dbh + write(fates_log(),*) 'elongf_leaf: ',elongf_leaf + write(fates_log(),*) 'elongf_fnrt: ',elongf_fnrt + write(fates_log(),*) 'elongf_stem: ',elongf_stem write(fates_log(),*) 'dCleaf_dd: ',target_dcdd(leaf_id) write(fates_log(),*) 'dCfnrt_dd: ',target_dcdd(fnrt_id) write(fates_log(),*) 'dCstore_dd: ',target_dcdd(store_id) @@ -1454,7 +1494,13 @@ subroutine CNPStatureGrowth(this,c_gain, n_gain, p_gain, & call bbgw_allom(dbh_tp1,ipft,bgw_c_target_tp1) call bdead_allom(agw_c_target_tp1,bgw_c_target_tp1, sapw_c_target_tp1, ipft, struct_c_target_tp1) call bstore_allom(dbh_tp1,ipft,canopy_trim,store_c_target_tp1) - + + ! Correct the targets based on the elongation factors + leaf_c_target_tp1 = elongf_leaf * leaf_c_target_tp1 + fnrt_c_target_tp1 = elongf_fnrt * fnrt_c_target_tp1 + sapw_c_target_tp1 = elongf_stem * sapw_c_target_tp1 + struct_c_target_tp1 = elongf_stem * struct_c_target_tp1 + write(fates_log(),*) 'leaf_c: ',leafc_tp1, leaf_c_target_tp1,leafc_tp1-leaf_c_target_tp1 write(fates_log(),*) 'fnrt_c: ',fnrtc_tp1, fnrt_c_target_tp1,fnrtc_tp1- fnrt_c_target_tp1 write(fates_log(),*) 'sapw_c: ',sapwc_tp1, sapw_c_target_tp1 ,sapwc_tp1- sapw_c_target_tp1 @@ -1710,15 +1756,19 @@ function GetNutrientTargetCNP(this,element_id,organ_id,stoich_mode) result(targe real(r8) :: leaf_c_target,fnrt_c_target real(r8) :: sapw_c_target,agw_c_target real(r8) :: bgw_c_target,struct_c_target + real(r8) :: elongf_leaf + real(r8) :: elongf_fnrt + real(r8) :: elongf_stem + - - - dbh => this%bc_inout(acnp_bc_inout_id_dbh)%rval canopy_trim = this%bc_in(acnp_bc_in_id_ctrim)%rval ipft = this%bc_in(acnp_bc_in_id_pft)%ival + elongf_leaf = this%bc_in(acnp_bc_in_id_efleaf)%rval + elongf_fnrt = this%bc_in(acnp_bc_in_id_effnrt)%rval + elongf_stem = this%bc_in(acnp_bc_in_id_efstem)%rval i_cvar = prt_global%sp_organ_map(organ_id,carbon12_element) - + ! Storage of nutrients are assumed to have different compartments than ! for carbon, and thus their targets are not associated with a tissue ! but is more represented as a fraction of the maximum amount of nutrient @@ -1733,6 +1783,12 @@ function GetNutrientTargetCNP(this,element_id,organ_id,stoich_mode) result(targe call bbgw_allom(dbh,ipft,bgw_c_target) call bdead_allom(agw_c_target,bgw_c_target, sapw_c_target, ipft, struct_c_target) + ! Correct the targets based on the elongation factors + leaf_c_target = elongf_leaf * leaf_c_target + fnrt_c_target = elongf_fnrt * fnrt_c_target + sapw_c_target = elongf_stem * sapw_c_target + struct_c_target = elongf_stem * struct_c_target + ! Target for storage is a fraction of the sum target of all ! non-reproductive organs @@ -2102,7 +2158,10 @@ function AllomCNPGrowthDeriv(l_state_array,l_state_mask,cbalance,intgr_params) r real(r8) :: total_dcdd_target ! target total (not reproductive) biomass derivative wrt d, (kgC/cm) real(r8) :: repro_fraction ! fraction of carbon balance directed towards reproduction (kgC/kgC) real(r8) :: total_dcostdd ! carbon cost for non-reproductive pools per unit increment of dbh - + + real(r8) :: elongf_leaf ! Leaf elongation factor (0-1) + real(r8) :: elongf_fnrt ! Fine-root "elongation factor" (0-1) + real(r8) :: elongf_stem ! Stem "elongation factor" (0-1) associate( dbh => l_state_array(dbh_id), & leaf_c => l_state_array(leaf_id), & @@ -2122,6 +2181,9 @@ function AllomCNPGrowthDeriv(l_state_array,l_state_mask,cbalance,intgr_params) r canopy_trim = intgr_params(acnp_bc_in_id_ctrim) ipft = int(intgr_params(acnp_bc_in_id_pft)) + elongf_leaf = intgr_params(acnp_bc_in_id_efleaf) + elongf_fnrt = intgr_params(acnp_bc_in_id_effnrt) + elongf_stem = intgr_params(acnp_bc_in_id_efstem) call bleaf(dbh,ipft,canopy_trim,leaf_c_target,leaf_dcdd_target) call bfineroot(dbh,ipft,canopy_trim,fnrt_c_target,fnrt_dcdd_target) @@ -2132,6 +2194,17 @@ function AllomCNPGrowthDeriv(l_state_array,l_state_mask,cbalance,intgr_params) r agw_dcdd_target, bgw_dcdd_target, sapw_dcdd_target, struct_dcdd_target) call bstore_allom(dbh,ipft,canopy_trim,store_c_target,store_dcdd_target) + ! Apply correction for partially deciduous plants. + leaf_c_target = elongf_leaf * leaf_c_target + fnrt_c_target = elongf_fnrt * fnrt_c_target + sapw_c_target = elongf_stem * sapw_c_target + struct_c_target = elongf_stem * struct_c_target + ! MLO - Need to double check that it is correct to apply factor to both stocks and growth + leaf_dcdd_target = elongf_leaf * leaf_dcdd_target + fnrt_dcdd_target = elongf_fnrt * fnrt_dcdd_target + sapw_dcdd_target = elongf_stem * sapw_dcdd_target + struct_dcdd_target = elongf_stem * struct_dcdd_target + if (mask_repro) then ! fraction of carbon going towards reproduction if (dbh <= prt_params%dbh_repro_threshold(ipft)) then @@ -2222,72 +2295,101 @@ end function AllomCNPGrowthDeriv ! ==================================================================================== - subroutine TargetAllometryCheck(bleaf,bfroot,bsap,bstore,bdead, & - bt_leaf,bt_froot,bt_sap,bt_store,bt_dead, & - grow_leaf,grow_froot,grow_sapw,grow_store) + subroutine TargetAllometryCheck(b0_leaf,b0_fnrt,b0_sapw,b0_store,b0_struct, & + bleaf,bfnrt,bsapw,bstore,bstruct, & + bt_leaf,bt_fnrt,bt_sapw,bt_store,bt_struct, & + carbon_balance, & + elongf_leaf,elongf_fnrt,elongf_stem,ipft,leaf_status, & + grow_leaf,grow_fnrt,grow_sapw,grow_store,grow_struct) ! Arguments - real(r8),intent(in) :: bleaf !actual - real(r8),intent(in) :: bfroot - real(r8),intent(in) :: bsap + real(r8),intent(in) :: b0_leaf !initial + real(r8),intent(in) :: b0_fnrt + real(r8),intent(in) :: b0_sapw + real(r8),intent(in) :: b0_store + real(r8),intent(in) :: b0_struct + real(r8),intent(in) :: bleaf !actual + real(r8),intent(in) :: bfnrt + real(r8),intent(in) :: bsapw real(r8),intent(in) :: bstore - real(r8),intent(in) :: bdead - real(r8),intent(in) :: bt_leaf !target - real(r8),intent(in) :: bt_froot - real(r8),intent(in) :: bt_sap + real(r8),intent(in) :: bstruct + real(r8),intent(in) :: bt_leaf !target + real(r8),intent(in) :: bt_fnrt + real(r8),intent(in) :: bt_sapw real(r8),intent(in) :: bt_store - real(r8),intent(in) :: bt_dead - logical,intent(out) :: grow_leaf !growth flag - logical,intent(out) :: grow_froot + real(r8),intent(in) :: bt_struct + real(r8),intent(in) :: carbon_balance !remaining carbon balance + real(r8),intent(in) :: elongf_leaf !elongation factors + real(r8),intent(in) :: elongf_fnrt + real(r8),intent(in) :: elongf_stem + integer,intent(in) :: ipft !Plant functional type + integer,intent(in) :: leaf_status !Phenology status + logical,intent(out) :: grow_leaf !growth flag + logical,intent(out) :: grow_fnrt logical,intent(out) :: grow_sapw logical,intent(out) :: grow_store - - if( (bt_leaf - bleaf)>calloc_abs_error) then - write(fates_log(),*) 'leaves are not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bleaf,bt_leaf - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( (bleaf - bt_leaf)>calloc_abs_error) then - ! leaf is above allometry, ignore - grow_leaf = .false. - else - grow_leaf = .true. - end if - - if( (bt_froot - bfroot)>calloc_abs_error) then - write(fates_log(),*) 'fineroots are not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bfroot, bt_froot - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bfroot-bt_froot)>calloc_abs_error ) then - grow_froot = .false. - else - grow_froot = .true. - end if - - if( (bt_sap - bsap)>calloc_abs_error) then - write(fates_log(),*) 'sapwood is not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bsap, bt_sap - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bsap-bt_sap)>calloc_abs_error ) then - grow_sapw = .false. + logical,intent(out) :: grow_struct + ! Local variables + logical :: fine_leaf + logical :: fine_fnrt + logical :: fine_sapw + logical :: fine_store + logical :: fine_struct + logical :: all_fine + ! Local constants + character(len= 3), parameter :: fmth = '(a)' + character(len=27), parameter :: fmtb = '(a,3(1x,es12.5,1x,a),1x,l1)' + character(len=13), parameter :: fmte = '(a,1x,es12.5)' + character(len=10), parameter :: fmti = '(a,1x,i12)' + + + ! First test whether or not each pool looks reasonable. + fine_leaf = (bt_leaf - bleaf ) <= calloc_abs_error + fine_fnrt = (bt_fnrt - bfnrt ) <= calloc_abs_error + fine_sapw = (bt_sapw - bsapw ) <= calloc_abs_error + fine_store = (bt_store - bstore ) <= calloc_abs_error + fine_struct = (bt_struct - bstruct) <= calloc_abs_error + all_fine = fine_leaf .and. fine_fnrt .and. fine_sapw .and. & + fine_store .and. fine_struct + + ! Decide whether or not to grow tissues (but only if all tissues look fine). + ! We grow only when biomass is less than target biomass (with tolerance). + if (all_fine) then + grow_leaf = ( bleaf - bt_leaf ) <= calloc_abs_error + grow_fnrt = ( bfnrt - bt_fnrt ) <= calloc_abs_error + grow_sapw = ( bsapw - bt_sapw ) <= calloc_abs_error + grow_store = ( bstore - bt_store ) <= calloc_abs_error + grow_struct = ( bstruct - bt_struct ) <= calloc_abs_error else - grow_sapw = .true. - end if - - if( (bt_store - bstore)>calloc_abs_error) then - write(fates_log(),*) 'storage is not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bstore,bt_store + ! If anything looks not fine, write a detailed report + write(fates_log(),fmt=fmth) '======' + write(fates_log(),fmt=fmth) ' At least one tissue is not on-allometry at the growth step' + write(fates_log(),fmt=fmth) '======' + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) ' Biomass and on-allometry test (''F'' means problem)' + write(fates_log(),fmt=fmth) '------' + write(fates_log(),fmt=fmth) ' Tissue | Initial | Current | Target | On-allometry' + write(fates_log(),fmt=fmtb) ' Leaf |',b0_leaf ,'|',bleaf ,'|',bt_leaf ,'|',fine_leaf + write(fates_log(),fmt=fmtb) ' Fine root |',b0_fnrt ,'|',bfnrt ,'|',bt_fnrt ,'|',fine_fnrt + write(fates_log(),fmt=fmtb) ' Sap wood |',b0_sapw ,'|',bsapw ,'|',bt_sapw ,'|',fine_sapw + write(fates_log(),fmt=fmtb) ' Storage |',b0_store ,'|',bstore ,'|',bt_store ,'|',fine_store + write(fates_log(),fmt=fmtb) ' Structural |',b0_struct ,'|',bstruct ,'|',bt_struct ,'|',fine_struct + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) ' Ancillary information' + write(fates_log(),fmt=fmth) '------' + write(fates_log(),fmt=fmti) ' PFT = ',ipft + write(fates_log(),fmt=fmti) ' leaf_status = ',leaf_status + write(fates_log(),fmt=fmte) ' elongf_leaf = ',elongf_leaf + write(fates_log(),fmt=fmte) ' elongf_fnrt = ',elongf_fnrt + write(fates_log(),fmt=fmte) ' elongf_stem = ',elongf_stem + write(fates_log(),fmt=fmte) ' carbon_balance = ',carbon_balance + write(fates_log(),fmt=fmte) ' calloc_abs_error = ',calloc_abs_error + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) '======' call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bstore-bt_store)>calloc_abs_error ) then - grow_store = .false. - else - grow_store = .true. end if - if( (bt_dead - bdead)>calloc_abs_error) then - write(fates_log(),*) 'structure not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bdead,bt_dead - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + return end subroutine TargetAllometryCheck diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index 5bdf624502..90375e8fb4 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -51,6 +51,10 @@ module PRTAllometricCarbonMod use PRTParametersMod , only : prt_params + use EDTypesMod , only : leaves_on + use EDTypesMod , only : leaves_off + use EDTypesMod , only : leaves_pshed + implicit none private @@ -90,10 +94,13 @@ module PRTAllometricCarbonMod integer, parameter :: num_bc_inout = 2 ! Number of in & output boundary conditions - integer, public, parameter :: ac_bc_in_id_pft = 1 ! Index for the PFT input BC - integer, public, parameter :: ac_bc_in_id_ctrim = 2 ! Index for the canopy trim function - integer, public, parameter :: ac_bc_in_id_lstat = 3 ! Leaf status (on or off) - integer, parameter :: num_bc_in = 3 ! Number of input boundary conditions + integer, public, parameter :: ac_bc_in_id_pft = 1 ! Index for the PFT input BC + integer, public, parameter :: ac_bc_in_id_ctrim = 2 ! Index for the canopy trim function + integer, public, parameter :: ac_bc_in_id_lstat = 3 ! Leaf status (on or off) + integer, public, parameter :: ac_bc_in_id_efleaf = 4 ! Elongation factor (leaves) + integer, public, parameter :: ac_bc_in_id_effnrt = 5 ! "Elongation factor" (fine roots) + integer, public, parameter :: ac_bc_in_id_efstem = 6 ! "Elongation factor" (stem) + integer, parameter :: num_bc_in = 6 ! Number of input boundary conditions ! THere are no purely output boundary conditions integer, parameter :: num_bc_out = 0 ! Number of purely output boundary condtions @@ -114,15 +121,24 @@ module PRTAllometricCarbonMod ! ------------------------------------------------------------------------------------- - type, public, extends(prt_vartypes) :: callom_prt_vartypes + type, public, extends(prt_vartypes) :: callom_prt_vartypes - contains + contains - procedure :: DailyPRT => DailyPRTAllometricCarbon - procedure :: FastPRT => FastPRTAllometricCarbon + procedure :: DailyPRT => DailyPRTAllometricCarbon + procedure :: FastPRT => FastPRTAllometricCarbon end type callom_prt_vartypes - + + type, public, extends(prt_vartypes) :: csimpler_allom_prt_vartypes + + contains + + procedure :: DailyPRT => DailyPRTAllometricCarbonSimpler + procedure :: FastPRT => FastPRTAllometricCarbonSimpler + + end type csimpler_allom_prt_vartypes + ! ------------------------------------------------------------------------------------ ! ! This next class is an extention of the base instance that maps state variables @@ -313,6 +329,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] @@ -338,6 +357,9 @@ subroutine DailyPRTAllometricCarbon(this) real(r8) :: repro_c0 ! "" real(r8) :: struct_c0 ! "" + logical :: is_hydecid_dormant ! Flag to signal that the cohort is drought deciduous and dormant + logical :: is_deciduous ! Flag to signal this is a deciduous PFT + logical :: grow_struct logical :: grow_leaf ! Are leaves at allometric target and should be grown? logical :: grow_fnrt ! Are fine-roots at allometric target and should be grown? @@ -358,6 +380,10 @@ subroutine DailyPRTAllometricCarbon(this) integer :: leaf_status ! are leaves on (2) or off (1) real(r8) :: leaf_age_flux ! carbon mass flux between leaf age classification pools + real(r8) :: elongf_leaf ! Leaf elongation factor + real(r8) :: elongf_fnrt ! Fine-root "elongation factor" + real(r8) :: elongf_stem ! Stem "elongation factor" + ! Integegrator variables c_pool is "mostly" carbon variables, it also includes ! dbh... @@ -402,12 +428,24 @@ subroutine DailyPRTAllometricCarbon(this) canopy_trim = this%bc_in(ac_bc_in_id_ctrim)%rval ipft = this%bc_in(ac_bc_in_id_pft)%ival leaf_status = this%bc_in(ac_bc_in_id_lstat)%ival + elongf_leaf = this%bc_in(ac_bc_in_id_efleaf)%rval + elongf_fnrt = this%bc_in(ac_bc_in_id_effnrt)%rval + elongf_stem = this%bc_in(ac_bc_in_id_efstem)%rval - intgr_params(:) = un_initialized - intgr_params(ac_bc_in_id_ctrim) = this%bc_in(ac_bc_in_id_ctrim)%rval - intgr_params(ac_bc_in_id_pft) = real(this%bc_in(ac_bc_in_id_pft)%ival) - - + intgr_params(:) = un_initialized + intgr_params(ac_bc_in_id_ctrim) = this%bc_in(ac_bc_in_id_ctrim)%rval + intgr_params(ac_bc_in_id_pft) = real(this%bc_in(ac_bc_in_id_pft)%ival,r8) + intgr_params(ac_bc_in_id_lstat) = real(this%bc_in(ac_bc_in_id_lstat)%ival,r8) + intgr_params(ac_bc_in_id_efleaf) = this%bc_in(ac_bc_in_id_efleaf)%rval + intgr_params(ac_bc_in_id_effnrt) = this%bc_in(ac_bc_in_id_effnrt)%rval + intgr_params(ac_bc_in_id_efstem) = this%bc_in(ac_bc_in_id_efstem)%rval + + + ! Set some logical flags to simplify "if" blocks + is_hydecid_dormant = ( prt_params%stress_decid(ipft) == 1 ) .and. & + any(leaf_status == [leaves_off,leaves_pshed] ) + is_deciduous = ( prt_params%stress_decid(ipft) == 1 ) .or. & + ( prt_params%season_decid(ipft) == 1 ) nleafage = prt_global%state_descriptor(leaf_c_id)%num_pos ! Number of leaf age class @@ -432,12 +470,12 @@ subroutine DailyPRTAllometricCarbon(this) store_c0 = store_c ! Set initial storage carbon repro_c0 = repro_c ! Set initial reproductive carbon struct_c0 = struct_c ! Set initial structural carbon - + ! ----------------------------------------------------------------------------------- ! II. Calculate target size of the biomass compartment for a given dbh. ! ----------------------------------------------------------------------------------- - + ! Target sapwood biomass according to allometry and trimming [kgC] call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) @@ -451,12 +489,8 @@ 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 - call bleaf(dbh,ipft,canopy_trim,target_leaf_c) - else - target_leaf_c = 0._r8 - end if - + call bleaf(dbh,ipft,canopy_trim,target_leaf_c) + ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) @@ -464,45 +498,85 @@ subroutine DailyPRTAllometricCarbon(this) call bstore_allom(dbh,ipft,canopy_trim,target_store_c) + + ! ----------------------------------------------------------------------------------- + ! II 1/2. Update target biomass based on the leaf elongation factor and the abscission + ! fraction for each non-leaf tissue. Elongation factor is binary for + ! cold-deciduous and original drought-deciduous, and always one for + ! evergreens. In case the plant is shedding leaves, we impose that any + ! positive carbon balance necessarily goes to storage, even if this causes + ! storage to go above allometry. + ! ----------------------------------------------------------------------------------- + if (is_hydecid_dormant) then + target_leaf_c = 0.0_r8 + target_fnrt_c = 0.0_r8 + target_sapw_c = 0.0_r8 + target_struct_c = 0.0_r8 + target_store_c = target_store_c + max(0.0_r8,carbon_balance) + else + target_leaf_c = elongf_leaf * target_leaf_c + target_fnrt_c = elongf_fnrt * target_fnrt_c + target_sapw_c = elongf_stem * target_sapw_c + target_struct_c = elongf_stem * target_struct_c + end if + + + ! ----------------------------------------------------------------------------------- ! III. Prioritize some amount of carbon to replace leaf/root turnover ! Make sure it isnt a negative payment, and either pay what is available ! or forcefully pay from storage. + ! MLO. Added a few conditions to decide what to do in case plants are deciduous. + ! Specifically, drought deciduous with leaves off should not replace fine + ! roots. They will be in negative carbon balance, and unlike cold deciduous, + ! the turnover rates will be high during the dry season (turnover is + ! temperature-dependent, but not moisture-dependent). Allocating carbon + ! to high-maintanence tissues will drain the storage with little benefit for + ! these plants. ! ----------------------------------------------------------------------------------- - - if( prt_params%evergreen(ipft) ==1 ) then + if ( is_hydecid_dormant ) then + ! Drought deciduous, dormant state. Set demands to both leaves and roots to zero. + leaf_c_demand = 0.0_r8 + fnrt_c_demand = 0.0_r8 + elseif ( is_deciduous ) then + ! Either cold deciduous plant, or drought deciduous with leaves on. Maintain roots. + leaf_c_demand = 0.0_r8 + fnrt_c_demand = max(0.0_r8, & + prt_params%leaf_stor_priority(ipft)*this%variables(fnrt_c_id)%turnover(icd)) + else + ! Evergreen PFT. Try to meet demands for both leaves and fine roots. + ! If this is not evergreen, this PFT isn't expected by FATES, and we assume + ! evergreen. 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 + fnrt_c_demand = max(0.0_r8, & + prt_params%leaf_stor_priority(ipft)*this%variables(fnrt_c_id)%turnover(icd)) end if - - 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 - ! We pay this even if we don't have the carbon + + if ( total_c_demand > nearzero ) then + + ! We pay this even if we don't have positive carbon balance. ! 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 +585,23 @@ 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 +611,28 @@ 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)) + 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_c_demand = leaf_c_demand + fnrt_c_demand - - if( (carbon_balance > nearzero ) .and. (total_c_demand>nearzero)) then + 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 @@ -562,40 +644,32 @@ subroutine DailyPRTAllometricCarbon(this) ! ----------------------------------------------------------------------------------- if( carbon_balance > nearzero ) then - leaf_below_target = max(target_leaf_c - sum(leaf_c(1:nleafage)),0.0_r8) - fnrt_below_target = max(target_fnrt_c - fnrt_c,0.0_r8) - sapw_below_target = max(target_sapw_c - sapw_c,0.0_r8) - store_below_target = max(target_store_c - store_c,0.0_r8) - + 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) + sapw_below_target = max(0.0_r8,target_sapw_c - sapw_c) + store_below_target = max(0.0_r8,target_store_c - store_c) + total_below_target = leaf_below_target + fnrt_below_target + & 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 - + if ( total_below_target > nearzero ) then + ! 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 @@ -606,9 +680,9 @@ subroutine DailyPRTAllometricCarbon(this) if( carbon_balance > nearzero ) then - struct_below_target = max(target_struct_c - struct_c ,0.0_r8) + struct_below_target = max(0.0_r8,target_struct_c - struct_c) - if ( struct_below_target > 0.0_r8) then + if ( struct_below_target > nearzero) then struct_c_flux = min(carbon_balance,struct_below_target) carbon_balance = carbon_balance - struct_c_flux @@ -617,7 +691,27 @@ subroutine DailyPRTAllometricCarbon(this) end if end if - + + + + ! ----------------------------------------------------------------------------------- + ! VII 1/2: If plant is semi-deciduous, there will be cases in which plant's carbon + ! balance is positive but plant is losing leaves, in which case the plant + ! should not invest in growth. Likewise, when plants are flushing but the + ! elongation factor is small, we may have need to "shed" excess storage by + ! placing it as carbon_balance again (not the neatest solution, but other- + ! wise the solver may fail to find a solution). + ! ----------------------------------------------------------------------------------- + select_stash_grow: select case (leaf_status) + case (leaves_off,leaves_pshed) + ! There is carbon balance, but plant is shedding leaves. We stash the carbon + ! to storage even if it makes their storage too large. + store_c_flux = carbon_balance + carbon_balance = carbon_balance - store_c_flux + store_c = store_c + store_c_flux + end select select_stash_grow + + ! ----------------------------------------------------------------------------------- ! VIII. If carbon is yet still available ... ! Our pools are now either on allometry or above (from fusion). @@ -632,31 +726,21 @@ subroutine DailyPRTAllometricCarbon(this) ! the plant has not been brought to be "on allometry", it thinks it has carbon ! left to allocate, and thus it must be on allometry when its not. ! ----------------------------------------------------------------------------------- - if_stature_growth: if( carbon_balance > calloc_abs_error ) then - + + ! This routine checks that actual carbon is not below that targets. It does ! allow actual pools to be above the target, and in these cases, it sends ! a false on the "grow_<>" flag, allowing the plant to grow into these pools. ! It also checks to make sure that structural biomass is not above the target. - if( (target_store_c - store_c)>calloc_abs_error) then - write(fates_log(),*) 'storage is not on-allometry at the growth step' - write(fates_log(),*) 'exiting' - write(fates_log(),*) 'cbal: ',carbon_balance - write(fates_log(),*) 'near-zero',nearzero - write(fates_log(),*) 'store_c: ',store_c - write(fates_log(),*) 'target c: ',target_store_c - write(fates_log(),*) 'store_c0:', store_c0 - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - - call TargetAllometryCheck(sum(leaf_c(1:nleafage)), fnrt_c, sapw_c, & - store_c, struct_c, & - target_leaf_c, target_fnrt_c, & - target_sapw_c, target_store_c, target_struct_c, & - grow_struct, grow_leaf, grow_fnrt, grow_sapw, grow_store) + call TargetAllometryCheck(sum(leaf_c0(1:nleafage)),fnrt_c0,sapw_c0,store_c0,struct_c0, & + sum(leaf_c(1:nleafage)), fnrt_c, sapw_c,store_c, struct_c, & + target_leaf_c, target_fnrt_c, target_sapw_c, & + target_store_c, target_struct_c, & + carbon_balance, & + elongf_leaf,elongf_fnrt,elongf_stem,ipft,leaf_status, & + grow_leaf, grow_fnrt, grow_sapw, grow_store, grow_struct) ! -------------------------------------------------------------------------------- ! The numerical integration of growth requires that the instantaneous state @@ -689,16 +773,27 @@ subroutine DailyPRTAllometricCarbon(this) c_pool(repro_c_id) = repro_c 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 + ! Only grow leaves if we are in a "leaf-on" status. For drought-deciduous, we + ! interrupt growth for all tissues when in dormant mode. + if (is_hydecid_dormant) then + c_mask(leaf_c_id) = .false. + c_mask(fnrt_c_id) = .false. + c_mask(sapw_c_id) = .false. + c_mask(struct_c_id) = .false. + else - c_mask(leaf_c_id) = .false. + 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(struct_c_id) = grow_struct + end if - c_mask(fnrt_c_id) = grow_fnrt - c_mask(sapw_c_id) = grow_sapw c_mask(store_c_id) = grow_store - c_mask(struct_c_id) = grow_struct c_mask(repro_c_id) = .true. ! Always calculate reproduction on growth c_mask(dbh_id) = .true. ! Always increment dbh on growth step @@ -732,6 +827,7 @@ subroutine DailyPRTAllometricCarbon(this) ! we remember the current step size as a good next guess. call CheckIntegratedAllometries(c_pool_out(dbh_id),ipft,canopy_trim, & + elongf_leaf, elongf_fnrt, elongf_stem, & c_pool_out(leaf_c_id), c_pool_out(fnrt_c_id), c_pool_out(sapw_c_id), & c_pool_out(store_c_id), c_pool_out(struct_c_id), & c_mask(leaf_c_id), c_mask(fnrt_c_id), c_mask(sapw_c_id), & @@ -830,7 +926,6 @@ subroutine DailyPRTAllometricCarbon(this) end if if_step_pass end do do_solve_check - end if if_stature_growth ! Track the net allocations and transport from this routine @@ -862,224 +957,856 @@ subroutine DailyPRTAllometricCarbon(this) end subroutine DailyPRTAllometricCarbon ! ===================================================================================== + - function AllomCGrowthDeriv(c_pools,c_mask,cbalance,intgr_params) result(dCdx) + ! ===================================================================================== - ! --------------------------------------------------------------------------------- - ! This function calculates the derivatives for the carbon pools - ! relative to the amount of carbon balance. This function is based completely - ! off of allometry, and assumes that there are no other species (ie nutrients) that - ! govern allocation. - ! --------------------------------------------------------------------------------- - - ! Arguments - real(r8),intent(in), dimension(:) :: c_pools ! Vector of carbon pools - ! dbh,leaf,root,sap,store,dead - logical,intent(in), dimension(:) :: c_mask ! logical mask of active pools - ! some may be turned off - real(r8),intent(in) :: cbalance ! The carbon balance of the - ! partial step (independant var) - - real(r8), intent(in),dimension(:) :: intgr_params ! Generic Array used to pass - ! parameters into this function - ! Return Value - ! Change in carbon (each pool) per change in total allocatable carbon (kgC/kgC) - real(r8),dimension(lbound(c_pools,dim=1):ubound(c_pools,dim=1)) :: dCdx - ! locals - integer :: ipft ! PFT index - real(r8) :: canopy_trim ! Canopy trimming function (boundary condition [0-1] - real(r8) :: ct_leaf ! target leaf biomass, dummy var (kgC) - real(r8) :: ct_fnrt ! target fine-root biomass, dummy var (kgC) - real(r8) :: ct_sap ! target sapwood biomass, dummy var (kgC) - real(r8) :: ct_agw ! target aboveground wood, dummy var (kgC) - real(r8) :: ct_bgw ! target belowground wood, dummy var (kgC) - real(r8) :: ct_store ! target storage, dummy var (kgC) - real(r8) :: ct_dead ! target structural biomas, dummy var (kgC) - real(r8) :: sapw_area ! dummy sapwood area - real(r8) :: ct_dleafdd ! target leaf biomass derivative wrt diameter, (kgC/cm) - real(r8) :: ct_dfnrtdd ! target fine-root biomass derivative wrt diameter, (kgC/cm) - real(r8) :: ct_dsapdd ! target sapwood biomass derivative wrt diameter, (kgC/cm) - real(r8) :: ct_dagwdd ! target AG wood biomass derivative wrt diameter, (kgC/cm) - real(r8) :: ct_dbgwdd ! target BG wood biomass derivative wrt diameter, (kgC/cm) - real(r8) :: ct_dstoredd ! target storage biomass derivative wrt diameter, (kgC/cm) - real(r8) :: ct_ddeaddd ! target structural biomass derivative wrt diameter, (kgC/cm) - real(r8) :: ct_dtotaldd ! target total (not reproductive) biomass derivative wrt diameter, (kgC/cm) - real(r8) :: repro_fraction ! fraction of carbon balance directed towards reproduction (kgC/kgC) + subroutine DailyPRTAllometricCarbonSimpler(this) + ! ----------------------------------------------------------------------------------- + ! + ! This is the main routine that handles allocation associated with the 1st + ! hypothesis; carbon only, and growth governed by allometry. + ! MLO. This is a slightly simpler alternative I am testing. + ! + ! This routine is explained in the technical documentation in detail. + ! + ! Some points: + ! 1) dbh, while not a PARTEH "state variable", is passed in from FATES (or other + ! model), is integrated along with the mass based state variables, and then + ! passed back to the ecosystem model. It is a "inout" style boundary condition. + ! + ! 2) It is assumed that both growth respiration, and maintenance respiration + ! costs have already been paid, and therefore the "carbon_balance" boundary + ! condition is the net carbon gained by the plant over the coarse of the day. + ! Think of "daily integrated NPP". + ! + ! 3) This routine will completely spend carbon_balance if it enters as a positive + ! value, or replace carbon balance (using storage) if it enters as a negative value. + ! + ! 4) It is assumed that the ecosystem model calling this routine has ensured that + ! the net amount of negative carbon is no greater than that which can be replaced + ! by storage. This routine will crash gracefully if that is not true. + ! + ! 5) Unlike the original sub-routine, here we do not distinguish carbon lost through + ! maintenance from long-term carbon "debt" (i.e. biomass below allometry). We simply + ! seek to bring the plant back to allometry in case there is carbon to do so. We only + ! maintain the priority of replenishing living tissues (plus storage) over structural + ! tissues (which should be sapwood turnover in any case). + ! + ! 6) If there is any carbon balance left after bringing living tissues to allometry, + ! we try to bring heartwood back on allometry. + ! + ! 7) Finally, if carbon is yet still available, it will grow all pools out concurrently + ! including some to reproduction. + ! + ! ---------------------------------------------------------------------------------- - associate( dbh => c_pools(dbh_id), & - cleaf => c_pools(leaf_c_id), & - cfnrt => c_pools(fnrt_c_id), & - csap => c_pools(sapw_c_id), & - cstore => c_pools(store_c_id), & - cdead => c_pools(struct_c_id), & - crepro => c_pools(repro_c_id), & ! Unused (memoryless) - mask_dbh => c_mask(dbh_id), & ! Unused (dbh always grows) - mask_leaf => c_mask(leaf_c_id), & - mask_fnrt => c_mask(fnrt_c_id), & - mask_sap => c_mask(sapw_c_id), & - mask_store => c_mask(store_c_id), & - mask_dead => c_mask(struct_c_id), & ! Unused (dead always grows) - mask_repro => c_mask(repro_c_id) ) - canopy_trim = intgr_params(ac_bc_in_id_ctrim) - ipft = int(intgr_params(ac_bc_in_id_pft)) - + ! The class is the only argument + class(csimpler_allom_prt_vartypes) :: this ! this class - call bleaf(dbh,ipft,canopy_trim,ct_leaf,ct_dleafdd) - call bfineroot(dbh,ipft,canopy_trim,ct_fnrt,ct_dfnrtdd) - call bsap_allom(dbh,ipft,canopy_trim,sapw_area,ct_sap,ct_dsapdd) + ! ----------------------------------------------------------------------------------- + ! These are local copies of the in/out boundary condition structure + ! ----------------------------------------------------------------------------------- - call bagw_allom(dbh,ipft,ct_agw,ct_dagwdd) - call bbgw_allom(dbh,ipft,ct_bgw,ct_dbgwdd) - call bdead_allom(ct_agw,ct_bgw, ct_sap, ipft, ct_dead, & - ct_dagwdd, ct_dbgwdd, ct_dsapdd, ct_ddeaddd) - call bstore_allom(dbh,ipft,canopy_trim,ct_store,ct_dstoredd) - - ! fraction of carbon going towards reproduction - if (dbh <= prt_params%dbh_repro_threshold(ipft)) then ! cap on leaf biomass - repro_fraction = prt_params%seed_alloc(ipft) - else - repro_fraction = prt_params%seed_alloc(ipft) + prt_params%seed_alloc_mature(ipft) - end if + real(r8),pointer :: dbh ! Diameter at breast height [cm] + ! this local will point to both in and out bc's + real(r8),pointer :: carbon_balance ! Daily carbon balance for this cohort [kgC] - dCdx = 0.0_r8 + real(r8) :: canopy_trim ! The canopy trimming function [0-1] + integer :: ipft ! Plant Functional Type index - ct_dtotaldd = ct_ddeaddd - if (mask_leaf) ct_dtotaldd = ct_dtotaldd + ct_dleafdd - if (mask_fnrt) ct_dtotaldd = ct_dtotaldd + ct_dfnrtdd - if (mask_sap) ct_dtotaldd = ct_dtotaldd + ct_dsapdd - if (mask_store) ct_dtotaldd = ct_dtotaldd + ct_dstoredd - ! It is possible that with some asymptotic, or hard - ! capped allometries, that all growth rates reach zero. - ! In this case, if there is carbon, give it to reproduction + real(r8) :: target_leaf_c ! target leaf carbon [kgC] + real(r8) :: target_fnrt_c ! target fine-root carbon [kgC] + real(r8) :: target_sapw_c ! target sapwood carbon [kgC] + real(r8) :: target_store_c ! target storage carbon [kgC] + real(r8) :: target_agw_c ! target above ground carbon in woody tissues [kgC] + real(r8) :: target_bgw_c ! target below ground carbon in woody tissues [kgC] + real(r8) :: target_struct_c ! target structural carbon [kgC] - if(ct_dtotaldd<=tiny(ct_dtotaldd))then + real(r8) :: sapw_area ! dummy var, x-section area of sapwood [m2] - dCdx(struct_c_id) = 0.0_r8 - dCdx(dbh_id) = 0.0_r8 - dCdx(leaf_c_id) = 0.0_r8 - dCdx(fnrt_c_id) = 0.0_r8 - dCdx(sapw_c_id) = 0.0_r8 - dCdx(store_c_id) = 0.0_r8 - dCdx(repro_c_id) = 1.0_r8 + real(r8) :: leaf_below_target ! fineroot biomass below target amount [kgC] + real(r8) :: fnrt_below_target ! fineroot biomass below target amount [kgC] + real(r8) :: sapw_below_target ! sapwood biomass below target amount [kgC] + real(r8) :: store_below_target ! storage biomass below target amount [kgC] + real(r8) :: struct_below_target ! dead (structural) biomass below target amount [kgC] + real(r8) :: total_below_target ! total biomass below the allometric target [kgC] - else + real(r8) :: allocation_factor ! allocation factor (relative to demand) to + ! reconstruct tissues - dCdx(struct_c_id) = (ct_ddeaddd/ct_dtotaldd)*(1.0_r8-repro_fraction) - dCdx(dbh_id) = (1.0_r8/ct_dtotaldd)*(1.0_r8-repro_fraction) - - if (mask_leaf) then - dCdx(leaf_c_id) = (ct_dleafdd/ct_dtotaldd)*(1.0_r8-repro_fraction) - else - dCdx(leaf_c_id) = 0.0_r8 - end if - - if (mask_fnrt) then - dCdx(fnrt_c_id) = (ct_dfnrtdd/ct_dtotaldd)*(1.0_r8-repro_fraction) - else - dCdx(fnrt_c_id) = 0.0_r8 - end if - - if (mask_sap) then - dCdx(sapw_c_id) = (ct_dsapdd/ct_dtotaldd)*(1.0_r8-repro_fraction) - else - dCdx(sapw_c_id) = 0.0_r8 - end if - - if (mask_store) then - dCdx(store_c_id) = (ct_dstoredd/ct_dtotaldd)*(1.0_r8-repro_fraction) - else - dCdx(store_c_id) = 0.0_r8 - end if - - dCdx(repro_c_id) = repro_fraction + real(r8) :: flux_adj ! adjustment made to growth flux term to minimize error [kgC] - end if - - end associate + logical :: step_pass ! Did the integration step pass? - return + real(r8) :: leaf_c_flux ! Transfer into leaves at various stages [kgC] + real(r8) :: fnrt_c_flux ! Transfer into fine-roots at various stages [kgC] + real(r8) :: sapw_c_flux ! Transfer into sapwood at various stages [kgC] + real(r8) :: store_c_flux ! Transfer into storage at various stages [kgC] + real(r8) :: repro_c_flux ! Transfer into reproduction at the final stage [kgC] + real(r8) :: struct_c_flux ! Transfer into structure at various stages [kgC] + + real(r8),dimension(max_nleafage) :: leaf_c0 + + ! Initial value of carbon used to determine net flux + real(r8) :: fnrt_c0 ! during this routine + real(r8) :: sapw_c0 ! "" + real(r8) :: store_c0 ! "" + real(r8) :: repro_c0 ! "" + real(r8) :: struct_c0 ! "" + + logical :: is_hydecid_dormant ! Flag to signal that the cohort is drought deciduous and dormant + logical :: is_deciduous ! Flag to signal this is a deciduous PFT + + logical :: grow_struct + logical :: grow_leaf ! Are leaves at allometric target and should be grown? + logical :: grow_fnrt ! Are fine-roots at allometric target and should be grown? + logical :: grow_sapw ! Is sapwood at allometric target and should be grown? + logical :: grow_store ! Is storage at allometric target and should be grown? + + ! integrator variables + real(r8) :: deltaC ! trial value for substep + integer :: ierr ! error flag for allometric growth step + integer :: nsteps ! number of sub-steps + integer :: istep ! current substep index + real(r8) :: totalC ! total carbon allocated over alometric growth step + real(r8) :: hite_out ! dummy height variable + + integer :: i_var ! index for iterating state variables + integer :: i_age ! index for iterating leaf ages + integer :: nleafage ! number of leaf age classifications + integer :: leaf_status ! are leaves on or off? + real(r8) :: leaf_age_flux ! carbon mass flux between leaf age classification pools + + real(r8) :: elongf_leaf ! Leaf elongation factor + real(r8) :: elongf_fnrt ! Fine-root "elongation factor" + real(r8) :: elongf_stem ! Stem "elongation factor" + + + ! Integrator variables c_pool are "mostly" carbon variables, but c_pool also includes + ! dbh... + ! ----------------------------------------------------------------------------------- + + real(r8),dimension(n_integration_vars) :: c_pool ! Vector of carbon pools passed to integrator + real(r8),dimension(n_integration_vars) :: c_pool_out ! Vector of carbon pools passed back from integrator + logical,dimension(n_integration_vars) :: c_mask ! Mask of active pools during integration + + integer , parameter :: max_substeps = 300 ! Maximum allowable iterations + + real(r8), parameter :: max_trunc_error = 1.0_r8 ! Maximum allowable truncation error + + integer, parameter :: ODESolve = 2 ! 1=RKF45, 2=Euler + + integer, parameter :: iexp_leaf = 1 ! index 1 is the expanding (i.e. youngest) + ! leaf age class, and therefore + ! all new allocation goes into that pool + character(len= 9), parameter :: fmti = '(a,1x,i5)' + character(len=13), parameter :: fmt0 = '(a,1x,es12.5)' + character(len=19), parameter :: fmth = '(a,1x,a5,3(1x,a12))' + character(len=22), parameter :: fmtg = '(a,5x,l1,3(1x,es12.5))' + + real(r8) :: intgr_params(num_bc_in) ! The boundary conditions to this routine, + ! are pressed into an array that is also + ! passed to the integrators + + associate( & + leaf_c => this%variables(leaf_c_id)%val, & + fnrt_c => this%variables(fnrt_c_id)%val(icd), & + sapw_c => this%variables(sapw_c_id)%val(icd), & + store_c => this%variables(store_c_id)%val(icd), & + repro_c => this%variables(repro_c_id)%val(icd), & + struct_c => this%variables(struct_c_id)%val(icd)) + + + ! ----------------------------------------------------------------------------------- + ! 0. + ! Copy the boundary conditions into readable local variables. + ! We don't use pointers for bc's that ar "in" only, only "in-out" and "out" + ! ----------------------------------------------------------------------------------- + + dbh => this%bc_inout(ac_bc_inout_id_dbh)%rval + carbon_balance => this%bc_inout(ac_bc_inout_id_netdc)%rval + + canopy_trim = this%bc_in(ac_bc_in_id_ctrim)%rval + ipft = this%bc_in(ac_bc_in_id_pft)%ival + leaf_status = this%bc_in(ac_bc_in_id_lstat)%ival + elongf_leaf = this%bc_in(ac_bc_in_id_efleaf)%rval + elongf_fnrt = this%bc_in(ac_bc_in_id_effnrt)%rval + elongf_stem = this%bc_in(ac_bc_in_id_efstem)%rval + + intgr_params(:) = un_initialized + intgr_params(ac_bc_in_id_ctrim) = this%bc_in(ac_bc_in_id_ctrim)%rval + intgr_params(ac_bc_in_id_pft) = real(this%bc_in(ac_bc_in_id_pft)%ival,r8) + intgr_params(ac_bc_in_id_lstat) = real(this%bc_in(ac_bc_in_id_lstat)%ival,r8) + intgr_params(ac_bc_in_id_efleaf) = this%bc_in(ac_bc_in_id_efleaf)%rval + intgr_params(ac_bc_in_id_effnrt) = this%bc_in(ac_bc_in_id_effnrt)%rval + intgr_params(ac_bc_in_id_efstem) = this%bc_in(ac_bc_in_id_efstem)%rval + + ! Set some logical flags to simplify "if" blocks + is_hydecid_dormant = ( prt_params%stress_decid(ipft) == 1 ) .and. & + any( leaf_status == [leaves_off,leaves_pshed] ) + is_deciduous = ( prt_params%stress_decid(ipft) == 1 ) .or. & + ( prt_params%season_decid(ipft) == 1 ) + + + nleafage = prt_global%state_descriptor(leaf_c_id)%num_pos ! Number of leaf age class + + ! ----------------------------------------------------------------------------------- + ! Call the routine that advances leaves in age. + ! This will move a portion of the leaf mass in each + ! age bin, to the next bin. This will not handle movement + ! of mass from the oldest bin into the litter pool, that is something else. + ! ----------------------------------------------------------------------------------- + + call this%AgeLeaves(ipft,sec_per_day) + + ! ----------------------------------------------------------------------------------- + ! I. Remember the values for the state variables at the beginning of this + ! routines. We will then use that to determine their net allocation and reactive + ! transport flux "%net_alloc" at the end. + ! ----------------------------------------------------------------------------------- + + leaf_c0(1:nleafage) = leaf_c(1:nleafage) ! Set initial leaf carbon + fnrt_c0 = fnrt_c ! Set initial fine-root carbon + sapw_c0 = sapw_c ! Set initial sapwood carbon + store_c0 = store_c ! Set initial storage carbon + repro_c0 = repro_c ! Set initial reproductive carbon + struct_c0 = struct_c ! Set initial structural carbon + + + ! ----------------------------------------------------------------------------------- + ! II. Calculate target size of the biomass compartment for a given dbh. + ! ----------------------------------------------------------------------------------- + ! Target sapwood biomass according to allometry and trimming [kgC] + call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) + + ! Target total above ground biomass in woody/fibrous tissues [kgC] + call bagw_allom(dbh,ipft,target_agw_c) + + ! Target total below ground biomass in woody/fibrous tissues [kgC] + call bbgw_allom(dbh,ipft,target_bgw_c) + + ! Target total dead (structrual) biomass [kgC] + call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) + + ! Target leaf biomass according to allometry and trimming + call bleaf(dbh,ipft,canopy_trim,target_leaf_c) + + ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] + call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) + + ! Target storage carbon [kgC,kgC/cm] + call bstore_allom(dbh,ipft,canopy_trim,target_store_c) + + + + ! ----------------------------------------------------------------------------------- + ! II 1/2. Update target biomass based on the leaf elongation factor and the abscission + ! fraction for each non-leaf tissue. Elongation factor is binary for + ! cold-deciduous and original drought-deciduous, and always one for + ! evergreens. In case the plant is shedding leaves, we impose that any + ! positive carbon balance necessarily goes to storage, even if this causes + ! storage to go above allometry. + ! ----------------------------------------------------------------------------------- + if (is_hydecid_dormant) then + target_leaf_c = 0.0_r8 + target_fnrt_c = 0.0_r8 + target_sapw_c = 0.0_r8 + target_struct_c = 0.0_r8 + target_store_c = target_store_c + max(0.0_r8,carbon_balance) + else + target_leaf_c = elongf_leaf * target_leaf_c + target_fnrt_c = elongf_fnrt * target_fnrt_c + target_sapw_c = elongf_stem * target_sapw_c + target_struct_c = elongf_stem * target_struct_c + end if + + + ! ----------------------------------------------------------------------------------- + ! III. If carbon is available, bring all the pools as close to the allometry + ! as possible. This also includes the storage pool, even though carbon may + ! be drawn from the storage. + ! ----------------------------------------------------------------------------------- + + ! Identify living organs (plus storage) that are under target. Priority is given + ! to the organs (or storage) that are the most depleted, without a pre-determined + ! sequence. + 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 ) + sapw_below_target = max( 0.0_r8, target_sapw_c - sapw_c ) + store_below_target = max( 0.0_r8, target_store_c - store_c ) + struct_below_target = max( 0.0_r8, target_struct_c - struct_c ) + total_below_target = leaf_below_target + fnrt_below_target + sapw_below_target + & + store_below_target + struct_below_target + + replenish_allom_check: if ( total_below_target > nearzero ) then + ! Available carbon for transfer is the sum of stored carbon and the daily + ! carbon balance. + allocation_factor = min(1.0_r8, (store_c + carbon_balance) / total_below_target ) + + ! Scale flux so pools can be replenished simultaneously. + 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 + struct_c_flux = struct_below_target * allocation_factor + + ! Replenish pools + 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 + struct_c = struct_c + struct_c_flux + + ! Update carbon balance + carbon_balance = carbon_balance - & + ( leaf_c_flux + fnrt_c_flux + sapw_c_flux + & + store_c_flux + struct_c_flux ) + + end if replenish_allom_check + + + ! ----------------------------------------------------------------------------------- + ! IV. If carbon balance is negative, reduce the storage pool. Otherwise, try to + ! fill the storage pool before growing. + ! ----------------------------------------------------------------------------------- + update_storage: if ( carbon_balance < 0.0_r8 ) then + ! If carbon balance is negative, store_c will be depleted. + store_c_flux = carbon_balance + store_c = store_c + store_c_flux + + ! After this operation, carbon_balance should be zero. + carbon_balance = carbon_balance - store_c_flux + + else + ! Non-negative carbon balance. Use left over carbon to fill the storage + ! pool and try to bring it to allometry before trying to grow in size. + store_below_target = max(0.0_r8,target_store_c - store_c) + store_c_flux = min( store_below_target,carbon_balance) + store_c = store_c + store_c_flux + + ! Update carbon balance + carbon_balance = carbon_balance - store_c_flux + + end if update_storage + + + ! ----------------------------------------------------------------------------------- + ! V. If carbon is yet still available ... + ! Our pools are now either on allometry or above (from fusion). + ! We we can increment those pools at or below, + ! including structure and reproduction according to their rates + ! Use an adaptive euler integration. If the error is not nominal, + ! the carbon balance sub-step (deltaC) will be halved and tried again + ! + ! Note that we compare against calloc_abs_error here because it is possible + ! that all the carbon was effectively used up, but a miniscule amount + ! remains due to numerical precision (ie -20 or so), so even though + ! the plant has not been brought to be "on allometry", it thinks it has carbon + ! left to allocate, and thus it must be on allometry when its not. + ! ----------------------------------------------------------------------------------- + if_stature_growth: if ( carbon_balance > calloc_abs_error ) then + + ! This routine checks that actual carbon is not below that targets. It does + ! allow actual pools to be above the target, and in these cases, it sends + ! a false on the "grow_<>" flag, allowing the plant to grow into these pools. + ! It also checks to make sure that structural biomass is not above the target. + call TargetAllometryCheck(sum(leaf_c0(1:nleafage)),fnrt_c0,sapw_c0,store_c0,struct_c0, & + sum(leaf_c(1:nleafage)), fnrt_c, sapw_c,store_c, struct_c, & + target_leaf_c, target_fnrt_c, target_sapw_c, & + target_store_c, target_struct_c, & + carbon_balance, & + elongf_leaf,elongf_fnrt,elongf_stem,ipft,leaf_status, & + grow_leaf, grow_fnrt, grow_sapw, grow_store, grow_struct) + + ! -------------------------------------------------------------------------------- + ! The numerical integration of growth requires that the instantaneous state + ! variables are passed in as an array. We call it "c_pool". + ! + ! Initialize the adaptive integrator arrays and flags + ! -------------------------------------------------------------------------------- + + ierr = 1 + totalC = carbon_balance + nsteps = 0 + + c_pool(:) = 0.0_r8 ! Zero state variable array + c_mask(:) = .false. ! This mask tells the integrator + ! which indices are active. Its possible + ! that due to fusion, or previous numerical + ! truncation errors, that one of these pools + ! may be larger than its target! We check + ! this, and if true, then we flag that + ! pool to be ignored. c_mask(i) = .false. + ! For grasses, since they don't grow very + ! large and thus won't accumulate such large + ! errors, we always mask as true. + + c_pool(leaf_c_id) = sum(leaf_c(1:nleafage)) + c_pool(fnrt_c_id) = fnrt_c + c_pool(sapw_c_id) = sapw_c + c_pool(store_c_id) = store_c + c_pool(struct_c_id) = struct_c + c_pool(repro_c_id) = repro_c + c_pool(dbh_id) = dbh + + ! Only grow leaves if we are in a "leaf-on" status. For drought-deciduous, we + ! interrupt growth for all tissues when in dormant mode. + if (is_hydecid_dormant) then + c_mask(leaf_c_id) = .false. + c_mask(fnrt_c_id) = .false. + c_mask(sapw_c_id) = .false. + c_mask(struct_c_id) = .false. + + else + 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(struct_c_id) = grow_struct + + end if + c_mask(store_c_id) = grow_store + c_mask(repro_c_id) = .true. ! Always calculate reproduction on growth + c_mask(dbh_id) = .true. ! Always increment dbh on growth step + + + ! When using the Euler method, we keep things simple. We always try + ! to make the first integration step to span the entirety of the integration + ! window for the independent variable (available carbon) + + if(ODESolve == 2) then + this%ode_opt_step = totalC + end if + + do_solve_check: do while( ierr .ne. 0 ) + + deltaC = min(totalC,this%ode_opt_step) + sel_ODESolve: select case (ODESolve) + case (1) + call RKF45(AllomCGrowthDeriv,c_pool,c_mask,deltaC,totalC, & + max_trunc_error,intgr_params,c_pool_out,this%ode_opt_step,step_pass) + + case (2) + call Euler(AllomCGrowthDeriv,c_pool,c_mask,deltaC,totalC,intgr_params,c_pool_out) + ! step_pass = .true. + + ! When integrating along the allometric curve, we have the luxury of perfect + ! hindsite. Ie, after we have made our step, we can see if the amount + ! of each carbon we have matches the target associated with the new dbh. + ! The following call evaluates how close we are to the allometically defined + ! targets. If we are too far (governed by max_trunc_error), then we + ! pass back the pass/fail flag (step_pass) as false. If false, then + ! we halve the step-size, and then retry. If that step was fine, then + ! we remember the current step size as a good next guess. + + call CheckIntegratedAllometries(c_pool_out(dbh_id),ipft,canopy_trim, & + elongf_leaf, elongf_fnrt, elongf_stem, & + c_pool_out(leaf_c_id), c_pool_out(fnrt_c_id), c_pool_out(sapw_c_id), & + c_pool_out(store_c_id), c_pool_out(struct_c_id), & + c_mask(leaf_c_id), c_mask(fnrt_c_id), c_mask(sapw_c_id), & + c_mask(store_c_id),c_mask(struct_c_id), max_trunc_error, step_pass) + + if(step_pass) then + this%ode_opt_step = deltaC + else + this%ode_opt_step = 0.5*deltaC + end if + case default + write(fates_log(),*) 'An integrator was chosen that does not exist' + write(fates_log(),*) 'ODESolve = ',ODESolve + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select sel_ODESolve + + nsteps = nsteps + 1 + + if (step_pass) then ! If true, then step is accepted + totalC = totalC - deltaC + c_pool(:) = c_pool_out(:) + end if + + if(nsteps > max_substeps ) then + write(fates_log(),fmt=*) '---~---' + write(fates_log(),fmt=*) 'Plant Growth Integrator could not find' + write(fates_log(),fmt=*) 'a solution in less than ',max_substeps,' tries.' + write(fates_log(),fmt=*) 'Aborting!' + write(fates_log(),fmt=*) '---~---' + write(fates_log(),fmt=fmti) 'Leaf status =',leaf_status + write(fates_log(),fmt=fmt0) 'Carbon_balance =',carbon_balance + write(fates_log(),fmt=fmt0) 'Elongf_leaf =',elongf_leaf + write(fates_log(),fmt=fmt0) 'Elongf_fnrt =',elongf_fnrt + write(fates_log(),fmt=fmt0) 'Elongf_stem =',elongf_stem + write(fates_log(),fmt=fmt0) 'deltaC =',deltaC + write(fates_log(),fmt=fmt0) 'totalC =',totalC + write(fates_log(),fmt=fmth) ' Tissue |', ' Grow',' Current',' Target' ,' Deficit' + write(fates_log(),fmt=fmtg) ' Leaf |', grow_leaf , sum(leaf_c(:)),target_leaf_c , target_leaf_c - sum(leaf_c(:)) + write(fates_log(),fmt=fmtg) ' Fine root |', grow_fnrt , fnrt_c,target_fnrt_c , target_fnrt_c - fnrt_c + write(fates_log(),fmt=fmtg) ' Sapwood |', grow_sapw , sapw_c,target_sapw_c , target_sapw_c - sapw_c + write(fates_log(),fmt=fmtg) ' Storage |', grow_store , store_c,target_store_c , target_store_c - store_c + write(fates_log(),fmt=fmtg) ' Structural |', grow_struct , struct_c,target_struct_c, target_struct_c - struct_c + write(fates_log(),fmt=*) '---~---' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! + ! TotalC should eventually be whittled down to near zero. + ! The solvers are not perfect, so we can't expect it to be perfectly zero. + ! Note that calloc_abs_error is 1e-9, which is really small (1 microgram of carbon) + ! yet also six orders of magnitude greater than typical rounding errors (~1e-15). + + ! At that point, update the actual states + ! -------------------------------------------------------------------------------- + if_step_pass: if( (totalC < calloc_abs_error) .and. (step_pass) )then + + ierr = 0 + leaf_c_flux = c_pool(leaf_c_id) - sum(leaf_c(1:nleafage)) + fnrt_c_flux = c_pool(fnrt_c_id) - fnrt_c + sapw_c_flux = c_pool(sapw_c_id) - sapw_c + store_c_flux = c_pool(store_c_id) - store_c + struct_c_flux = c_pool(struct_c_id) - struct_c + repro_c_flux = c_pool(repro_c_id) - repro_c + + ! Make an adjustment to flux partitions to make it match remaining c balance + flux_adj = carbon_balance/(leaf_c_flux+fnrt_c_flux+sapw_c_flux + & + store_c_flux+struct_c_flux+repro_c_flux) + + + leaf_c_flux = leaf_c_flux*flux_adj + fnrt_c_flux = fnrt_c_flux*flux_adj + sapw_c_flux = sapw_c_flux*flux_adj + store_c_flux = store_c_flux*flux_adj + struct_c_flux = struct_c_flux*flux_adj + repro_c_flux = repro_c_flux*flux_adj + + 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 + + carbon_balance = carbon_balance - struct_c_flux + struct_c = struct_c + struct_c_flux + + carbon_balance = carbon_balance - repro_c_flux + repro_c = repro_c + repro_c_flux + + dbh = c_pool(dbh_id) + + if( abs(carbon_balance)>calloc_abs_error ) then + write(fates_log(),*) 'carbon conservation error while integrating pools' + write(fates_log(),*) 'along alometric curve' + write(fates_log(),*) 'carbon_balance = ',carbon_balance,totalC + write(fates_log(),*) 'exiting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + end if if_step_pass + + end do do_solve_check + + end if if_stature_growth + + ! Track the net allocations and transport from this routine + ! (the AgeLeaves() routine handled tracking allocation through aging) + + this%variables(leaf_c_id)%net_alloc(icd) = & + this%variables(leaf_c_id)%net_alloc(icd) + (leaf_c(icd) - leaf_c0(icd)) + + this%variables(fnrt_c_id)%net_alloc(icd) = & + this%variables(fnrt_c_id)%net_alloc(icd) + (fnrt_c - fnrt_c0) + + this%variables(sapw_c_id)%net_alloc(icd) = & + this%variables(sapw_c_id)%net_alloc(icd) + (sapw_c - sapw_c0) + + this%variables(store_c_id)%net_alloc(icd) = & + this%variables(store_c_id)%net_alloc(icd) + (store_c - store_c0) + + this%variables(repro_c_id)%net_alloc(icd) = & + this%variables(repro_c_id)%net_alloc(icd) + (repro_c - repro_c0) + + this%variables(struct_c_id)%net_alloc(icd) = & + this%variables(struct_c_id)%net_alloc(icd) + (struct_c - struct_c0) + + + + end associate + + return + end subroutine DailyPRTAllometricCarbonSimpler + + ! ===================================================================================== + + + function AllomCGrowthDeriv(c_pools,c_mask,cbalance,intgr_params) result(dCdx) + + ! --------------------------------------------------------------------------------- + ! This function calculates the derivatives for the carbon pools + ! relative to the amount of carbon balance. This function is based completely + ! off of allometry, and assumes that there are no other species (ie nutrients) that + ! govern allocation. + ! --------------------------------------------------------------------------------- + + ! Arguments + real(r8),intent(in), dimension(:) :: c_pools ! Vector of carbon pools + ! dbh,leaf,root,sap,store,dead + logical,intent(in), dimension(:) :: c_mask ! logical mask of active pools + ! some may be turned off + real(r8),intent(in) :: cbalance ! The carbon balance of the + ! partial step (independant var) + + real(r8), intent(in),dimension(:) :: intgr_params ! Generic Array used to pass + ! parameters into this function + + + ! Return Value + ! Change in carbon (each pool) per change in total allocatable carbon (kgC/kgC) + real(r8),dimension(lbound(c_pools,dim=1):ubound(c_pools,dim=1)) :: dCdx + + ! locals + integer :: ipft ! PFT index + real(r8) :: canopy_trim ! Canopy trimming function (boundary condition) [0-1] + real(r8) :: elongf_leaf ! Leaf elongation factor (boundary condition) [0-1] + real(r8) :: elongf_fnrt ! Fine-root "elongation factor" (boundary condition) [0-1] + real(r8) :: elongf_stem ! Stem "elongation factor" (boundary condition) [0-1] + real(r8) :: ct_leaf ! target leaf biomass, dummy var (kgC) + real(r8) :: ct_fnrt ! target fine-root biomass, dummy var (kgC) + real(r8) :: ct_sap ! target sapwood biomass, dummy var (kgC) + real(r8) :: ct_agw ! target aboveground wood, dummy var (kgC) + real(r8) :: ct_bgw ! target belowground wood, dummy var (kgC) + real(r8) :: ct_store ! target storage, dummy var (kgC) + real(r8) :: ct_dead ! target structural biomas, dummy var (kgC) + real(r8) :: sapw_area ! dummy sapwood area + real(r8) :: ct_dleafdd ! target leaf biomass derivative wrt diameter, (kgC/cm) + real(r8) :: ct_dfnrtdd ! target fine-root biomass derivative wrt diameter, (kgC/cm) + real(r8) :: ct_dsapdd ! target sapwood biomass derivative wrt diameter, (kgC/cm) + real(r8) :: ct_dagwdd ! target AG wood biomass derivative wrt diameter, (kgC/cm) + real(r8) :: ct_dbgwdd ! target BG wood biomass derivative wrt diameter, (kgC/cm) + real(r8) :: ct_dstoredd ! target storage biomass derivative wrt diameter, (kgC/cm) + real(r8) :: ct_ddeaddd ! target structural biomass derivative wrt diameter, (kgC/cm) + real(r8) :: ct_dtotaldd ! target total (not reproductive) biomass derivative wrt diameter, (kgC/cm) + real(r8) :: repro_fraction ! fraction of carbon balance directed towards reproduction (kgC/kgC) + + + associate( dbh => c_pools(dbh_id), & + cleaf => c_pools(leaf_c_id), & + cfnrt => c_pools(fnrt_c_id), & + csap => c_pools(sapw_c_id), & + cstore => c_pools(store_c_id), & + cdead => c_pools(struct_c_id), & + crepro => c_pools(repro_c_id), & ! Unused (memoryless) + mask_dbh => c_mask(dbh_id), & ! Unused (dbh always grows) + mask_leaf => c_mask(leaf_c_id), & + mask_fnrt => c_mask(fnrt_c_id), & + mask_sap => c_mask(sapw_c_id), & + mask_store => c_mask(store_c_id), & + mask_dead => c_mask(struct_c_id), & ! Unused (dead always grows) + mask_repro => c_mask(repro_c_id) ) + + canopy_trim = intgr_params(ac_bc_in_id_ctrim) + ipft = int(intgr_params(ac_bc_in_id_pft)) + elongf_leaf = intgr_params(ac_bc_in_id_efleaf) + elongf_fnrt = intgr_params(ac_bc_in_id_effnrt) + elongf_stem = intgr_params(ac_bc_in_id_efstem) + + call bleaf(dbh,ipft,canopy_trim,ct_leaf,ct_dleafdd) + call bfineroot(dbh,ipft,canopy_trim,ct_fnrt,ct_dfnrtdd) + call bsap_allom(dbh,ipft,canopy_trim,sapw_area,ct_sap,ct_dsapdd) + + call bagw_allom(dbh,ipft,ct_agw,ct_dagwdd) + call bbgw_allom(dbh,ipft,ct_bgw,ct_dbgwdd) + call bdead_allom(ct_agw,ct_bgw, ct_sap, ipft, ct_dead, & + ct_dagwdd, ct_dbgwdd, ct_dsapdd, ct_ddeaddd) + call bstore_allom(dbh,ipft,canopy_trim,ct_store,ct_dstoredd) + + ! Apply elongation factor correction to targets + ct_leaf = elongf_leaf * ct_leaf + ct_fnrt = elongf_fnrt * ct_fnrt + ct_sap = elongf_stem * ct_sap + ct_dead = elongf_stem * ct_dead + ! MLO - Need to double check that it is correct to multiply derivatives too. + ct_dleafdd = elongf_leaf * ct_dleafdd + ct_dfnrtdd = elongf_fnrt * ct_dfnrtdd + ct_dsapdd = elongf_stem * ct_dsapdd + ct_ddeaddd = elongf_stem * ct_ddeaddd + + ! fraction of carbon going towards reproduction + if (dbh <= prt_params%dbh_repro_threshold(ipft)) then ! cap on leaf biomass + repro_fraction = prt_params%seed_alloc(ipft) + else + repro_fraction = prt_params%seed_alloc(ipft) + prt_params%seed_alloc_mature(ipft) + end if + + dCdx = 0.0_r8 + + ct_dtotaldd = ct_ddeaddd + if (mask_leaf) ct_dtotaldd = ct_dtotaldd + ct_dleafdd + if (mask_fnrt) ct_dtotaldd = ct_dtotaldd + ct_dfnrtdd + if (mask_sap) ct_dtotaldd = ct_dtotaldd + ct_dsapdd + if (mask_store) ct_dtotaldd = ct_dtotaldd + ct_dstoredd + + ! It is possible that with some asymptotic, or hard + ! capped allometries, that all growth rates reach zero. + ! In this case, if there is carbon, give it to reproduction + + if(ct_dtotaldd<=tiny(ct_dtotaldd))then + + dCdx(struct_c_id) = 0.0_r8 + dCdx(dbh_id) = 0.0_r8 + dCdx(leaf_c_id) = 0.0_r8 + dCdx(fnrt_c_id) = 0.0_r8 + dCdx(sapw_c_id) = 0.0_r8 + dCdx(store_c_id) = 0.0_r8 + dCdx(repro_c_id) = 1.0_r8 + + else + + dCdx(struct_c_id) = (ct_ddeaddd/ct_dtotaldd)*(1.0_r8-repro_fraction) + dCdx(dbh_id) = (1.0_r8/ct_dtotaldd)*(1.0_r8-repro_fraction) + + if (mask_leaf) then + dCdx(leaf_c_id) = (ct_dleafdd/ct_dtotaldd)*(1.0_r8-repro_fraction) + else + dCdx(leaf_c_id) = 0.0_r8 + end if + + if (mask_fnrt) then + dCdx(fnrt_c_id) = (ct_dfnrtdd/ct_dtotaldd)*(1.0_r8-repro_fraction) + else + dCdx(fnrt_c_id) = 0.0_r8 + end if + + if (mask_sap) then + dCdx(sapw_c_id) = (ct_dsapdd/ct_dtotaldd)*(1.0_r8-repro_fraction) + else + dCdx(sapw_c_id) = 0.0_r8 + end if + + if (mask_store) then + dCdx(store_c_id) = (ct_dstoredd/ct_dtotaldd)*(1.0_r8-repro_fraction) + else + dCdx(store_c_id) = 0.0_r8 + end if + + dCdx(repro_c_id) = repro_fraction + + end if + + end associate + + return end function AllomCGrowthDeriv ! ==================================================================================== - subroutine TargetAllometryCheck(bleaf,bfroot,bsap,bstore,bdead, & - bt_leaf,bt_froot,bt_sap,bt_store,bt_dead, & - grow_dead,grow_leaf,grow_froot,grow_sapw,grow_store) + subroutine TargetAllometryCheck(b0_leaf,b0_fnrt,b0_sapw,b0_store,b0_struct, & + bleaf,bfnrt,bsapw,bstore,bstruct, & + bt_leaf,bt_fnrt,bt_sapw,bt_store,bt_struct, & + carbon_balance, & + elongf_leaf,elongf_fnrt,elongf_stem,ipft,leaf_status, & + grow_leaf,grow_fnrt,grow_sapw,grow_store,grow_struct) ! Arguments - real(r8),intent(in) :: bleaf !actual - real(r8),intent(in) :: bfroot - real(r8),intent(in) :: bsap + real(r8),intent(in) :: b0_leaf !initial + real(r8),intent(in) :: b0_fnrt + real(r8),intent(in) :: b0_sapw + real(r8),intent(in) :: b0_store + real(r8),intent(in) :: b0_struct + real(r8),intent(in) :: bleaf !actual + real(r8),intent(in) :: bfnrt + real(r8),intent(in) :: bsapw real(r8),intent(in) :: bstore - real(r8),intent(in) :: bdead - real(r8),intent(in) :: bt_leaf !target - real(r8),intent(in) :: bt_froot - real(r8),intent(in) :: bt_sap + real(r8),intent(in) :: bstruct + real(r8),intent(in) :: bt_leaf !target + real(r8),intent(in) :: bt_fnrt + real(r8),intent(in) :: bt_sapw real(r8),intent(in) :: bt_store - real(r8),intent(in) :: bt_dead - logical,intent(out) :: grow_leaf !growth flag - logical,intent(out) :: grow_froot + real(r8),intent(in) :: bt_struct + real(r8),intent(in) :: carbon_balance !remaining carbon balance + real(r8),intent(in) :: elongf_leaf !elongation factors + real(r8),intent(in) :: elongf_fnrt + real(r8),intent(in) :: elongf_stem + integer,intent(in) :: ipft !Plant functional type + integer,intent(in) :: leaf_status !Phenology status + logical,intent(out) :: grow_leaf !growth flag + logical,intent(out) :: grow_fnrt logical,intent(out) :: grow_sapw logical,intent(out) :: grow_store - logical,intent(out) :: grow_dead - - if( (bt_leaf - bleaf)>calloc_abs_error) then - write(fates_log(),*) 'leaves are not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bleaf,bt_leaf - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( (bleaf - bt_leaf)>calloc_abs_error) then - ! leaf is above allometry, ignore - grow_leaf = .false. - else - grow_leaf = .true. - end if - - if( (bt_froot - bfroot)>calloc_abs_error) then - write(fates_log(),*) 'fineroots are not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bfroot, bt_froot - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bfroot-bt_froot)>calloc_abs_error ) then - grow_froot = .false. - else - grow_froot = .true. - end if - - if( (bt_sap - bsap)>calloc_abs_error) then - write(fates_log(),*) 'sapwood is not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bsap, bt_sap - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bsap-bt_sap)>calloc_abs_error ) then - grow_sapw = .false. - else - grow_sapw = .true. - end if - - if( (bt_store - bstore)>calloc_abs_error) then - write(fates_log(),*) 'storage is not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bstore,bt_store - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bstore-bt_store)>calloc_abs_error ) then - grow_store = .false. + logical,intent(out) :: grow_struct + ! Local variables + logical :: fine_leaf + logical :: fine_fnrt + logical :: fine_sapw + logical :: fine_store + logical :: fine_struct + logical :: all_fine + ! Local constants + character(len= 3), parameter :: fmth = '(a)' + character(len=27), parameter :: fmtb = '(a,3(1x,es12.5,1x,a),1x,l1)' + character(len=13), parameter :: fmte = '(a,1x,es12.5)' + character(len=10), parameter :: fmti = '(a,1x,i12)' + + + ! First test whether or not each pool looks reasonable. + fine_leaf = (bt_leaf - bleaf ) <= calloc_abs_error + fine_fnrt = (bt_fnrt - bfnrt ) <= calloc_abs_error + fine_sapw = (bt_sapw - bsapw ) <= calloc_abs_error + fine_store = (bt_store - bstore ) <= calloc_abs_error + fine_struct = (bt_struct - bstruct) <= calloc_abs_error + all_fine = fine_leaf .and. fine_fnrt .and. fine_sapw .and. & + fine_store .and. fine_struct + + ! Decide whether or not to grow tissues (but only if all tissues look fine). + ! We grow only when biomass is less than target biomass (with tolerance). + if (all_fine) then + grow_leaf = ( bleaf - bt_leaf ) <= calloc_abs_error + grow_fnrt = ( bfnrt - bt_fnrt ) <= calloc_abs_error + grow_sapw = ( bsapw - bt_sapw ) <= calloc_abs_error + grow_store = ( bstore - bt_store ) <= calloc_abs_error + grow_struct = ( bstruct - bt_struct ) <= calloc_abs_error else - grow_store = .true. - end if - - if( (bt_dead - bdead)>calloc_abs_error) then - write(fates_log(),*) 'structure not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bdead,bt_dead + ! If anything looks not fine, write a detailed report + write(fates_log(),fmt=fmth) '======' + write(fates_log(),fmt=fmth) ' At least one tissue is not on-allometry at the growth step' + write(fates_log(),fmt=fmth) '======' + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) ' Biomass and on-allometry test (''F'' means problem)' + write(fates_log(),fmt=fmth) '------' + write(fates_log(),fmt=fmth) ' Tissue | Initial | Current | Target | On-allometry' + write(fates_log(),fmt=fmtb) ' Leaf |',b0_leaf ,'|',bleaf ,'|',bt_leaf ,'|',fine_leaf + write(fates_log(),fmt=fmtb) ' Fine root |',b0_fnrt ,'|',bfnrt ,'|',bt_fnrt ,'|',fine_fnrt + write(fates_log(),fmt=fmtb) ' Sap wood |',b0_sapw ,'|',bsapw ,'|',bt_sapw ,'|',fine_sapw + write(fates_log(),fmt=fmtb) ' Storage |',b0_store ,'|',bstore ,'|',bt_store ,'|',fine_store + write(fates_log(),fmt=fmtb) ' Structural |',b0_struct ,'|',bstruct ,'|',bt_struct ,'|',fine_struct + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) ' Ancillary information' + write(fates_log(),fmt=fmth) '------' + write(fates_log(),fmt=fmti) ' PFT = ',ipft + write(fates_log(),fmt=fmti) ' leaf_status = ',leaf_status + write(fates_log(),fmt=fmte) ' elongf_leaf = ',elongf_leaf + write(fates_log(),fmt=fmte) ' elongf_fnrt = ',elongf_fnrt + write(fates_log(),fmt=fmte) ' elongf_stem = ',elongf_stem + write(fates_log(),fmt=fmte) ' carbon_balance = ',carbon_balance + write(fates_log(),fmt=fmte) ' calloc_abs_error = ',calloc_abs_error + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) '======' call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( (bdead-bt_dead)> calloc_abs_error) then - grow_dead = .false. - else - grow_dead = .true. end if - return end subroutine TargetAllometryCheck @@ -1100,5 +1827,20 @@ subroutine FastPRTAllometricCarbon(this) end subroutine FastPRTAllometricCarbon + ! ===================================================================================== + + subroutine FastPRTAllometricCarbonSimpler(this) + + implicit none + class(csimpler_allom_prt_vartypes) :: this ! this class + + ! This routine does nothing, because in the carbon only allometric RT model + ! we currently don't have any fast-timestep processes + ! Think of this as a stub. + + + return + end subroutine FastPRTAllometricCarbonSimpler + end module PRTAllometricCarbonMod diff --git a/parteh/PRTGenericMod.F90 b/parteh/PRTGenericMod.F90 index 3dab9563a3..4d686d8ac6 100644 --- a/parteh/PRTGenericMod.F90 +++ b/parteh/PRTGenericMod.F90 @@ -66,6 +66,7 @@ module PRTGenericMod ! These should each have their own module ! ------------------------------------------------------------------------------------- + integer, parameter, public :: prt_csimpler_allom_hyp = 0 integer, parameter, public :: prt_carbon_allom_hyp = 1 integer, parameter, public :: prt_cnp_flex_allom_hyp = 2 ! Still under development diff --git a/parteh/PRTLossFluxesMod.F90 b/parteh/PRTLossFluxesMod.F90 index 13b09b2e37..faba7b50a0 100644 --- a/parteh/PRTLossFluxesMod.F90 +++ b/parteh/PRTLossFluxesMod.F90 @@ -109,8 +109,9 @@ subroutine PRTPhenologyFlush(prt, ipft, organ_id, c_store_transfer_frac) ! those parameters and clauses need to be added !if(organ_id .ne. leaf_organ) then - if(organ_id .ne. leaf_organ .AND. prt_params%woody(ipft) == itrue) then - write(fates_log(),*) 'Deciduous drop and re-flushing only allowed in leaves' + if( organ_id /= leaf_organ .and. organ_id /= fnrt_organ .AND. & + prt_params%woody(ipft) == itrue ) then + write(fates_log(),*) 'Deciduous drop and re-flushing only allowed in leaves and fine roots' write(fates_log(),*) ' leaf_organ: ',leaf_organ write(fates_log(),*) ' organ: ',organ_id write(fates_log(),*) 'Exiting' @@ -424,8 +425,9 @@ subroutine PRTDeciduousTurnover(prt,ipft,organ_id,mass_fraction) ! those parameters and clauses need to be added !if(organ_id .ne. leaf_organ) then - if(organ_id .ne. leaf_organ .AND. prt_params%woody(ipft) == itrue) then - write(fates_log(),*) 'Deciduous drop and re-flushing only allowed in leaves' + if( organ_id /= leaf_organ .and. organ_id /= fnrt_organ .AND. & + prt_params%woody(ipft) == itrue) then + write(fates_log(),*) 'Deciduous drop and re-flushing only allowed in leaves and fine roots' write(fates_log(),*) ' leaf_organ: ',leaf_organ write(fates_log(),*) ' organ: ',organ_id write(fates_log(),*) 'Exiting' diff --git a/parteh/PRTParametersMod.F90 b/parteh/PRTParametersMod.F90 index 04a0f5dda0..36e648c9ab 100644 --- a/parteh/PRTParametersMod.F90 +++ b/parteh/PRTParametersMod.F90 @@ -17,7 +17,10 @@ module PRTParametersMod integer, allocatable :: season_decid(:) ! Is the plant seasonally deciduous (1=yes, 0=no) integer, allocatable :: evergreen(:) ! Is the plant an evergreen (1=yes, 0=no) - + ! Drop fraction for tissues other than leaves (PFT-dependent) + real(r8), allocatable :: phen_fnrt_drop_fraction(:) ! Abscission fraction of fine roots + real(r8), allocatable :: phen_stem_drop_fraction(:) ! Abscission fraction of stems + ! Growth and Turnover Parameters real(r8), allocatable :: senleaf_long_fdrought(:) ! Multiplication factor for leaf longevity of senescent ! leaves during drought( 1.0 indicates no change) @@ -100,6 +103,8 @@ module PRTParametersMod real(r8), allocatable :: c2b(:) ! Carbon to biomass multiplier [kg/kgC] real(r8), allocatable :: wood_density(:) ! wood density g cm^-3 ... + + ! MLO - Shouldn't this be an integer? real(r8), allocatable :: woody(:) ! Does the plant have wood? (1=yes, 0=no) real(r8), allocatable :: crown(:) ! fraction of the height of the plant ! that is occupied by crown diff --git a/parteh/PRTParamsFATESMod.F90 b/parteh/PRTParamsFATESMod.F90 index 208ff848fb..ba0d9cafad 100644 --- a/parteh/PRTParamsFATESMod.F90 +++ b/parteh/PRTParamsFATESMod.F90 @@ -17,7 +17,9 @@ module PRTInitParamsFatesMod use FatesGlobals, only : fates_log use shr_log_mod, only : errMsg => shr_log_errMsg use EDPftvarcon, only : EDPftvarcon_inst - use PRTGenericMod, only : prt_cnp_flex_allom_hyp,prt_carbon_allom_hyp + use PRTGenericMod, only : prt_cnp_flex_allom_hyp + use PRTGenericMod, only : prt_carbon_allom_hyp + use PRTGenericMod, only : prt_csimpler_allom_hyp use FatesAllometryMod , only : h_allom use FatesAllometryMod , only : h2d_allom use FatesAllometryMod , only : bagw_allom @@ -163,6 +165,14 @@ subroutine PRTRegisterPFT(fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_phen_fnrt_drop_fraction' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_phen_stem_drop_fraction' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_fnrt_prof_a' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -411,6 +421,14 @@ subroutine PRTReceivePFT(fates_params) call ArrayNint(tmpreal,prt_params%evergreen) deallocate(tmpreal) + name = 'fates_phen_fnrt_drop_fraction' + call fates_params%RetreiveParameterAllocate(name=name, & + data=prt_params%phen_fnrt_drop_fraction) + + name = 'fates_phen_stem_drop_fraction' + call fates_params%RetreiveParameterAllocate(name=name, & + data=prt_params%phen_stem_drop_fraction) + name = 'fates_leaf_slamax' call fates_params%RetreiveParameterAllocate(name=name, & data=prt_params%slamax) @@ -840,7 +858,8 @@ subroutine FatesReportPFTParams(is_master) logical, intent(in) :: is_master ! Only log if this is the master proc logical, parameter :: debug_report = .false. - character(len=32),parameter :: fmt0 = '(a,100(F12.4,1X))' + character(len=15),parameter :: fmti = '(a,100(I12,1X))' + character(len=17),parameter :: fmt0 = '(a,100(F12.4,1X))' integer :: npft,ipft @@ -856,9 +875,11 @@ subroutine FatesReportPFTParams(is_master) end if write(fates_log(),*) '----------- FATES PARTEH Parameters -----------------' - write(fates_log(),fmt0) 'stress_decid = ',prt_params%stress_decid - write(fates_log(),fmt0) 'season_decid = ',prt_params%season_decid - write(fates_log(),fmt0) 'evergreen = ',prt_params%evergreen + write(fates_log(),fmti) 'stress_decid = ',prt_params%stress_decid + write(fates_log(),fmti) 'season_decid = ',prt_params%season_decid + write(fates_log(),fmti) 'evergreen = ',prt_params%evergreen + write(fates_log(),fmt0) 'phen_fnrt_drop_fraction = ',prt_params%phen_fnrt_drop_fraction + write(fates_log(),fmt0) 'phen_stem_drop_fraction = ',prt_params%phen_stem_drop_fraction write(fates_log(),fmt0) 'wood_density = ',prt_params%wood_density write(fates_log(),fmt0) 'dbh max height = ',prt_params%allom_dbh_maxheight write(fates_log(),fmt0) 'dbh mature = ',prt_params%dbh_repro_threshold @@ -918,7 +939,7 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'turnover_carb_retrans = ',prt_params%turnover_carb_retrans write(fates_log(),fmt0) 'turnover_nitr_retrans = ',prt_params%turnover_nitr_retrans write(fates_log(),fmt0) 'turnover_phos_retrans = ',prt_params%turnover_phos_retrans - write(fates_log(),fmt0) 'organ_id = ',prt_params%organ_id + write(fates_log(),fmti) 'organ_id = ',prt_params%organ_id write(fates_log(),fmt0) 'nitr_store_ratio = ',prt_params%nitr_store_ratio write(fates_log(),fmt0) 'phos_store_ratio = ',prt_params%phos_store_ratio write(fates_log(),*) '-------------------------------------------------' @@ -1014,12 +1035,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_csimpler_allom_hyp,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' @@ -1036,7 +1057,7 @@ subroutine PRTCheckParams(is_master) end if end do - end if + end select pftloop: do ipft = 1,npft @@ -1128,9 +1149,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_csimpler_allom_hyp,prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp) do i = 1,norgans io = prt_params%organ_id(i) @@ -1165,10 +1185,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 +1256,8 @@ 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 +1278,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_csimpler_allom_hyp,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 +1289,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 +1330,7 @@ subroutine PRTCheckParams(is_master) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end select ! Check turnover time-scales