diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 6d9f8cbcb5..1acc2a7835 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 @@ -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/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 54f68036cf..d82ca05867 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_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,259 +981,130 @@ 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 + !---~--- + ! 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)) / & @@ -1245,38 +1119,14 @@ subroutine phenology( currentSite, bc_in ) 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 !---~--- @@ -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/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/EDInitMod.F90 b/main/EDInitMod.F90 index 625e65fabc..c4895cbe49 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,20 +267,15 @@ 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%liqvol_memory(:,:) = nan @@ -394,24 +394,24 @@ 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) :: 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 +425,107 @@ 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 + 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)%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 +1205,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 +1238,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..79e567cef5 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -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,37 @@ 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) :: 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..72554cec91 --- /dev/null +++ b/main/FatesCumulativeMemoryMod.F90 @@ -0,0 +1,252 @@ +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 : area_inv + 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 : 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 + 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 + + + + ! The soil memory variables are defined for each PFT + pft_memory_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_memory_loop + + + end subroutine UpdateMemoryMoisture +end module FatesCumulativeMemoryMod diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 86d3bbeddb..537ab84d5f 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,12 @@ 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_meanliqvol_si_pft integer :: ih_meansmp_si_pft integer :: ih_elong_factor_si_pft @@ -2472,12 +2467,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 +2520,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 +3238,13 @@ 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_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 +3321,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) @@ -6290,41 +6277,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 +7017,43 @@ 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_complx, ivar=ivar, initialize=initialize_variables, & - index=ih_dleafoff_si_pft) - - call this%set_history_var(vname='FATES_DAYSINCE_DROUGHTLEAFON_PF', & - units='days', & - long='PFT-level days elapsed since drought leaf flush', & - use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & - upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & - index=ih_dleafon_si_pft) + 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_MEANLIQVOL_DROUGHTPHEN_PF', & units='m3 m-3', & diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index b5fe1acfe9..2f4632c1a7 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 @@ -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) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index f2c78051a1..0a86af6af4 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -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 @@ -225,11 +217,14 @@ module FatesRestartInterfaceMod integer :: ir_litter_moisture_pa_nfsc ! 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_liqvolmem_siwmft integer :: ir_smpmem_siwmft @@ -683,34 +678,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 +686,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 ) @@ -1304,25 +1267,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, & @@ -2104,16 +2079,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,11 +2137,14 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_agesinceanthrodist_pa => this%rvars(ir_agesinceanthrodist_pa)%r81d, & rio_nocomp_pft_label_pa => this%rvars(ir_nocomp_pft_label_pa)%int1d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & - rio_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_liqvolmem_siwmft => this%rvars(ir_liqvolmem_siwmft)%r81d, & rio_smpmem_siwmft => this%rvars(ir_smpmem_siwmft)%r81d, & @@ -2324,14 +2294,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 @@ -2728,14 +2703,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 @@ -3102,16 +3069,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,11 +3127,14 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_agesinceanthrodist_pa => this%rvars(ir_agesinceanthrodist_pa)%r81d, & rio_nocomp_pft_label_pa => this%rvars(ir_nocomp_pft_label_pa)%int1d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & - rio_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_liqvolmem_siwmft => this%rvars(ir_liqvolmem_siwmft)%r81d, & rio_smpmem_siwmft => this%rvars(ir_smpmem_siwmft)%r81d, & @@ -3309,11 +3271,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) @@ -3759,16 +3724,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/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 55d1d0d41c..855295af0e 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" ; 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" ; + 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,32 +1462,52 @@ 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_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_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_evergreen = 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0 ; - 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_gddtemp = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; + + fates_phen_gddthresh_a = -68, -68, -68, -68, -68, -68, -68, -68, -68, -68, + -68, -68, -68, -68 ; + + fates_phen_gddthresh_b = 638, 638, 638, 638, 638, 638, 638, 638, 638, 638, + 638, 638, 638, 638 ; + + 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_mindaysoff = 100, 100, 90, 100, 100, 90, 100, 100, 90, 100, + 90, 90, 100, 100 ; + + fates_phen_leaf_habit = 1, 1, 2, 1, 3, 2, 1, 3, 2, 1, 2, 2, 3, 3 ; + fates_phen_mindaysoff = 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100 ; + fates_phen_mindayson = 90, 90, 90, 90, 100, 90, 90, 100, 90, 90, 90, 90, + 100, 100 ; + 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_season_decid = 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0 ; + 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_phen_stress_decid = 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1 ; - 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 +1837,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..e4e62d8e24 100644 --- a/parteh/PRTParametersMod.F90 +++ b/parteh/PRTParametersMod.F90 @@ -11,23 +11,42 @@ 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 + + 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 +60,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..2742040573 100644 --- a/parteh/PRTParamsFATESMod.F90 +++ b/parteh/PRTParamsFATESMod.F90 @@ -4,7 +4,7 @@ module PRTInitParamsFatesMod ! the CLM/ELM module system. use FatesConstantsMod, only : r8 => fates_r8 - use FatesConstantsMod, only : itrue,ifalse + use FatesConstantsMod, only : itrue use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : years_per_day use FatesInterfaceTypesMod, only : hlm_parteh_mode @@ -32,7 +32,8 @@ 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 : isemi_stress_decid ! ! !PUBLIC TYPES: @@ -153,23 +154,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 +198,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 +477,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' - 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' + name = 'fates_phen_leaf_habit' 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 +492,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 +1033,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 +1132,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 +1156,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,10 +1172,6 @@ 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? integer :: nerror ! Count number of errors. If this is not ! zero by theend of the subroutine, stop @@ -1137,7 +1181,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 +1262,10 @@ 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 + if (prt_params%phen_leaf_habit(ipft) == isemi_stress_decid) then 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 @@ -1260,7 +1276,7 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) ' the dry threshold. Positive = soil water content [m3/m3],' 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 +1291,7 @@ 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(),*) '---~---' @@ -1286,7 +1302,7 @@ subroutine PRTCheckParams(is_master) end if ! For all deciduous PFTs, check that abscission fractions are all bounded. - if (prt_params%evergreen(ipft) == ifalse) then + if (prt_params%phen_leaf_habit(ipft) /= ievergreen) then ! 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 +1311,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,7 +1342,7 @@ 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(),*) '' @@ -1793,20 +1809,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 +1873,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..1f4446887c 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 @@ -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