diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index f7d4782f15..4e3ea3a29b 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -689,7 +689,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr,bc_in) !!allocate(copyc%l2fr_ema) ! Note, no need to give a starter value here, ! that will be taken care of in copy() - !!call copyc%l2fr_ema%InitRMean(ema_60day) + !!call copyc%l2fr_ema%InitRSumm(ema_60day) ! Initialize the PARTEH object and point to the ! correct boundary condition fields @@ -1157,7 +1157,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) !!allocate(copyc%l2fr_ema) ! Note, no need to give a starter value here, ! that will be taken care of in copy() - !!call copyc%l2fr_ema%InitRMean(ema_60day) + !!call copyc%l2fr_ema%InitRSumm(ema_60day) ! Initialize the PARTEH object and point to the ! correct boundary condition fields @@ -1172,7 +1172,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! (keep as an example) ! Initialize running means !allocate(copyc%tveg_lpa) - !call copyc%tveg_lpa%InitRMean(ema_lpa,& + !call copyc%tveg_lpa%InitRSumm(ema_lpa,& ! init_value=currentPatch%tveg_lpa%GetMean()) call currentCohort%Copy(copyc) !makes an identical copy... diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 6d9f8cbcb5..3e69a0913a 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -39,6 +39,7 @@ Module EDCohortDynamicsMod use PRTGenericMod , only : num_elements use FatesConstantsMod , only : leaves_off use FatesConstantsMod , only : leaves_shedding + use FatesConstantsMod , only : ihard_season_decid use FatesConstantsMod , only : ihard_stress_decid use FatesConstantsMod , only : isemi_stress_decid use EDParamsMod , only : ED_val_cohort_age_fusion_tol @@ -173,7 +174,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, height, coage, dbh, ! (Keeping as an example) !! allocate(newCohort%tveg_lpa) - !! call newCohort%tveg_lpa%InitRMean(ema_lpa,init_value=patchptr%tveg_lpa%GetMean()) + !! call newCohort%tveg_lpa%InitRSumm(ema_lpa,init_value=patchptr%tveg_lpa%GetMean()) if (hlm_use_planthydro .eq. itrue) then @@ -1389,10 +1390,10 @@ subroutine DamageRecovery(csite,cpatch,ccohort,newly_recovered) !--- Set some logical flags to simplify "if" blocks is_hydecid_dormant = & - any(prt_params%stress_decid(ipft) == [ihard_stress_decid,isemi_stress_decid] ) & + any(prt_params%phen_leaf_habit(ipft) == [ihard_stress_decid,isemi_stress_decid] ) & .and. any(ccohort%status_coh == [leaves_off,leaves_shedding] ) is_sedecid_dormant = & - ( prt_params%season_decid(ipft) == itrue ) & + ( prt_params%phen_leaf_habit(ipft) == ihard_season_decid ) & .and. any(ccohort%status_coh == [leaves_off,leaves_shedding] ) ! If plants are drought deciduous and are losing or lost all leaves, they cannot diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 053b1ca3bc..312b01bb83 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -16,6 +16,7 @@ module EDMortalityFunctionsMod use FatesConstantsMod , only : cstarvation_model_lin use FatesConstantsMod , only : cstarvation_model_exp use FatesConstantsMod , only : nearzero + use FatesConstantsMod , only : ihard_season_decid use FatesConstantsMod , only : ihard_stress_decid use FatesConstantsMod , only : isemi_stress_decid use FatesConstantsMod , only : leaves_off @@ -112,11 +113,9 @@ subroutine mortality_rates( cohort_in,bc_in, btran_ft, mean_temp, & ! the future we could accelerate senescence to avoid mortality. Note that both drought ! deciduous and cold deciduous are considered here to be consistent with the idea that ! plants without leaves cannot die of hydraulic failure. - is_decid_dormant = & ! - ( prt_params%stress_decid(cohort_in%pft) == ihard_stress_decid .or. & ! Drought deciduous - prt_params%stress_decid(cohort_in%pft) == isemi_stress_decid .or. & ! Semi-deciduous - prt_params%season_decid(cohort_in%pft) == itrue ) .and. & ! Cold deciduous - ( cohort_in%status_coh == leaves_off ) ! ! Fully abscised + is_decid_dormant = & ! + any ( prt_params%phen_leaf_habit(cohort_in%pft) == [ihard_season_decid,ihard_stress_decid,isemi_stress_decid]) .and. & ! Deciduous + ( cohort_in%status_coh == leaves_off ) ! ! Fully abscised ! Size Dependent Senescence ! rate (r) and inflection point (ip) define the increase in mortality rate with dbh diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index b1fc9af66d..872b3d9039 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -106,9 +106,9 @@ module EDPatchDynamicsMod use EDParamsMod, only : logging_event_code use EDParamsMod, only : logging_export_frac use EDParamsMod, only : maxpatches_by_landuse - use FatesRunningMeanMod, only : ema_sdlng_mdd - use FatesRunningMeanMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par, ema_sdlng2sap_par - use FatesRunningMeanMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm + use FatesRunningSummMod, only : ema_sdlng_mdd + use FatesRunningSummMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par, ema_sdlng2sap_par + use FatesRunningSummMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm use FatesRadiationMemMod, only : num_swb ! CIME globals @@ -808,7 +808,7 @@ subroutine spawn_patches( currentSite, bc_in) ! (Keeping as an example) ! Allocate running mean functions !allocate(nc%tveg_lpa) - !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=newPatch%tveg_lpa%GetMean()) + !call nc%tveg_lpa%InitRSumm(ema_lpa,init_value=newPatch%tveg_lpa%GetMean()) call nc%ZeroValues() @@ -1701,7 +1701,7 @@ subroutine split_patch(currentSite, currentPatch, new_patch, fraction_to_keep, a ! (Keeping as an example) ! Allocate running mean functions !allocate(nc%tveg_lpa) - !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) + !call nc%tveg_lpa%InitRSumm(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) call nc%ZeroValues() @@ -3183,21 +3183,26 @@ subroutine fuse_2_patches(csite, dp, rp) call endrun(msg=errMsg(sourcefile, __LINE__)) endif - ! Weighted mean of the running means - call rp%tveg24%FuseRMean(dp%tveg24,rp%area*inv_sum_area) - call rp%tveg_lpa%FuseRMean(dp%tveg_lpa,rp%area*inv_sum_area) + ! Weighted mean of the running summaries + call rp%tveg24%FuseRSumm(dp%tveg24,rp%area*inv_sum_area) + call rp%tveg_lpa%FuseRSumm(dp%tveg_lpa,rp%area*inv_sum_area) + + do pft = 1,numpft + call rp%btran24_ft(pft)%p%FuseRSumm(dp%btran24_ft(pft)%p,rp%area*inv_sum_area) + enddo + if ( hlm_regeneration_model == TRS_regeneration ) then - call rp%seedling_layer_par24%FuseRMean(dp%seedling_layer_par24,rp%area*inv_sum_area) - call rp%sdlng_mort_par%FuseRMean(dp%sdlng_mort_par,rp%area*inv_sum_area) - call rp%sdlng2sap_par%FuseRMean(dp%sdlng2sap_par,rp%area*inv_sum_area) + call rp%seedling_layer_par24%FuseRSumm(dp%seedling_layer_par24,rp%area*inv_sum_area) + call rp%sdlng_mort_par%FuseRSumm(dp%sdlng_mort_par,rp%area*inv_sum_area) + call rp%sdlng2sap_par%FuseRSumm(dp%sdlng2sap_par,rp%area*inv_sum_area) do pft = 1,numpft - call rp%sdlng_emerg_smp(pft)%p%FuseRMean(dp%sdlng_emerg_smp(pft)%p,rp%area*inv_sum_area) - call rp%sdlng_mdd(pft)%p%FuseRMean(dp%sdlng_mdd(pft)%p,rp%area*inv_sum_area) + call rp%sdlng_emerg_smp(pft)%p%FuseRSumm(dp%sdlng_emerg_smp(pft)%p,rp%area*inv_sum_area) + call rp%sdlng_mdd(pft)%p%FuseRSumm(dp%sdlng_mdd(pft)%p,rp%area*inv_sum_area) enddo end if - call rp%tveg_longterm%FuseRMean(dp%tveg_longterm,rp%area*inv_sum_area) + call rp%tveg_longterm%FuseRSumm(dp%tveg_longterm,rp%area*inv_sum_area) rp%livegrass = (dp%livegrass*dp%area + rp%livegrass*rp%area) * inv_sum_area rp%ros_front = (dp%ros_front*dp%area + rp%ros_front*rp%area) * inv_sum_area @@ -3811,6 +3816,10 @@ subroutine CopyPatchMeansTimers(dp, rp) call rp%tveg_lpa%CopyFromDonor(dp%tveg_lpa) call rp%tveg_longterm%CopyFromDonor(dp%tveg_longterm) + do ipft = 1,numpft + call rp%btran24_ft(ipft)%p%CopyFromDonor(dp%btran24_ft(ipft)%p) + end do + if ( hlm_regeneration_model == TRS_regeneration ) then call rp%seedling_layer_par24%CopyFromDonor(dp%seedling_layer_par24) call rp%sdlng_mort_par%CopyFromDonor(dp%sdlng_mort_par) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 54f68036cf..5ed7ed12fb 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -67,6 +67,8 @@ module EDPhysiologyMod use EDParamsMod , only : nclmax use EDTypesMod , only : AREA,AREA_INV use FatesConstantsMod , only : leaves_shedding + use FatesConstantsMod , only : ievergreen + use FatesConstantsMod , only : ihard_season_decid use FatesConstantsMod , only : ihard_stress_decid use FatesConstantsMod , only : isemi_stress_decid use EDParamsMod , only : nlevleaf @@ -82,15 +84,16 @@ module EDPhysiologyMod use PRTGenericMod , only : element_list use PRTGenericMod , only : element_pos use EDTypesMod , only : site_fluxdiags_type - use EDTypesMod , only : phen_cstat_nevercold - use EDTypesMod , only : phen_cstat_iscold - use EDTypesMod , only : phen_cstat_notcold + use EDTypesMod , only : phen_estat_evergreen + use EDTypesMod , only : phen_cstat_timeoff + use EDTypesMod , only : phen_cstat_tempoff + use EDTypesMod , only : phen_cstat_tempon + use EDTypesMod , only : phen_cstat_timeon use EDTypesMod , only : phen_dstat_timeoff 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 : phen_dstat_pshed use EDTypesMod , only : init_recruit_trim use shr_log_mod , only : errMsg => shr_log_errMsg use FatesGlobals , only : fates_log @@ -172,8 +175,6 @@ module EDPhysiologyMod integer :: istat ! return status code character(len=255) :: smsg ! Message string for deallocation errors - - integer, parameter :: dleafon_drycheck = 100 ! Drought deciduous leaves max days on check parameter real(r8), parameter :: decid_leaf_long_max = 1.0_r8 ! Maximum leaf lifespan for ! deciduous PFTs [years] @@ -196,10 +197,6 @@ module EDPhysiologyMod ! a residual amount of leaves, which may create ! computational problems. The current threshold ! is the same used in ED-2.2. - - real(r8), parameter :: smp_lwr_bound = -1000000._r8 ! Imposed soil matric potential lower bound for - ! frozen or excessively dry soils, used when - ! computing water stress. ! ============================================================================ contains @@ -477,7 +474,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(1:numpft), bc_in, currentPatch) + call SeedGermination(litt, currentSite%phen_status(1:numpft), bc_in, currentPatch) ! Send fluxes from newly created litter into the litter pools ! This litter flux is from non-disturbance inducing mortality, as well @@ -766,20 +763,21 @@ subroutine trim_canopy( currentSite ) sla_levleaf = min(sla_max,prt_params%slatop(ipft)/nscaler_levleaf) ! Find the realised leaf lifespan, depending on the leaf phenology. - if (prt_params%season_decid(ipft) == itrue) then + select case (prt_params%phen_leaf_habit(ipft)) + case (ihard_season_decid) ! Cold-deciduous costs. Assume time-span to be 1 year to be consistent ! with FATES default pft_leaf_lifespan = decid_leaf_long_max - elseif (any(prt_params%stress_decid(ipft) == [ihard_stress_decid,isemi_stress_decid]) )then + case (ihard_stress_decid,isemi_stress_decid) ! Drought-decidous costs. Assume time-span to be the least between ! 1 year and the life span provided by the parameter file. pft_leaf_lifespan = & min(decid_leaf_long_max,leaf_long) - else !evergreen costs + case (ievergreen) !evergreen costs pft_leaf_lifespan = leaf_long - end if + end select ! Leaf cost at leaf level z (kgC m-2 year-1) accounting for sla profile ! (Convert from SLA in m2g-1 to m2kg-1) @@ -906,11 +904,6 @@ subroutine phenology( currentSite, bc_in ) ! ! !USES: use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - use EDParamsMod, only : ED_val_phen_a, ED_val_phen_b, ED_val_phen_c - use EDParamsMod, only : ED_val_phen_chiltemp - use EDParamsMod, only : ED_val_phen_mindayson - use EDParamsMod, only : ED_val_phen_ncolddayslim - use EDParamsMod, only : ED_val_phen_coldtemp use EDBtranMod, only : check_layer_water ! ! !ARGUMENTS: @@ -920,25 +913,16 @@ subroutine phenology( currentSite, bc_in ) ! ! !LOCAL VARIABLES: - type(fates_patch_type),pointer :: cpatch integer :: model_day_int ! integer model day 1 - inf 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 :: ipft ! plant functional type index - integer :: j ! Soil layer index - real(r8) :: mean_10day_liqvol ! mean soil liquid volume over last 10 days [m3/m3] + real(r8) :: mean_10day_btran ! mean transpitation wetness factor (btran) over last 10 days. real(r8) :: mean_10day_smp ! mean soil matric potential over last 10 days [mm] - 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, - 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. - integer :: nlevroot ! Number of rooting levels to consider + integer :: gddend ! end of counting period for growing degree days. real(r8) :: temp_in_C ! daily averaged temperature in celsius real(r8) :: temp_wgt ! canopy area weighting factor for daily average ! vegetation temperature calculation @@ -949,6 +933,27 @@ subroutine phenology( currentSite, bc_in ) ! lifespan and the maximum lifespan of drought ! deciduous (see parameter decid_leaf_long_max ! at the beginning of this file). + + ! The three parameters below are for the growing degree days + ! threshold function, which is modulated by the number of + ! chilling days (NCD): GDD_Thresh = a + b * exp (c * NCD) + real(r8) :: phen_gddthresh_a ! - a parameter (intercept) + real(r8) :: phen_gddthresh_b ! - b parameter (multiplier) + real(r8) :: phen_gddthresh_c ! - c parameter (exponent) + + real(r8) :: phen_gddtemp ! (Lower) Temperature threshold for accumulating + ! growing degree days (used for leaf flushing) + real(r8) :: phen_chilltemp ! (Upper) Temperature threshold for accumulating + ! number of chilling days (used for leaf flushing) + real(r8) :: phen_coldtemp ! (Upper) daily average temperature threshold below which the + ! day is considered cold and will count towards leaf + ! abscission + real(r8) :: phen_mindayson ! Minimum number of days that plants must remain with leaves + ! flushed before they can abscise leaves again. + real(r8) :: phen_mindaysoff ! Minimum number of days that plants must remain leafless + ! before theu can flush leaves again. + real(r8) :: phen_ncolddayslim ! Minimum number of recent days below the temperature threshold + ! that initiates temperature-driven leaf abscission. real(r8) :: phen_drought_threshold ! For drought hard-deciduous, this is the threshold ! below which plants will abscise leaves, and ! above which plants will flush leaves. For semi- @@ -966,8 +971,6 @@ subroutine phenology( currentSite, bc_in ) ! the values are soil matric potential [mm]. ! Ignored for hard-deciduous and evergreen ! plants. - real(r8) :: phen_doff_time ! Minimum number of days that plants must remain - ! leafless before flushing leaves again. ! Logical tests to make code more readable logical :: smoist_below_threshold ! Is soil moisture below threshold? @@ -978,305 +981,152 @@ subroutine phenology( currentSite, bc_in ) 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? + logical :: growing_season ! Is this the growing season for cold-deciduous? + logical :: gdd_above_threshold ! Has the growing degree day threshold been exceeded? + logical :: colddays_above_threshold ! Has there been many cold days recently? - - ! This is the integer model day. The first day of the simulation is 1, and it - ! continues monotonically, indefinitely - ! Advance it. (this should be a global, no reason - ! for site level, but we don't have global scalars in the - ! restart file) - currentSite%phen_model_date = currentSite%phen_model_date + 1 + ! This is the elapsed number of days since the very beginning of the simulation time model_day_int = currentSite%phen_model_date - ! 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 - temp_in_C = 0._r8 - temp_wgt = 0._r8 - cpatch => CurrentSite%oldest_patch - do while(associated(cpatch)) - temp_in_C = temp_in_C + cpatch%tveg24%GetMean()*cpatch%total_canopy_area - temp_wgt = temp_wgt + cpatch%total_canopy_area - cpatch => cpatch%younger - end do - if(temp_wgt>nearzero)then - temp_in_C = temp_in_C/temp_wgt - tfrz - else - ! If there is no canopy area, we use the veg temperature - ! of the first patch, which is the forcing air temperature - ! as defined in CLM/ELM. The forcing air temperature - ! should be the same among all patches. (Although - ! it is unlikely there are more than 1 in this scenario) - temp_in_C = CurrentSite%oldest_patch%tveg24%GetMean() - tfrz - end if + ! Retrieve average canopy temperature of the current day + temp_in_C = currentSite%vegtemp_memory(1) + !-----------------Cold Phenology--------------------! !Zero growing degree and chilling day counters if (currentSite%lat > 0)then - ncdstart = 270 !Northern Hemisphere begining November - gddstart = 1 !Northern Hemisphere begining January + ncdstart = 270 ! Northern Hemisphere begining November + gddstart = 1 ! Northern Hemisphere begining January + gddend = 180 ! Northern Hemisphere ending June else - ncdstart = 120 !Southern Hemisphere beginning May - gddstart = 181 !Northern Hemisphere begining July + ncdstart = 120 ! Southern Hemisphere beginning May + gddstart = 181 ! Southern Hemisphere begining July + gddend = 365 ! Southern Hemisphere ending December endif + ! Define growing season based on start and end dates. Check for cases in which the growing season + ! includes 31 December and 1 January. + if (gddend >= gddstart) then + growing_season = hlm_day_of_year >= gddstart .and. hlm_day_of_year <= gddend + else + growing_season = hlm_day_of_year >= gddstart .or. hlm_day_of_year <= gddend + end if + + ! Count the number of chilling days over a seasonal window. ! For comparing against GDD, we start calculating chilling ! in the late autumn. ! This value is used to determine the GDD exceedance threshold if (hlm_day_of_year == ncdstart)then - currentSite%nchilldays = 0 - endif - - !Accumulate growing/chilling days after start of counting period - if (temp_in_C < ED_val_phen_chiltemp)then - currentSite%nchilldays = currentSite%nchilldays + 1 + currentSite%nchilldays(1:numpft) = 0 endif - - !GDD accumulation function, which also depends on chilling days. - ! -68 + 638 * (-0.001 * ncd) - gdd_threshold = ED_val_phen_a + ED_val_phen_b*exp(ED_val_phen_c*real(currentSite%nchilldays,r8)) - - !Accumulate temperature of last 10 days. - currentSite%vegtemp_memory(2:num_vegtemp_mem) = currentSite%vegtemp_memory(1:num_vegtemp_mem-1) - currentSite%vegtemp_memory(1) = temp_in_C - - !count number of days for leaves off - ncolddays = 0 - do i_tmem = 1,num_vegtemp_mem - if (currentSite%vegtemp_memory(i_tmem) < ED_val_phen_coldtemp)then - ncolddays = ncolddays + 1 - endif - enddo - - ! Here is where we do the GDD accumulation calculation ! ! reset GDD on set dates if (hlm_day_of_year == gddstart)then - currentSite%grow_deg_days = 0._r8 - endif - ! - ! accumulate the GDD using daily mean temperatures - ! Don't accumulate GDD during the growing season (that wouldn't make sense) - if (temp_in_C .gt. 0._r8 .and. currentSite%cstatus == phen_cstat_iscold) then - currentSite%grow_deg_days = currentSite%grow_deg_days + temp_in_C + currentSite%grow_deg_days(1:numpft) = 0._r8 endif - !this logic is to prevent GDD accumulating after the leaves have fallen and before the - ! beginnning of the accumulation period, to prevend erroneous autumn leaf flushing. - if(model_day_int> ndays_per_year)then !only do this after the first year to prevent odd behaviour - - if(currentSite%lat .gt. 0.0_r8)then !Northern Hemisphere - ! In the north, don't accumulate when we are past the leaf fall date. - ! Accumulation starts on day 1 of year in NH. - ! The 180 is to prevent going into an 'always off' state after initialization - if( model_day_int .gt. currentSite%cleafoffdate.and.hlm_day_of_year.gt.180)then ! - currentSite%grow_deg_days = 0._r8 - endif - else !Southern Hemisphere - ! In the South, don't accumulate after the leaf off date, and before the start of - ! the accumulation phase (day 181). - if(model_day_int .gt. currentSite%cleafoffdate.and.hlm_day_of_year.lt.gddstart) then! - currentSite%grow_deg_days = 0._r8 - endif - endif - endif !year1 - - ! Calculate the number of days since the leaves last came on - ! and off. If this is the beginning of the simulation, that day might - ! not had occured yet, so set it to last year to get things rolling - - if (model_day_int < currentSite%cleafoffdate) then - currentSite%cndaysleafoff = model_day_int - (currentSite%cleafoffdate - ndays_per_year) - else - currentSite%cndaysleafoff = model_day_int - currentSite%cleafoffdate - end if - - if (model_day_int < currentSite%cleafondate) then - currentSite%cndaysleafon = model_day_int - (currentSite%cleafondate - ndays_per_year) - else - currentSite%cndaysleafon = model_day_int - currentSite%cleafondate - end if - - - - !LEAF ON: COLD DECIDUOUS. Needs to - !1) have exceeded the growing degree day threshold - !2) The leaves should not be on already - !3) There should have been at least one chilling day in the counting period. - ! this prevents tropical or warm climate plants that are "cold-deciduous" - ! from ever re-flushing after they have reached their maximum age (thus - ! preventing them from competing - - if ( any(currentSite%cstatus == [phen_cstat_iscold,phen_cstat_nevercold]) .and. & - (currentSite%grow_deg_days > gdd_threshold) .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 - 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 + !---~--- + ! Loop through every PFT to assign the elongation factor. + ! Add PFT loop to account for different parameters and PFT-specific rooting depth profiles. + !---~--- + pft_elong_loop: do ipft=1,numpft + ! Copy values to a local variable to make code more legible. + phen_gddthresh_a = prt_params%phen_gddthresh_a (ipft) + phen_gddthresh_b = prt_params%phen_gddthresh_b (ipft) + phen_gddthresh_c = prt_params%phen_gddthresh_c (ipft) + phen_gddtemp = prt_params%phen_gddtemp (ipft) + phen_chilltemp = prt_params%phen_chilltemp (ipft) + phen_coldtemp = prt_params%phen_coldtemp (ipft) + phen_mindayson = prt_params%phen_mindayson (ipft) + phen_mindaysoff = prt_params%phen_mindaysoff (ipft) + phen_ncolddayslim = prt_params%phen_ncolddayslim (ipft) + phen_drought_threshold = prt_params%phen_drought_threshold(ipft) + phen_moist_threshold = prt_params%phen_moist_threshold (ipft) + ! PFT leaf lifespan in days. This is the shortest between the leaf longevity + ! (defined as a PFT parameter) and the maximum canopy leaf life span allowed + ! for drought deciduous (local parameter). The sum term accounts for the + ! total leaf life span of this cohort. + ! Note we only use canopy leaf lifespan here and assume that understory cohorts + ! would behave the same as canopy cohorts with regards to phenology. + ndays_pft_leaf_lifespan = & + nint(ndays_per_year*min(decid_leaf_long_max,sum(prt_params%leaf_long(ipft,:)))) - !LEAF OFF: COLD THRESHOLD - !Needs to: - !1) have exceeded the number of cold days threshold - !2) have exceeded the minimum leafon time. - !3) The leaves should not be off already - !4) The day of simulation should be larger than the counting period. - - - if ( (currentSite%cstatus == phen_cstat_notcold) .and. & - (model_day_int > num_vegtemp_mem) .and. & - (ncolddays > ED_val_phen_ncolddayslim) .and. & - (currentSite%cndaysleafon > ED_val_phen_mindayson) )then - - currentSite%grow_deg_days = 0._r8 ! The equations for Botta et al - ! are for calculations of - ! first flush, but if we dont - ! 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%cndaysleafoff = 0 - if ( debug ) write(fates_log(),*) 'leaves off' - endif + ! Accumulate growing/chilling days after start of counting period + if (temp_in_C < phen_chilltemp) then + currentSite%nchilldays(ipft) = currentSite%nchilldays(ipft) + 1 + end if - ! LEAF OFF: COLD LIFESPAN THRESHOLD - ! NOTE: Some areas of the planet will never generate a cold day - ! and thus %nchilldays will never go from zero to 1. The following logic - ! when coupled with this fact will essentially prevent cold-deciduous - ! plants from re-emerging in areas without at least some cold days - - if( (currentSite%cstatus == phen_cstat_notcold) .and. & - (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 - currentSite%cndaysleafoff = 0 - - if ( debug ) write(fates_log(),*) 'leaves off' - endif + ! GDD accumulation function, which also depends on chilling days. + ! -68 + 638 * (-0.001 * ncd) + gdd_threshold = phen_gddthresh_a + & + phen_gddthresh_b * exp(phen_gddthresh_c*real(currentSite%nchilldays(ipft),r8)) + + ! Count number of days for leaves off + ncolddays = count(currentSite%vegtemp_memory(1:num_vegtemp_mem) < phen_coldtemp) + + ! Accumulate the GDD using daily mean temperatures. For now, we only accumulate GDD + ! when plants are leafless and when this is the growing season. This may change in + ! future pull requests. + if ( temp_in_C > phen_gddtemp .and. & ! Warm day + any( currentSite%phen_status(ipft) == [phen_cstat_tempoff,phen_cstat_timeoff] ) .and. & ! PFT is dormant + growing_season) then ! Growing season + currentSite%grow_deg_days(ipft) = currentSite%grow_deg_days(ipft) + temp_in_C - phen_gddtemp + end if - ! Loop through every PFT to assign the elongation factor. - ! Add PFT look to account for different PFT rooting depth profiles. - pft_elong_loop: do ipft=1,numpft - ! Copy values to a local variable to make code more legible. - phen_drought_threshold = prt_params%phen_drought_threshold(ipft) - phen_moist_threshold = prt_params%phen_moist_threshold (ipft) - phen_doff_time = prt_params%phen_doff_time (ipft) + !---~--- + ! We update values relevant for phenology controls, but we do not update leaf + ! phenology until the memory variables are fully populated. + !---~--- + case_past_spinup: select case(prt_params%phen_leaf_habit(ipft)) + case (ihard_season_decid) + if (model_day_int <= num_vegtemp_mem) cycle pft_elong_loop + case (ihard_stress_decid,isemi_stress_decid) + if (model_day_int <= numWaterMem ) cycle pft_elong_loop + end select case_past_spinup + !---~--- - ! 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,ipft) = currentSite%liqvol_memory(i_wmem-1,ipft) - currentSite%smp_memory (i_wmem,ipft) = currentSite%smp_memory (i_wmem-1,ipft) - end do - ! Find the rooting depth distribution for PFT - call set_root_fraction( currentSite%rootfrac_scr, ipft, 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 - ! 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,ipft) = sum( bc_in%h2o_liqvol_sl (2:nlevroot) * & - currentSite%rootfrac_scr(2:nlevroot) ) / & - rootfrac_notop - currentSite%smp_memory (1,ipft) = 0._r8 - do j = 2,nlevroot - if(check_layer_water(bc_in%h2o_liqvol_sl(j),bc_in%tempk_sl(j)) ) then - currentSite%smp_memory (1,ipft) = currentSite%smp_memory (1,ipft) + & - bc_in%smp_sl (j) * & - currentSite%rootfrac_scr(j) / & - rootfrac_notop - else - ! Nominal extreme suction for frozen or unreasonably dry soil - currentSite%smp_memory (1,ipft) = currentSite%smp_memory (1,ipft) + & - smp_lwr_bound * & - currentSite%rootfrac_scr(j) / & - rootfrac_notop - end if - end do - - ! 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,ipft)) / & - real(numWaterMem,r8) - mean_10day_smp = sum(currentSite%smp_memory (1:numWaterMem,ipft)) / & - real(numWaterMem,r8) + !---~--- + ! Calculate the mean transpiration wetness factor (btran) and matric potential + ! (mm) over the last 10 days + !---~--- + mean_10day_btran = sum(currentSite%btran_memory(1:numWaterMem,ipft)) / & + real(numWaterMem,r8) + mean_10day_smp = sum(currentSite%smp_memory (1:numWaterMem,ipft)) / & + real(numWaterMem,r8) ! Compare the moisture with the threshold. if ( phen_drought_threshold >= 0. ) then - ! Liquid volume in reference layer (m3/m3) - smoist_below_threshold = mean_10day_liqvol < phen_drought_threshold + ! Transpiration wetness factor (btran) + smoist_below_threshold = mean_10day_btran < phen_drought_threshold else ! Soil matric potential in reference layer (mm) smoist_below_threshold = mean_10day_smp < phen_drought_threshold end if - ! 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(ipft)) then - currentSite%dndaysleafoff(ipft) = model_day_int - (currentSite%dleafoffdate(ipft)-ndays_per_year) - else - currentSite%dndaysleafoff(ipft) = model_day_int - currentSite%dleafoffdate(ipft) - end if - if (model_day_int < currentSite%dleafondate(ipft)) then - currentSite%dndaysleafon(ipft) = model_day_int - (currentSite%dleafondate(ipft)-ndays_per_year) - else - currentSite%dndaysleafon(ipft) = model_day_int - currentSite%dleafondate(ipft) - end if - ! Elongation factor from the previous step. elongf_prev = currentSite%elong_factor(ipft) - ! PFT leaf lifespan in days. This is the shortest between the leaf longevity - ! (defined as a PFT parameter) and the maximum canopy leaf life span allowed - ! for drought deciduous (local parameter). The sum term accounts for the - ! total leaf life span of this cohort. - ! Note we only use canopy leaf lifespan here and assume that understory cohorts - ! would behave the same as canopy cohorts with regards to phenology. - ndays_pft_leaf_lifespan = & - nint(ndays_per_year*min(decid_leaf_long_max,sum(prt_params%leaf_long(ipft,:)))) - - !---~--- - ! Find elongation factors by comparing the moisture with the thresholds. For each - ! tissue --- leaves, fine roots, and stems (sapwood+heartwood) --- elongation factor + ! Find elongation factors by comparing the temperature and moisture memories with thresholds. For each + ! tissue---leaves, fine roots, and stems (sapwood+heartwood)---, elongation factor ! is the maximum fraction of biomass (relative to maximum biomass given allometry) ! that can be allocated to each tissue due to phenology. In this select case, we ! define the elongation factor based on the PFT-specific phenology strategy of this @@ -1289,17 +1139,174 @@ subroutine phenology( currentSite, bc_in ) ! the leaf biomass will be capped at 40% of the biomass the cohort would have if ! it were in well-watered conditions. !---~--- - case_drought_phen: select case (prt_params%stress_decid(ipft)) + case_update_phen: select case(prt_params%phen_leaf_habit(ipft)) + case (ievergreen) + !---~--- + ! EVERGREEN + ! + ! Trivial case, nothing happens, enforce that leaves must remain flushed. + ! This case could be suppressed from here too. + !---~--- + currentSite%phen_status (ipft) = phen_estat_evergreen + currentSite%elong_factor(ipft) = 1.0_r8 + case (ihard_season_decid) + !---~--- + ! COLD (SEASON) DECIDUOUS + ! + ! The decision on whether to abscise (shed) or flush leaves is in principle + ! defined by the number of cold days and the growing degree days. However, we must + ! also account the time since last abscission or flushing event, to avoid excessive + ! "flickering" of the leaf elongation factor, and to ensure plants do not remain + ! leafless if temperatures fail to exceed threshold or if temperature never drops + ! below the threshold. + !---~--- + + + + !---~--- + ! Save some conditions in logical variables to simplify code below + !---~--- + ! Leaves have been "on" for longer than the minimum number of days. + exceed_min_on_period = & + any( currentSite%phen_status(ipft) == [phen_cstat_tempon,phen_cstat_timeon] ) .and. & + ( currentSite%ndaysleafon(ipft) > phen_mindayson ) + ! Leaves have been "off" for longer than the minimum number of days. + exceed_min_off_period = & + any( currentSite%phen_status(ipft) == [phen_cstat_tempoff,phen_cstat_timeoff] ) .and. & + ( currentSite%ndaysleafoff(ipft) > phen_mindaysoff ) + ! Leaves have been "on" for longer than the leaf lifetime. Note that this is slightly different + ! than the original code. This now uses the same logic as the drought deciduous, leaves should not + ! last longer than their nominal leaf lifespan (unless the nominal leaf lifespan is more than one + ! year. + prolonged_on_period = & + any( currentSite%phen_status(ipft) == [phen_cstat_tempon,phen_cstat_timeon] ) .and. & + ( currentSite%ndaysleafoff(ipft) > ndays_pft_leaf_lifespan) + ! Last flushing was a very long time ago. For now we only force flushing when the last + ! abscission event was thermally driven, so it remains consistent with the requirement of + ! having cold nights for cold-deciduous plants to exist. This may be revised. + last_flush_long_ago = & + ( currentSite%phen_status(ipft) == phen_cstat_tempoff ) .and. & + ( currentSite%ndaysleafoff(ipft) > phen_mindaysoff ) .and. & + ( currentSite%ndaysleafon(ipft) >= ndays_per_year-dd_offon_toler ) .and. & + ( currentSite%ndaysleafon(ipft) <= ndays_per_year+dd_offon_toler ) + ! Check whether GDD has exceeded the threshold for flushing or not. Currently we + ! also check that there has been at least one chilling day in the counting period. + ! This prevents tropical or warm climate plants that are "cold-deciduous" from ever + ! re-flushing after they have reached their maximum age (thus preventing them from + ! competing). + ! MLO - Is the latter a desirable feature? I wonder if this blocks seasonal + ! deciduous plants in areas where they exist but winters are mild, like + ! the temperate west coast of most continents. + gdd_above_threshold = & + ( currentSite%grow_deg_days(ipft) > gdd_threshold ) .and. & + ( currentSite%nchilldays(ipft) >= 1 ) + ! Check whether there have been multiple cold days recently (except when it is + ! very early in the simulation and the memory count has not been fully populated). + colddays_above_threshold = ncolddays > phen_ncolddayslim + !---~--- + + + + !---~--- + ! Revision of the conditions, added an if/elseif/else structure to ensure only + ! up to one change occurs at any given time. + !---~--- + cold_temp_ifelse: if ( exceed_min_off_period .and. gdd_above_threshold ) then + !---~--- + ! LEAF ON: COLD DECIDUOUS, THERMAL THRESHOLD EXCEEDED + ! + ! Here we check that the following conditions were met, allowing + ! leaves to flush: + ! 1) PFT has exceeded the growing degree day threshold + ! 2) PFT is in dormant state (leafless) + ! 3) There was at least one chilling day in the counting period. This + ! last check prevents tropical or warm climate plants that are "cold-deciduous" + ! from ever re-flushing after they have reached their maximum age (thus preventing + ! them from competing). + ! MLO - Revise the third case? I imagine this may prevent deciduous trees to + ! establish in mild temperature climates, where they may still occur + ! (e.g., coastal California) + !---~--- + currentSite%phen_status (ipft) = phen_cstat_tempon ! Set to not-cold status (leaves can come on) + currentSite%leafondate (ipft) = model_day_int ! Record flushing date + currentSite%ndaysleafon (ipft) = 0 ! Reset time since flushing + currentSite%grow_deg_days(ipft) = 0._r8 ! zero GDD for the rest of the year until counting season begins. + if ( debug ) write(fates_log(),*) 'leaves on' + elseif ( last_flush_long_ago ) then + !---~--- + ! LEAF ON: COLD DECIDUOUS, TIME OFF EXCEEDED + ! + ! This did not exist in previous versions, but it mimics what is done in drought deciduous. + ! If the plants have been dormant for a very long time and the conditions never allowed for + ! them to start flushing leaves, we force a bud burst, otherwise they can remain dormant for + ! ever. This will likley lead to high mortality, but this is fine, because this PFT's thermal + ! requirements are not suitable for this location. + !---~--- + currentSite%phen_status (ipft) = phen_cstat_timeon ! Set to not-cold status (leaves can come on) + currentSite%leafondate (ipft) = model_day_int ! Record flushing date + currentSite%ndaysleafon (ipft) = 0 ! Reset time since flushing + currentSite%grow_deg_days(ipft) = 0._r8 ! zero GDD for the rest of the year until counting season begins. + if ( debug ) write(fates_log(),*) 'leaves on' + + elseif ( prolonged_on_period ) then + !---~--- + ! LEAF OFF: COLD DECIDUOUS, LIFESPAN THRESHOLD + ! + ! Abscise leaves once they exceed their nominal leaf lifespan or 1 year, + ! whichever happens first. This will prevent deciduous plants to remain permanently green. + !---~--- + currentSite%grow_deg_days(ipft) = 0._r8 ! Reset GDD + currentSite%phen_status (ipft) = phen_cstat_timeoff ! Alter status to imply that this site + ! is never really cold enough for + ! this cold deciduous PFT + currentSite%leafoffdate (ipft) = model_day_int ! Record abscission date + currentSite%ndaysleafoff (ipft) = 0 ! Reset time since abscission + + if ( debug ) write(fates_log(),*) 'leaves off' + + elseif ( exceed_min_on_period .and. colddays_above_threshold ) then + !---~--- + ! LEAF OFF: COLD DECIDUOUS, THERMAL THRESHOLD EXCEEDED + ! + ! Here we check that the following conditions are met, allowing + ! leaf abscission: + ! 1) The PFT has leaves on + ! 2) PFT has exceeded the number of cold days threshold + ! 3) PFT had leaves for more time than the minimum time required. + !---~--- + currentSite%grow_deg_days(ipft) = 0._r8 ! Reset GDD + currentSite%phen_status (ipft) = phen_cstat_tempoff ! Alter status of site to 'leaves off' + currentSite%leafoffdate (ipft) = model_day_int ! Record abscission date + currentSite%ndaysleafoff (ipft) = 0 ! Reset time since abscission + + if ( debug ) write(fates_log(),*) 'leaves off' + end if cold_temp_ifelse + !---~--- + + + + !---~--- + ! Assign elongation factors for cold-drought deciduous PFTs, which will be used + ! to define the cohort status. These are "hard" cold-deciduous, so the value should + ! be either zero or one. + !---~--- + select case (currentSite%phen_status(ipft)) + case (phen_cstat_tempoff,phen_cstat_timeoff) + currentSite%elong_factor(ipft) = 0.0_r8 + case (phen_cstat_tempon,phen_cstat_timeon) + currentSite%elong_factor(ipft) = 1.0_r8 + end select + !---~--- + case (ihard_stress_decid) !---~--- + ! HYDRO (STRESS) DECIDUOUS, "HARD" + ! ! Default ("hard") drought deciduous phenology. The decision on whether to ! abscise (shed) or flush leaves is in principle defined by the soil moisture ! in the rooting zone. However, we must also account the time since last ! abscission or flushing event, to avoid excessive "flickering" of the leaf ! elongation factor if soil moisture is right at the threshold. - ! - ! (MLO thought: maybe we should define moisture equivalents of GDD and chilling - ! days to simplify the cases a bit...) !---~--- @@ -1308,92 +1315,88 @@ subroutine phenology( currentSite, bc_in ) !---~--- ! Leaves have been "on" for longer than the minimum number of days. exceed_min_on_period = & - any( currentSite%dstatus(ipft) == [phen_dstat_timeon,phen_dstat_moiston] ) .and. & - (currentSite%dndaysleafon(ipft) > dleafon_drycheck) + any( currentSite%phen_status(ipft) == [phen_dstat_timeon,phen_dstat_moiston] ) .and. & + ( currentSite%ndaysleafon(ipft) > phen_mindayson ) ! Leaves have been "off" for longer than the minimum number of days. exceed_min_off_period = & - ( currentSite%dstatus(ipft) == phen_dstat_timeoff ) .and. & - ( currentSite%dndaysleafoff(ipft) > min_daysoff_dforcedflush ) + ( currentSite%phen_status(ipft) == phen_dstat_timeoff ) .and. & + ( currentSite%ndaysleafoff(ipft) > min_daysoff_dforcedflush ) ! Leaves have been "on" for longer than the leaf lifetime. prolonged_on_period = & - any( currentSite%dstatus(ipft) == [phen_dstat_timeon,phen_dstat_moiston] ) .and. & - ( currentSite%dndaysleafon(ipft) > ndays_pft_leaf_lifespan ) + any( currentSite%phen_status(ipft) == [phen_dstat_timeon,phen_dstat_moiston] ) .and. & + ( currentSite%ndaysleafon(ipft) > ndays_pft_leaf_lifespan ) ! Leaves have been "off" for a sufficiently long time and the last flushing ! was about one year ago (+/- tolerance). prolonged_off_period = & - any( currentSite%dstatus(ipft) == [phen_dstat_timeoff,phen_dstat_moistoff] ) .and. & - ( currentSite%dndaysleafoff(ipft) > phen_doff_time ) .and. & - ( currentSite%dndaysleafon(ipft) >= ndays_per_year-dd_offon_toler ) .and. & - ( currentSite%dndaysleafon(ipft) <= ndays_per_year+dd_offon_toler ) + any( currentSite%phen_status(ipft) == [phen_dstat_timeoff,phen_dstat_moistoff] ) .and. & + ( currentSite%ndaysleafoff(ipft) > phen_mindaysoff ) .and. & + ( currentSite%ndaysleafon(ipft) >= ndays_per_year-dd_offon_toler ) .and. & + ( currentSite%ndaysleafon(ipft) <= ndays_per_year+dd_offon_toler ) ! Last flushing was a very long time ago. last_flush_long_ago = & - ( currentSite%dstatus(ipft) == phen_dstat_moistoff ) .and. & - ( currentSite%dndaysleafon(ipft) > ndays_per_year+dd_offon_toler ) + ( currentSite%phen_status(ipft) == phen_dstat_moistoff ) .and. & + ( currentSite%ndaysleafon(ipft) > ndays_per_year+dd_offon_toler ) !---~--- !---~--- - ! Revision of the conditions, added an if/elseif/else structure to ensure only - ! up to one change occurs at any given time. Also, prevent changes until the - ! soil moisture memory is populated (the outer if check). + ! Revision of the conditions, added an if/elseif/else structure to ensure only + ! up to one change occurs at any given time. !---~--- - past_spinup_ifelse: if (model_day_int > numWaterMem) then - drought_smoist_ifelse: if ( 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(ipft) = phen_dstat_moiston ! set status to leaf-on - currentSite%dleafondate(ipft) = model_day_int ! save the model day we start flushing - currentSite%dndaysleafon(ipft) = 0 - currentSite%elong_factor(ipft) = 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(ipft) = phen_dstat_timeon ! force budburst! - currentSite%dleafondate(ipft) = model_day_int ! record leaf on date - currentSite%dndaysleafon(ipft) = 0 - currentSite%elong_factor(ipft) = 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(ipft) = phen_dstat_timeon ! force budburst! - currentSite%dleafondate(ipft) = model_day_int ! record leaf on date - currentSite%dndaysleafon(ipft) = 0 - currentSite%elong_factor(ipft) = 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(ipft) = phen_dstat_timeoff !alter status of site to 'leaves off' - currentSite%dleafoffdate(ipft) = model_day_int !record leaf on date - currentSite%dndaysleafoff(ipft) = 0 - currentSite%elong_factor(ipft) = 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(ipft) = phen_dstat_moistoff ! alter status of site to 'leaves off' - currentSite%dleafoffdate(ipft) = model_day_int ! record leaf on date - currentSite%dndaysleafoff(ipft) = 0 - currentSite%elong_factor(ipft) = 0. - end if drought_smoist_ifelse - end if past_spinup_ifelse + drought_smoist_ifelse: if ( 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%phen_status(ipft) = phen_dstat_moiston ! set status to leaf-on + currentSite%leafondate(ipft) = model_day_int ! save the model day we start flushing + currentSite%ndaysleafon(ipft) = 0 + currentSite%elong_factor(ipft) = 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%phen_status(ipft) = phen_dstat_timeon ! force budburst! + currentSite%leafondate(ipft) = model_day_int ! record leaf on date + currentSite%ndaysleafon(ipft) = 0 + currentSite%elong_factor(ipft) = 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%phen_status(ipft) = phen_dstat_timeon ! force budburst! + currentSite%leafondate(ipft) = model_day_int ! record leaf on date + currentSite%ndaysleafon(ipft) = 0 + currentSite%elong_factor(ipft) = 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%phen_status(ipft) = phen_dstat_timeoff !alter status of site to 'leaves off' + currentSite%leafoffdate(ipft) = model_day_int !record leaf on date + currentSite%ndaysleafoff(ipft) = 0 + currentSite%elong_factor(ipft) = 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%phen_status(ipft) = phen_dstat_moistoff ! alter status of site to 'leaves off' + currentSite%leafoffdate(ipft) = model_day_int ! record leaf on date + currentSite%ndaysleafoff(ipft) = 0 + currentSite%elong_factor(ipft) = 0. + end if drought_smoist_ifelse !---~--- @@ -1412,7 +1415,7 @@ subroutine phenology( currentSite, bc_in ) !---~--- if (phen_drought_threshold >= 0.) then elongf_1st = elongf_min + (1.0_r8 - elongf_min ) * & - ( mean_10day_liqvol - phen_drought_threshold ) / & + ( mean_10day_btran - phen_drought_threshold ) / & ( phen_moist_threshold - phen_drought_threshold ) else elongf_1st = elongf_min + (1.0_r8 - elongf_min ) * & @@ -1429,47 +1432,43 @@ subroutine phenology( currentSite, bc_in ) !---~--- ! Leaves have been flushing for a short period of time. recent_flush = elongf_prev >= elongf_min .and. & - ( currentSite%dndaysleafon(ipft) <= dleafon_drycheck ) + ( currentSite%ndaysleafon(ipft) <= phen_mindayson ) ! Leaves have been abscissing for a short period of time. recent_abscission = elongf_prev < elongf_min .and. & - ( currentSite%dndaysleafoff(ipft) <= min_daysoff_dforcedflush ) + ( currentSite%ndaysleafoff(ipft) <= min_daysoff_dforcedflush ) ! Leaves have been flushing for longer than their time span. prolonged_on_period = all( [elongf_prev,elongf_1st] >= elongf_min ) .and. & - ( currentSite%dndaysleafon(ipft) > ndays_pft_leaf_lifespan ) + ( currentSite%ndaysleafon(ipft) > ndays_pft_leaf_lifespan ) ! It's been a long time since the plants had flushed their leaves. last_flush_long_ago = all( [elongf_prev,elongf_1st] < elongf_min ) .and. & - ( currentSite%dndaysleafon(ipft) > ndays_per_year+dd_offon_toler ) + ( currentSite%ndaysleafon(ipft) > ndays_per_year+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(ipft) = elongf_prev - - elseif ( prolonged_on_period ) then + drought_gradual_ifelse: if ( prolonged_on_period ) then ! Leaves have been on for too long and exceeded leaf lifespan. Force abscission currentSite%elong_factor(ipft) = 0.0_r8 ! Force full budburst - currentSite%dstatus(ipft) = phen_dstat_timeoff ! Flag that this has been forced - currentSite%dleafoffdate(ipft) = model_day_int ! Record leaf off date - currentSite%dndaysleafoff(ipft) = 0 ! Reset clock + currentSite%phen_status(ipft) = phen_dstat_timeoff ! Flag that this has been forced + currentSite%leafoffdate(ipft) = model_day_int ! Record leaf off date + currentSite%ndaysleafoff(ipft) = 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(ipft) = elongf_min ! Force minimum budburst - currentSite%dstatus(ipft) = phen_dstat_timeon ! Flag that this has been forced - currentSite%dleafondate(ipft) = model_day_int ! Record leaf on date - currentSite%dndaysleafon(ipft) = 0 ! Reset clock + currentSite%phen_status(ipft) = phen_dstat_timeon ! Flag that this has been forced + currentSite%leafondate(ipft) = model_day_int ! Record leaf on date + currentSite%ndaysleafon(ipft) = 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(ipft) = elongf_prev ! Elongation factor cannot decrease - currentSite%dstatus(ipft) = phen_dstat_timeon ! Flag that this has been forced + currentSite%phen_status(ipft) = 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(ipft) = 0.0_r8 ! Elongation factor must remain 0. - currentSite%dstatus(ipft) = phen_dstat_timeoff ! Flag that this has been forced + currentSite%phen_status(ipft) = 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. @@ -1477,9 +1476,9 @@ subroutine phenology( currentSite, bc_in ) if (elongf_prev >= elongf_min ) then ! This is the first day moisture fell below minimum. Flag change of status. - currentSite%dstatus(ipft) = phen_dstat_moistoff ! Flag that this has not been forced - currentSite%dleafoffdate(ipft) = model_day_int ! Record leaf off date - currentSite%dndaysleafoff(ipft) = 0 ! Reset clock + currentSite%phen_status(ipft) = phen_dstat_moistoff ! Flag that this has not been forced + currentSite%leafoffdate(ipft) = model_day_int ! Record leaf off date + currentSite%ndaysleafoff(ipft) = 0 ! Reset clock end if else @@ -1489,40 +1488,15 @@ subroutine phenology( currentSite, bc_in ) if (elongf_prev < elongf_min ) then ! This is the first day moisture allows leaves to exist. Flag change of status. - currentSite%dstatus(ipft) = phen_dstat_moiston ! Flag that this has not been forced - currentSite%dleafondate(ipft) = model_day_int ! Record leaf on date - currentSite%dndaysleafon(ipft) = 0 ! Reset clock + currentSite%phen_status(ipft) = phen_dstat_moiston ! Flag that this has not been forced + currentSite%leafondate(ipft) = model_day_int ! Record leaf on date + currentSite%ndaysleafon(ipft) = 0 ! Reset clock elseif (elongf_1st < elongf_prev) then - currentSite%dstatus(ipft) = phen_dstat_pshed ! Flag partial shedding, + currentSite%phen_status(ipft) = phen_dstat_pshed ! Flag partial shedding, ! but do not reset the clock end if end if drought_gradual_ifelse - - - case default - ! Neither hard deciduous or semi-deciduous. For now we treat this as synonym - ! of non-drought deciduous. In the future we may consider other drought deciduous - ! strategies (e.g., abscission driven by moisture, flushing driven by photo- - ! period). - currentSite%dstatus(ipft) = phen_dstat_moiston - - ! Assign elongation factors for non-drought deciduous PFTs, which will be used - ! to define the cohort status. - case_cold_phen: select case(prt_params%season_decid(ipft)) - case (ifalse) - ! Evergreen, ensure that elongation factor is always one. - currentSite%elong_factor(ipft) = 1.0_r8 - case (itrue) - ! Cold-deciduous. Define elongation factor based on cold status - select case (currentSite%cstatus) - case (phen_cstat_nevercold,phen_cstat_iscold) - currentSite%elong_factor(ipft) = 0.0_r8 - case (phen_cstat_notcold) - currentSite%elong_factor(ipft) = 1.0_r8 - end select - end select case_cold_phen - - end select case_drought_phen + end select case_update_phen end do pft_elong_loop @@ -1611,31 +1585,36 @@ subroutine phenology_leafonoff(currentSite) ! 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. - is_time_block: if (prt_params%season_decid(ipft) == itrue) then ! Cold deciduous + is_time_block: select case (prt_params%phen_leaf_habit(ipft)) + case (ihard_season_decid) ! Cold deciduous ! 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 + is_flushing_time = & + any( currentSite%phen_status(ipft) == [phen_cstat_tempon,phen_cstat_timeon]) .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. + is_shedding_time = & + any(currentSite%phen_status(ipft) == [phen_cstat_tempoff,phen_cstat_timeoff]) .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 (any(prt_params%stress_decid(ipft) == [ihard_stress_decid,isemi_stress_decid]) ) then ! Drought deciduous + case (ihard_stress_decid,isemi_stress_decid) ! Drought deciduous ! 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_shedding] ) + is_flushing_time = & + any( currentSite%phen_status(ipft) == [phen_dstat_moiston,phen_dstat_timeon] ) .and. & ! Leaf flushing time (moisture or time) + any( currentCohort%status_coh == [leaves_off,leaves_shedding] ) ! 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_shedding] ) - else + is_shedding_time = & + any( currentSite%phen_status(ipft) == [phen_dstat_moistoff,phen_dstat_timeoff,phen_dstat_pshed] ) .and. & + any( currentCohort%status_coh == [leaves_on,leaves_shedding] ) + case (ievergreen) ! This PFT is not deciduous. is_flushing_time = .false. is_shedding_time = .false. - end if is_time_block + end select is_time_block @@ -2319,7 +2298,7 @@ subroutine SeedDecay( litt , currentPatch, bc_in ) end subroutine SeedDecay ! ============================================================================ - subroutine SeedGermination( litt, cold_stat, drought_stat, bc_in, currentPatch ) + subroutine SeedGermination( litt, phen_stat, bc_in, currentPatch ) ! ! !DESCRIPTION: ! Flux from seed bank into the seedling pool @@ -2329,8 +2308,7 @@ subroutine SeedGermination( litt, cold_stat, drought_stat, bc_in, currentPatch ) ! ! !ARGUMENTS type(litter_type) :: litt - integer , intent(in) :: cold_stat ! Is the site in cold leaf-off status? - integer, dimension(numpft), intent(in) :: drought_stat ! Is the site in drought leaf-off status? + integer, dimension(numpft), intent(in) :: phen_stat ! Phenological status type(bc_in_type), intent(in) :: bc_in type(fates_patch_type), intent(in) :: currentPatch ! @@ -2421,20 +2399,18 @@ subroutine SeedGermination( litt, cold_stat, drought_stat, bc_in, currentPatch ) litt%seed_germ_in(pft) = litt%seed(pft) * seedling_emerg_rate end if if_tfs_or_def - - !set the germination only under the growing season...c.xu - - if ((prt_params%season_decid(pft) == itrue ) .and. & - (any(cold_stat == [phen_cstat_nevercold,phen_cstat_iscold]))) then - ! no germination for all PFTs when cold - litt%seed_germ_in(pft) = 0.0_r8 - endif - ! Drought deciduous, halt germination when status is shedding, even leaves are not - ! completely abscissed. MLO - select case (prt_params%stress_decid(pft)) + select case (prt_params%phen_leaf_habit(pft)) + case (ihard_season_decid) + !set the germination only under the growing season...c.xu + if (any(phen_stat(pft) == [phen_cstat_tempoff,phen_cstat_timeoff])) then + ! no germination for all PFTs when cold + litt%seed_germ_in(pft) = 0.0_r8 + end if case (ihard_stress_decid,isemi_stress_decid) - if (any(drought_stat(pft) == [phen_dstat_timeoff,phen_dstat_moistoff,phen_dstat_pshed])) then + ! Drought deciduous, halt germination when status is shedding, even leaves are not + ! completely abscissed. MLO + if (any(phen_stat(pft) == [phen_dstat_timeoff,phen_dstat_moistoff,phen_dstat_pshed])) then litt%seed_germ_in(pft) = 0.0_r8 end if end select @@ -2542,23 +2518,24 @@ subroutine recruitment(currentSite, currentPatch, bc_in) efstem_coh = 1.0_r8 leaf_status = leaves_on - ! 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 - efleaf_coh = 0.0_r8 - effnrt_coh = 1.0_r8 - fnrt_drop_fraction - efstem_coh = 1.0_r8 - stem_drop_fraction - leaf_status = leaves_off - 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. - select case (prt_params%stress_decid(ft)) + ! look for cases in which leaves should be off + select case (prt_params%phen_leaf_habit(ft)) + case (ihard_season_decid) + select case(currentSite%phen_status(ft)) + case (phen_cstat_tempoff, phen_cstat_timeoff) + ! If the plant is seasonally (cold) deciduous, and the site status is flagged + ! as "cold", then set the cohort's status to leaves_off. + efleaf_coh = 0.0_r8 + effnrt_coh = 1.0_r8 - fnrt_drop_fraction + efstem_coh = 1.0_r8 - stem_drop_fraction + leaf_status = leaves_off + end select case (ihard_stress_decid, isemi_stress_decid) + ! 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. efleaf_coh = currentSite%elong_factor(ft) effnrt_coh = 1.0_r8 - (1.0_r8 - efleaf_coh)*fnrt_drop_fraction efstem_coh = 1.0_r8 - (1.0_r8 - efleaf_coh)*stem_drop_fraction diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 9b3b9ef919..9ae53b3025 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -12,7 +12,7 @@ module FatesPatchMod use FatesUtilsMod, only : check_hlm_list use FatesUtilsMod, only : check_var_real use FatesCohortMod, only : fates_cohort_type - use FatesRunningMeanMod, only : rmean_type, rmean_arr_type + use FatesRunningSummMod, only : rsumm_type, rsumm_arr_type use FatesLitterMod, only : litter_type use FatesFuelMod, only : fuel_type use PRTGenericMod, only : num_elements @@ -24,9 +24,9 @@ module FatesPatchMod use EDParamsMod, only : nlevleaf, nclmax, maxpft use FatesConstantsMod, only : n_dbh_bins, n_dist_types use FatesConstantsMod, only : t_water_freeze_k_1atm - use FatesRunningMeanMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm - use FatesRunningMeanMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par - use FatesRunningMeanMod, only : ema_sdlng2sap_par, ema_sdlng_mdd + use FatesRunningSummMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm + use FatesRunningSummMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par + use FatesRunningSummMod, only : ema_sdlng2sap_par, ema_sdlng_mdd use TwoStreamMLPEMod, only : twostream_type use FatesRadiationMemMod, only : num_swb use FatesRadiationMemMod, only : num_rad_stream_types @@ -75,23 +75,24 @@ module FatesPatchMod !--------------------------------------------------------------------------- ! RUNNING MEANS - !class(rmean_type), pointer :: t2m ! place-holder for 2m air temperature (variable window-size) - class(rmean_type), pointer :: tveg24 ! 24-hour mean vegetation temperature [K] - class(rmean_type), pointer :: tveg_lpa ! running mean of vegetation temperature at the + !class(rsumm_type), pointer :: t2m ! place-holder for 2m air temperature (variable window-size) + class(rsumm_type), pointer :: tveg24 ! 24-hour summary vegetation temperature [K] + class(rsumm_type), pointer :: tveg_lpa ! running mean of vegetation temperature at the ! leaf photosynthesis acclimation timescale [K] - class(rmean_type), pointer :: tveg_longterm ! long-term running mean of vegetation temperature at the + class(rsumm_type), pointer :: tveg_longterm ! long-term running mean of vegetation temperature at the ! leaf photosynthesis acclimation timescale [K] (i.e T_home) - class(rmean_type), pointer :: seedling_layer_par24 ! 24-hour mean of photosynthetically active radiation at seedling layer [W/m2] - class(rmean_arr_type), pointer :: sdlng_emerg_smp(:) ! running mean of soil matric potential at the seedling + class(rsumm_arr_type), pointer :: btran24_ft(:) ! 24-hour summary of transpiration wetness factor (aka btran) + class(rsumm_type), pointer :: seedling_layer_par24 ! 24-hour summary of photosynthetically active radiation at seedling layer [W/m2] + class(rsumm_arr_type), pointer :: sdlng_emerg_smp(:) ! running mean of soil matric potential at the seedling ! rooting depth at the H2O seedling emergence timescale (see sdlng_emerg_h2o_timescale parameter) - class(rmean_type), pointer :: sdlng_mort_par ! running mean of photosythetically active radiation + class(rsumm_type), pointer :: sdlng_mort_par ! running mean of photosythetically active radiation ! at the seedling layer and at the par-based seedling ! mortality timescale (sdlng_mort_par_timescale) - class(rmean_arr_type), pointer :: sdlng_mdd(:) ! running mean of moisture deficit days + class(rsumm_arr_type), pointer :: sdlng_mdd(:) ! running mean of moisture deficit days ! at the seedling layer and at the mdd-based seedling ! mortality timescale (sdlng_mdd_timescale) ! (sdlng2sap_par_timescale) - class(rmean_type), pointer :: sdlng2sap_par ! running mean of photosythetically active radiation + class(rsumm_type), pointer :: sdlng2sap_par ! running mean of photosythetically active radiation ! at the seedling layer and at the par-based seedling ! to sapling transition timescale ! (sdlng2sap_par_timescale) @@ -614,18 +615,28 @@ subroutine InitRunningMeans(this, current_tod, regeneration_model, numpft) ! Until bc's are pointed to by sites give veg a default temp [K] real(r8), parameter :: temp_init_veg = 15._r8 + t_water_freeze_k_1atm real(r8), parameter :: init_seedling_par = 5.0_r8 ! arbitrary initialization for seedling layer [MJ m-2 d-1] - real(r8), parameter :: init_seedling_smp = -26652.0_r8 ! abitrary initialization of smp [mm] + real(r8), parameter :: init_seedling_smp = -26652.0_r8 ! arbitrary initialization of smp [mm] + real(r8), parameter :: init_btran = 1.0_r8 ! arbitrary initial value for btran [fraction] integer :: pft ! pft looping index allocate(this%tveg24) allocate(this%tveg_lpa) allocate(this%tveg_longterm) + allocate(this%btran24_ft(numpft)) - ! set initial values for running means - call this%tveg24%InitRMean(fixed_24hr, init_value=temp_init_veg, & + ! set initial values for running summaries + call this%tveg24%InitRSumm(fixed_24hr, init_value=temp_init_veg, & init_offset=real(current_tod, r8)) - call this%tveg_lpa%InitRmean(ema_lpa, init_value=temp_init_veg) - call this%tveg_longterm%InitRmean(ema_longterm, init_value=temp_init_veg) + call this%tveg_lpa%InitRSumm(ema_lpa, init_value=temp_init_veg) + call this%tveg_longterm%InitRSumm(ema_longterm, init_value=temp_init_veg) + + do pft = 1,numpft + allocate(this%btran24_ft(pft)%p) + + call this%btran24_ft(pft)%p%InitRSumm(fixed_24hr,init_value=init_btran, & + init_offset=real(current_tod, r8)) + end do + if (regeneration_model == TRS_regeneration) then allocate(this%seedling_layer_par24) @@ -634,20 +645,20 @@ subroutine InitRunningMeans(this, current_tod, regeneration_model, numpft) allocate(this%sdlng_mort_par) allocate(this%sdlng2sap_par) - call this%seedling_layer_par24%InitRMean(fixed_24hr, & + call this%seedling_layer_par24%InitRSumm(fixed_24hr, & init_value=init_seedling_par, init_offset=real(current_tod, r8)) - call this%sdlng_mort_par%InitRMean(ema_sdlng_mort_par, & + call this%sdlng_mort_par%InitRSumm(ema_sdlng_mort_par, & init_value=temp_init_veg) - call this%sdlng2sap_par%InitRMean(ema_sdlng2sap_par, & + call this%sdlng2sap_par%InitRSumm(ema_sdlng2sap_par, & init_value=init_seedling_par) do pft = 1,numpft allocate(this%sdlng_mdd(pft)%p) allocate(this%sdlng_emerg_smp(pft)%p) - call this%sdlng_mdd(pft)%p%InitRMean(ema_sdlng_mdd, & + call this%sdlng_mdd(pft)%p%InitRSumm(ema_sdlng_mdd, & init_value=0.0_r8) - call this%sdlng_emerg_smp(pft)%p%InitRMean(ema_sdlng_emerg_h2o, & + call this%sdlng_emerg_smp(pft)%p%InitRSumm(ema_sdlng_emerg_h2o, & init_value=init_seedling_smp) end do end if @@ -911,7 +922,7 @@ subroutine FreeMemory(this, regeneration_model, numpft) call endrun(msg=errMsg(sourcefile, __LINE__)) endif - ! deallocate running means + ! deallocate running summaries deallocate(this%tveg24, stat=istat, errmsg=smsg) if (istat/=0) then write(fates_log(),*) 'dealloc011: fail on deallocate(this%tveg24):'//trim(smsg) @@ -928,6 +939,12 @@ subroutine FreeMemory(this, regeneration_model, numpft) call endrun(msg=errMsg(sourcefile, __LINE__)) endif + do pft = 1, numpft + deallocate(this%btran24_ft(pft)%p) + end do + deallocate(this%btran24_ft) + + if (regeneration_model == TRS_regeneration) then deallocate(this%seedling_layer_par24) deallocate(this%sdlng_mort_par) diff --git a/functional_unit_testing/parteh/PartehDriver.py b/functional_unit_testing/parteh/PartehDriver.py index 88ca296255..c8319865ae 100644 --- a/functional_unit_testing/parteh/PartehDriver.py +++ b/functional_unit_testing/parteh/PartehDriver.py @@ -176,35 +176,9 @@ def main(): for iplnt in range(num_plants): ipft = use_pfts[iplnt] - evergreen = np.int(fates_params['evergreen'].data[ipft]) - cold_deciduous = np.int(fates_params['season_decid'].data[ipft]) - stress_deciduous = np.int(fates_params['stress_decid'].data[ipft]) - if(evergreen==1): - if(cold_deciduous==1): - print("Poorly defined phenology mode 0") - exit(2) - if(stress_deciduous==1): - print("Poorly defined phenology mode 1") - exit(2) - phen_type.append(1) - elif(cold_deciduous==1): - if(evergreen==1): - print("Poorly defined phenology mode 2") - exit(2) - if(stress_deciduous==1): - print("Poorly defined phenology mode 3") - exit(2) - phen_type.append(2) - elif(stress_deciduous==1): - if(evergreen==1): - print("Poorly defined phenology mode 4") - exit(2) - if(cold_deciduous==1): - print("Poorly defined phenology mode 5") - exit(2) - phen_type.append(3) - else: - print("Unknown phenology mode ? {} {} {}".format(evergreen,cold_deciduous,stress_deciduous)) + phen_leaf_habit = np.int(fates_params['phen_leaf_habit'].data[ipft]) + if(phen_leaf_habit < 1 or phen_leaf_habit > 4): + print("Unknown phenology mode ? {}".format(phen_leaf_habit)) exit(2) diff --git a/functional_unit_testing/parteh/parteh_controls_phenevents_v2.xml b/functional_unit_testing/parteh/parteh_controls_phenevents_v2.xml index 18cf824c62..9d6d25f6eb 100644 --- a/functional_unit_testing/parteh/parteh_controls_phenevents_v2.xml +++ b/functional_unit_testing/parteh/parteh_controls_phenevents_v2.xml @@ -57,9 +57,7 @@ 1 , 1 , 2 , 2 , 2 - 1 , 0 , 1 , 0 , 0 - 0 , 1 , 0 , 1 , 1 - 0 , 0 , 0 , 0 , 0 + 1 , 2 , 1 , 2 , 2 0.2 , 0.2 , 0.2 , 0.2 , 0.2 0.2 , 0.2, 0.2, 0.2, 0.2 30.0 , 30.0 , 30.0, 30.0 , 30.0 diff --git a/functional_unit_testing/parteh/parteh_controls_smoketests.xml b/functional_unit_testing/parteh/parteh_controls_smoketests.xml index d7675c7276..662be047d1 100644 --- a/functional_unit_testing/parteh/parteh_controls_smoketests.xml +++ b/functional_unit_testing/parteh/parteh_controls_smoketests.xml @@ -57,9 +57,7 @@ 1 , 2 , 2 , 2 , 2 - 1 , 1 , 1 , 1 , 1 - 0 , 0 , 0 , 0 , 0 - 0 , 0 , 0 , 0 , 0 + 1 , 1 , 1 , 1 , 1 0.2 , 0.2 , 0.2 , 0.2 , 0.2 0.2 , 0.2, 0.2, 0.2, 0.2 30.0 , 30.0 , 30.0, 30.0 , 30.0 diff --git a/functional_unit_testing/parteh/parteh_controls_variable_netc.xml b/functional_unit_testing/parteh/parteh_controls_variable_netc.xml index 47be70426d..156ef5e47d 100644 --- a/functional_unit_testing/parteh/parteh_controls_variable_netc.xml +++ b/functional_unit_testing/parteh/parteh_controls_variable_netc.xml @@ -55,7 +55,7 @@ 1 , 2 , 2 - 1 , 1, 1 + 1 , 1, 1 0.2 , 0.2, 0.2 0.2 , 0.2 , 0.2 30.0 , 30.0 , 30.0 diff --git a/main/CMakeLists.txt b/main/CMakeLists.txt index 8db1d06331..5c52174d1c 100644 --- a/main/CMakeLists.txt +++ b/main/CMakeLists.txt @@ -22,7 +22,7 @@ list(APPEND fates_sources EDParamsMod.F90 EDTypesMod.F90 FatesHydraulicsMemMod.F90 - FatesRunningMeanMod.F90 + FatesRunningSummMod.F90 FatesParameterDerivedMod.F90 FatesSizeAgeTypeIndicesMod.F90 FatesIntegratorsMod.F90 diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 625e65fabc..7f6a5c299b 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -11,6 +11,7 @@ module EDInitMod use FatesConstantsMod , only : primaryland use FatesConstantsMod , only : nearzero use FatesConstantsMod , only : rsnbl_math_prec + use FatesConstantsMod , only : ndays_per_year use EDTypesMod , only : min_patch_area_forced use FatesConstantsMod , only : n_landuse_cats use FatesConstantsMod , only : is_crop @@ -22,6 +23,7 @@ module EDInitMod use FatesInterfaceTypesMod , only : hlm_is_restart use FatesInterfaceTypesMod , only : hlm_current_tod use FatesInterfaceTypesMod , only : hlm_regeneration_model + use FatesInterfaceTypesMod , only : hlm_model_day use EDPftvarcon , only : EDPftvarcon_inst use PRTParametersMod , only : prt_params use EDCohortDynamicsMod , only : create_cohort, fuse_cohorts @@ -41,16 +43,19 @@ module EDInitMod use EDTypesMod , only : init_spread_inventory use FatesConstantsMod , only : leaves_on use FatesConstantsMod , only : leaves_off + use FatesConstantsMod , only : ievergreen + use FatesConstantsMod , only : ihard_season_decid use FatesConstantsMod , only : ihard_stress_decid use FatesConstantsMod , only : isemi_stress_decid use PRTGenericMod , only : num_elements use PRTGenericMod , only : element_list - use EDTypesMod , only : phen_cstat_nevercold - use EDTypesMod , only : phen_cstat_iscold + use EDTypesMod , only : phen_cstat_timeoff + use EDTypesMod , only : phen_cstat_tempoff use EDTypesMod , only : phen_dstat_timeoff use EDTypesMod , only : phen_dstat_moistoff - use EDTypesMod , only : phen_cstat_notcold + use EDTypesMod , only : phen_cstat_tempon use EDTypesMod , only : phen_dstat_moiston + use EDTypesMod , only : phen_estat_evergreen use FatesInterfaceTypesMod , only : bc_in_type,bc_out_type use FatesInterfaceTypesMod , only : hlm_use_planthydro use FatesInterfaceTypesMod , only : hlm_use_inventory_init @@ -262,22 +267,18 @@ 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%grow_deg_days = nan ! growing degree days + site_in%phen_status(:) = fates_unset_int ! are leaves in this pixel on or off? + 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 (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%nchilldays(:) = fates_unset_int ! number of chilling days + site_in%ncolddays(:) = fates_unset_int ! number of cold days + site_in%leafondate(:) = fates_unset_int ! doy of last leaf flushing event + site_in%leafoffdate(:) = fates_unset_int ! doy of last leaf abscission event + site_in%ndaysleafon(:) = fates_unset_int ! days since last leaf flushing + site_in%ndaysleafoff(: ) = fates_unset_int ! days since last leaf abscission site_in%elong_factor(:) = nan ! Elongation factor (0 - full abscission; 1 - fully flushed) + site_in%btran_memory(:,:) = nan site_in%liqvol_memory(:,:) = nan site_in%smp_memory(:,:) = nan site_in%vegtemp_memory(:) = nan ! record of last 10 days temperature for senescence model. @@ -394,24 +395,25 @@ subroutine set_site_properties( nsites, sites,bc_in ) ! ! !LOCAL VARIABLES: integer :: s - integer :: cstat ! cold status phenology flag - real(r8) :: GDD - integer :: dstat ! drought status phenology flag - 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 - integer :: fates_pft ! used in fixed biogeog mode + real(r8) :: gdd_1st ! First guess for growing degree days + real(r8) :: btranmem_1st ! First guess for transpiration wetness factor memory + real(r8) :: liqvolmem_1st ! First guess for soil water content memory + real(r8) :: smpmem_1st ! First guess for soil matric potential memory + real(r8) :: tempmem_1st ! First guess for temperature memory + real(r8) :: elong_factor_1st ! First guess for elongation factor (fraction of leaves + ! flushed relative to the maximum given plant size) + ! 0 - completely leafless + ! 0.5 - 50% of leaves flushed + ! 1.0 - 100% of leaves flushed + integer :: leafon ! DOY for last leaf flushing event, initial guess + integer :: leafoff ! DOY for last leaf abscission event, initial guess + integer :: ndleafon ! Elapsed days since last flushing, initial guess + integer :: ndleafoff ! Elapsed days since last abscission, initial guess + integer :: ft ! PFT loop index + integer :: model_day_int ! Model integration date + real(r8) :: sumarea ! area of PFTs in nocomp mode. + integer :: hlm_pft ! used in fixed biogeog mode + integer :: fates_pft ! used in fixed biogeog mode integer :: i_landusetype real(r8) :: temp_vec(numpft) ! temporary vector integer :: i_pftcount @@ -425,45 +427,109 @@ subroutine set_site_properties( nsites, sites,bc_in ) if ( hlm_is_restart == ifalse ) then - GDD = 30.0_r8 - cleafon = 100 - cleafoff = 300 - cndleafon = 0 - cndleafoff = 0 - cstat = phen_cstat_notcold ! Leaves are on - dstat = phen_dstat_moiston ! Leaves are on - dleafoff = 300 - dleafon = 100 - dndleafon = 0 - dndleafoff = 0 - liqvolmem = 0.5_r8 - smpmem = 0._r8 - elong_factor = 1._r8 + ! For now, we impose that deciduous plants start the simulation fully flushed, + ! but with the possibility of switching status if needed. + gdd_1st = 30.0_r8 ! Assume sufficiently warm conditions + btranmem_1st = 1._r8 ! Assume well watered conditions + liqvolmem_1st = 0.5_r8 ! Assume well watered conditions + smpmem_1st = 0._r8 ! Assume well watered conditions + tempmem_1st = 15.0_r8 ! Assume sufficiently warm conditions + elong_factor_1st = 1._r8 ! Assume leaves fully flushed. + + ! Set the model integration date since the beginning of time. + model_day_int = floor(hlm_model_day) do s = 1,nsites - sites(s)%nchilldays = 0 - sites(s)%ncolddays = 0 ! recalculated in phenology - ! immediately, so yes this - ! is memory-less, but needed - ! for first value in history file - sites(s)%phen_model_date = 0 - sites(s)%cleafondate = cleafon - hlm_day_of_year - sites(s)%cleafoffdate = cleafoff - hlm_day_of_year - sites(s)%cndaysleafon = cndleafon - sites(s)%cndaysleafoff = cndleafoff - sites(s)%dleafoffdate (1:numpft) = dleafoff - hlm_day_of_year - sites(s)%dleafondate (1:numpft) = dleafon - hlm_day_of_year - sites(s)%dndaysleafon (1:numpft) = dndleafon - sites(s)%dndaysleafoff(1:numpft) = dndleafoff - sites(s)%grow_deg_days = GDD - - 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(1:numpft) = dstat - sites(s)%elong_factor(1:numpft) = elong_factor + ! Define the last flushing and abscission dates based on the hemisphere. This + ! is to avoid first guess dates that are too far from optimal in the Southern + ! Hemisphere, which could make deciduous plants spend multiple years until + ! synchronising with the local climate. For now, both cold- and drought-deciduous + ! dates default to leaf flushing during spring, and to leaf abscission during + ! autumn. This assumes winter dry season, which is common in many tropical + ! regions. This is still not ideal for some tropical areas and for Mediterranean + ! climates, in which peak greenness occurs during autumn or spring, respectively. + ! In the future, we could consider implementing some map that allows assigning + ! the expected dates for each location. + if (sites(s)%lat > 0._r8) then + leafon = 100 - hlm_day_of_year + leafoff = 300 - hlm_day_of_year + else + leafon = 280 - hlm_day_of_year + leafoff = 120 - hlm_day_of_year + end if + + + ! Make sure the initial dates for last leaf abscission and leaf flushing + ! events precede the the beginning of the simulation. This check used to occur + ! inside the integrator, but it is better to ensure this during initialisation, + ! so we reduce the number of logical tests during time steps. + if ( leafon > model_day_int ) then + leafon = leafon - ndays_per_year + end if + if ( leafoff > model_day_int ) then + leafoff = leafoff - ndays_per_year + end if + + ! Initialise the number of days since last phenological transitions as the + ! negative of the last transition dates. This is consistent with the original + ! calculation, but different from the original cold start that assumed that + ! both elapsed dates were + ndleafon = model_day_int - leafon + ndleafoff = model_day_int - leafoff + + + !---~--- + ! Initialise site-level variables that are not PFT specific or that can be + ! set for all PFTs without considering the leaf phenological habit. Some of + ! the memory and cumulative variables are currently used only for leaf + ! phenology, but they may be useful for other modules, so they are set for + ! all PFTs. + ! Some variables such as ncolddays are recalculated every day (memoryless) + ! but they need an initial value so they are defined for the first history + ! file. + !---~--- + sites(s)%phen_model_date = model_day_int + sites(s)%vegtemp_memory(1:num_vegtemp_mem) = tempmem_1st + sites(s)%btran_memory (1:numWaterMem,1:numpft) = btranmem_1st + sites(s)%liqvol_memory (1:numWaterMem,1:numpft) = liqvolmem_1st + sites(s)%smp_memory (1:numWaterMem,1:numpft) = smpmem_1st + sites(s)%grow_deg_days (1:numpft) = gdd_1st + sites(s)%nchilldays (1:numpft) = 0 + sites(s)%ncolddays (1:numpft) = 0 + sites(s)%elong_factor (1:numpft) = elong_factor_1st + + + + !---~--- + ! Fill in variables that do require checks for leaf phenology habit. + !---~--- + do_phen_pft_init: do ft = 1, numpft + case_leaf_habit: select case (prt_params%phen_leaf_habit(ft)) + case (ievergreen) + ! Evergreens. Phenological status is not used, set it accordingly. + sites(s)%phen_status (ft) = phen_estat_evergreen + case (ihard_season_decid) + ! Cold deciduous. Initialise both status and time since last flushing + ! and abscission days. We assume the initial status to be with leaves + ! fully flushed. + sites(s)%phen_status (ft) = phen_cstat_tempon + sites(s)%leafondate (ft) = leafon + sites(s)%leafoffdate (ft) = leafoff + sites(s)%ndaysleafon (ft) = ndleafon + sites(s)%ndaysleafoff(ft) = ndleafoff + case (ihard_stress_decid,isemi_stress_decid) + ! Drought deciduous. Initialise both status and time since last flushing + ! and abscission days. We assume the initial status to be with leaves + ! fully flushed. + sites(s)%phen_status (ft) = phen_dstat_moiston + sites(s)%leafondate (ft) = leafon + sites(s)%leafoffdate (ft) = leafoff + sites(s)%ndaysleafon (ft) = ndleafon + sites(s)%ndaysleafoff(ft) = ndleafoff + end select case_leaf_habit + end do do_phen_pft_init + !---~--- + sites(s)%NF = 0.0_r8 sites(s)%NF_successful = 0.0_r8 @@ -1143,15 +1209,24 @@ subroutine init_cohorts(site_in, patch_in, bc_in) efstem_coh = 1.0_r8 leaf_status = leaves_on else - ! use built-in phenology - if (prt_params%season_decid(pft) == itrue .and. & - any(site_in%cstatus == [phen_cstat_nevercold, phen_cstat_iscold])) then - ! Cold deciduous, off season, assume complete abscission - efleaf_coh = 0.0_r8 - effnrt_coh = 1.0_r8 - fnrt_drop_fraction - efstem_coh = 1.0_r8 - stem_drop_fraction - leaf_status = leaves_off - else if (any(prt_params%stress_decid(pft) == [ihard_stress_decid, isemi_stress_decid])) then + ! use built-in phenology + phen_select: select case (prt_params%phen_leaf_habit(pft)) + case (ihard_season_decid) + select case (site_in%phen_status(pft)) + case (phen_cstat_tempoff, phen_cstat_timeoff) + ! Cold deciduous, off season, assume complete abscission + efleaf_coh = 0.0_r8 + effnrt_coh = 1.0_r8 - fnrt_drop_fraction + efstem_coh = 1.0_r8 - stem_drop_fraction + leaf_status = leaves_off + case default + ! Cold deciduous, growing season, assume leaves fully flushed + efleaf_coh = 1.0_r8 + effnrt_coh = 1.0_r8 + efstem_coh = 1.0_r8 + leaf_status = leaves_on + end select + case (ihard_stress_decid, isemi_stress_decid) ! 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 @@ -1167,14 +1242,13 @@ subroutine init_cohorts(site_in, patch_in, bc_in) else leaf_status = leaves_off end if - else - ! Evergreens, or deciduous during growing season - ! Assume leaves fully flushed - efleaf_coh = 1.0_r8 - effnrt_coh = 1.0_r8 - efstem_coh = 1.0_r8 - leaf_status = leaves_on - end if + case (ievergreen) + ! Evergreens, assume leaves fully flushed + efleaf_coh = 1.0_r8 + effnrt_coh = 1.0_r8 + efstem_coh = 1.0_r8 + leaf_status = leaves_on + end select phen_select end if if_spmode ! If positive EDPftvarcon_inst%initd is interpreted as initial recruit density. diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index bfba6d7d56..fc6989ff32 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -106,6 +106,7 @@ module EDMainMod use EDPftvarcon, only : EDPftvarcon_inst use FatesHistoryInterfaceMod, only : fates_hist use FatesLandUseChangeMod, only: FatesGrazing + use FatesCumulativeMemoryMod, only : UpdateCumulativeMemoryVars ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -168,6 +169,10 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out) end do call currentSite%flux_diags%ZeroFluxDiags() + ! Call a routine that will compute cumulative variables and "memory" averages for + ! a suite of variables. These variables are mostly used for leaf phenology, but they + ! may be useful for other components (disturbances, mortality, management). + call UpdateCumulativeMemoryVars(currentSite,bc_in) ! Call a routine that simply identifies if logging should occur ! This is limited to a global event until more structured event handling is enabled @@ -553,7 +558,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(ft) == [phen_dstat_moiston,phen_dstat_timeon])) then + if(any(currentSite%phen_status(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 b90b547dad..10567c1e86 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -52,13 +52,6 @@ module EDParamsMod real(r8),protected, public :: ED_val_cwd_fcel ! Cellulose fraction for CWD real(r8),protected, public :: ED_val_cwd_flig ! Lignin fraction of coarse woody debris real(r8),protected, public :: maintresp_nonleaf_baserate ! Base maintenance respiration rate for plant tissues - real(r8),protected, public :: ED_val_phen_a ! GDD accumulation function, intercept parameter: gdd_thesh = a + b exp(c*ncd) - real(r8),protected, public :: ED_val_phen_b ! GDD accumulation function, multiplier parameter: gdd_thesh = a + b exp(c*ncd) - real(r8),protected, public :: ED_val_phen_c ! GDD accumulation function, exponent parameter: gdd_thesh = a + b exp(c*ncd) - real(r8),protected, public :: ED_val_phen_chiltemp ! chilling day counting threshold for vegetation - real(r8),protected, public :: ED_val_phen_mindayson ! day threshold compared against days since leaves became on-allometry - real(r8),protected, public :: ED_val_phen_ncolddayslim ! day threshold exceedance for temperature leaf-drop - real(r8),protected, public :: ED_val_phen_coldtemp ! vegetation temperature exceedance that flags a cold-day for leaf-drop real(r8),protected, public :: ED_val_cohort_size_fusion_tol ! minimum fraction in difference in dbh between cohorts real(r8),protected, public :: ED_val_cohort_age_fusion_tol ! minimum fraction in differece in cohort age between cohorts real(r8),protected, public :: ED_val_patch_fusion_tol ! minimum fraction in difference in profiles between patches @@ -141,13 +134,6 @@ module EDParamsMod character(len=param_string_length),parameter,public :: ED_name_cwd_fcel= "fates_frag_cwd_fcel" character(len=param_string_length),parameter,public :: ED_name_cwd_flig= "fates_frag_cwd_flig" character(len=param_string_length),parameter,public :: fates_name_maintresp_nonleaf_baserate= "fates_maintresp_nonleaf_baserate" - character(len=param_string_length),parameter,public :: ED_name_phen_a= "fates_phen_gddthresh_a" - character(len=param_string_length),parameter,public :: ED_name_phen_b= "fates_phen_gddthresh_b" - character(len=param_string_length),parameter,public :: ED_name_phen_c= "fates_phen_gddthresh_c" - character(len=param_string_length),parameter,public :: ED_name_phen_chiltemp= "fates_phen_chilltemp" - character(len=param_string_length),parameter,public :: ED_name_phen_mindayson= "fates_phen_mindayson" - character(len=param_string_length),parameter,public :: ED_name_phen_ncolddayslim= "fates_phen_ncolddayslim" - character(len=param_string_length),parameter,public :: ED_name_phen_coldtemp= "fates_phen_coldtemp" character(len=param_string_length),parameter,public :: ED_name_cohort_size_fusion_tol= "fates_cohort_size_fusion_tol" character(len=param_string_length),parameter,public :: ED_name_cohort_age_fusion_tol = "fates_cohort_age_fusion_tol" character(len=param_string_length),parameter,public :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol" @@ -310,13 +296,6 @@ subroutine FatesParamsInit() ED_val_cwd_fcel = nan ED_val_cwd_flig = nan maintresp_nonleaf_baserate = nan - ED_val_phen_a = nan - ED_val_phen_b = nan - ED_val_phen_c = nan - ED_val_phen_chiltemp = nan - ED_val_phen_mindayson = nan - ED_val_phen_ncolddayslim = nan - ED_val_phen_coldtemp = nan ED_val_cohort_size_fusion_tol = nan ED_val_cohort_age_fusion_tol = nan ED_val_patch_fusion_tol = nan @@ -422,27 +401,6 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=fates_name_maintresp_nonleaf_baserate, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=ED_name_phen_a, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - - call fates_params%RegisterParameter(name=ED_name_phen_b, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - - call fates_params%RegisterParameter(name=ED_name_phen_c, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - - call fates_params%RegisterParameter(name=ED_name_phen_chiltemp, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - - call fates_params%RegisterParameter(name=ED_name_phen_mindayson, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - - call fates_params%RegisterParameter(name=ED_name_phen_ncolddayslim, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - - call fates_params%RegisterParameter(name=ED_name_phen_coldtemp, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=ED_name_cohort_size_fusion_tol, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -632,27 +590,6 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetrieveParameter(name=fates_name_maintresp_nonleaf_baserate, & data=maintresp_nonleaf_baserate) - call fates_params%RetrieveParameter(name=ED_name_phen_a, & - data=ED_val_phen_a) - - call fates_params%RetrieveParameter(name=ED_name_phen_b, & - data=ED_val_phen_b) - - call fates_params%RetrieveParameter(name=ED_name_phen_c, & - data=ED_val_phen_c) - - call fates_params%RetrieveParameter(name=ED_name_phen_chiltemp, & - data=ED_val_phen_chiltemp) - - call fates_params%RetrieveParameter(name=ED_name_phen_mindayson, & - data=ED_val_phen_mindayson) - - call fates_params%RetrieveParameter(name=ED_name_phen_ncolddayslim, & - data=ED_val_phen_ncolddayslim) - - call fates_params%RetrieveParameter(name=ED_name_phen_coldtemp, & - data=ED_val_phen_coldtemp) - call fates_params%RetrieveParameter(name=ED_name_cohort_size_fusion_tol, & data=ED_val_cohort_size_fusion_tol) @@ -832,13 +769,6 @@ subroutine FatesReportParams(is_master) write(fates_log(),fmt0) 'ED_val_cwd_fcel = ',ED_val_cwd_fcel write(fates_log(),fmt0) 'ED_val_cwd_flig = ',ED_val_cwd_flig write(fates_log(),fmt0) 'fates_maintresp_nonleaf_baserate = ', maintresp_nonleaf_baserate - write(fates_log(),fmt0) 'ED_val_phen_a = ',ED_val_phen_a - write(fates_log(),fmt0) 'ED_val_phen_b = ',ED_val_phen_b - write(fates_log(),fmt0) 'ED_val_phen_c = ',ED_val_phen_c - write(fates_log(),fmt0) 'ED_val_phen_chiltemp = ',ED_val_phen_chiltemp - write(fates_log(),fmt0) 'ED_val_phen_mindayson = ',ED_val_phen_mindayson - write(fates_log(),fmt0) 'ED_val_phen_ncolddayslim = ',ED_val_phen_ncolddayslim - write(fates_log(),fmt0) 'ED_val_phen_coldtemp = ',ED_val_phen_coldtemp write(fates_log(),fmt0) 'ED_val_cohort_size_fusion_tol = ',ED_val_cohort_size_fusion_tol write(fates_log(),fmt0) 'ED_val_cohort_age_fusion_tol = ',ED_val_cohort_age_fusion_tol write(fates_log(),fmt0) 'ED_val_patch_fusion_tol = ',ED_val_patch_fusion_tol diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 319fc5ff51..835ffab36f 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -26,6 +26,7 @@ module EDPftvarcon use FatesInterfaceTypesMod, only : hlm_nitrogen_spec, hlm_phosphorus_spec use FatesInterfaceTypesMod, only : hlm_parteh_mode use FatesInterfaceTypesMod, only : hlm_nu_com + use FatesConstantsMod , only : ievergreen use FatesConstantsMod , only : prescribed_p_uptake use FatesConstantsMod , only : prescribed_n_uptake use FatesConstantsMod , only : coupled_p_uptake @@ -1673,9 +1674,8 @@ subroutine FatesCheckParams(is_master) ! This subroutine performs logical checks on user supplied parameters. It cross ! compares various parameters and will fail if they don't make sense. ! Examples: - ! A tree can not be defined as both evergreen and deciduous. A woody plant - ! cannot have a structural biomass allometry intercept of 0, and a non-woody - ! plant (grass) can't have a non-zero intercept... + ! A woody plant cannot have a structural biomass allometry intercept of 0, and a + ! non-woody plant (grass) can't have a non-zero intercept... ! ----------------------------------------------------------------------------------- use FatesConstantsMod , only : fates_check_param_set use FatesConstantsMod , only : itrue, ifalse @@ -1962,7 +1962,7 @@ subroutine FatesCheckParams(is_master) ! Check if the fraction of storage used for flushing deciduous trees ! is greater than zero, and less than or equal to 1. - if (prt_params%evergreen(ipft) == ifalse) then + if (prt_params%phen_leaf_habit(ipft) /= ievergreen) then if ( ( EDPftvarcon_inst%phenflush_fraction(ipft) < nearzero ) .or. & ( EDPFtvarcon_inst%phenflush_fraction(ipft) > 1 ) ) then @@ -1970,7 +1970,8 @@ subroutine FatesCheckParams(is_master) write(fates_log(),*) ' on bud-burst. If phenflush_fraction is not greater than 0' write(fates_log(),*) ' it will not be able to put out any leaves. Plants need leaves.' write(fates_log(),*) ' PFT#: ',ipft - write(fates_log(),*) ' evergreen flag: (should be 0):',int(prt_params%evergreen(ipft)) + write(fates_log(),*) ' phen_leaf_habit: (evergreen should be ',ievergreen,'):', & + int(prt_params%phen_leaf_habit(ipft)) write(fates_log(),*) ' phenflush_fraction: ', EDPFtvarcon_inst%phenflush_fraction(ipft) write(fates_log(),*) ' Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index c21cdd6fe1..535a7f841d 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -25,7 +25,7 @@ module EDTypesMod use FatesLitterMod, only : litter_type use FatesLitterMod, only : ncwd use FatesConstantsMod, only : days_per_year - use FatesRunningMeanMod, only : rmean_type,rmean_arr_type + use FatesRunningSummMod, only : rsumm_type,rsumm_arr_type use FatesConstantsMod, only : fates_unset_r8 use FatesInterfaceTypesMod,only : bc_in_type use FatesInterfaceTypesMod,only : bc_out_type @@ -93,18 +93,24 @@ module EDTypesMod ! BIOLOGY/BIOGEOCHEMISTRY integer , parameter, public :: num_vegtemp_mem = 10 ! Window of time over which we track temp for cold sensecence (days) - ! Phenology status flag definitions (cold type is cstat, dry type is dstat) + ! Phenology status flag definitions (cold type is cstat, dry type is dstat). + ! MLO - Now that the phenology status is shared across both cold and drought deciduous, make sure these indices are unique + ! The default phen_estat_evergreen is a dummy value reserved for evergreen PFTs. - integer, parameter, public :: phen_cstat_nevercold = 0 ! This (location/plant) has not experienced a cold period over a large number - ! of days, leaves are dropped and flagged as non-cold region - 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_estat_evergreen = 0 ! This PFT is evergreen. - 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) + ! Season (cold) deciduous flags + integer, parameter, public :: phen_cstat_timeoff = 1 ! - Leaves off due to time, leaves exceeded life span without cold days + integer, parameter, public :: phen_cstat_tempoff = 2 ! - Leaves off after reaching multiple cold days + integer, parameter, public :: phen_cstat_tempon = 3 ! - Leaves on after exceeding the growing degree day threshold + integer, parameter, public :: phen_cstat_timeon = 4 ! - Leaves on after extended cold period, forced flushing + + ! Stress (drought) deciduous flags + integer, parameter, public :: phen_dstat_timeoff = 5 ! Leaves off due to time, leaves exceeded life span without desiccation + integer, parameter, public :: phen_dstat_moistoff = 6 ! Leaves off after moisture availability fell below the threshold + integer, parameter, public :: phen_dstat_moiston = 7 ! Leaves on after moisture availability increased above the threshold + integer, parameter, public :: phen_dstat_timeon = 8 ! Leaves on after extended dry period, forced flushing + integer, parameter, public :: phen_dstat_pshed = 9 ! Leaves partially abscissing due to drying conditions ! PATCH FUSION real(r8), parameter, public :: force_patchfuse_min_biomass = 0.005_r8 ! min biomass (kg / m2 patch area) below which to force-fuse patches @@ -414,40 +420,38 @@ module EDTypesMod ! PHENOLOGY - real(r8) :: grow_deg_days ! Phenology growing degree days real(r8) :: snow_depth ! site-level snow depth (used for ELAI/TLAI calcs) - integer :: cstatus ! are leaves in this pixel on or off for cold decid - ! 0 = this site has not experienced a cold period over at least + real(r8), dimension(maxpft) :: grow_deg_days ! Phenology growing degree days + integer, dimension(maxpft) :: phen_status ! Phenology status of this PFT. This is shared by all + ! phenology habits: + ! 0 = this is an evergreen PFT + ! 1 = this site has not experienced a cold period over at least ! 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(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 :: 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 + ! 2 = this site is in a cold-state where leaves should have fallen + ! 3 = this site is in a warm-state where leaves are allowed to flush + ! 4 = leaves off due to time exceedance + ! 5 = leaves off due to moisture avail + ! 6 = leaves on due to moisture avail + ! 7 = leaves on due to time exceedance + ! 8 = leaves partially on (ED2-like phenology) + integer, dimension(maxpft) :: nchilldays ! Number of chilling days (for Botta's GDD trheshold calculation) + integer, dimension(maxpft) :: ncolddays ! Number of cold days (for temperature-driven leaf abscission) + real(r8), dimension(num_vegtemp_mem) :: vegtemp_memory ! record of last 10 days temperature for senescence model. deg C + integer, dimension(maxpft) :: leafondate ! Model date (day integer) of last leaf flushing event + integer, dimension(maxpft) :: leafoffdate ! Model date (day integer) of last leaf abscission event + integer, dimension(maxpft) :: ndaysleafon ! Number of days since last leaf flushing event + integer, dimension(maxpft) :: ndaysleafoff ! Number of days since last leaf abscission event + real(r8), dimension(maxpft) :: elong_factor ! Elongation factor (ED2-like phenology). This is zero when leaves are ! completely off, and one when they are completely flushed. integer :: phen_model_date ! current model date (day integer) ! this date stays continuous when ! in runs that are restarted, regardless of ! the conditions of restart - 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) + real(r8), dimension(numWaterMem,maxpft) :: btran_memory ! last 10 days of transpiration wetness factor (drought phenology) + real(r8), dimension(numWaterMem,maxpft) :: liqvol_memory ! last 10 days of soil liquid water volume (drought phenology) + real(r8), dimension(numWaterMem,maxpft) :: smp_memory ! last 10 days of soil matric potential (drought phenology) ! FIRE diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 0c69199f4c..4d535f9de4 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -72,21 +72,37 @@ module FatesConstantsMod ! but is shedding them (partial shedding). This plant ! should not allocate carbon towards growth or ! reproduction. -integer, parameter, public :: ihard_stress_decid = 1 ! If the PFT is stress (drought) deciduous, - ! this flag is used to tell that the PFT - ! is a "hard" deciduous (i.e., the plant - ! has only two statuses, the plant either - ! sheds all leaves when it's time, or seeks - ! to flush the leaves back to allometry - ! when conditions improve. -integer, parameter, public :: isemi_stress_decid = 2 ! If the PFT is stress (drought) deciduous, - ! this flag is used to tell that the PFT - ! is a semi-deciduous (i.e., the plant - ! can downregulate the amount of leaves - ! relative to the allometry based on - ! soil moisture conditions. It can still - ! shed all leaves if conditions are very - ! dry. +integer, parameter, public :: ievergreen = 1 ! Flag that indicates that the leaf phenology + ! habit for a PFT is evergreen. This means + ! that seasonal environmental conditions do + ! not directly impact leaf biomass, although + ! the total leaf biomass can fall below + ! allometry if the plant's NPP is negative + ! and causes a significant depletion of the + ! storage pool. +integer, parameter, public :: ihard_season_decid = 2 ! Flag that indicates that the leaf phenology + ! habit for a PFT is a "hard" season (cold) + ! deciduous. This means that the plant + ! has only two statuses, the plant either + ! abscises all leaves when conditions + ! deteriorate, or flushes leaves to bring it + ! back to allometry when conditions improve. +integer, parameter, public :: ihard_stress_decid = 3 ! Flag that indicates that the leaf phenology + ! habit for a PFT is a "hard" stress + ! (drought) deciduous. This means that the + ! plant has only two statuses, the plant + ! either abscises all leaves when conditions + ! deteriorate, or flushes leaves to bring it + ! back to allometry when conditions improve +integer, parameter, public :: isemi_stress_decid = 4 ! Flag that indicates that the leaf phenology + ! habit for a PFT is a stress (hydro) + ! semi-deciduous. This means that the plant + ! can partially abscise or flush leaves + ! based on water availability, and + ! conditions. It can still abscise all leaves + ! when conditions are very dry, and flush all + ! leaves back to allometry when water is not + ! limiting. 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 diff --git a/main/FatesCumulativeMemoryMod.F90 b/main/FatesCumulativeMemoryMod.F90 new file mode 100644 index 0000000000..405403d898 --- /dev/null +++ b/main/FatesCumulativeMemoryMod.F90 @@ -0,0 +1,292 @@ +module FatesCumulativeMemoryMod + !---~--- + ! This module contains procedures that update multiple cumulative/memory variables + ! that drive leaf phenology, and potentially mortality, disturbances, and management. + ! These sub-routines are related to FatesRunningMeanMod, however, they need a separate + ! module because many of these procedures require access to EDTypesMod, which in turn + ! depends on FatesRunningMeanMod. + !---~--- + + + + use EDBtranMod , only : check_layer_water + use EDTypesMod , only : ed_site_type + use EDTypesMod , only : num_vegtemp_mem + use EDTypesMod , only : numWaterMem + use FatesAllometryMod , only : set_root_fraction + use FatesConstantsMod , only : ievergreen + use FatesConstantsMod , only : ihard_season_decid + use FatesConstantsMod , only : ihard_stress_decid + use FatesConstantsMod , only : isemi_stress_decid + use FatesConstantsMod , only : ndays_per_year + use FatesConstantsMod , only : nearzero + use FatesConstantsMod , only : nocomp_bareground + use FatesConstantsMod , only : r8 => fates_r8 + use FatesConstantsMod , only : tfrz => t_water_freeze_k_1atm + use FatesInterfaceTypesMod, only : bc_in_type + use FatesInterfaceTypesMod, only : hlm_model_day + use FatesInterfaceTypesMod, only : numpft + use FatesPatchMod , only : fates_patch_type + use PRTParametersMod , only : prt_params + + implicit none + private + + !---~--- + ! The main sub-routine that updates all cumulative and memory-related variables. + ! This is the only sub-routine that needs to be public, the specific procedures that + ! update time, temperature and moisture memories are called by this sub-routine. + !---~--- + public :: UpdateCumulativeMemoryVars + + !---~--- + ! Imposed soil matric potential lower bound for frozen or excessively dry soils, + ! used when computing water stress. + !---~--- + real(r8), public, parameter :: smp_lwr_bound = -1000000._r8 + +contains + + subroutine UpdateCumulativeMemoryVars(currentSite,bc_in) + !---~--- + ! This sub-routine updates multiple variables that retain cumulative or "memory" + ! related information. These are mostly used for phenology, but some of the variables + ! can be useful for disturbances, mortality and management. + !---~--- + + ! Arguments + type(ed_site_type), intent(inout), target :: currentSite + type(bc_in_type), intent(in) :: bc_in + + + ! Update elapsed time + call UpdatePhenologyDate(currentSite) + + ! Update temperature-related variables. Currently this is just the temperature + ! memory, but it will eventually contain growing degree days and number of + ! chilling days. + call UpdateCumulativeThermal(currentSite,bc_in) + + ! Update moisture-related memory variables. + call UpdateMemoryMoisture(currentSite,bc_in) + + end subroutine UpdateCumulativeMemoryVars + + + + + + + subroutine UpdatePhenologyDate(currentSite) + !---~--- + ! This sub-routine simply updates the number of elapsed days in this simulation. + ! The first day of the simulation is 1, and it continues monotonically, + ! indefinitely. + !---~--- + + ! Arguments + type(ed_site_type), intent(inout), target :: currentSite + ! Local variables + integer :: ipft ! PFT index + + !---~--- + ! Advance elapsed time, by using the host land model variable. The only reason + ! this is a site variable instead of a global variable is that we need to save this + ! information to the restart file, and we do not have global scalars in the restart + ! file. This value starts at zero and increases indefinitely. + !---~--- + currentSite%phen_model_date = floor(hlm_model_day) + + + ! Update the number of days since last flushing and abscission events. This is + ! done for deciduous PFTs only. + do ipft = 1, numpft + select case (prt_params%phen_leaf_habit(ipft)) + case (ihard_season_decid,ihard_stress_decid,isemi_stress_decid) + !---~--- + ! Update the number of days since last flushing and abscission events, by + ! using the current date and the dates of the last events. Note that we no + ! longer need to check whether this is the beginning of the simulation, + ! because the cold start initialisation already fixes the dates of the last + ! event to be prior to the beginning of the simulation. + !---~--- + currentSite%ndaysleafoff(ipft) = & + currentSite%phen_model_date - currentSite%leafoffdate(ipft) + currentSite%ndaysleafon (ipft) = & + currentSite%phen_model_date - currentSite%leafondate (ipft) + end select + end do + !---~--- + + end subroutine UpdatePhenologyDate + + + + + + + subroutine UpdateCumulativeThermal(currentSite,bc_in) + !---~--- + ! Currently this sub-routine updates the canopy temperature memory variable only. + ! In the future, this will also host updates to growing degree days and number of + ! chilling days. + !---~--- + + ! Arguments + type(ed_site_type), intent(inout), target :: currentSite + type(bc_in_type), intent(in) :: bc_in + + ! Local variables + type(fates_patch_type), pointer :: cpatch ! Current patch + real(r8) :: temp_in_C ! Daily averaged temperature in degC + real(r8) :: temp_wgt ! Canopy area weighting factor for daily average + ! vegetation temperature calculation + + ! Find the average temperature for the previous 24 hours + temp_in_C = 0._r8 + temp_wgt = 0._r8 + cpatch => CurrentSite%oldest_patch + do while( associated(cpatch) ) + temp_in_C = temp_in_C + cpatch%tveg24%GetMean()*cpatch%total_canopy_area + temp_wgt = temp_wgt + cpatch%total_canopy_area + cpatch => cpatch%younger + end do + if ( temp_wgt > nearzero ) then + temp_in_C = temp_in_C/temp_wgt - tfrz + else + ! If there is no canopy area, we use the vegetation temperature + ! of the first patch, which is the forcing air temperature + ! as defined in CLM/ELM. The forcing air temperature + ! should be the same across all patches. (Although + ! it is unlikely there are more than 1 in this scenario) + temp_in_C = CurrentSite%oldest_patch%tveg24%GetMean() - tfrz + end if + + + ! Record temperature for the last 10 days to the memory variable. + currentSite%vegtemp_memory(2:num_vegtemp_mem) = currentSite%vegtemp_memory(1:num_vegtemp_mem-1) + currentSite%vegtemp_memory(1) = temp_in_C + + + end subroutine UpdateCumulativeThermal + + + + + + + subroutine UpdateMemoryMoisture(currentSite,bc_in) + !---~--- + ! This sub-routine updates the moisture-related memory variables. These variables + ! are used for leaf phenology, and may become useful for other modules too. + !---~--- + + ! Arguments + type(ed_site_type), intent(inout), target :: currentSite + type(bc_in_type), intent(in) :: bc_in + + ! Local variables + type(fates_patch_type), pointer :: cpatch ! Current patch + integer :: ipft ! Plant Functional Type index + integer :: i_wmem ! Loop counter for water memory days + integer :: j ! Soil layer index + integer :: nlevroot ! Number of rooting levels to consider + real(r8) :: rootfrac_notop ! Total rooting fraction excluding the top soil layer + real(r8) :: min_btran ! Minimum transpiration wetness factor + real(r8) :: site_veg_area ! Fraction of the site area that is not bare + + + ! Transfer the transpiration wetness factor memory (we always track the last 10 + ! days). We shift the memory from days to the previous day, and make room for + ! current day + pft_btransfer_loop: do ipft=1,numpft + do i_wmem = numWaterMem,2,-1 + currentSite%btran_memory (i_wmem,ipft) = currentSite%btran_memory (i_wmem-1,ipft) + end do + end do pft_btransfer_loop + + + ! We now loop through all the patches to update the current-day minimum btran + currentSite%btran_memory (1,:) = 0._r8 + site_veg_area = 0._r8 + cpatch => CurrentSite%oldest_patch + patch_loop: do while( associated(cpatch) ) + ! Bypass bare patches + if_not_bare: if (cpatch%nocomp_pft_label /= nocomp_bareground) then + pft_btranupdate_loop: do ipft=1,numpft + currentSite%btran_memory(1,ipft) = currentSite%btran_memory (1,ipft) + & + cpatch%btran24_ft(ipft)%p%GetMin() * cpatch%area + end do pft_btranupdate_loop + + site_veg_area = site_veg_area + cpatch%area + end if if_not_bare + + cpatch => cpatch%younger + end do patch_loop + + ! Check if there is any vegetated patch in this site. + if (site_veg_area > nearzero) then + ! Normalise site-level btran by the vegetated area + currentSite%btran_memory(1,:) = currentSite%btran_memory(1,:) / site_veg_area + else + ! Dummy btran for non-vegetated area. This shouldn't be used by anything. + currentSite%btran_memory(1,:) = 0.5_r8 + end if + + + ! The soil memory variables are defined for each PFT + pft_soilmem_loop: do ipft=1,numpft + + ! Update soil moisture information memory (we always track the last 10 days). + ! We shift the memory from days to the previous day, and make room for current day + do i_wmem = numWaterMem,2,-1 + currentSite%liqvol_memory(i_wmem,ipft) = currentSite%liqvol_memory(i_wmem-1,ipft) + currentSite%smp_memory (i_wmem,ipft) = currentSite%smp_memory (i_wmem-1,ipft) + end do + + ! Find the rooting depth distribution for PFT + call set_root_fraction( currentSite%rootfrac_scr, ipft, 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 + + ! Set the memory to be the weighted average of the liquid volumetric soil water, + ! using the root fraction of each layer (except the topmost one) as the weighting + ! factor. + currentSite%liqvol_memory(1,ipft) = & + sum( bc_in%h2o_liqvol_sl(2:nlevroot) * currentSite%rootfrac_scr(2:nlevroot) ) / & + rootfrac_notop + + ! For the soil matric potential, we must check whether or not to include the layer + ! based on the layer temperature. + currentSite%smp_memory (1,ipft) = 0._r8 + root_loop: do j = 2,nlevroot + if ( check_layer_water(bc_in%h2o_liqvol_sl(j),bc_in%tempk_sl(j)) ) then + currentSite%smp_memory (1,ipft) = & + currentSite%smp_memory (1,ipft) + & + bc_in%smp_sl(j) * currentSite%rootfrac_scr(j) / rootfrac_notop + else + ! Nominal extreme suction for frozen or unreasonably dry soil + currentSite%smp_memory (1,ipft) = & + currentSite%smp_memory (1,ipft) + & + smp_lwr_bound * currentSite%rootfrac_scr(j) / rootfrac_notop + end if + end do root_loop + end do pft_soilmem_loop + + + end subroutine UpdateMemoryMoisture +end module FatesCumulativeMemoryMod diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 86d3bbeddb..09f5363d71 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -56,7 +56,6 @@ module FatesHistoryInterfaceMod use FatesInterfaceTypesMod , only : hlm_parteh_mode use FatesInterfaceTypesMod , only : hlm_use_sp use EDParamsMod , only : ED_val_comp_excln - use EDParamsMod , only : ED_val_phen_coldtemp use EDParamsMod , only : nlevleaf use EDParamsMod , only : ED_val_history_height_bin_edges use EDParamsMod , only : ED_val_history_ageclass_bin_edges @@ -438,13 +437,6 @@ module FatesHistoryInterfaceMod integer :: ih_h2oveg_hydro_err_si integer :: ih_lai_si integer :: ih_elai_si - - integer :: ih_site_cstatus_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_nesterov_fire_danger_si integer :: ih_fire_nignitions_si @@ -627,9 +619,14 @@ module FatesHistoryInterfaceMod integer :: ih_crownarea_si_cnlf 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_site_phen_status_si_pft + integer :: ih_gdd_si_pft + integer :: ih_site_nchilldays_si_pft + integer :: ih_site_ncolddays_si_pft + integer :: ih_leafoff_si_pft + integer :: ih_leafon_si_pft + integer :: ih_meanbtran24_si_pft + integer :: ih_minbtran24_si_pft integer :: ih_meanliqvol_si_pft integer :: ih_meansmp_si_pft integer :: ih_elong_factor_si_pft @@ -2472,12 +2469,6 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites,bc_in) hio_canopy_mortality_carbonflux_si => this%hvars(ih_canopy_mortality_carbonflux_si)%r81d, & hio_ustory_mortality_carbonflux_si => this%hvars(ih_understory_mortality_carbonflux_si)%r81d, & hio_woodproduct_si => this%hvars(ih_woodproduct_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_site_cstatus_si => this%hvars(ih_site_cstatus_si)%r81d, & - hio_cleafoff_si => this%hvars(ih_cleafoff_si)%r81d, & - hio_cleafon_si => this%hvars(ih_cleafon_si)%r81d, & hio_cbal_err_fates_si => this%hvars(ih_cbal_err_fates_si)%r81d, & hio_tveg24 => this%hvars(ih_tveg24_si)%r81d, & hio_tlongterm => this%hvars(ih_tlongterm_si)%r81d, & @@ -2531,20 +2522,6 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites,bc_in) ! Canopy spread index (0-1) hio_canopy_spread_si(io_si) = sites(s)%spread - ! 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) - hio_site_ncolddays_si(io_si) = real(sites(s)%ncolddays,r8) - - ! Growing degree-days - hio_gdd_si(io_si) = sites(s)%grow_deg_days - - ! Model days elapsed since leaf on/off for cold-deciduous - hio_cleafoff_si(io_si) = real(sites(s)%phen_model_date - sites(s)%cleafoffdate,r8) - hio_cleafon_si(io_si) = real(sites(s)%phen_model_date - sites(s)%cleafondate,r8) - ! track total wood product accumulation at the site level hio_woodproduct_si(io_si) = sites(s)%resources_management%trunk_product_site & * AREA_INV @@ -3263,9 +3240,15 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) hio_crownarea_cl => this%hvars(ih_crownarea_cl)%r82d) ! Break up associates for NAG compilers - associate( 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, & + associate( & + hio_gdd_si_pft => this%hvars(ih_gdd_si_pft)%r82d, & + hio_site_ncolddays_si_pft => this%hvars(ih_site_ncolddays_si_pft)%r82d, & + hio_site_nchilldays_si_pft => this%hvars(ih_site_nchilldays_si_pft)%r82d, & + hio_site_phen_status_si_pft => this%hvars(ih_site_phen_status_si_pft)%r82d, & + hio_leafoff_si_pft => this%hvars(ih_leafoff_si_pft)%r82d, & + hio_leafon_si_pft => this%hvars(ih_leafon_si_pft)%r82d, & + hio_meanbtran24_si_pft => this%hvars(ih_meanbtran24_si_pft)%r82d, & + hio_minbtran24_si_pft => this%hvars(ih_minbtran24_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, & @@ -3342,14 +3325,22 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) end do - ! Update drought deciduous information (now separated by PFT). + ! Update deciduous information (now separated by PFT). do ft = 1,numpft - ! Update the site-PFT status for drought deciduous - hio_site_dstatus_si_pft(io_si,ft) = real(sites(s)%dstatus(ft),r8) - ! Model days elapsed since leaf off/on for drought deciduous - hio_dleafoff_si_pft(io_si,ft) = real(sites(s)%dndaysleafon (ft),r8) - hio_dleafon_si_pft(io_si,ft) = real(sites(s)%dndaysleafoff(ft),r8) + ! Update the site phenological status for this PFT + hio_site_phen_status_si_pft(io_si,ft) = real(sites(s)%phen_status(ft),r8) + + ! Number of chill days and cold days + hio_site_nchilldays_si_pft(io_si,ft) = real(sites(s)%nchilldays(ft),r8) + hio_site_ncolddays_si_pft(io_si,ft) = real(sites(s)%ncolddays(ft),r8) + + ! Growing degree-days + hio_gdd_si_pft(io_si,ft) = sites(s)%grow_deg_days(ft) + + ! Model days elapsed since leaf on/off for deciduous + hio_leafoff_si_pft(io_si,ft) = real(sites(s)%ndaysleafoff(ft),r8) + hio_leafon_si_pft(io_si,ft) = real(sites(s)%ndaysleafon (ft),r8) ! Leaf elongation factor (0 means fully abscissed, 1 means fully flushed). hio_elong_factor_si_pft(io_si,ft) = sites(s)%elong_factor(ft) @@ -3422,6 +3413,13 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) cpatch%frac_burnt * cpatch%area * AREA_INV / sec_per_day endif + ! 24hr mean btran + hio_meanbtran24_si_pft(io_si,ft) = hio_meanbtran24_si_pft(io_si,ft) + & + cpatch%btran24_ft(ft)%p%GetMean() * cpatch%area * AREA_INV + + ! 24hr minimum btran + hio_minbtran24_si_pft(io_si,ft) = hio_minbtran24_si_pft(io_si,ft) + & + cpatch%btran24_ft(ft)%p%GetMin() * cpatch%area * AREA_INV end do ! loop through cohorts on patch @@ -6290,41 +6288,6 @@ subroutine define_history_vars(this, initialize_variables) upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & index=ih_ca_weighted_height_si) - call this%set_history_var(vname='FATES_COLD_STATUS', units='', & - long='site-level cold status, 0=not cold-dec, 1=too cold for leaves, 2=not too cold', & - use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & - upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & - index=ih_site_cstatus_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', & - upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, index=ih_gdd_si) - - call this%set_history_var(vname='FATES_NCHILLDAYS', units = 'days', & - long='site-level number of chill days', use_default='active', & - avgflag='A', vtype=site_r8, hlms='CLM:ALM', & - upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & - index=ih_site_nchilldays_si) - - call this%set_history_var(vname='FATES_NCOLDDAYS', units = 'days', & - long='site-level number of cold days', use_default='active', & - avgflag='A', vtype=site_r8, hlms='CLM:ALM', & - upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & - index=ih_site_ncolddays_si) - - call this%set_history_var(vname='FATES_DAYSINCE_COLDLEAFOFF', & - units='days', long='site-level days elapsed since cold leaf drop', & - use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & - upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & - index=ih_cleafoff_si) - - call this%set_history_var(vname='FATES_DAYSINCE_COLDLEAFON', & - units='days', long='site-level days elapsed since cold leaf flush', & - use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & - upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & - index=ih_cleafon_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', & @@ -7065,26 +7028,57 @@ subroutine define_history_vars(this, initialize_variables) upfreq=group_dyna_complx, 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=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & - index=ih_site_dstatus_si_pft) + !MLO - Phenology variables are now defined for each PFT. - call this%set_history_var(vname='FATES_DAYSINCE_DROUGHTLEAFOFF_PF', & - units='days', long='PFT-level days elapsed since drought leaf drop', & + call this%set_history_var(vname='FATES_PHEN_STATUS_PF', units='', & + long='Phenology status: 0 evergreen; 1 never cold; 2 cold state; 3 warm state; 4 forced abscission;'// & + ' 5 dry state; 6 moist state; 7 forced flush; 8 partial abscission', & use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_site_phen_status_si_pft) + + call this%set_history_var(vname='FATES_GDD_PF', units='degree_Celsius', & + long='PFT-level growing degree days', use_default='active', & + avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, index=ih_gdd_si_pft) + + call this%set_history_var(vname='FATES_NCHILLDAYS_PF', units = 'days', & + long='PFT-level number of chill days', use_default='active', & + avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_site_nchilldays_si_pft) + + call this%set_history_var(vname='FATES_NCOLDDAYS_PF', units = 'days', & + long='PFT-level number of cold days', use_default='active', & + avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_site_ncolddays_si_pft) + + call this%set_history_var(vname='FATES_DAYSINCE_LEAFOFF_PF', & + units='days', long='PFT-level days elapsed since last leaf abscission event', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_leafoff_si_pft) + + call this%set_history_var(vname='FATES_DAYSINCE_LEAFON_PF', & + units='days', long='PFT-level days elapsed since last leaf flushing event', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_leafon_si_pft) + + call this%set_history_var(vname='FATES_MEANBTRAN24_PF', & + units='1', & + long='PFT-level 24hr mean transpiration wetness factor (btran)', & + use_default='inactive', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & - index=ih_dleafoff_si_pft) + index=ih_meanbtran24_si_pft) - call this%set_history_var(vname='FATES_DAYSINCE_DROUGHTLEAFON_PF', & - units='days', & - long='PFT-level days elapsed since drought leaf flush', & + call this%set_history_var(vname='FATES_MINBTRAN24_PF', & + units='1', & + long='PFT-level 24hr minimum transpiration wetness factor (btran)', & use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & - index=ih_dleafon_si_pft) + index=ih_minbtran24_si_pft) call this%set_history_var(vname='FATES_MEANLIQVOL_DROUGHTPHEN_PF', & units='m3 m-3', & diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 555fa526c1..614f722938 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -99,15 +99,15 @@ module FatesInterfaceMod use PRTInitParamsFatesMod , only : PRTCheckParams, PRTDerivedParams use PRTAllometricCarbonMod , only : InitPRTGlobalAllometricCarbon use PRTAllometricCNPMod , only : InitPRTGlobalAllometricCNP - use FatesRunningMeanMod , only : ema_24hr - use FatesRunningMeanMod , only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par - use FatesRunningMeanMod , only : ema_sdlng_mdd, ema_sdlng2sap_par - use FatesRunningMeanMod , only : fixed_24hr - use FatesRunningMeanMod , only : ema_lpa - use FatesRunningMeanMod , only : ema_longterm - use FatesRunningMeanMod , only : ema_60day - use FatesRunningMeanMod , only : moving_ema_window - use FatesRunningMeanMod , only : fixed_window + use FatesRunningSummMod , only : ema_24hr + use FatesRunningSummMod , only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par + use FatesRunningSummMod , only : ema_sdlng_mdd, ema_sdlng2sap_par + use FatesRunningSummMod , only : fixed_24hr + use FatesRunningSummMod , only : ema_lpa + use FatesRunningSummMod , only : ema_longterm + use FatesRunningSummMod , only : ema_60day + use FatesRunningSummMod , only : moving_ema_window + use FatesRunningSummMod , only : fixed_window use FatesHistoryInterfaceMod , only : fates_hist use FatesHydraulicsMemMod , only : nshell use FatesHydraulicsMemMod , only : nlevsoi_hyd_max @@ -530,6 +530,7 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats, allocate(bc_in%tgcm_pa(maxpatch_total)) allocate(bc_in%t_soisno_sl(nlevsoil_in)) + ! Canopy Radiation bc_in%coszen = nan allocate(bc_in%fcansno_pa(maxpatch_total)) @@ -1085,8 +1086,8 @@ subroutine InitTimeAveragingGlobals() !allocate(ema_60day) !call ema_60day%define(prt_params%fnrt_adapt_tscl*sec_per_day,sec_per_day,moving_ema_window) - !class(rmean_arr_type), pointer :: ema_fnrt_tscale(:) - !rmean_arr_type + !class(rsumm_arr_type), pointer :: ema_fnrt_tscale(:) + !rsumm_arr_type return @@ -2233,7 +2234,7 @@ end subroutine FatesReportParameters subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) ! In this routine, we update any FATES buffers where - ! we calculate running means. It is assumed that this buffer is updated + ! we calculate running summaries. It is assumed that this buffer is updated ! on the model time-step. type(ed_site_type), intent(inout) :: sites(:) @@ -2265,9 +2266,15 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) nocomp_bare: if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then - call cpatch%tveg24%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) - call cpatch%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) - call cpatch%tveg_longterm%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) + call cpatch%tveg24%UpdateRSumm(bc_in(s)%t_veg_pa(ifp)) + call cpatch%tveg_lpa%UpdateRSumm(bc_in(s)%t_veg_pa(ifp)) + call cpatch%tveg_longterm%UpdateRSumm(bc_in(s)%t_veg_pa(ifp)) + + ! Update btran + do pft = 1, numpft + call cpatch%btran24_ft(pft)%p%UpdateRSumm(cpatch%btran_ft(pft)) + end do + ! Update the seedling layer par running means if ( hlm_regeneration_model == TRS_regeneration ) then @@ -2285,9 +2292,9 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) new_seedling_layer_par = seedling_par_high*par_high_frac + seedling_par_low*par_low_frac - call cpatch%seedling_layer_par24%UpdateRMean(new_seedling_layer_par) - call cpatch%sdlng_mort_par%UpdateRMean(new_seedling_layer_par) - call cpatch%sdlng2sap_par%UpdateRMean(new_seedling_layer_par) + call cpatch%seedling_layer_par24%UpdateRSumm(new_seedling_layer_par) + call cpatch%sdlng_mort_par%UpdateRSumm(new_seedling_layer_par) + call cpatch%sdlng2sap_par%UpdateRSumm(new_seedling_layer_par) do pft = 1,numpft @@ -2307,8 +2314,8 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) endif ! Update the seedling layer smp and mdd running means - call cpatch%sdlng_emerg_smp(pft)%p%UpdateRMean(new_seedling_layer_smp) - call cpatch%sdlng_mdd(pft)%p%UpdateRMean(new_seedling_mdd) + call cpatch%sdlng_emerg_smp(pft)%p%UpdateRSumm(new_seedling_layer_smp) + call cpatch%sdlng_mdd(pft)%p%UpdateRSumm(new_seedling_mdd) enddo !end pft loop @@ -2316,7 +2323,7 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) ccohort => cpatch%tallest do while (associated(ccohort)) - ! call ccohort%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) + ! call ccohort%tveg_lpa%UpdateRSumm(bc_in(s)%t_veg_pa(ifp)) if(.not.ccohort%isnew)then ! [kgC/plant/yr] -> [gC/m2/s] site_npp = site_npp + ccohort%npp_acc_hold * ccohort%n*area_inv * & diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index b5fe1acfe9..ec38a80678 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -28,7 +28,6 @@ module FatesInventoryInitMod ! FATES GLOBALS use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : pi_const - use FatesConstantsMod, only : itrue use FatesConstantsMod, only : nearzero use FatesGlobals , only : endrun => fates_endrun use FatesGlobals , only : fates_log @@ -45,12 +44,14 @@ module FatesInventoryInitMod use EDTypesMod , only : area use FatesConstantsMod, only : leaves_on use FatesConstantsMod, only : leaves_off + use FatesConstantsMod, only : ievergreen + use FatesConstantsMod, only : ihard_season_decid use FatesConstantsMod, only : ihard_stress_decid use FatesConstantsMod, only : isemi_stress_decid use PRTGenericMod , only : num_elements use PRTGenericMod , only : element_list - use EDTypesMod , only : phen_cstat_nevercold - use EDTypesMod , only : phen_cstat_iscold + use EDTypesMod , only : phen_cstat_timeoff + use EDTypesMod , only : phen_cstat_tempoff use EDTypesMod , only : phen_dstat_timeoff use EDTypesMod , only : phen_dstat_moistoff use PRTParametersMod , only : prt_params @@ -71,7 +72,7 @@ module FatesInventoryInitMod use PRTGenericMod, only : phosphorus_element use PRTGenericMod, only : SetState use FatesConstantsMod, only : primaryland - use FatesRunningMeanMod, only : ema_lpa + use FatesRunningSummMod, only : ema_lpa use PRTGenericMod, only : StorageNutrientTarget use FatesConstantsMod, only : fates_unset_int use EDCanopyStructureMod, only : canopy_summarization, canopy_structure @@ -970,17 +971,25 @@ subroutine set_inventory_cohort_type1(csite,bc_in,css_file_unit,npatches, & fnrt_drop_fraction = prt_params%phen_fnrt_drop_fraction(temp_cohort%pft) stem_drop_fraction = prt_params%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 - ! Cold deciduous and season is for leaves off. Set leaf status and - ! elongation factors accordingly - temp_cohort%efleaf_coh = 0.0_r8 - temp_cohort%effnrt_coh = 1._r8 - fnrt_drop_fraction - temp_cohort%efstem_coh = 1._r8 - stem_drop_fraction - - temp_cohort%status_coh = leaves_off + phen_select: select case (prt_params%phen_leaf_habit(temp_cohort%pft)) + case (ihard_season_decid) + select case(csite%phen_status(temp_cohort%pft)) + case (phen_cstat_tempoff,phen_cstat_timeoff) + ! Cold deciduous and season is for leaves off. Set leaf status and + ! elongation factors accordingly + temp_cohort%efleaf_coh = 0.0_r8 + temp_cohort%effnrt_coh = 1._r8 - fnrt_drop_fraction + temp_cohort%efstem_coh = 1._r8 - stem_drop_fraction + temp_cohort%status_coh = leaves_off + case default + ! Cold deciduous during the growing season. Assume tissues are fully flushed. + temp_cohort%efleaf_coh = 1.0_r8 + temp_cohort%effnrt_coh = 1.0_r8 + temp_cohort%efstem_coh = 1.0_r8 + temp_cohort%status_coh = leaves_on + end select - elseif ( any(prt_params%stress_decid(temp_cohort%pft) == [ihard_stress_decid,isemi_stress_decid])) then + case (ihard_stress_decid,isemi_stress_decid) ! 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. @@ -1002,14 +1011,13 @@ subroutine set_inventory_cohort_type1(csite,bc_in,css_file_unit,npatches, & ! Leaves are off (abscissing). temp_cohort%status_coh = leaves_off end if - else - ! Evergreen, or deciduous PFT during the growing season. Assume tissues are fully flushed. + case (ievergreen) + ! Evergreen. Assume tissues are fully flushed. temp_cohort%efleaf_coh = 1.0_r8 temp_cohort%effnrt_coh = 1.0_r8 temp_cohort%efstem_coh = 1.0_r8 - temp_cohort%status_coh = leaves_on - end if + end select phen_select call bagw_allom(temp_cohort%dbh,temp_cohort%pft, & temp_cohort%crowndamage, temp_cohort%efstem_coh, c_agw) @@ -1040,7 +1048,7 @@ subroutine set_inventory_cohort_type1(csite,bc_in,css_file_unit,npatches, & ! (Keeping as an example) ! Allocate running mean functions !allocate(temp_cohort%tveg_lpa) - !call temp_cohort%tveg_lpa%InitRMean(ema_lpa,init_value=cpatch%tveg_lpa%GetMean()) + !call temp_cohort%tveg_lpa%InitRSumm(ema_lpa,init_value=cpatch%tveg_lpa%GetMean()) do el = 1,num_elements diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index f2c78051a1..bf30afc641 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -51,8 +51,8 @@ module FatesRestartInterfaceMod use EDParamsMod, only : nlevleaf use PRTGenericMod, only : prt_global use PRTGenericMod, only : num_elements - use FatesRunningMeanMod, only : rmean_type - use FatesRunningMeanMod, only : ema_lpa + use FatesRunningSummMod, only : rsumm_type + use FatesRunningSummMod, only : ema_lpa use FatesRadiationMemMod, only : num_swb,norman_solver,twostr_solver use TwoStreamMLPEMod, only : normalized_upper_boundary use FatesConstantsMod, only : n_term_mort_types @@ -96,16 +96,8 @@ module FatesRestartInterfaceMod ! Indices to the restart variable object integer :: ir_npatch_si - integer :: ir_cd_status_si - integer :: ir_nchill_days_si - integer :: ir_ncold_days_si - integer :: ir_cleafondate_si - integer :: ir_cleafoffdate_si - integer :: ir_cndaysleafon_si - integer :: ir_cndaysleafoff_si integer :: ir_phenmodeldate_si integer :: ir_fireweather_index_si - integer :: ir_gdd_si integer :: ir_min_allowed_landuse_fraction_si integer :: ir_landuse_vector_gt_min_si integer :: ir_area_bareground_si @@ -224,13 +216,19 @@ module FatesRestartInterfaceMod integer :: ir_scorch_ht_pa_pft integer :: ir_litter_moisture_pa_nfsc + integer :: ir_btran24_pa_pft + ! Site level - integer :: ir_dd_status_sift - integer :: ir_dleafondate_sift - integer :: ir_dleafoffdate_sift - integer :: ir_dndaysleafon_sift - integer :: ir_dndaysleafoff_sift + integer :: ir_phen_status_sift + integer :: ir_nchill_days_sift + integer :: ir_ncold_days_sift + integer :: ir_leafondate_sift + integer :: ir_leafoffdate_sift + integer :: ir_ndaysleafon_sift + integer :: ir_ndaysleafoff_sift + integer :: ir_gdd_sift integer :: ir_elong_factor_sift + integer :: ir_btranmem_siwmft integer :: ir_liqvolmem_siwmft integer :: ir_smpmem_siwmft integer :: ir_recl2fr_sipfcl @@ -391,9 +389,9 @@ module FatesRestartInterfaceMod procedure, private :: GetCohortRealVector procedure, private :: SetCohortRealVector procedure, private :: RegisterCohortVector - procedure, private :: DefineRMeanRestartVar - procedure, private :: GetRMeanRestartVar - procedure, private :: SetRMeanRestartVar + procedure, private :: DefineRSummRestartVar + procedure, private :: GetRSummRestartVar + procedure, private :: SetRSummRestartVar end type fates_restart_interface_type @@ -683,34 +681,6 @@ subroutine define_restart_vars(this, initialize_variables) long_name='Total number of FATES patches per column', units='none', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_npatch_si ) - call this%set_restart_var(vname='fates_cold_dec_status', vtype=site_int, & - 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_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 ) - - call this%set_restart_var(vname='fates_cold_days', vtype=site_int, & - long_name='cold day counter', units='unitless', flushval = flushinvalid, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_ncold_days_si ) - - call this%set_restart_var(vname='fates_cold_leafondate', vtype=site_int, & - long_name='the model day of last cold leaf on', units='absolute integer day', flushval = flushinvalid, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_cleafondate_si ) - - call this%set_restart_var(vname='fates_cold_leafoffdate', vtype=site_int, & - long_name='the model day last cold leaf off', units='absolute integer day', flushval = flushinvalid, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_cleafoffdate_si ) - - call this%set_restart_var(vname='fates_cold_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_cold_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_phen_model_date', vtype=site_int, & long_name='integer model day used for phen timing', units='absolute integer day', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_phenmodeldate_si ) @@ -719,10 +689,6 @@ subroutine define_restart_vars(this, initialize_variables) long_name='a nesterov index accumulator', units='unitless', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fireweather_index_si ) - call this%set_restart_var(vname='fates_gdd_site', vtype=site_r8, & - long_name='growing degree days at each site', units='degC days', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_gdd_si ) - call this%set_restart_var(vname='fates_min_allowed_landuse_fraction_site', vtype=site_r8, & long_name='minimum allowed land use fraction at each site', units='fraction', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_min_allowed_landuse_fraction_si ) @@ -1055,7 +1021,10 @@ subroutine define_restart_vars(this, initialize_variables) long_name='scorch height', units='m', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_litter_moisture_pa_nfsc) end if - + + call this%DefineRSummRestartVar(vname='fates_btran24_pa_pft',vtype=cohort_r8, & + long_name='24-hour patch transpiration wetness factor', & + units='1', initialize=initialize_variables,ivar=ivar, index = ir_btran24_pa_pft) call this%RegisterCohortVector(symbol_base='fates_year_net_up', vtype=cohort_r8, & long_name_base='yearly net uptake at leaf layers', & @@ -1304,25 +1273,37 @@ subroutine define_restart_vars(this, initialize_variables) ! site x time level vars ! - 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_phen_status', vtype=site_int, & + long_name='phenology status flag', units='unitless', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_phen_status_sift ) + + 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_sift ) + + call this%set_restart_var(vname='fates_cold_days', vtype=site_int, & + long_name='cold day counter', units='unitless', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_ncold_days_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_phen_leafondate', vtype=site_int, & + long_name='the model day of last leaf flushing event', units='absolute integer day', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_leafondate_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_phen_leafoffdate', vtype=site_int, & + long_name='the model day last leaf abscission event', units='absolute integer day', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_leafoffdate_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_phen_ndaysleafon', vtype=site_int, & + long_name='number of days since last leaf flushing event', units='days', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_ndaysleafon_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_phen_ndaysleafoff', vtype=site_int, & + long_name='number of days since last leaf abscission event', units='days', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_ndaysleafoff_sift ) + + call this%set_restart_var(vname='fates_gdd_site', vtype=site_r8, & + long_name='growing degree days at each site', units='degC days', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_gdd_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, & @@ -1333,13 +1314,18 @@ subroutine define_restart_vars(this, initialize_variables) units='-', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_recl2fr_sipfcl) + call this%set_restart_var(vname='fates_btran_memory', vtype=cohort_r8, & + long_name='last 10 days of minimum transpiration wetness factor, by site x day-index x PFT', & + units='1', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_btranmem_siwmft ) + 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', & + long_name='last 10 days of volumetric soil water, by site x day-index x PFT', & units='m3/m3', flushval = flushzero, & 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', & + long_name='last 10 days of soil matric potential, by site x day-index x PFT', & units='mm', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_smpmem_siwmft ) @@ -1564,48 +1550,49 @@ subroutine define_restart_vars(this, initialize_variables) units='kg/m2/yr', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_emanpp_si) - call this%DefineRMeanRestartVar(vname='fates_tveg24patch',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_tveg24patch',vtype=cohort_r8, & long_name='24-hour patch veg temp', & units='K', initialize=initialize_variables,ivar=ivar, index = ir_tveg24_pa) - call this%DefineRMeanRestartVar(vname='fates_disturbance_rates',vtype=cohort_r8, & + call this%set_restart_var(vname='fates_disturbance_rates',vtype=cohort_r8, & long_name='disturbance rates by donor land-use type, receiver land-use type, and disturbance type', & - units='1/day', initialize=initialize_variables,ivar=ivar, index = ir_disturbance_rates_siluludi) + units='1/day', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_disturbance_rates_siluludi) if ( hlm_regeneration_model == TRS_regeneration ) then - call this%DefineRMeanRestartVar(vname='fates_seedling_layer_par24',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_seedling_layer_par24',vtype=cohort_r8, & long_name='24-hour seedling layer PAR', & units='W m2-1', initialize=initialize_variables,ivar=ivar, index = ir_seedling_layer_par24_pa) - call this%DefineRMeanRestartVar(vname='fates_sdlng_emerg_smp',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_sdlng_emerg_smp',vtype=cohort_r8, & long_name='seedling layer PAR on the seedling emergence timescale', & units='mm suction', initialize=initialize_variables,ivar=ivar, index = ir_sdlng_emerg_smp_pa) - call this%DefineRMeanRestartVar(vname='fates_sdlng_mort_par',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_sdlng_mort_par',vtype=cohort_r8, & long_name='seedling layer PAR on the seedling mortality timescale', & units='W m2-1', initialize=initialize_variables,ivar=ivar, index = ir_sdlng_mort_par_pa) - call this%DefineRMeanRestartVar(vname='fates_sdlng2sap_par',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_sdlng2sap_par',vtype=cohort_r8, & long_name='seedling layer PAR on the seedling to sapling transition timescale', & units='W m2-1', initialize=initialize_variables,ivar=ivar, index = ir_sdlng2sap_par_pa) - call this%DefineRMeanRestartVar(vname='fates_sdlng_mdd',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_sdlng_mdd',vtype=cohort_r8, & long_name='seedling moisture deficit days', & units='mm days', initialize=initialize_variables,ivar=ivar, index = ir_sdlng_mdd_pa) end if - call this%DefineRMeanRestartVar(vname='fates_tveglpapatch',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_tveglpapatch',vtype=cohort_r8, & long_name='running average (EMA) of patch veg temp for photo acclim', & units='K', initialize=initialize_variables,ivar=ivar, index = ir_tveglpa_pa) - call this%DefineRMeanRestartVar(vname='fates_tveglongtermpatch',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_tveglongtermpatch',vtype=cohort_r8, & long_name='long-term (T_home) running average (EMA) of patch veg temp for photo acclim', & units='K', initialize=initialize_variables,ivar=ivar, index = ir_tveglongterm_pa) ! (Keeping as an example) - !call this%DefineRMeanRestartVar(vname='fates_tveglpacohort',vtype=cohort_r8, & + !call this%DefineRSummRestartVar(vname='fates_tveglpacohort',vtype=cohort_r8, & ! long_name='running average (EMA) of cohort veg temp for photo acclim', & ! units='K', initialize=initialize_variables,ivar=ivar, index = ir_tveglpa_co) @@ -1624,7 +1611,7 @@ end subroutine define_restart_vars ! ===================================================================================== - subroutine DefineRMeanRestartVar(this,vname,vtype,long_name,units,initialize,ivar,index) + subroutine DefineRSummRestartVar(this,vname,vtype,long_name,units,initialize,ivar,index) class(fates_restart_interface_type) :: this character(len=*),intent(in) :: vname @@ -1647,6 +1634,26 @@ subroutine DefineRMeanRestartVar(this,vname,vtype,long_name,units,initialize,iva units=units, flushval = flushzero, & hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = dummy_index ) + call this%set_restart_var(vname= trim(vname)//'_cminimum', vtype=vtype, & + long_name=long_name//' current minimum', & + units=units, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = index ) + + call this%set_restart_var(vname= trim(vname)//'_lminimum', vtype=vtype, & + long_name=long_name//' latest minimum', & + units=units, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = dummy_index ) + + call this%set_restart_var(vname= trim(vname)//'_cmaximum', vtype=vtype, & + long_name=long_name//' current maximum', & + units=units, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = index ) + + call this%set_restart_var(vname= trim(vname)//'_lmaximum', vtype=vtype, & + long_name=long_name//' latest maximum', & + units=units, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = dummy_index ) + call this%set_restart_var(vname= trim(vname)//'_cindex', vtype=vtype, & long_name=long_name//' index', & units='index', flushval = flushzero, & @@ -1654,15 +1661,15 @@ subroutine DefineRMeanRestartVar(this,vname,vtype,long_name,units,initialize,iva return - end subroutine DefineRMeanRestartVar + end subroutine DefineRSummRestartVar ! ===================================================================================== - subroutine GetRMeanRestartVar(this, rmean_var, ir_var_index, position_index) + subroutine GetRSummRestartVar(this, rsumm_var, ir_var_index, position_index) class(fates_restart_interface_type) , intent(inout) :: this - class(rmean_type), intent(inout) :: rmean_var + class(rsumm_type), intent(inout) :: rsumm_var integer,intent(in) :: ir_var_index integer,intent(in) :: position_index @@ -1671,21 +1678,29 @@ subroutine GetRMeanRestartVar(this, rmean_var, ir_var_index, position_index) integer :: ir_pos_var ! global variable index - rmean_var%c_mean = this%rvars(ir_var_index)%r81d(position_index) + rsumm_var%c_mean = this%rvars(ir_var_index)%r81d(position_index) + + rsumm_var%l_mean = this%rvars(ir_var_index+1)%r81d(position_index) + + rsumm_var%c_minimum = this%rvars(ir_var_index+2)%r81d(position_index) - rmean_var%l_mean = this%rvars(ir_var_index+1)%r81d(position_index) + rsumm_var%l_minimum = this%rvars(ir_var_index+3)%r81d(position_index) - rmean_var%c_index = nint(this%rvars(ir_var_index+2)%r81d(position_index)) + rsumm_var%c_maximum = this%rvars(ir_var_index+4)%r81d(position_index) + + rsumm_var%l_maximum = this%rvars(ir_var_index+5)%r81d(position_index) + + rsumm_var%c_index = nint(this%rvars(ir_var_index+6)%r81d(position_index)) return - end subroutine GetRMeanRestartVar + end subroutine GetRSummRestartVar ! ======================================================================================= - subroutine SetRMeanRestartVar(this, rmean_var, ir_var_index, position_index) + subroutine SetRSummRestartVar(this, rsumm_var, ir_var_index, position_index) class(fates_restart_interface_type) , intent(inout) :: this - class(rmean_type), intent(inout) :: rmean_var + class(rsumm_type), intent(inout) :: rsumm_var integer,intent(in) :: ir_var_index integer,intent(in) :: position_index @@ -1693,14 +1708,22 @@ subroutine SetRMeanRestartVar(this, rmean_var, ir_var_index, position_index) integer :: i_pos ! vector position loop index integer :: ir_pos_var ! global variable index - this%rvars(ir_var_index)%r81d(position_index) = rmean_var%c_mean - - this%rvars(ir_var_index+1)%r81d(position_index) = rmean_var%l_mean - - this%rvars(ir_var_index+2)%r81d(position_index) = real(rmean_var%c_index,r8) + this%rvars(ir_var_index)%r81d(position_index) = rsumm_var%c_mean + + this%rvars(ir_var_index+1)%r81d(position_index) = rsumm_var%l_mean + + this%rvars(ir_var_index+2)%r81d(position_index) = rsumm_var%c_minimum + + this%rvars(ir_var_index+3)%r81d(position_index) = rsumm_var%l_minimum + + this%rvars(ir_var_index+4)%r81d(position_index) = rsumm_var%c_maximum + + this%rvars(ir_var_index+5)%r81d(position_index) = rsumm_var%l_maximum + + this%rvars(ir_var_index+6)%r81d(position_index) = real(rsumm_var%c_index,r8) return - end subroutine SetRMeanRestartVar + end subroutine SetRSummRestartVar ! ===================================================================================== @@ -2104,16 +2127,8 @@ 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_nchill_days_si => this%rvars(ir_nchill_days_si)%int1d, & - rio_ncold_days_si => this%rvars(ir_ncold_days_si)%int1d, & - rio_cleafondate_si => this%rvars(ir_cleafondate_si)%int1d, & - rio_cleafoffdate_si => this%rvars(ir_cleafoffdate_si)%int1d, & - rio_cndaysleafon_si => this%rvars(ir_cndaysleafon_si)%int1d, & - rio_cndaysleafoff_si => this%rvars(ir_cndaysleafoff_si)%int1d, & rio_phenmodeldate_si => this%rvars(ir_phenmodeldate_si)%int1d, & rio_fireweather_index_si => this%rvars(ir_fireweather_index_si)%r81d, & - rio_gdd_si => this%rvars(ir_gdd_si)%r81d, & rio_min_allowed_landuse_fraction_si => this%rvars(ir_min_allowed_landuse_fraction_si)%r81d, & rio_landuse_vector_gt_min_si => this%rvars(ir_landuse_vector_gt_min_si)%int1d, & rio_area_bareground_si => this%rvars(ir_area_bareground_si)%r81d, & @@ -2170,12 +2185,16 @@ 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_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_phen_status_sift => this%rvars(ir_phen_status_sift)%int1d, & + rio_nchill_days_sift => this%rvars(ir_nchill_days_sift)%int1d, & + rio_ncold_days_sift => this%rvars(ir_ncold_days_sift)%int1d, & + rio_leafondate_sift => this%rvars(ir_leafondate_sift)%int1d, & + rio_leafoffdate_sift => this%rvars(ir_leafoffdate_sift)%int1d, & + rio_ndaysleafon_sift => this%rvars(ir_ndaysleafon_sift)%int1d, & + rio_ndaysleafoff_sift => this%rvars(ir_ndaysleafoff_sift)%int1d, & + rio_gdd_sift => this%rvars(ir_gdd_sift)%r81d, & rio_elong_factor_sift => this%rvars(ir_elong_factor_sift)%r81d, & + rio_btranmem_siwmft => this%rvars(ir_btranmem_siwmft)%r81d, & rio_liqvolmem_siwmft => this%rvars(ir_liqvolmem_siwmft)%r81d, & rio_smpmem_siwmft => this%rvars(ir_smpmem_siwmft)%r81d, & rio_recl2fr_sipfcl => this%rvars(ir_recl2fr_sipfcl)%r81d, & @@ -2324,14 +2343,19 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_fmortcflux_cano_sipft(io_idx_si_pft) = sites(s)%fmort_carbonflux_canopy(i_pft) rio_fmortcflux_usto_sipft(io_idx_si_pft) = sites(s)%fmort_carbonflux_ustory(i_pft) rio_imortcflux_sipft(io_idx_si_pft) = sites(s)%imort_carbonflux(i_pft) - 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_phen_status_sift(io_idx_si_pft) = sites(s)%phen_status(i_pft) + rio_nchill_days_sift(io_idx_si_pft) = sites(s)%nchilldays(i_pft) + rio_ncold_days_sift(io_idx_si_pft) = sites(s)%ncolddays(i_pft) + rio_leafondate_sift(io_idx_si_pft) = sites(s)%leafondate(i_pft) + rio_leafoffdate_sift(io_idx_si_pft) = sites(s)%leafoffdate(i_pft) + rio_ndaysleafon_sift(io_idx_si_pft) = sites(s)%ndaysleafon(i_pft) + rio_ndaysleafoff_sift(io_idx_si_pft) = sites(s)%ndaysleafoff(i_pft) + rio_gdd_sift(io_idx_si_pft) = sites(s)%grow_deg_days(i_pft) rio_elong_factor_sift(io_idx_si_pft) = sites(s)%elong_factor(i_pft) rio_seed_in_sift(io_idx_si_pft) = sites(s)%seed_in(i_pft) - rio_seed_out_sift(io_idx_si_pft) = sites(s)%seed_out(i_pft) + rio_seed_out_sift(io_idx_si_pft) = sites(s)%seed_out(i_pft) + io_idx_si_pft = io_idx_si_pft + 1 end do @@ -2533,7 +2557,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) endif ! (Keeping as an example) - ! call this%SetRMeanRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) + ! call this%SetRSummRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) io_idx_co = io_idx_co + 1 @@ -2551,19 +2575,25 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_nocomp_pft_label_pa(io_idx_co_1st)= cpatch%nocomp_pft_label rio_area_pa(io_idx_co_1st) = cpatch%area - ! Patch level running means - call this%SetRMeanRestartVar(cpatch%tveg24, ir_tveg24_pa, io_idx_co_1st) - call this%SetRMeanRestartVar(cpatch%tveg_lpa, ir_tveglpa_pa, io_idx_co_1st) - call this%SetRMeanRestartVar(cpatch%tveg_longterm, ir_tveglongterm_pa, io_idx_co_1st) + ! Patch level running summaries + call this%SetRSummRestartVar(cpatch%tveg24, ir_tveg24_pa, io_idx_co_1st) + call this%SetRSummRestartVar(cpatch%tveg_lpa, ir_tveglpa_pa, io_idx_co_1st) + call this%SetRSummRestartVar(cpatch%tveg_longterm, ir_tveglongterm_pa, io_idx_co_1st) + + io_idx_pa_pft = io_idx_co_1st + do i_pft = 1, numpft + call this%SetRSummRestartVar(cpatch%btran24_ft(i_pft)%p, ir_btran24_pa_pft,io_idx_pa_pft) + io_idx_pa_pft = io_idx_pa_pft + 1 + end do if ( hlm_regeneration_model == TRS_regeneration ) then - call this%SetRMeanRestartVar(cpatch%seedling_layer_par24, ir_seedling_layer_par24_pa, io_idx_co_1st) - call this%SetRMeanRestartVar(cpatch%sdlng_mort_par, ir_sdlng_mort_par_pa,io_idx_co_1st) - call this%SetRMeanRestartVar(cpatch%sdlng2sap_par, ir_sdlng2sap_par_pa,io_idx_co_1st) + call this%SetRSummRestartVar(cpatch%seedling_layer_par24, ir_seedling_layer_par24_pa, io_idx_co_1st) + call this%SetRSummRestartVar(cpatch%sdlng_mort_par, ir_sdlng_mort_par_pa,io_idx_co_1st) + call this%SetRSummRestartVar(cpatch%sdlng2sap_par, ir_sdlng2sap_par_pa,io_idx_co_1st) io_idx_pa_pft = io_idx_co_1st do i_pft = 1, numpft - call this%SetRMeanRestartVar(cpatch%sdlng_mdd(i_pft)%p, ir_sdlng_mdd_pa,io_idx_pa_pft) - call this%SetRMeanRestartVar(cpatch%sdlng_emerg_smp(i_pft)%p, ir_sdlng_emerg_smp_pa,io_idx_pa_pft) + call this%SetRSummRestartVar(cpatch%sdlng_mdd(i_pft)%p, ir_sdlng_mdd_pa,io_idx_pa_pft) + call this%SetRSummRestartVar(cpatch%sdlng_emerg_smp(i_pft)%p, ir_sdlng_emerg_smp_pa,io_idx_pa_pft) io_idx_pa_pft = io_idx_pa_pft + 1 enddo end if @@ -2728,14 +2758,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_fmortcarea_cano_si(io_idx_si) = sites(s)%fmort_crownarea_canopy rio_fmortcarea_usto_si(io_idx_si) = sites(s)%fmort_crownarea_ustory - rio_cd_status_si(io_idx_si) = sites(s)%cstatus - rio_nchill_days_si(io_idx_si) = sites(s)%nchilldays - rio_ncold_days_si(io_idx_si) = sites(s)%ncolddays - rio_cleafondate_si(io_idx_si) = sites(s)%cleafondate - rio_cleafoffdate_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 - rio_gdd_si(io_idx_si) = sites(s)%grow_deg_days rio_phenmodeldate_si(io_idx_si) = sites(s)%phen_model_date rio_solar_zenith_angle(io_idx_si) = sites(s)%coszen @@ -2762,6 +2784,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) do i = 1,numWaterMem ! numWaterMem currently 10 do i_pft=1,numpft + rio_btranmem_siwmft ( io_idx_si_wmem ) = sites(s)%btran_memory(i,i_pft) 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 @@ -2960,7 +2983,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in, bc_out) ! (Keeping as an example) ! Allocate running mean functions !allocate(new_cohort%tveg_lpa) - !call new_cohort%tveg_lpa%InitRMean(ema_lpa) + !call new_cohort%tveg_lpa%InitRSumm(ema_lpa) ! Update the previous @@ -3102,16 +3125,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: i_term_type ! loop counter for termination type associate( rio_npatch_si => this%rvars(ir_npatch_si)%int1d, & - rio_cd_status_si => this%rvars(ir_cd_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_cleafondate_si => this%rvars(ir_cleafondate_si)%int1d, & - rio_cleafoffdate_si => this%rvars(ir_cleafoffdate_si)%int1d, & - rio_cndaysleafon_si => this%rvars(ir_cndaysleafon_si)%int1d, & - rio_cndaysleafoff_si => this%rvars(ir_cndaysleafoff_si)%int1d, & rio_phenmodeldate_si => this%rvars(ir_phenmodeldate_si)%int1d, & rio_fireweather_index_si => this%rvars(ir_fireweather_index_si)%r81d, & - rio_gdd_si => this%rvars(ir_gdd_si)%r81d, & rio_min_allowed_landuse_fraction_si => this%rvars(ir_min_allowed_landuse_fraction_si)%r81d, & rio_landuse_vector_gt_min_si => this%rvars(ir_landuse_vector_gt_min_si)%int1d, & rio_area_bareground_si => this%rvars(ir_area_bareground_si)%r81d, & @@ -3168,12 +3183,16 @@ 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_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_phen_status_sift => this%rvars(ir_phen_status_sift)%int1d, & + rio_nchill_days_sift => this%rvars(ir_nchill_days_sift)%int1d, & + rio_ncold_days_sift => this%rvars(ir_ncold_days_sift)%int1d, & + rio_leafondate_sift => this%rvars(ir_leafondate_sift)%int1d, & + rio_leafoffdate_sift => this%rvars(ir_leafoffdate_sift)%int1d, & + rio_ndaysleafon_sift => this%rvars(ir_ndaysleafon_sift)%int1d, & + rio_ndaysleafoff_sift => this%rvars(ir_ndaysleafoff_sift)%int1d, & + rio_gdd_sift => this%rvars(ir_gdd_sift)%r81d, & rio_elong_factor_sift => this%rvars(ir_elong_factor_sift)%r81d, & + rio_btranmem_siwmft => this%rvars(ir_btranmem_siwmft)%r81d, & rio_liqvolmem_siwmft => this%rvars(ir_liqvolmem_siwmft)%r81d, & rio_smpmem_siwmft => this%rvars(ir_smpmem_siwmft)%r81d, & rio_recl2fr_sipfcl => this%rvars(ir_recl2fr_sipfcl)%r81d, & @@ -3309,11 +3328,14 @@ subroutine get_restart_vectors(this, nc, nsites, sites) sites(s)%fmort_carbonflux_canopy(i_pft) = rio_fmortcflux_cano_sipft(io_idx_si_pft) sites(s)%fmort_carbonflux_ustory(i_pft) = rio_fmortcflux_usto_sipft(io_idx_si_pft) sites(s)%imort_carbonflux(i_pft) = rio_imortcflux_sipft(io_idx_si_pft) - 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)%phen_status(i_pft) = rio_phen_status_sift(io_idx_si_pft) + sites(s)%nchilldays(i_pft) = rio_nchill_days_sift(io_idx_si_pft) + sites(s)%ncolddays(i_pft) = rio_ncold_days_sift(io_idx_si_pft) + sites(s)%leafondate(i_pft) = rio_leafondate_sift(io_idx_si_pft) + sites(s)%leafoffdate(i_pft) = rio_leafoffdate_sift(io_idx_si_pft) + sites(s)%ndaysleafon(i_pft) = rio_ndaysleafon_sift(io_idx_si_pft) + sites(s)%ndaysleafoff(i_pft) = rio_ndaysleafoff_sift(io_idx_si_pft) + sites(s)%grow_deg_days(i_pft) = rio_gdd_sift(io_idx_si_pft) sites(s)%elong_factor(i_pft) = rio_elong_factor_sift(io_idx_si_pft) sites(s)%seed_in(i_pft) = rio_seed_in_sift(io_idx_si_pft) sites(s)%seed_out(i_pft) = rio_seed_out_sift(io_idx_si_pft) @@ -3502,7 +3524,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end if ! (Keeping as an example) - !call this%GetRMeanRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) + !call this%GetRSummRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) ccohort%c_area = this%rvars(ir_c_area_co)%r81d(io_idx_co) ccohort%treelai = this%rvars(ir_treelai_co)%r81d(io_idx_co) @@ -3537,18 +3559,25 @@ subroutine get_restart_vectors(this, nc, nsites, sites) cpatch%ncl_p = this%rvars(ir_nclp_pa)%int1d(io_idx_co_1st) cpatch%zstar = this%rvars(ir_zstar_pa)%r81d(io_idx_co_1st) - call this%GetRMeanRestartVar(cpatch%tveg24, ir_tveg24_pa, io_idx_co_1st) - call this%GetRMeanRestartVar(cpatch%tveg_lpa, ir_tveglpa_pa, io_idx_co_1st) - call this%GetRMeanRestartVar(cpatch%tveg_longterm, ir_tveglongterm_pa, io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%tveg24, ir_tveg24_pa, io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%tveg_lpa, ir_tveglpa_pa, io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%tveg_longterm, ir_tveglongterm_pa, io_idx_co_1st) + + ! Daily summary for btran + io_idx_pa_pft = io_idx_co_1st + do pft = 1, numpft + call this%GetRSummRestartVar(cpatch%btran24_ft(pft)%p, ir_btran24_pa_pft, io_idx_pa_pft) + io_idx_pa_pft = io_idx_pa_pft + 1 + enddo if ( hlm_regeneration_model == TRS_regeneration ) then - call this%GetRMeanRestartVar(cpatch%seedling_layer_par24, ir_seedling_layer_par24_pa, io_idx_co_1st) - call this%GetRMeanRestartVar(cpatch%sdlng_mort_par, ir_sdlng_mort_par_pa,io_idx_co_1st) - call this%GetRMeanRestartVar(cpatch%sdlng2sap_par, ir_sdlng2sap_par_pa,io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%seedling_layer_par24, ir_seedling_layer_par24_pa, io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%sdlng_mort_par, ir_sdlng_mort_par_pa,io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%sdlng2sap_par, ir_sdlng2sap_par_pa,io_idx_co_1st) io_idx_pa_pft = io_idx_co_1st do pft = 1, numpft - call this%GetRMeanRestartVar(cpatch%sdlng_mdd(pft)%p, ir_sdlng_mdd_pa,io_idx_pa_pft) - call this%GetRMeanRestartVar(cpatch%sdlng_emerg_smp(pft)%p, ir_sdlng_emerg_smp_pa,io_idx_pa_pft) + call this%GetRSummRestartVar(cpatch%sdlng_mdd(pft)%p, ir_sdlng_mdd_pa,io_idx_pa_pft) + call this%GetRSummRestartVar(cpatch%sdlng_emerg_smp(pft)%p, ir_sdlng_emerg_smp_pa,io_idx_pa_pft) io_idx_pa_pft = io_idx_pa_pft + 1 enddo end if @@ -3675,6 +3704,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) do i = 1,numWaterMem do i_pft=1,numpft + sites(s)%btran_memory (i,i_pft) = rio_btranmem_siwmft ( io_idx_si_wmem ) 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 @@ -3759,16 +3789,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) sites(s)%demotion_carbonflux = rio_democflux_si(io_idx_si) sites(s)%promotion_carbonflux = rio_promcflux_si(io_idx_si) - ! Site level phenology status flags - - sites(s)%cstatus = rio_cd_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_cleafondate_si(io_idx_si) - sites(s)%cleafoffdate = rio_cleafoffdate_si(io_idx_si) - sites(s)%cndaysleafon = rio_cndaysleafon_si(io_idx_si) - sites(s)%cndaysleafoff = rio_cndaysleafoff_si(io_idx_si) - sites(s)%grow_deg_days = rio_gdd_si(io_idx_si) sites(s)%phen_model_date= rio_phenmodeldate_si(io_idx_si) sites(s)%coszen = rio_solar_zenith_angle(io_idx_si) diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningSummMod.F90 similarity index 55% rename from main/FatesRunningMeanMod.F90 rename to main/FatesRunningSummMod.F90 index 721030b974..cd75e6f7b2 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningSummMod.F90 @@ -1,8 +1,9 @@ -module FatesRunningMeanMod +module FatesRunningSummMod use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : fates_huge use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) use shr_log_mod, only : errMsg => shr_log_errMsg use FatesGlobals, only : endrun => fates_endrun @@ -29,7 +30,7 @@ module FatesRunningMeanMod ! memory period is, and if it should be zero'd ! These are globally defined on the proc. - type, public :: rmean_def_type + type, public :: rsumm_def_type real(r8) :: mem_period ! The total integration period (s) real(r8) :: up_period ! The period between updates (s) @@ -40,13 +41,13 @@ module FatesRunningMeanMod procedure :: define - end type rmean_def_type + end type rsumm_def_type ! This holds the time varying information for the mean ! which is instantiated on sites, patches, and cohorts - type, public :: rmean_type + type, public :: rsumm_type real(r8) :: c_mean ! The current mean value, if this ! is a moving window, it is the mean. @@ -59,6 +60,22 @@ module FatesRunningMeanMod ! as c_mean for moving windows, and for fixed windows ! it is the mean value when the time integration window ! last completed. + + real(r8) :: c_minimum ! The current minimum value. Depending on the time, + ! this will be a partial minimum as it needs to + ! compare with all to fully determine the minimum + + real(r8) :: l_minimum ! The latest reportable minimum value. This is the + ! minimum value for the last window period that was + ! fully assessed. + + real(r8) :: c_maximum ! The current maximum value. Depending on the time, + ! this will be a partial minimum as it needs to + ! compare with all to fully determine the maximum + + real(r8) :: l_maximum ! The latest reportable maximum value. This is the + ! maximum value for the last window period that was + ! fully assessed. integer :: c_index ! The number of values that have ! been added to the mean so far @@ -67,17 +84,19 @@ module FatesRunningMeanMod ! This points to the global structure that ! defines the nature of this mean/avg - type(rmean_def_type), pointer :: def_type + type(rsumm_def_type), pointer :: def_type contains procedure :: GetMean - procedure :: InitRMean - procedure :: UpdateRMean - procedure :: FuseRMean + procedure :: GetMin + procedure :: GetMax + procedure :: InitRSumm + procedure :: UpdateRSumm + procedure :: FuseRSumm procedure :: CopyFromDonor - end type rmean_type + end type rsumm_type logical, parameter :: debug = .true. @@ -87,37 +106,38 @@ module FatesRunningMeanMod ! Define the time methods that we want to have available to us - class(rmean_def_type), public, pointer :: ema_24hr ! Exponential moving average - 24hr window - class(rmean_def_type), public, pointer :: fixed_24hr ! Fixed, 24-hour window - class(rmean_def_type), public, pointer :: ema_lpa ! Exponential moving average - leaf photo acclimation - class(rmean_def_type), public, pointer :: ema_longterm ! Exponential moving average - long-term leaf photo acclimation - class(rmean_def_type), public, pointer :: ema_60day ! Exponential moving average, 60 day + class(rsumm_def_type), public, pointer :: ema_24hr ! Exponential moving average - 24hr window + class(rsumm_def_type), public, pointer :: fixed_24hr ! Fixed, 24-hour window + class(rsumm_def_type), public, pointer :: ema_lpa ! Exponential moving average - leaf photo acclimation + class(rsumm_def_type), public, pointer :: ema_longterm ! Exponential moving average - long-term leaf photo acclimation + class(rsumm_def_type), public, pointer :: ema_60day ! Exponential moving average, 60 day ! Updated daily - class(rmean_def_type), public, pointer :: ema_storemem ! EMA used for smoothing N/C and P/C storage - class(rmean_def_type), public, pointer :: ema_sdlng_mort_par ! EMA for seedling mort from light stress - class(rmean_def_type), public, pointer :: ema_sdlng2sap_par ! EMA for seedling to sapling transition rates + class(rsumm_def_type), public, pointer :: ema_storemem ! EMA used for smoothing N/C and P/C storage + class(rsumm_def_type), public, pointer :: ema_sdlng_mort_par ! EMA for seedling mort from light stress + class(rsumm_def_type), public, pointer :: ema_sdlng2sap_par ! EMA for seedling to sapling transition rates ! based in par - class(rmean_def_type), public, pointer :: ema_sdlng_emerg_h2o ! EMA for moisture-based seedling emergence - class(rmean_def_type), public, pointer :: ema_sdlng_mdd ! EMA for seedling moisture deficit days + class(rsumm_def_type), public, pointer :: ema_sdlng_emerg_h2o ! EMA for moisture-based seedling emergence + class(rsumm_def_type), public, pointer :: ema_sdlng_mdd ! EMA for seedling moisture deficit days ! If we want to have different running mean specs based on ! pft or other types of constants - type, public :: rmean_arr_type - class(rmean_type), pointer :: p - end type rmean_arr_type + type, public :: rsumm_arr_type + class(rsumm_type), pointer :: p + end type rsumm_arr_type contains subroutine define(this,mem_period,up_period,method) + ! This sub-routine sets parameters and method for this running summary structure - class(rmean_def_type) :: this + class(rsumm_def_type) :: this - real(r8),intent(in) :: mem_period - real(r8),intent(in) :: up_period - integer,intent(in) :: method + real(r8),intent(in) :: mem_period ! The total integration period (s) + real(r8),intent(in) :: up_period ! The period between updates (s) + integer,intent(in) :: method ! Is this a fixed or moving window? ! Check the memory and update periods if(debug) then @@ -143,42 +163,84 @@ end subroutine define ! ===================================================================================== function GetMean(this) + ! Function to retrieve the running average for this variable. - class(rmean_type) :: this + class(rsumm_type) :: this real(r8) :: GetMean - if(this%def_type%method .eq. moving_ema_window) then + select case(this%def_type%method) + case (moving_ema_window) if(this%c_index == 0 .and. debug) then write(fates_log(), *) 'attempting to get a running mean from a variable' write(fates_log(), *) 'that has not been given a value yet' call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end select GetMean = this%l_mean end function GetMean + + ! ===================================================================================== + function GetMin(this) + ! Function to retrieve the minimum value in this time window for this variable. + + class(rsumm_type) :: this + real(r8) :: GetMin + + select case(this%def_type%method) + case (moving_ema_window) + write(fates_log(), *) ' Attempting to get the minimum value from a variable' + write(fates_log(), *) 'that is integrated using exponential moving averages.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + case (fixed_window) + GetMin = this%l_minimum + end select + + + end function GetMin + ! ===================================================================================== - subroutine InitRMean(this,rmean_def,init_value,init_offset) + function GetMax(this) + ! Function to retrieve the maximum value in this time window for this variable. - class(rmean_type) :: this - type(rmean_def_type), target :: rmean_def - real(r8),intent(in),optional :: init_value - real(r8),intent(in),optional :: init_offset + class(rsumm_type) :: this + real(r8) :: GetMax + + select case(this%def_type%method) + case (moving_ema_window) + write(fates_log(), *) ' Attempting to get the maximum value from a variable' + write(fates_log(), *) 'that is integrated using exponential moving averages.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + case (fixed_window) + GetMax = this%l_maximum + end select + + end function GetMax + + ! ===================================================================================== + + subroutine InitRSumm(this,rsumm_def,init_value,init_offset) + ! Sub-routine to initialize the summary structure variables. ! If the initialization happens part-way through a fixed averaging window ! we need to account for this. The current method moves the position ! index to match the offset, and then assumes that the init_value provided ! was a constant over the offset period. - ! If the first value is offset, such that the we are a portion of the - ! way through the window, we need to account for this. + class(rsumm_type) :: this + type(rsumm_def_type), target :: rsumm_def + real(r8),intent(in),optional :: init_value + real(r8),intent(in),optional :: init_offset ! Point to the averaging type - this%def_type => rmean_def + this%def_type => rsumm_def - if(this%def_type%method .eq. fixed_window) then + select case(this%def_type%method) + case (fixed_window) if(debug) then if(.not.(present(init_offset).and.present(init_value)) )then @@ -189,11 +251,11 @@ subroutine InitRMean(this,rmean_def,init_value,init_offset) end if ! Check to see if the offset is an even increment of the update frequency - if( abs(real(nint(init_offset/rmean_def%up_period),r8)-(init_offset/rmean_def%up_period)) > nearzero ) then + if( abs(real(nint(init_offset/rsumm_def%up_period),r8)-(init_offset/rsumm_def%up_period)) > nearzero ) then write(fates_log(), *) 'when initializing a temporal mean on a fixed window' - write(fates_log(), *) 'the time offset must be an inrement of the update frequency' + write(fates_log(), *) 'the time offset must be an increment of the update frequency' write(fates_log(), *) 'offset: ',init_offset - write(fates_log(), *) 'up freq: ',rmean_def%up_period + write(fates_log(), *) 'up freq: ',rsumm_def%up_period call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -204,17 +266,26 @@ subroutine InitRMean(this,rmean_def,init_value,init_offset) end if - this%c_index = modulo(nint(init_offset/rmean_def%up_period),rmean_def%n_mem) - this%c_mean = real(this%c_index,r8)/real(rmean_def%n_mem,r8)*init_value + + ! If the first value is offset, such that the we are a portion of the + ! way through the window, we need to account for this. + this%c_index = modulo(nint(init_offset/rsumm_def%up_period),rsumm_def%n_mem) + this%c_mean = real(this%c_index,r8)/real(rsumm_def%n_mem,r8)*init_value this%l_mean = init_value - elseif(this%def_type%method .eq. moving_ema_window) then - - if(present(init_value))then + ! Minima and maxima are not weighted so they can both assigned the + ! initial value directly. + this%c_minimum = init_value + this%l_minimum = init_value + this%c_maximum = init_value + this%l_maximum = init_value + + case (moving_ema_window) + if (present(init_value)) then this%c_mean = init_value this%l_mean = init_value if(present(init_offset))then - this%c_index = min(nint(init_offset/rmean_def%up_period),rmean_def%n_mem) + this%c_index = min(nint(init_offset/rsumm_def%up_period),rsumm_def%n_mem) else this%c_index = 1 end if @@ -223,20 +294,26 @@ subroutine InitRMean(this,rmean_def,init_value,init_offset) this%l_mean = nan this%c_index = 0 end if - - end if - + + ! We do not track minima and maxima when handling exponential moving averages. + this%c_minimum = nan + this%l_minimum = nan + this%c_maximum = nan + this%l_maximum = nan + end select + return - end subroutine InitRMean + end subroutine InitRSumm ! ===================================================================================== subroutine CopyFromDonor(this, donor) + ! Sub-routine to transfer running summary information from one structure to another - class(rmean_type) :: this - class(rmean_type),intent(in) :: donor + class(rsumm_type) :: this ! Receptor structure + class(rsumm_type),intent(in) :: donor ! Donor structure if( .not.associated(this%def_type)) then write(fates_log(), *) 'Attempted to copy over running mean' @@ -246,11 +323,15 @@ subroutine CopyFromDonor(this, donor) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - this%c_mean = donor%c_mean - this%l_mean = donor%l_mean - this%c_index = donor%c_index + this%c_mean = donor%c_mean + this%l_mean = donor%l_mean + this%c_minimum = donor%c_minimum + this%l_minimum = donor%l_minimum + this%c_maximum = donor%c_maximum + this%l_maximum = donor%l_maximum + this%c_index = donor%c_index + - return end subroutine CopyFromDonor @@ -258,26 +339,28 @@ end subroutine CopyFromDonor ! ===================================================================================== - subroutine UpdateRMean(this, new_value) + subroutine UpdateRSumm(this, new_value) + ! Sub-routine that updates the running summary metrics. - class(rmean_type) :: this - real(r8) :: new_value ! The newest value added to the running mean - real(r8) :: wgt + class(rsumm_type) :: this + real(r8), intent(in) :: new_value ! The newest value added to the running mean + real(r8) :: wgt - if(this%def_type%method.eq.moving_ema_window) then + select case (this%def_type%method) + case (moving_ema_window) this%c_index = min(this%def_type%n_mem,this%c_index + 1) - + if(this%c_index==1) then this%c_mean = new_value else wgt = 1._r8/real(this%c_index,r8) this%c_mean = this%c_mean*(1._r8-wgt) + wgt*new_value end if - + this%l_mean = this%c_mean - - else + + case (fixed_window) ! If the last time we updated we had hit the ! end of the averaging memory period, and @@ -287,23 +370,29 @@ subroutine UpdateRMean(this, new_value) this%c_index = this%c_index + 1 wgt = this%def_type%up_period/this%def_type%mem_period this%c_mean = this%c_mean + new_value*wgt + this%c_minimum = min(this%c_minimum,new_value) + this%c_maximum = max(this%c_maximum,new_value) if(this%c_index == this%def_type%n_mem) then - this%l_mean = this%c_mean - this%c_mean = 0._r8 + this%l_mean = this%c_mean + this%c_mean = 0._r8 + this%l_minimum = this%c_minimum + this%c_minimum = fates_huge + this%l_maximum = this%c_maximum + this%c_maximum = -fates_huge this%c_index = 0 end if - - end if + end select return - end subroutine UpdateRmean + end subroutine UpdateRSumm ! ===================================================================================== - subroutine FuseRMean(this,donor,recip_wgt) - + subroutine FuseRSumm(this,donor,recip_wgt) + ! Sub-routine that fuses data from two running summary structures. + ! ! Rules for fusion: ! If both entities have valid means already, then you simply use the ! weight provided to combine them. @@ -313,18 +402,14 @@ subroutine FuseRMean(this,donor,recip_wgt) ! both. - class(rmean_type) :: this - class(rmean_type), pointer :: donor + class(rsumm_type) :: this + class(rsumm_type), pointer :: donor real(r8),intent(in) :: recip_wgt ! Weighting factor for recipient (0-1) - ! Unecessary - !if(this%def_type%n_mem .ne. donor%def_type%n_mem) then - ! write(fates_log(), *) 'memory size is somehow different during fusion' - ! write(fates_log(), *) 'of two running mean variables: '!,this%name,donor%name - ! call endrun(msg=errMsg(sourcefile, __LINE__)) - !end if - if(this%def_type%method .eq. fixed_window ) then + + select case (this%def_type%method) + case (fixed_window) if (this%c_index .ne. donor%c_index) then write(fates_log(), *) 'trying to fuse two fixed-window averages' write(fates_log(), *) 'that are at different points in the window?' @@ -333,21 +418,36 @@ subroutine FuseRMean(this,donor,recip_wgt) write(fates_log(), *) 'c_index', this%c_index, donor%c_index call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end select ! This last logic clause just simply prevents us from doing math ! on uninitialized values. If both are unintiailized, then ! leave the result as uninitialized - if( .not. (donor%c_index==0) ) then - - if(this%c_index==0) then - this%c_mean = donor%c_mean - this%l_mean = donor%l_mean - this%c_index = donor%c_index + if ( donor%c_index /= 0 ) then + + if ( this%c_index == 0 ) then + this%c_mean = donor%c_mean + this%l_mean = donor%l_mean + this%c_minimum = donor%c_minimum + this%l_minimum = donor%l_minimum + this%c_maximum = donor%c_maximum + this%l_maximum = donor%l_maximum + this%c_index = donor%c_index else ! Take the weighted mean between the two this%c_mean = this%c_mean*recip_wgt + donor%c_mean*(1._r8-recip_wgt) this%l_mean = this%l_mean*recip_wgt + donor%l_mean*(1._r8-recip_wgt) + + ! Minima and maxima are fused only for fixed window, otherwise they are + ! not computed and must remain nan + select case (this%def_type%method) + case (fixed_window) + this%c_minimum = this%c_minimum*recip_wgt + donor%c_minimum*(1._r8-recip_wgt) + this%l_minimum = this%l_minimum*recip_wgt + donor%l_minimum*(1._r8-recip_wgt) + this%c_maximum = this%c_maximum*recip_wgt + donor%c_maximum*(1._r8-recip_wgt) + this%l_maximum = this%l_maximum*recip_wgt + donor%l_maximum*(1._r8-recip_wgt) + end select + ! Update the index to the larger of the two this%c_index = max(this%c_index,donor%c_index) end if @@ -355,7 +455,7 @@ subroutine FuseRMean(this,donor,recip_wgt) end if return - end subroutine FuseRMean + end subroutine FuseRSumm -end module FatesRunningMeanMod +end module FatesRunningSummMod diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 55d1d0d41c..17377b3964 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -486,36 +486,54 @@ variables: double fates_nonhydro_smpso(fates_pft) ; fates_nonhydro_smpso:units = "mm" ; fates_nonhydro_smpso:long_name = "Soil water potential at full stomatal opening" ; + double fates_phen_chilltemp(fates_pft) ; + fates_phen_chilltemp:units = "degrees C" ; + fates_phen_chilltemp:long_name = "(Upper) Temperature threshold for accumulating number of chilling days (used for leaf flushing)" ; double fates_phen_cold_size_threshold(fates_pft) ; fates_phen_cold_size_threshold:units = "cm" ; fates_phen_cold_size_threshold:long_name = "the dbh size above which will lead to phenology-related stem and leaf drop" ; + double fates_phen_coldtemp(fates_pft) ; + fates_phen_coldtemp:units = "degrees C" ; + fates_phen_coldtemp:long_name = "(Upper) daily average temperature threshold below which the day is considered cold and will count towards leaf abscission" ; double fates_phen_drought_threshold(fates_pft) ; - fates_phen_drought_threshold:units = "m3/m3 or mm" ; - fates_phen_drought_threshold:long_name = "threshold for drought phenology (or lower threshold for semi-deciduous PFTs); 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_evergreen(fates_pft) ; - fates_phen_evergreen:units = "logical flag" ; - fates_phen_evergreen:long_name = "Binary flag for evergreen leaf habit" ; + fates_phen_drought_threshold:units = "1 or mm" ; + fates_phen_drought_threshold:long_name = "threshold for drought phenology (or lower threshold for semi-deciduous PFTs); the quantity depends on the sign: if positive, the threshold is water availability (btran units). If negative, the threshold is soil matric potentical (mm)" ; double fates_phen_flush_fraction(fates_pft) ; fates_phen_flush_fraction:units = "fraction" ; fates_phen_flush_fraction:long_name = "Upon bud-burst, the maximum fraction of storage carbon used for flushing leaves" ; 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_gddtemp(fates_pft) ; + fates_phen_gddtemp:units = "degrees C" ; + fates_phen_gddtemp:long_name = "(Lower) Temperature threshold for accumulating growing degree days (used for leaf flushing)" ; + double fates_phen_gddthresh_a(fates_pft) ; + fates_phen_gddthresh_a:units = "degrees C" ; + fates_phen_gddthresh_a:long_name = "Intercept (a) of the growing degree days (GDD) threshold modulated by number of chilling days (NCD), GDD = a + b * exp(c*NCD)" ; + double fates_phen_gddthresh_b(fates_pft) ; + fates_phen_gddthresh_b:units = "degrees C" ; + fates_phen_gddthresh_b:long_name = "Slope (b) of the growing degree days (GDD) threshold modulated by number of chilling days (NCD), GDD = a + b * exp(c*NCD)" ; + double fates_phen_gddthresh_c(fates_pft) ; + fates_phen_gddthresh_c:units = "1/days" ; + fates_phen_gddthresh_c:long_name = "Exponent slope (c) of the growing degree days (GDD) threshold modulated by number of chilling days (NCD), GDD = a + b * exp(c*NCD)" ; + double fates_phen_leaf_habit(fates_pft) ; + fates_phen_leaf_habit:units = "flag" ; + fates_phen_leaf_habit:long_name = "Flag for leaf phenology habit. 1 - evergreen; 2 - season (cold) deciduous; 3 - stress (hydro) deciduous; 4 - stress (hydro) semi-deciduous" ; double fates_phen_mindaysoff(fates_pft) ; fates_phen_mindaysoff:units = "days" ; - fates_phen_mindaysoff:long_name = "day threshold compared against days since leaves abscised (shed)" ; + fates_phen_mindaysoff:long_name = "Minimum number of days that plants must remain leafless before flushing leaves again" ; + double fates_phen_mindayson(fates_pft) ; + fates_phen_mindayson:units = "days" ; + fates_phen_mindayson:long_name = "Minimum number of days that plants must remain with leaves flushed before they can abscise leaves again" ; double fates_phen_moist_threshold(fates_pft) ; - fates_phen_moist_threshold:units = "m3/m3 or mm" ; - fates_phen_moist_threshold:long_name = "upper threshold for drought phenology (only for drought semi-deciduous PFTs); 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_season_decid(fates_pft) ; - fates_phen_season_decid:units = "logical flag" ; - fates_phen_season_decid:long_name = "Binary flag for seasonal-deciduous leaf habit" ; + fates_phen_moist_threshold:units = "1 or mm" ; + fates_phen_moist_threshold:long_name = "upper threshold for drought phenology (only for drought semi-deciduous PFTs); the quantity depends on the sign: if positive, the threshold is water availability (btran units). If negative, the threshold is soil matric potentical (mm)" ; + double fates_phen_ncolddayslim(fates_pft) ; + fates_phen_ncolddayslim:units = "days" ; + fates_phen_ncolddayslim:long_name = "Minimum number of recent days below the temperature threshold that initiates temperature-driven leaf abscission" ; double fates_phen_stem_drop_fraction(fates_pft) ; fates_phen_stem_drop_fraction:units = "fraction" ; fates_phen_stem_drop_fraction:long_name = "fraction of stems to drop for non-woody species during drought/cold" ; - double fates_phen_stress_decid(fates_pft) ; - fates_phen_stress_decid:units = "logical flag" ; - fates_phen_stress_decid:long_name = "Flag for stress/drought-deciduous leaf habit. 0 - not stress deciduous; 1 - default drought deciduous (two target states only, fully flushed or fully abscised); 2 - semi-deciduous" ; double fates_prescribed_npp_canopy(fates_pft) ; fates_prescribed_npp_canopy:units = "kgC / m^2 / yr" ; fates_prescribed_npp_canopy:long_name = "NPP per unit crown area of canopy trees for prescribed physiology mode" ; @@ -873,27 +891,6 @@ variables: double fates_patch_fusion_tol ; fates_patch_fusion_tol:units = "unitless" ; fates_patch_fusion_tol:long_name = "minimum fraction in difference in profiles between patches" ; - double fates_phen_chilltemp ; - fates_phen_chilltemp:units = "degrees C" ; - fates_phen_chilltemp:long_name = "chilling day counting threshold for vegetation" ; - double fates_phen_coldtemp ; - fates_phen_coldtemp:units = "degrees C" ; - fates_phen_coldtemp:long_name = "vegetation temperature exceedance that flags a cold-day for leaf-drop" ; - double fates_phen_gddthresh_a ; - fates_phen_gddthresh_a:units = "none" ; - fates_phen_gddthresh_a:long_name = "GDD accumulation function, intercept parameter: gdd_thesh = a + b exp(c*ncd)" ; - double fates_phen_gddthresh_b ; - fates_phen_gddthresh_b:units = "none" ; - fates_phen_gddthresh_b:long_name = "GDD accumulation function, multiplier parameter: gdd_thesh = a + b exp(c*ncd)" ; - double fates_phen_gddthresh_c ; - fates_phen_gddthresh_c:units = "none" ; - fates_phen_gddthresh_c:long_name = "GDD accumulation function, exponent parameter: gdd_thesh = a + b exp(c*ncd)" ; - 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_ncolddayslim ; - fates_phen_ncolddayslim:units = "days" ; - fates_phen_ncolddayslim:long_name = "day threshold exceedance for temperature leaf-drop" ; double fates_q10_froz ; fates_q10_froz:units = "unitless" ; fates_q10_froz:long_name = "Q10 for frozen-soil respiration rates" ; @@ -1465,31 +1462,46 @@ data: fates_nonhydro_smpso = -66000, -66000, -66000, -66000, -66000, -66000, -66000, -66000, -66000, -66000, -66000, -66000, -66000, -66000 ; + fates_phen_chilltemp = 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5 ; + fates_phen_cold_size_threshold = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; - fates_phen_drought_threshold = -152957.4, -152957.4, -152957.4, -152957.4, - -152957.4, -152957.4, -152957.4, -152957.4, -152957.4, -152957.4, - -152957.4, -152957.4, -152957.4, -152957.4 ; + fates_phen_coldtemp = 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, + 7.5, 7.5, 7.5, 7.5 ; - fates_phen_evergreen = 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0 ; + fates_phen_drought_threshold = 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, + 0.2, 0.2, 0.2, 0.2, 0.2 ; fates_phen_flush_fraction = _, _, 0.5, _, 0.5, 0.5, _, 0.5, 0.5, _, 0.5, 0.5, 0.5, 0.5 ; fates_phen_fnrt_drop_fraction = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; - fates_phen_mindaysoff = 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, - 100, 100, 100, 100 ; + fates_phen_gddtemp = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; - fates_phen_moist_threshold = -122365.9, -122365.9, -122365.9, -122365.9, - -122365.9, -122365.9, -122365.9, -122365.9, -122365.9, -122365.9, - -122365.9, -122365.9, -122365.9, -122365.9 ; + fates_phen_gddthresh_a = -68, -68, -68, -68, -68, -68, -68, -68, -68, -68, + -68, -68, -68, -68 ; - fates_phen_season_decid = 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0 ; + fates_phen_gddthresh_b = 638, 638, 638, 638, 638, 638, 638, 638, 638, 638, + 638, 638, 638, 638 ; - fates_phen_stem_drop_fraction = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; + fates_phen_gddthresh_c = -0.01, -0.01, -0.01, -0.01, -0.01, -0.01, -0.01, + -0.01, -0.01, -0.01, -0.01, -0.01, -0.01, -0.01 ; + + fates_phen_leaf_habit = 1, 1, 2, 1, 3, 2, 1, 3, 2, 1, 2, 2, 3, 3 ; + + fates_phen_mindaysoff = 100, 100, 90, 100, 100, 90, 100, 100, 90, 100, + 90, 90, 100, 100 ; - fates_phen_stress_decid = 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1 ; + fates_phen_mindayson = 90, 90, 90, 90, 100, 90, 90, 100, 90, 90, 90, 90, + 100, 100 ; + + fates_phen_moist_threshold = 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, + 0.4, 0.4, 0.4, 0.4, 0.4 ; + + fates_phen_ncolddayslim = 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5 ; + + fates_phen_stem_drop_fraction = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; fates_prescribed_npp_canopy = 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4 ; @@ -1820,20 +1832,6 @@ data: fates_patch_fusion_tol = 0.05 ; - fates_phen_chilltemp = 5 ; - - fates_phen_coldtemp = 7.5 ; - - fates_phen_gddthresh_a = -68 ; - - fates_phen_gddthresh_b = 638 ; - - fates_phen_gddthresh_c = -0.01 ; - - fates_phen_mindayson = 90 ; - - fates_phen_ncolddayslim = 5 ; - fates_q10_froz = 1.5 ; fates_q10_mr = 1.5 ; diff --git a/parteh/PRTAllometricCNPMod.F90 b/parteh/PRTAllometricCNPMod.F90 index da049161a9..47d6d95443 100644 --- a/parteh/PRTAllometricCNPMod.F90 +++ b/parteh/PRTAllometricCNPMod.F90 @@ -52,7 +52,7 @@ module PRTAllometricCNPMod use FatesIntegratorsMod , only : Euler use FatesConstantsMod , only : calloc_abs_error use FatesConstantsMod , only : nearzero - use FatesConstantsMod , only : itrue + use FatesConstantsMod , only : ievergreen use FatesConstantsMod , only : fates_unset_r8 use FatesConstantsMod , only : fates_unset_int use FatesConstantsMod , only : sec_per_day @@ -1039,7 +1039,7 @@ subroutine CNPPrioritizedReplacement(this,c_gain, n_gain, p_gain, target_c) ! Also, dont allocate to replace turnover if this is not evergreen ! (this prevents accidental re-flushing on the day they drop) if( ( any(leaf_status == [leaves_off,leaves_shedding]) .or. & - (prt_params%evergreen(ipft) /= itrue) ) & + (prt_params%phen_leaf_habit(ipft) /= ievergreen) ) & .and. (i_org == leaf_organ)) cycle ! The priority code associated with this organ diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index 500140b2c8..c5e46f783e 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -60,6 +60,7 @@ module PRTAllometricCarbonMod use FatesConstantsMod , only : leaves_on use FatesConstantsMod , only : leaves_off use FatesConstantsMod , only : leaves_shedding + use FatesConstantsMod , only : ihard_season_decid use FatesConstantsMod , only : ihard_stress_decid use FatesConstantsMod , only : isemi_stress_decid @@ -433,10 +434,10 @@ subroutine DailyPRTAllometricCarbon(this,phase) elongf_fnrt = this%bc_in(ac_bc_in_id_effnrt)%rval elongf_stem = this%bc_in(ac_bc_in_id_efstem)%rval !--- Set some logical flags to simplify "if" blocks - is_hydecid_dormant = any(prt_params%stress_decid(ipft) == [ihard_stress_decid,isemi_stress_decid] ) & + is_hydecid_dormant = any( prt_params%phen_leaf_habit(ipft) == [ihard_stress_decid,isemi_stress_decid] ) & .and. any(leaf_status == [leaves_off,leaves_shedding] ) - is_deciduous = any(prt_params%stress_decid(ipft) == [ihard_stress_decid,isemi_stress_decid] ) & - .or. ( prt_params%season_decid(ipft) == itrue ) + is_deciduous = & + any( prt_params%phen_leaf_habit(ipft) == [ihard_season_decid,ihard_stress_decid,isemi_stress_decid] ) nleafage = prt_global%state_descriptor(leaf_c_id)%num_pos ! Number of leaf age class diff --git a/parteh/PRTLossFluxesMod.F90 b/parteh/PRTLossFluxesMod.F90 index a3be88b044..6cca20182f 100644 --- a/parteh/PRTLossFluxesMod.F90 +++ b/parteh/PRTLossFluxesMod.F90 @@ -24,6 +24,7 @@ module PRTLossFluxesMod use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : calloc_abs_error use FatesConstantsMod, only : itrue + use FatesConstantsMod, only : ievergreen use FatesGlobals , only : endrun => fates_endrun use FatesGlobals , only : fates_log use shr_log_mod , only : errMsg => shr_log_errMsg @@ -756,7 +757,7 @@ subroutine MaintTurnoverSimpleRetranslocation(prt,ipft,icanlayer,is_drought) ! Only evergreens have maintenance turnover (must also change trimming logic ! if we want to change this) ! ------------------------------------------------------------------------------------- - if ( leaf_long > nearzero .and. prt_params%evergreen(ipft)==itrue ) then + if ( leaf_long > nearzero .and. prt_params%phen_leaf_habit(ipft) == ievergreen ) then if(is_drought) then base_turnover(leaf_organ) = years_per_day / & diff --git a/parteh/PRTParametersMod.F90 b/parteh/PRTParametersMod.F90 index a4d6ddb3af..c1d588a743 100644 --- a/parteh/PRTParametersMod.F90 +++ b/parteh/PRTParametersMod.F90 @@ -11,23 +11,39 @@ module PRTParametersMod type,public :: prt_param_type - ! The following three PFT classes - ! are mutually exclusive - ! MLO: perhaps we should replace these three parameters with a single - ! parameter (phenology(:)) that is assigned different indices? - integer, allocatable :: stress_decid(:) ! Is the plant stress deciduous? - ! 0 - No - ! 1 - Drought "hard" deciduous (i.e., PFT - ! sheds leaves all at once when stressed) - ! 2 - Drought semi-deciduous (i.e., PFT - ! sheds leaves gradually as drought - ! conditions deteriorate) - 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 + integer, allocatable :: phen_leaf_habit(:) ! Leaf phenological habit? Current options include the following: + ! (actual values defined in FatesConstantsMod.F90) + ! - ievergreen - evergreen + ! - ihard_season_decid - obligate cold deciduous (i.e., + ! leaves will abscise and flush every winter) + ! - ihard_stress_decid - obligate drought deciduous (i.e., + ! leaves will abscie and flush at least once a year) + ! - isemi_stress_decid - drought semi-deciduous (i.e., + ! partial abscission and flushing are allowed). + + ! 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 + + ! The three parameters below are for the growing degree days + ! threshold function, which is modulated by the number of + ! chilling days (NCD): GDD_Thresh = a + b * exp (c * NCD) + real(r8), allocatable :: phen_gddthresh_a(:) ! - a parameter (intercept) + real(r8), allocatable :: phen_gddthresh_b(:) ! - b parameter (multiplier) + real(r8), allocatable :: phen_gddthresh_c(:) ! - c parameter (exponent) + + real(r8), allocatable :: phen_gddtemp(:) ! (Lower) Temperature threshold for accumulating + ! growing degree days (used for leaf flushing) + real(r8), allocatable :: phen_chilltemp(:) ! (Upper) Temperature threshold for accumulating + ! number of chilling days (used for leaf flushing) + real(r8), allocatable :: phen_coldtemp(:) ! (Upper) daily average temperature threshold below which the + ! day is considered cold and will count towards leaf abscission + real(r8), allocatable :: phen_mindayson(:) ! Minimum number of days that plants must remain with leaves + ! flushed before they can abscise leaves again. + real(r8), allocatable :: phen_mindaysoff(:) ! Minimum number of days that plants must remain leafless + ! before flushing leaves again. + real(r8), allocatable :: phen_ncolddayslim(:) ! Minimum number of recent days below the temperature threshold + ! that initiates temperature-driven leaf abscission. real(r8), allocatable :: phen_drought_threshold(:) ! For obligate (hard) drought deciduous, this is the threshold ! below which plants will abscise leaves, and ! above which plants will flush leaves. For semi-deciduous @@ -41,8 +57,6 @@ module PRTParametersMod ! are soil volumetric water content [m3/m3]. If ! negative, the values are soil matric potential [mm]. ! Ignored for non-deciduous plants. - real(r8), allocatable :: phen_doff_time(:) ! Minimum number of days that plants must remain leafless before - ! flushing leaves again. ! Growth and Turnover Parameters diff --git a/parteh/PRTParamsFATESMod.F90 b/parteh/PRTParamsFATESMod.F90 index 15eddb91b1..4876b02da6 100644 --- a/parteh/PRTParamsFATESMod.F90 +++ b/parteh/PRTParamsFATESMod.F90 @@ -4,10 +4,12 @@ module PRTInitParamsFatesMod ! the CLM/ELM module system. use FatesConstantsMod, only : r8 => fates_r8 - use FatesConstantsMod, only : itrue,ifalse + use FatesConstantsMod, only : ifalse + use FatesConstantsMod, only : itrue use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : years_per_day use FatesInterfaceTypesMod, only : hlm_parteh_mode + use FatesInterfaceTypesMod, only : hlm_use_planthydro use PRTParametersMod, only : prt_params use PRTGenericMod, only : num_organ_types use PRTGenericMod, only : leaf_organ, fnrt_organ, store_organ @@ -32,7 +34,10 @@ module PRTInitParamsFatesMod use FatesAllometryMod, only : set_root_fraction use PRTGenericMod, only : StorageNutrientTarget use EDTypesMod, only : init_recruit_trim - use FatesConstantsMod, only : ihard_stress_decid, isemi_stress_decid + use FatesConstantsMod, only : ievergreen + use FatesConstantsMod, only : ihard_season_decid + use FatesConstantsMod, only : ihard_stress_decid + use FatesConstantsMod, only : isemi_stress_decid ! ! !PUBLIC TYPES: @@ -153,23 +158,43 @@ subroutine PRTRegisterPFT(fates_params) character(len=param_string_length) :: name - name = 'fates_phen_stress_decid' + name = 'fates_phen_leaf_habit' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_phen_season_decid' + 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) + dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_phen_evergreen' + 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) + dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_phen_stem_drop_fraction' + name = 'fates_phen_gddthresh_a' 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' + name = 'fates_phen_gddthresh_b' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_phen_gddthresh_c' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_phen_gddtemp' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_phen_chilltemp' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_phen_coldtemp' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_phen_mindayson' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -177,6 +202,10 @@ 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_ncolddayslim' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_phen_drought_threshold' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -452,25 +481,11 @@ subroutine PRTReceivePFT(fates_params) real(r8), allocatable :: tmpreal(:) ! Temporary variable to hold floats ! that are converted to ints - name = 'fates_phen_stress_decid' - call fates_params%RetrieveParameterAllocate(name=name, & - data=tmpreal) - allocate(prt_params%stress_decid(size(tmpreal,dim=1))) - call ArrayNint(tmpreal,prt_params%stress_decid) - deallocate(tmpreal) - - name = 'fates_phen_season_decid' + name = 'fates_phen_leaf_habit' call fates_params%RetrieveParameterAllocate(name=name, & data=tmpreal) - allocate(prt_params%season_decid(size(tmpreal,dim=1))) - call ArrayNint(tmpreal,prt_params%season_decid) - deallocate(tmpreal) - - name = 'fates_phen_evergreen' - call fates_params%RetrieveParameterAllocate(name=name, & - data=tmpreal) - allocate(prt_params%evergreen(size(tmpreal,dim=1))) - call ArrayNint(tmpreal,prt_params%evergreen) + allocate(prt_params%phen_leaf_habit(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,prt_params%phen_leaf_habit) deallocate(tmpreal) name = 'fates_phen_stem_drop_fraction' @@ -481,9 +496,41 @@ subroutine PRTReceivePFT(fates_params) call fates_params%RetrieveParameterAllocate(name=name, & data=prt_params%phen_fnrt_drop_fraction) + name = 'fates_phen_gddthresh_a' + call fates_params%RetrieveParameterAllocate(name=name, & + data=prt_params%phen_gddthresh_a) + + name = 'fates_phen_gddthresh_b' + call fates_params%RetrieveParameterAllocate(name=name, & + data=prt_params%phen_gddthresh_b) + + name = 'fates_phen_gddthresh_c' + call fates_params%RetrieveParameterAllocate(name=name, & + data=prt_params%phen_gddthresh_c) + + name = 'fates_phen_gddtemp' + call fates_params%RetrieveParameterAllocate(name=name, & + data=prt_params%phen_gddtemp) + + name = 'fates_phen_chilltemp' + call fates_params%RetrieveParameterAllocate(name=name, & + data=prt_params%phen_chilltemp) + + name = 'fates_phen_coldtemp' + call fates_params%RetrieveParameterAllocate(name=name, & + data=prt_params%phen_coldtemp) + + name = 'fates_phen_mindayson' + call fates_params%RetrieveParameterAllocate(name=name, & + data=prt_params%phen_mindayson) + name = 'fates_phen_mindaysoff' call fates_params%RetrieveParameterAllocate(name=name, & - data=prt_params%phen_doff_time) + data=prt_params%phen_mindaysoff) + + name = 'fates_phen_ncolddayslim' + call fates_params%RetrieveParameterAllocate(name=name, & + data=prt_params%phen_ncolddayslim) name = 'fates_phen_drought_threshold' call fates_params%RetrieveParameterAllocate(name=name, & @@ -990,12 +1037,18 @@ subroutine FatesReportPFTParams(is_master) end if write(fates_log(),*) '----------- FATES PARTEH Parameters -----------------' - 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(),fmti) 'phen_leaf_habit = ',prt_params%phen_leaf_habit 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) 'phen_doff_time = ',prt_params%phen_doff_time + write(fates_log(),fmt0) 'phen_gddthresh_a = ',prt_params%phen_gddthresh_a + write(fates_log(),fmt0) 'phen_gddthresh_b = ',prt_params%phen_gddthresh_b + write(fates_log(),fmt0) 'phen_gddthresh_c = ',prt_params%phen_gddthresh_c + write(fates_log(),fmt0) 'phen_gddtemp = ',prt_params%phen_gddtemp + write(fates_log(),fmt0) 'phen_chilltemp = ',prt_params%phen_chilltemp + write(fates_log(),fmt0) 'phen_coldtemp = ',prt_params%phen_coldtemp + write(fates_log(),fmt0) 'phen_mindayson = ',prt_params%phen_mindayson + write(fates_log(),fmt0) 'phen_mindaysoff = ',prt_params%phen_mindaysoff + write(fates_log(),fmt0) 'phen_ncolddayslim = ',prt_params%phen_ncolddayslim write(fates_log(),fmt0) 'phen_drought_threshold = ',prt_params%phen_drought_threshold write(fates_log(),fmt0) 'phen_moist_threshold = ',prt_params%phen_moist_threshold write(fates_log(),fmt0) 'wood_density = ',prt_params%wood_density @@ -1083,7 +1136,7 @@ subroutine PRTDerivedParams() integer :: i, io ! generic loop index and organ loop index norgans = size(prt_params%organ_id,1) - npft = size(prt_params%evergreen,1) + npft = size(prt_params%phen_leaf_habit,1) ! Set the reverse lookup map for organs to the parameter file index allocate(prt_params%organ_param_id(num_organ_types)) @@ -1107,9 +1160,8 @@ subroutine PRTCheckParams(is_master) ! This subroutine performs logical checks on user supplied parameters. It cross ! compares various parameters and will fail if they don't make sense. ! Examples: - ! A tree can not be defined as both evergreen and deciduous. A woody plant - ! cannot have a structural biomass allometry intercept of 0, and a non-woody - ! plant (grass) can't have a non-zero intercept... + ! A woody plant cannot have a structural biomass allometry intercept of 0, and a + ! non-woody plant (grass) can't have a non-zero intercept... ! ----------------------------------------------------------------------------------- @@ -1124,11 +1176,10 @@ subroutine PRTCheckParams(is_master) integer :: iage ! leaf age class index integer :: norgans ! size of the plant organ dimension integer :: i, io ! generic loop index and organ loop index - logical :: is_evergreen ! Is the PFT evergreen - logical :: is_season_decid ! Is the PFT cold-deciduous? - logical :: is_stress_decid ! Is the PFT drought-deciduous? - logical :: is_semi_decid ! Is the PFT drought semi-deciduous? logical :: is_hmode_fine ! Did the height allometry pass the check? + logical :: risk_hf_mort_dd ! Is there a risk for drought deciduous to die of + ! hydraulic failure before abscising leaves, + ! due to parameter settings? integer :: nerror ! Count number of errors. If this is not ! zero by theend of the subroutine, stop ! the run. @@ -1137,7 +1188,7 @@ subroutine PRTCheckParams(is_master) - npft = size(prt_params%evergreen,1) + npft = size(prt_params%phen_leaf_habit,1) ! Prior to performing checks copy grperc to the ! organ dimensioned version @@ -1218,38 +1269,28 @@ subroutine PRTCheckParams(is_master) pftloop: do ipft = 1,npft - ! Check to see if evergreen, deciduous flags are mutually exclusive - ! By the way, if these are mutually exclusive, shouldn't we define a - ! single prt_params%leaf_phenology and a list of codes for the different - ! types (i.e., ievergreen, iseason_decid, istress_hard, istress_semi, etc.)? - ! ---------------------------------------------------------------------------------- - is_evergreen = prt_params%evergreen(ipft) == itrue - is_season_decid = prt_params%season_decid(ipft) == itrue - is_stress_decid = any(prt_params%stress_decid(ipft) == [ihard_stress_decid,isemi_stress_decid]) - is_semi_decid = prt_params%stress_decid(ipft) == isemi_stress_decid - - if ( ( is_evergreen .and. is_season_decid ) .or. & - ( is_evergreen .and. is_stress_decid ) .or. & - ( is_season_decid .and. is_stress_decid ) ) then - - write(fates_log(),*) '---~---' - write(fates_log(),*) 'PFT # ',ipft,' must be defined as having one of three' - write(fates_log(),*) 'phenology habits, ie, only one of the flags below should' - write(fates_log(),*) 'be different than ',ifalse - write(fates_log(),*) 'stress_decid: ',prt_params%stress_decid(ipft) - write(fates_log(),*) 'season_decid: ',prt_params%season_decid(ipft) - write(fates_log(),*) 'evergreen: ',prt_params%evergreen(ipft) - write(fates_log(),*) '---~---' - write(fates_log(),*) '' - write(fates_log(),*) '' - nerror = nerror + 1 - end if - - ! When using the the drought semi-deciduous 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). - if (is_semi_decid) then + select_phen_thresh: select case (prt_params%phen_leaf_habit(ipft)) + case (ihard_stress_decid) + if ( prt_params%phen_drought_threshold(ipft) > 1._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(),*) '---~---' + write(fates_log(),*) ' When using drought-deciduous phenology and a positive parameter,' + write(fates_log(),*) ' (water availability [btran units]) the threshold must be less' + write(fates_log(),*) ' than or equal to 1.' + write(fates_log(),*) ' PFT = ',ipft + write(fates_log(),*) ' phen_leaf_habit = ',prt_params%phen_leaf_habit(ipft) + write(fates_log(),*) ' fates_phen_drought_threshold = ',prt_params%phen_drought_threshold(ipft) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 + end if + case (isemi_stress_decid) if ( prt_params%phen_drought_threshold(ipft)*prt_params%phen_moist_threshold(ipft) < 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 @@ -1257,10 +1298,10 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) '---~---' write(fates_log(),*) ' When using drought semi-deciduous phenology,' write(fates_log(),*) ' the moist threshold must have the same sign as' - write(fates_log(),*) ' the dry threshold. Positive = soil water content [m3/m3],' + write(fates_log(),*) ' the dry threshold. Positive = water availability [btran units],' write(fates_log(),*) ' Negative = soil matric potential [mm].' write(fates_log(),*) ' PFT = ',ipft - write(fates_log(),*) ' Stress_decid = ',prt_params%stress_decid(ipft) + write(fates_log(),*) ' phen_leaf_habit = ',prt_params%phen_leaf_habit(ipft) write(fates_log(),*) ' fates_phen_drought_threshold = ',prt_params%phen_drought_threshold(ipft) write(fates_log(),*) ' fates_phen_moist_threshold = ',prt_params%phen_moist_threshold (ipft) write(fates_log(),*) '---~---' @@ -1275,7 +1316,24 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) ' By greater we mean more positive or less negative, and' write(fates_log(),*) ' they cannot be the identical.' write(fates_log(),*) ' PFT = ',ipft - write(fates_log(),*) ' Stress_decid = ',prt_params%stress_decid(ipft) + write(fates_log(),*) ' phen_leaf_habit = ',prt_params%phen_leaf_habit(ipft) + write(fates_log(),*) ' fates_phen_drought_threshold = ',prt_params%phen_drought_threshold(ipft) + write(fates_log(),*) ' fates_phen_moist_threshold = ',prt_params%phen_moist_threshold (ipft) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 + + elseif ( any([prt_params%phen_drought_threshold(ipft),prt_params%phen_moist_threshold(ipft)] > 1._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(),*) '---~---' + write(fates_log(),*) ' When using drought semi-deciduous phenology and positive parameters,' + write(fates_log(),*) ' (water availability [btran units]) both parameters must be less' + write(fates_log(),*) ' than or equal to 1.' + write(fates_log(),*) ' PFT = ',ipft + write(fates_log(),*) ' phen_leaf_habit = ',prt_params%phen_leaf_habit(ipft) write(fates_log(),*) ' fates_phen_drought_threshold = ',prt_params%phen_drought_threshold(ipft) write(fates_log(),*) ' fates_phen_moist_threshold = ',prt_params%phen_moist_threshold (ipft) write(fates_log(),*) '---~---' @@ -1283,10 +1341,11 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) '' nerror = nerror + 1 end if - end if + end select select_phen_thresh ! For all deciduous PFTs, check that abscission fractions are all bounded. - if (prt_params%evergreen(ipft) == ifalse) then + select_drop_frac: select case (prt_params%phen_leaf_habit(ipft)) + case (ihard_season_decid,ihard_stress_decid,isemi_stress_decid) ! Check if the fraction of fine roots to be actively abscised relative to leaf abscission ! is bounded between 0 and 1 (exactly 0 and 1 are acceptable). if ( ( prt_params%phen_fnrt_drop_fraction(ipft) < 0.0_r8 ) .or. & @@ -1295,7 +1354,7 @@ subroutine PRTCheckParams(is_master) 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: (should be 0):',prt_params%evergreen(ipft) + write(fates_log(),*) ' phen_leaf_habit: ',prt_params%phen_leaf_habit(ipft) write(fates_log(),*) ' phen_fnrt_drop_fraction: ', prt_params%phen_fnrt_drop_fraction(ipft) write(fates_log(),*) '---~---' write(fates_log(),*) '' @@ -1326,16 +1385,101 @@ subroutine PRTCheckParams(is_master) 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 - write(fates_log(),*) ' evergreen flag: (should be 0):',prt_params%evergreen(ipft) + write(fates_log(),*) ' phen_leaf_habit: ',prt_params%phen_leaf_habit(ipft) write(fates_log(),*) ' phen_stem_drop_fraction: ', prt_params%phen_stem_drop_fraction(ipft) write(fates_log(),*) '---~---' write(fates_log(),*) '' write(fates_log(),*) '' nerror = nerror + 1 end if - end if - - + end select select_drop_frac + + + ! When using the the drought deciduous phenology, it is now deprecated to set + ! the parameters as negative, because it makes it really hard to coordinate with + ! the hydraulic failure mortality parameter. If the values are positive, we compare + ! the parameter against the hydraulic failure mortality thresholds and warn the user + ! if their parameters may lead to hydraulic failure mortality before leaf abscission. + ! For now these will be warnings only, but they may be escalated to errors in the + ! future. + select_phen_hfmort: select case (prt_params%phen_leaf_habit(ipft)) + case (ihard_stress_decid,isemi_stress_decid) + if_phen_deprecated: if (prt_params%phen_drought_threshold(ipft) < 0._r8) then + ! Negative parameter, warn the user this is now deprecated. + write(fates_log(),*) "---~---" + write(fates_log(),*) " WARNING!" + write(fates_log(),*) "---~---" + write(fates_log(),*) " The drought phenology threshold is negative (i.e., based" + write(fates_log(),*) " on soil matric potential). This setting makes coordination" + write(fates_log(),*) " with hydraulic failure mortality hard, potentially result-" + write(fates_log(),*) " ing in excessive hydraulic failure mortality for the" + write(fates_log(),*) " following PFT. This is still allowed for back compatibility," + write(fates_log(),*) " but it is deprecated. We recommend using positive values," + write(fates_log(),*) " which sets phenology in ""btran"" units." + write(fates_log(),*) " " + write(fates_log(),*) " PFT index: ",ipft + write(fates_log(),*) " phen_drought_threshold: ",prt_params%phen_drought_threshold(ipft) + select case (prt_params%phen_leaf_habit(ipft)) + case (isemi_stress_decid) + write(fates_log(),*) " phen_moist_threshold: ",prt_params%phen_moist_threshold(ipft) + end select + select case (hlm_use_planthydro) + case (ifalse) + write(fates_log(),*) " hf_sm_threshold: ",EDPftvarcon_inst%hf_sm_threshold(ipft) + case (itrue) + write(fates_log(),*) " hf_flc_threshold: ",EDPftvarcon_inst%hf_flc_threshold(ipft) + end select + write(fates_log(),*) " " + write(fates_log(),*) "---~---" + write(fates_log(),*) "" + write(fates_log(),*) "" + else + ! Check for potential parameter-driven hydraulic failure mortality for drought deciduous + ! plants. + select case (hlm_use_planthydro) + case (ifalse) + risk_hf_mort_dd = & + prt_params%phen_drought_threshold(ipft) < EDPftvarcon_inst%hf_sm_threshold(ipft) + case (itrue) + risk_hf_mort_dd = & + prt_params%phen_drought_threshold(ipft) < (1._r8 - EDPftvarcon_inst%hf_flc_threshold(ipft)) + end select + + if_phen_hfmort: if (risk_hf_mort_dd) then + ! Positive parameter, but threshold is lower than the hydraulic failure threshold. + write(fates_log(),*) "---~---" + write(fates_log(),*) " WARNING!" + write(fates_log(),*) "---~---" + write(fates_log(),*) " The drought phenology threshold requires drier conditions than" + write(fates_log(),*) " the hydraulic failure mortality. This will likely cause" + write(fates_log(),*) " deciduous trees to die before they can abscise leaves. If this" + write(fates_log(),*) " is not the intended behaviour, revise the parameters below." + write(fates_log(),*) " " + write(fates_log(),*) " PFT index: ",ipft + write(fates_log(),*) " phen_drought_threshold: ",prt_params%phen_drought_threshold(ipft) + select case (prt_params%phen_leaf_habit(ipft)) + case (isemi_stress_decid) + write(fates_log(),*) " phen_moist_threshold: ",prt_params%phen_moist_threshold(ipft) + end select + select case (hlm_use_planthydro) + case (ifalse) + write(fates_log(),*) " hf_sm_threshold: ",EDPftvarcon_inst%hf_sm_threshold(ipft) + case (itrue) + write(fates_log(),*) " hf_flc_threshold: ",EDPftvarcon_inst%hf_flc_threshold(ipft) + write(fates_log(),*) " " + write(fates_log(),*) " Note: ""hf_flc_threshold"" is fraction of LOSS of conductivity," + write(fates_log(),*) " so the higher hydraulic the number, the less likely the plant" + write(fates_log(),*) " will experience failure mortality. Parameter" + write(fates_log(),*) " ""fates_phen_drought_threshold"" is defined in ""btran""" + write(fates_log(),*) " units, which is equivalent to 1-flc." + end select + write(fates_log(),*) " " + write(fates_log(),*) "---~---" + write(fates_log(),*) "" + write(fates_log(),*) "" + end if if_phen_hfmort + end if if_phen_deprecated + end select select_phen_hfmort ! Check to see if mature and base seed allocation is greater than 1 @@ -1793,20 +1937,19 @@ subroutine PRTCheckParams(is_master) nerror = nerror + 1 end if - else - if (prt_params%evergreen(ipft) .eq. itrue) then - write(fates_log(),*) "---~---" - write(fates_log(),*) 'You specified zero leaf turnover: ' - write(fates_log(),*) 'ipft: ',ipft,' iage: ',iage - write(fates_log(),*) 'leaf_long(ipft,iage): ',prt_params%leaf_long(ipft,iage) - write(fates_log(),*) 'yet this is an evergreen PFT, and it only makes sense' - write(fates_log(),*) 'that an evergreen would have leaf maintenance turnover' - write(fates_log(),*) 'disable this error if you are ok with this' - write(fates_log(),*) "---~---" - write(fates_log(),*) '' - write(fates_log(),*) '' - nerror = nerror + 1 - end if + elseif (prt_params%phen_leaf_habit(ipft) == ievergreen) then + write(fates_log(),*) "---~---" + write(fates_log(),*) 'You specified zero leaf turnover: ' + write(fates_log(),*) 'ipft: ',ipft,' iage: ',iage + write(fates_log(),*) 'phen_leaf_habit: ',prt_params%phen_leaf_habit(ipft) + write(fates_log(),*) 'leaf_long(ipft,iage): ',prt_params%leaf_long(ipft,iage) + write(fates_log(),*) 'yet this is an evergreen PFT, and it only makes sense' + write(fates_log(),*) 'that an evergreen would have leaf maintenance turnover' + write(fates_log(),*) 'disable this error if you are ok with this' + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end do @@ -1858,21 +2001,19 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) '' nerror = nerror + 1 end if - - else - if (prt_params%evergreen(ipft) .eq. itrue) then - write(fates_log(),*) "---~---" - write(fates_log(),*) 'You specified zero root turnover: ' - write(fates_log(),*) 'ipft: ',ipft - write(fates_log(),*) 'root_long(ipft): ',prt_params%root_long(ipft) - write(fates_log(),*) 'yet this is an evergreen PFT, and it only makes sense' - write(fates_log(),*) 'that an evergreen would have root maintenance turnover' - write(fates_log(),*) 'disable this error if you are ok with this' - write(fates_log(),*) "---~---" - write(fates_log(),*) '' - write(fates_log(),*) '' - nerror = nerror + 1 - end if + elseif (prt_params%phen_leaf_habit(ipft) == ievergreen) then + write(fates_log(),*) "---~---" + write(fates_log(),*) 'You specified zero root turnover: ' + write(fates_log(),*) 'ipft: ',ipft + write(fates_log(),*) 'phen_leaf_habit: ',prt_params%phen_leaf_habit(ipft) + write(fates_log(),*) 'root_long(ipft): ',prt_params%root_long(ipft) + write(fates_log(),*) 'yet this is an evergreen PFT, and it only makes sense' + write(fates_log(),*) 'that an evergreen would have root maintenance turnover' + write(fates_log(),*) 'disable this error if you are ok with this' + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if ! Check Branch turnover doesn't exceed one day diff --git a/testing/testing_shr/FatesFactoryMod.F90 b/testing/testing_shr/FatesFactoryMod.F90 index aede3dbab9..f5533f4776 100644 --- a/testing/testing_shr/FatesFactoryMod.F90 +++ b/testing/testing_shr/FatesFactoryMod.F90 @@ -2,7 +2,8 @@ module FatesFactoryMod use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : leaves_on, leaves_off - use FatesConstantsMod, only : itrue + use FatesConstantsMod, only : ievergreen + use FatesConstantsMod, only : ihard_season_decid use FatesConstantsMod, only : ihard_stress_decid use FatesConstantsMod, only : isemi_stress_decid use FatesConstantsMod, only : primaryland @@ -21,8 +22,8 @@ module FatesFactoryMod use EDParamsMod, only : nclmax use EDParamsMod, only : photo_temp_acclim_timescale use EDParamsMod, only : photo_temp_acclim_thome_time - use FatesRunningMeanMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm - use FatesRunningMeanMod, only : moving_ema_window, fixed_window + use FatesRunningSummMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm + use FatesRunningSummMod, only : moving_ema_window, fixed_window use EDCohortDynamicsMod, only : InitPRTObject use PRTParametersMod, only : prt_params use PRTGenericMod, only : element_pos @@ -307,27 +308,34 @@ subroutine CohortFactory(cohort, pft, can_lai, dbh, number, crown_damage, status end if ! set leaf elongation factors - if (prt_params%season_decid(pft) == itrue .and. status_local == leaves_off) then - elongf_leaf = 0.0_r8 - elongf_fnrt = 1.0_r8 - prt_params%phen_fnrt_drop_fraction(pft) - elongf_stem = 1.0_r8 - prt_params%phen_stem_drop_fraction(pft) - - else if (any(prt_params%stress_decid(pft) == [ihard_stress_decid, isemi_stress_decid])) then + phen_select: select case (prt_params%phen_leaf_habit(pft)) + case (ihard_season_decid) + if (status_local == leaves_off) then + elongf_leaf = 0.0_r8 + elongf_fnrt = 1.0_r8 - prt_params%phen_fnrt_drop_fraction(pft) + elongf_stem = 1.0_r8 - prt_params%phen_stem_drop_fraction(pft) + else + elongf_leaf = 1.0_r8 + elongf_fnrt = 1.0_r8 + elongf_stem = 1.0_r8 + end if + + case (ihard_stress_decid, isemi_stress_decid) elongf_leaf = elong_fact_local elongf_fnrt = 1.0_r8 - (1.0_r8 - elongf_leaf)*prt_params%phen_fnrt_drop_fraction(pft) elongf_stem = 1.0_r8 - (1.0_r8 - elongf_leaf)*prt_params%phen_stem_drop_fraction(pft) - - if (elongf_leaf > 0.0_r8) then + + if (elongf_leaf > 0.0_r8) then status_local = leaves_on else status_local = leaves_off end if - - else + + case (ievergreen) elongf_leaf = 1.0_r8 - elongf_fnrt = 1.0_r8 - elongf_stem = 1.0_r8 - end if + elongf_fnrt = 1.0_r8 + elongf_stem = 1.0_r8 + end select phen_select ! calculate allometric properties