diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index f7d4782f15..4e3ea3a29b 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -689,7 +689,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr,bc_in) !!allocate(copyc%l2fr_ema) ! Note, no need to give a starter value here, ! that will be taken care of in copy() - !!call copyc%l2fr_ema%InitRMean(ema_60day) + !!call copyc%l2fr_ema%InitRSumm(ema_60day) ! Initialize the PARTEH object and point to the ! correct boundary condition fields @@ -1157,7 +1157,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) !!allocate(copyc%l2fr_ema) ! Note, no need to give a starter value here, ! that will be taken care of in copy() - !!call copyc%l2fr_ema%InitRMean(ema_60day) + !!call copyc%l2fr_ema%InitRSumm(ema_60day) ! Initialize the PARTEH object and point to the ! correct boundary condition fields @@ -1172,7 +1172,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! (keep as an example) ! Initialize running means !allocate(copyc%tveg_lpa) - !call copyc%tveg_lpa%InitRMean(ema_lpa,& + !call copyc%tveg_lpa%InitRSumm(ema_lpa,& ! init_value=currentPatch%tveg_lpa%GetMean()) call currentCohort%Copy(copyc) !makes an identical copy... diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 6d9f8cbcb5..39c4174b74 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -173,7 +173,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, height, coage, dbh, ! (Keeping as an example) !! allocate(newCohort%tveg_lpa) - !! call newCohort%tveg_lpa%InitRMean(ema_lpa,init_value=patchptr%tveg_lpa%GetMean()) + !! call newCohort%tveg_lpa%InitRSumm(ema_lpa,init_value=patchptr%tveg_lpa%GetMean()) if (hlm_use_planthydro .eq. itrue) then diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index b1fc9af66d..872b3d9039 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -106,9 +106,9 @@ module EDPatchDynamicsMod use EDParamsMod, only : logging_event_code use EDParamsMod, only : logging_export_frac use EDParamsMod, only : maxpatches_by_landuse - use FatesRunningMeanMod, only : ema_sdlng_mdd - use FatesRunningMeanMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par, ema_sdlng2sap_par - use FatesRunningMeanMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm + use FatesRunningSummMod, only : ema_sdlng_mdd + use FatesRunningSummMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par, ema_sdlng2sap_par + use FatesRunningSummMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm use FatesRadiationMemMod, only : num_swb ! CIME globals @@ -808,7 +808,7 @@ subroutine spawn_patches( currentSite, bc_in) ! (Keeping as an example) ! Allocate running mean functions !allocate(nc%tveg_lpa) - !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=newPatch%tveg_lpa%GetMean()) + !call nc%tveg_lpa%InitRSumm(ema_lpa,init_value=newPatch%tveg_lpa%GetMean()) call nc%ZeroValues() @@ -1701,7 +1701,7 @@ subroutine split_patch(currentSite, currentPatch, new_patch, fraction_to_keep, a ! (Keeping as an example) ! Allocate running mean functions !allocate(nc%tveg_lpa) - !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) + !call nc%tveg_lpa%InitRSumm(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) call nc%ZeroValues() @@ -3183,21 +3183,26 @@ subroutine fuse_2_patches(csite, dp, rp) call endrun(msg=errMsg(sourcefile, __LINE__)) endif - ! Weighted mean of the running means - call rp%tveg24%FuseRMean(dp%tveg24,rp%area*inv_sum_area) - call rp%tveg_lpa%FuseRMean(dp%tveg_lpa,rp%area*inv_sum_area) + ! Weighted mean of the running summaries + call rp%tveg24%FuseRSumm(dp%tveg24,rp%area*inv_sum_area) + call rp%tveg_lpa%FuseRSumm(dp%tveg_lpa,rp%area*inv_sum_area) + + do pft = 1,numpft + call rp%btran24_ft(pft)%p%FuseRSumm(dp%btran24_ft(pft)%p,rp%area*inv_sum_area) + enddo + if ( hlm_regeneration_model == TRS_regeneration ) then - call rp%seedling_layer_par24%FuseRMean(dp%seedling_layer_par24,rp%area*inv_sum_area) - call rp%sdlng_mort_par%FuseRMean(dp%sdlng_mort_par,rp%area*inv_sum_area) - call rp%sdlng2sap_par%FuseRMean(dp%sdlng2sap_par,rp%area*inv_sum_area) + call rp%seedling_layer_par24%FuseRSumm(dp%seedling_layer_par24,rp%area*inv_sum_area) + call rp%sdlng_mort_par%FuseRSumm(dp%sdlng_mort_par,rp%area*inv_sum_area) + call rp%sdlng2sap_par%FuseRSumm(dp%sdlng2sap_par,rp%area*inv_sum_area) do pft = 1,numpft - call rp%sdlng_emerg_smp(pft)%p%FuseRMean(dp%sdlng_emerg_smp(pft)%p,rp%area*inv_sum_area) - call rp%sdlng_mdd(pft)%p%FuseRMean(dp%sdlng_mdd(pft)%p,rp%area*inv_sum_area) + call rp%sdlng_emerg_smp(pft)%p%FuseRSumm(dp%sdlng_emerg_smp(pft)%p,rp%area*inv_sum_area) + call rp%sdlng_mdd(pft)%p%FuseRSumm(dp%sdlng_mdd(pft)%p,rp%area*inv_sum_area) enddo end if - call rp%tveg_longterm%FuseRMean(dp%tveg_longterm,rp%area*inv_sum_area) + call rp%tveg_longterm%FuseRSumm(dp%tveg_longterm,rp%area*inv_sum_area) rp%livegrass = (dp%livegrass*dp%area + rp%livegrass*rp%area) * inv_sum_area rp%ros_front = (dp%ros_front*dp%area + rp%ros_front*rp%area) * inv_sum_area @@ -3811,6 +3816,10 @@ subroutine CopyPatchMeansTimers(dp, rp) call rp%tveg_lpa%CopyFromDonor(dp%tveg_lpa) call rp%tveg_longterm%CopyFromDonor(dp%tveg_longterm) + do ipft = 1,numpft + call rp%btran24_ft(ipft)%p%CopyFromDonor(dp%btran24_ft(ipft)%p) + end do + if ( hlm_regeneration_model == TRS_regeneration ) then call rp%seedling_layer_par24%CopyFromDonor(dp%seedling_layer_par24) call rp%sdlng_mort_par%CopyFromDonor(dp%sdlng_mort_par) diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 9b3b9ef919..9ae53b3025 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -12,7 +12,7 @@ module FatesPatchMod use FatesUtilsMod, only : check_hlm_list use FatesUtilsMod, only : check_var_real use FatesCohortMod, only : fates_cohort_type - use FatesRunningMeanMod, only : rmean_type, rmean_arr_type + use FatesRunningSummMod, only : rsumm_type, rsumm_arr_type use FatesLitterMod, only : litter_type use FatesFuelMod, only : fuel_type use PRTGenericMod, only : num_elements @@ -24,9 +24,9 @@ module FatesPatchMod use EDParamsMod, only : nlevleaf, nclmax, maxpft use FatesConstantsMod, only : n_dbh_bins, n_dist_types use FatesConstantsMod, only : t_water_freeze_k_1atm - use FatesRunningMeanMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm - use FatesRunningMeanMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par - use FatesRunningMeanMod, only : ema_sdlng2sap_par, ema_sdlng_mdd + use FatesRunningSummMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm + use FatesRunningSummMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par + use FatesRunningSummMod, only : ema_sdlng2sap_par, ema_sdlng_mdd use TwoStreamMLPEMod, only : twostream_type use FatesRadiationMemMod, only : num_swb use FatesRadiationMemMod, only : num_rad_stream_types @@ -75,23 +75,24 @@ module FatesPatchMod !--------------------------------------------------------------------------- ! RUNNING MEANS - !class(rmean_type), pointer :: t2m ! place-holder for 2m air temperature (variable window-size) - class(rmean_type), pointer :: tveg24 ! 24-hour mean vegetation temperature [K] - class(rmean_type), pointer :: tveg_lpa ! running mean of vegetation temperature at the + !class(rsumm_type), pointer :: t2m ! place-holder for 2m air temperature (variable window-size) + class(rsumm_type), pointer :: tveg24 ! 24-hour summary vegetation temperature [K] + class(rsumm_type), pointer :: tveg_lpa ! running mean of vegetation temperature at the ! leaf photosynthesis acclimation timescale [K] - class(rmean_type), pointer :: tveg_longterm ! long-term running mean of vegetation temperature at the + class(rsumm_type), pointer :: tveg_longterm ! long-term running mean of vegetation temperature at the ! leaf photosynthesis acclimation timescale [K] (i.e T_home) - class(rmean_type), pointer :: seedling_layer_par24 ! 24-hour mean of photosynthetically active radiation at seedling layer [W/m2] - class(rmean_arr_type), pointer :: sdlng_emerg_smp(:) ! running mean of soil matric potential at the seedling + class(rsumm_arr_type), pointer :: btran24_ft(:) ! 24-hour summary of transpiration wetness factor (aka btran) + class(rsumm_type), pointer :: seedling_layer_par24 ! 24-hour summary of photosynthetically active radiation at seedling layer [W/m2] + class(rsumm_arr_type), pointer :: sdlng_emerg_smp(:) ! running mean of soil matric potential at the seedling ! rooting depth at the H2O seedling emergence timescale (see sdlng_emerg_h2o_timescale parameter) - class(rmean_type), pointer :: sdlng_mort_par ! running mean of photosythetically active radiation + class(rsumm_type), pointer :: sdlng_mort_par ! running mean of photosythetically active radiation ! at the seedling layer and at the par-based seedling ! mortality timescale (sdlng_mort_par_timescale) - class(rmean_arr_type), pointer :: sdlng_mdd(:) ! running mean of moisture deficit days + class(rsumm_arr_type), pointer :: sdlng_mdd(:) ! running mean of moisture deficit days ! at the seedling layer and at the mdd-based seedling ! mortality timescale (sdlng_mdd_timescale) ! (sdlng2sap_par_timescale) - class(rmean_type), pointer :: sdlng2sap_par ! running mean of photosythetically active radiation + class(rsumm_type), pointer :: sdlng2sap_par ! running mean of photosythetically active radiation ! at the seedling layer and at the par-based seedling ! to sapling transition timescale ! (sdlng2sap_par_timescale) @@ -614,18 +615,28 @@ subroutine InitRunningMeans(this, current_tod, regeneration_model, numpft) ! Until bc's are pointed to by sites give veg a default temp [K] real(r8), parameter :: temp_init_veg = 15._r8 + t_water_freeze_k_1atm real(r8), parameter :: init_seedling_par = 5.0_r8 ! arbitrary initialization for seedling layer [MJ m-2 d-1] - real(r8), parameter :: init_seedling_smp = -26652.0_r8 ! abitrary initialization of smp [mm] + real(r8), parameter :: init_seedling_smp = -26652.0_r8 ! arbitrary initialization of smp [mm] + real(r8), parameter :: init_btran = 1.0_r8 ! arbitrary initial value for btran [fraction] integer :: pft ! pft looping index allocate(this%tveg24) allocate(this%tveg_lpa) allocate(this%tveg_longterm) + allocate(this%btran24_ft(numpft)) - ! set initial values for running means - call this%tveg24%InitRMean(fixed_24hr, init_value=temp_init_veg, & + ! set initial values for running summaries + call this%tveg24%InitRSumm(fixed_24hr, init_value=temp_init_veg, & init_offset=real(current_tod, r8)) - call this%tveg_lpa%InitRmean(ema_lpa, init_value=temp_init_veg) - call this%tveg_longterm%InitRmean(ema_longterm, init_value=temp_init_veg) + call this%tveg_lpa%InitRSumm(ema_lpa, init_value=temp_init_veg) + call this%tveg_longterm%InitRSumm(ema_longterm, init_value=temp_init_veg) + + do pft = 1,numpft + allocate(this%btran24_ft(pft)%p) + + call this%btran24_ft(pft)%p%InitRSumm(fixed_24hr,init_value=init_btran, & + init_offset=real(current_tod, r8)) + end do + if (regeneration_model == TRS_regeneration) then allocate(this%seedling_layer_par24) @@ -634,20 +645,20 @@ subroutine InitRunningMeans(this, current_tod, regeneration_model, numpft) allocate(this%sdlng_mort_par) allocate(this%sdlng2sap_par) - call this%seedling_layer_par24%InitRMean(fixed_24hr, & + call this%seedling_layer_par24%InitRSumm(fixed_24hr, & init_value=init_seedling_par, init_offset=real(current_tod, r8)) - call this%sdlng_mort_par%InitRMean(ema_sdlng_mort_par, & + call this%sdlng_mort_par%InitRSumm(ema_sdlng_mort_par, & init_value=temp_init_veg) - call this%sdlng2sap_par%InitRMean(ema_sdlng2sap_par, & + call this%sdlng2sap_par%InitRSumm(ema_sdlng2sap_par, & init_value=init_seedling_par) do pft = 1,numpft allocate(this%sdlng_mdd(pft)%p) allocate(this%sdlng_emerg_smp(pft)%p) - call this%sdlng_mdd(pft)%p%InitRMean(ema_sdlng_mdd, & + call this%sdlng_mdd(pft)%p%InitRSumm(ema_sdlng_mdd, & init_value=0.0_r8) - call this%sdlng_emerg_smp(pft)%p%InitRMean(ema_sdlng_emerg_h2o, & + call this%sdlng_emerg_smp(pft)%p%InitRSumm(ema_sdlng_emerg_h2o, & init_value=init_seedling_smp) end do end if @@ -911,7 +922,7 @@ subroutine FreeMemory(this, regeneration_model, numpft) call endrun(msg=errMsg(sourcefile, __LINE__)) endif - ! deallocate running means + ! deallocate running summaries deallocate(this%tveg24, stat=istat, errmsg=smsg) if (istat/=0) then write(fates_log(),*) 'dealloc011: fail on deallocate(this%tveg24):'//trim(smsg) @@ -928,6 +939,12 @@ subroutine FreeMemory(this, regeneration_model, numpft) call endrun(msg=errMsg(sourcefile, __LINE__)) endif + do pft = 1, numpft + deallocate(this%btran24_ft(pft)%p) + end do + deallocate(this%btran24_ft) + + if (regeneration_model == TRS_regeneration) then deallocate(this%seedling_layer_par24) deallocate(this%sdlng_mort_par) diff --git a/main/CMakeLists.txt b/main/CMakeLists.txt index 8db1d06331..5c52174d1c 100644 --- a/main/CMakeLists.txt +++ b/main/CMakeLists.txt @@ -22,7 +22,7 @@ list(APPEND fates_sources EDParamsMod.F90 EDTypesMod.F90 FatesHydraulicsMemMod.F90 - FatesRunningMeanMod.F90 + FatesRunningSummMod.F90 FatesParameterDerivedMod.F90 FatesSizeAgeTypeIndicesMod.F90 FatesIntegratorsMod.F90 diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index c21cdd6fe1..d1c9078114 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -25,7 +25,7 @@ module EDTypesMod use FatesLitterMod, only : litter_type use FatesLitterMod, only : ncwd use FatesConstantsMod, only : days_per_year - use FatesRunningMeanMod, only : rmean_type,rmean_arr_type + use FatesRunningSummMod, only : rsumm_type,rsumm_arr_type use FatesConstantsMod, only : fates_unset_r8 use FatesInterfaceTypesMod,only : bc_in_type use FatesInterfaceTypesMod,only : bc_out_type diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 86d3bbeddb..89100df6b7 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -630,6 +630,8 @@ module FatesHistoryInterfaceMod integer :: ih_site_dstatus_si_pft integer :: ih_dleafoff_si_pft integer :: ih_dleafon_si_pft + integer :: ih_meanbtran24_si_pft + integer :: ih_minbtran24_si_pft integer :: ih_meanliqvol_si_pft integer :: ih_meansmp_si_pft integer :: ih_elong_factor_si_pft @@ -3266,6 +3268,8 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) 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, & + hio_meanbtran24_si_pft => this%hvars(ih_meanbtran24_si_pft)%r82d, & + hio_minbtran24_si_pft => this%hvars(ih_minbtran24_si_pft)%r82d, & hio_meanliqvol_si_pft => this%hvars(ih_meanliqvol_si_pft)%r82d, & hio_meansmp_si_pft => this%hvars(ih_meansmp_si_pft)%r82d, & hio_elong_factor_si_pft => this%hvars(ih_elong_factor_si_pft)%r82d, & @@ -3422,6 +3426,13 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) cpatch%frac_burnt * cpatch%area * AREA_INV / sec_per_day endif + ! 24hr mean btran + hio_meanbtran24_si_pft(io_si,ft) = hio_meanbtran24_si_pft(io_si,ft) + & + cpatch%btran24_ft(ft)%p%GetMean() * cpatch%area * AREA_INV + + ! 24hr minimum btran + hio_minbtran24_si_pft(io_si,ft) = hio_minbtran24_si_pft(io_si,ft) + & + cpatch%btran24_ft(ft)%p%GetMin() * cpatch%area * AREA_INV end do ! loop through cohorts on patch @@ -7086,6 +7097,20 @@ subroutine define_history_vars(this, initialize_variables) upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & index=ih_dleafon_si_pft) + call this%set_history_var(vname='FATES_MEANBTRAN24_PF', & + units='1', & + long='PFT-level 24hr mean transpiration wetness factor (btran)', & + use_default='inactive', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & + index=ih_meanbtran24_si_pft) + + call this%set_history_var(vname='FATES_MINBTRAN24_PF', & + units='1', & + long='PFT-level 24hr minimum transpiration wetness factor (btran)', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & + index=ih_minbtran24_si_pft) + call this%set_history_var(vname='FATES_MEANLIQVOL_DROUGHTPHEN_PF', & units='m3 m-3', & long='PFT-level mean liquid water volume for drought phenolgy', & diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 555fa526c1..614f722938 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -99,15 +99,15 @@ module FatesInterfaceMod use PRTInitParamsFatesMod , only : PRTCheckParams, PRTDerivedParams use PRTAllometricCarbonMod , only : InitPRTGlobalAllometricCarbon use PRTAllometricCNPMod , only : InitPRTGlobalAllometricCNP - use FatesRunningMeanMod , only : ema_24hr - use FatesRunningMeanMod , only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par - use FatesRunningMeanMod , only : ema_sdlng_mdd, ema_sdlng2sap_par - use FatesRunningMeanMod , only : fixed_24hr - use FatesRunningMeanMod , only : ema_lpa - use FatesRunningMeanMod , only : ema_longterm - use FatesRunningMeanMod , only : ema_60day - use FatesRunningMeanMod , only : moving_ema_window - use FatesRunningMeanMod , only : fixed_window + use FatesRunningSummMod , only : ema_24hr + use FatesRunningSummMod , only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par + use FatesRunningSummMod , only : ema_sdlng_mdd, ema_sdlng2sap_par + use FatesRunningSummMod , only : fixed_24hr + use FatesRunningSummMod , only : ema_lpa + use FatesRunningSummMod , only : ema_longterm + use FatesRunningSummMod , only : ema_60day + use FatesRunningSummMod , only : moving_ema_window + use FatesRunningSummMod , only : fixed_window use FatesHistoryInterfaceMod , only : fates_hist use FatesHydraulicsMemMod , only : nshell use FatesHydraulicsMemMod , only : nlevsoi_hyd_max @@ -530,6 +530,7 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats, allocate(bc_in%tgcm_pa(maxpatch_total)) allocate(bc_in%t_soisno_sl(nlevsoil_in)) + ! Canopy Radiation bc_in%coszen = nan allocate(bc_in%fcansno_pa(maxpatch_total)) @@ -1085,8 +1086,8 @@ subroutine InitTimeAveragingGlobals() !allocate(ema_60day) !call ema_60day%define(prt_params%fnrt_adapt_tscl*sec_per_day,sec_per_day,moving_ema_window) - !class(rmean_arr_type), pointer :: ema_fnrt_tscale(:) - !rmean_arr_type + !class(rsumm_arr_type), pointer :: ema_fnrt_tscale(:) + !rsumm_arr_type return @@ -2233,7 +2234,7 @@ end subroutine FatesReportParameters subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) ! In this routine, we update any FATES buffers where - ! we calculate running means. It is assumed that this buffer is updated + ! we calculate running summaries. It is assumed that this buffer is updated ! on the model time-step. type(ed_site_type), intent(inout) :: sites(:) @@ -2265,9 +2266,15 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) nocomp_bare: if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then - call cpatch%tveg24%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) - call cpatch%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) - call cpatch%tveg_longterm%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) + call cpatch%tveg24%UpdateRSumm(bc_in(s)%t_veg_pa(ifp)) + call cpatch%tveg_lpa%UpdateRSumm(bc_in(s)%t_veg_pa(ifp)) + call cpatch%tveg_longterm%UpdateRSumm(bc_in(s)%t_veg_pa(ifp)) + + ! Update btran + do pft = 1, numpft + call cpatch%btran24_ft(pft)%p%UpdateRSumm(cpatch%btran_ft(pft)) + end do + ! Update the seedling layer par running means if ( hlm_regeneration_model == TRS_regeneration ) then @@ -2285,9 +2292,9 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) new_seedling_layer_par = seedling_par_high*par_high_frac + seedling_par_low*par_low_frac - call cpatch%seedling_layer_par24%UpdateRMean(new_seedling_layer_par) - call cpatch%sdlng_mort_par%UpdateRMean(new_seedling_layer_par) - call cpatch%sdlng2sap_par%UpdateRMean(new_seedling_layer_par) + call cpatch%seedling_layer_par24%UpdateRSumm(new_seedling_layer_par) + call cpatch%sdlng_mort_par%UpdateRSumm(new_seedling_layer_par) + call cpatch%sdlng2sap_par%UpdateRSumm(new_seedling_layer_par) do pft = 1,numpft @@ -2307,8 +2314,8 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) endif ! Update the seedling layer smp and mdd running means - call cpatch%sdlng_emerg_smp(pft)%p%UpdateRMean(new_seedling_layer_smp) - call cpatch%sdlng_mdd(pft)%p%UpdateRMean(new_seedling_mdd) + call cpatch%sdlng_emerg_smp(pft)%p%UpdateRSumm(new_seedling_layer_smp) + call cpatch%sdlng_mdd(pft)%p%UpdateRSumm(new_seedling_mdd) enddo !end pft loop @@ -2316,7 +2323,7 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out) ccohort => cpatch%tallest do while (associated(ccohort)) - ! call ccohort%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) + ! call ccohort%tveg_lpa%UpdateRSumm(bc_in(s)%t_veg_pa(ifp)) if(.not.ccohort%isnew)then ! [kgC/plant/yr] -> [gC/m2/s] site_npp = site_npp + ccohort%npp_acc_hold * ccohort%n*area_inv * & diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index b5fe1acfe9..d710988023 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -71,7 +71,7 @@ module FatesInventoryInitMod use PRTGenericMod, only : phosphorus_element use PRTGenericMod, only : SetState use FatesConstantsMod, only : primaryland - use FatesRunningMeanMod, only : ema_lpa + use FatesRunningSummMod, only : ema_lpa use PRTGenericMod, only : StorageNutrientTarget use FatesConstantsMod, only : fates_unset_int use EDCanopyStructureMod, only : canopy_summarization, canopy_structure @@ -1040,7 +1040,7 @@ subroutine set_inventory_cohort_type1(csite,bc_in,css_file_unit,npatches, & ! (Keeping as an example) ! Allocate running mean functions !allocate(temp_cohort%tveg_lpa) - !call temp_cohort%tveg_lpa%InitRMean(ema_lpa,init_value=cpatch%tveg_lpa%GetMean()) + !call temp_cohort%tveg_lpa%InitRSumm(ema_lpa,init_value=cpatch%tveg_lpa%GetMean()) do el = 1,num_elements diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index f2c78051a1..80fe70e6f4 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -51,8 +51,8 @@ module FatesRestartInterfaceMod use EDParamsMod, only : nlevleaf use PRTGenericMod, only : prt_global use PRTGenericMod, only : num_elements - use FatesRunningMeanMod, only : rmean_type - use FatesRunningMeanMod, only : ema_lpa + use FatesRunningSummMod, only : rsumm_type + use FatesRunningSummMod, only : ema_lpa use FatesRadiationMemMod, only : num_swb,norman_solver,twostr_solver use TwoStreamMLPEMod, only : normalized_upper_boundary use FatesConstantsMod, only : n_term_mort_types @@ -224,6 +224,8 @@ module FatesRestartInterfaceMod integer :: ir_scorch_ht_pa_pft integer :: ir_litter_moisture_pa_nfsc + integer :: ir_btran24_pa_pft + ! Site level integer :: ir_dd_status_sift integer :: ir_dleafondate_sift @@ -391,9 +393,9 @@ module FatesRestartInterfaceMod procedure, private :: GetCohortRealVector procedure, private :: SetCohortRealVector procedure, private :: RegisterCohortVector - procedure, private :: DefineRMeanRestartVar - procedure, private :: GetRMeanRestartVar - procedure, private :: SetRMeanRestartVar + procedure, private :: DefineRSummRestartVar + procedure, private :: GetRSummRestartVar + procedure, private :: SetRSummRestartVar end type fates_restart_interface_type @@ -1055,7 +1057,10 @@ subroutine define_restart_vars(this, initialize_variables) long_name='scorch height', units='m', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_litter_moisture_pa_nfsc) end if - + + call this%DefineRSummRestartVar(vname='fates_btran24_pa_pft',vtype=cohort_r8, & + long_name='24-hour patch transpiration wetness factor', & + units='1', initialize=initialize_variables,ivar=ivar, index = ir_btran24_pa_pft) call this%RegisterCohortVector(symbol_base='fates_year_net_up', vtype=cohort_r8, & long_name_base='yearly net uptake at leaf layers', & @@ -1564,48 +1569,49 @@ subroutine define_restart_vars(this, initialize_variables) units='kg/m2/yr', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_emanpp_si) - call this%DefineRMeanRestartVar(vname='fates_tveg24patch',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_tveg24patch',vtype=cohort_r8, & long_name='24-hour patch veg temp', & units='K', initialize=initialize_variables,ivar=ivar, index = ir_tveg24_pa) - call this%DefineRMeanRestartVar(vname='fates_disturbance_rates',vtype=cohort_r8, & + call this%set_restart_var(vname='fates_disturbance_rates',vtype=cohort_r8, & long_name='disturbance rates by donor land-use type, receiver land-use type, and disturbance type', & - units='1/day', initialize=initialize_variables,ivar=ivar, index = ir_disturbance_rates_siluludi) + units='1/day', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_disturbance_rates_siluludi) if ( hlm_regeneration_model == TRS_regeneration ) then - call this%DefineRMeanRestartVar(vname='fates_seedling_layer_par24',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_seedling_layer_par24',vtype=cohort_r8, & long_name='24-hour seedling layer PAR', & units='W m2-1', initialize=initialize_variables,ivar=ivar, index = ir_seedling_layer_par24_pa) - call this%DefineRMeanRestartVar(vname='fates_sdlng_emerg_smp',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_sdlng_emerg_smp',vtype=cohort_r8, & long_name='seedling layer PAR on the seedling emergence timescale', & units='mm suction', initialize=initialize_variables,ivar=ivar, index = ir_sdlng_emerg_smp_pa) - call this%DefineRMeanRestartVar(vname='fates_sdlng_mort_par',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_sdlng_mort_par',vtype=cohort_r8, & long_name='seedling layer PAR on the seedling mortality timescale', & units='W m2-1', initialize=initialize_variables,ivar=ivar, index = ir_sdlng_mort_par_pa) - call this%DefineRMeanRestartVar(vname='fates_sdlng2sap_par',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_sdlng2sap_par',vtype=cohort_r8, & long_name='seedling layer PAR on the seedling to sapling transition timescale', & units='W m2-1', initialize=initialize_variables,ivar=ivar, index = ir_sdlng2sap_par_pa) - call this%DefineRMeanRestartVar(vname='fates_sdlng_mdd',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_sdlng_mdd',vtype=cohort_r8, & long_name='seedling moisture deficit days', & units='mm days', initialize=initialize_variables,ivar=ivar, index = ir_sdlng_mdd_pa) end if - call this%DefineRMeanRestartVar(vname='fates_tveglpapatch',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_tveglpapatch',vtype=cohort_r8, & long_name='running average (EMA) of patch veg temp for photo acclim', & units='K', initialize=initialize_variables,ivar=ivar, index = ir_tveglpa_pa) - call this%DefineRMeanRestartVar(vname='fates_tveglongtermpatch',vtype=cohort_r8, & + call this%DefineRSummRestartVar(vname='fates_tveglongtermpatch',vtype=cohort_r8, & long_name='long-term (T_home) running average (EMA) of patch veg temp for photo acclim', & units='K', initialize=initialize_variables,ivar=ivar, index = ir_tveglongterm_pa) ! (Keeping as an example) - !call this%DefineRMeanRestartVar(vname='fates_tveglpacohort',vtype=cohort_r8, & + !call this%DefineRSummRestartVar(vname='fates_tveglpacohort',vtype=cohort_r8, & ! long_name='running average (EMA) of cohort veg temp for photo acclim', & ! units='K', initialize=initialize_variables,ivar=ivar, index = ir_tveglpa_co) @@ -1624,7 +1630,7 @@ end subroutine define_restart_vars ! ===================================================================================== - subroutine DefineRMeanRestartVar(this,vname,vtype,long_name,units,initialize,ivar,index) + subroutine DefineRSummRestartVar(this,vname,vtype,long_name,units,initialize,ivar,index) class(fates_restart_interface_type) :: this character(len=*),intent(in) :: vname @@ -1647,6 +1653,26 @@ subroutine DefineRMeanRestartVar(this,vname,vtype,long_name,units,initialize,iva units=units, flushval = flushzero, & hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = dummy_index ) + call this%set_restart_var(vname= trim(vname)//'_cminimum', vtype=vtype, & + long_name=long_name//' current minimum', & + units=units, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = index ) + + call this%set_restart_var(vname= trim(vname)//'_lminimum', vtype=vtype, & + long_name=long_name//' latest minimum', & + units=units, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = dummy_index ) + + call this%set_restart_var(vname= trim(vname)//'_cmaximum', vtype=vtype, & + long_name=long_name//' current maximum', & + units=units, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = index ) + + call this%set_restart_var(vname= trim(vname)//'_lmaximum', vtype=vtype, & + long_name=long_name//' latest maximum', & + units=units, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = dummy_index ) + call this%set_restart_var(vname= trim(vname)//'_cindex', vtype=vtype, & long_name=long_name//' index', & units='index', flushval = flushzero, & @@ -1654,15 +1680,15 @@ subroutine DefineRMeanRestartVar(this,vname,vtype,long_name,units,initialize,iva return - end subroutine DefineRMeanRestartVar + end subroutine DefineRSummRestartVar ! ===================================================================================== - subroutine GetRMeanRestartVar(this, rmean_var, ir_var_index, position_index) + subroutine GetRSummRestartVar(this, rsumm_var, ir_var_index, position_index) class(fates_restart_interface_type) , intent(inout) :: this - class(rmean_type), intent(inout) :: rmean_var + class(rsumm_type), intent(inout) :: rsumm_var integer,intent(in) :: ir_var_index integer,intent(in) :: position_index @@ -1671,21 +1697,29 @@ subroutine GetRMeanRestartVar(this, rmean_var, ir_var_index, position_index) integer :: ir_pos_var ! global variable index - rmean_var%c_mean = this%rvars(ir_var_index)%r81d(position_index) + rsumm_var%c_mean = this%rvars(ir_var_index)%r81d(position_index) - rmean_var%l_mean = this%rvars(ir_var_index+1)%r81d(position_index) + rsumm_var%l_mean = this%rvars(ir_var_index+1)%r81d(position_index) - rmean_var%c_index = nint(this%rvars(ir_var_index+2)%r81d(position_index)) + rsumm_var%c_minimum = this%rvars(ir_var_index+2)%r81d(position_index) + + rsumm_var%l_minimum = this%rvars(ir_var_index+3)%r81d(position_index) + + rsumm_var%c_maximum = this%rvars(ir_var_index+4)%r81d(position_index) + + rsumm_var%l_maximum = this%rvars(ir_var_index+5)%r81d(position_index) + + rsumm_var%c_index = nint(this%rvars(ir_var_index+6)%r81d(position_index)) return - end subroutine GetRMeanRestartVar + end subroutine GetRSummRestartVar ! ======================================================================================= - subroutine SetRMeanRestartVar(this, rmean_var, ir_var_index, position_index) + subroutine SetRSummRestartVar(this, rsumm_var, ir_var_index, position_index) class(fates_restart_interface_type) , intent(inout) :: this - class(rmean_type), intent(inout) :: rmean_var + class(rsumm_type), intent(inout) :: rsumm_var integer,intent(in) :: ir_var_index integer,intent(in) :: position_index @@ -1693,14 +1727,22 @@ subroutine SetRMeanRestartVar(this, rmean_var, ir_var_index, position_index) integer :: i_pos ! vector position loop index integer :: ir_pos_var ! global variable index - this%rvars(ir_var_index)%r81d(position_index) = rmean_var%c_mean - - this%rvars(ir_var_index+1)%r81d(position_index) = rmean_var%l_mean - - this%rvars(ir_var_index+2)%r81d(position_index) = real(rmean_var%c_index,r8) + this%rvars(ir_var_index)%r81d(position_index) = rsumm_var%c_mean + + this%rvars(ir_var_index+1)%r81d(position_index) = rsumm_var%l_mean + + this%rvars(ir_var_index+2)%r81d(position_index) = rsumm_var%c_minimum + + this%rvars(ir_var_index+3)%r81d(position_index) = rsumm_var%l_minimum + + this%rvars(ir_var_index+4)%r81d(position_index) = rsumm_var%c_maximum + + this%rvars(ir_var_index+5)%r81d(position_index) = rsumm_var%l_maximum + + this%rvars(ir_var_index+6)%r81d(position_index) = real(rsumm_var%c_index,r8) return - end subroutine SetRMeanRestartVar + end subroutine SetRSummRestartVar ! ===================================================================================== @@ -2533,7 +2575,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) endif ! (Keeping as an example) - ! call this%SetRMeanRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) + ! call this%SetRSummRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) io_idx_co = io_idx_co + 1 @@ -2551,19 +2593,25 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_nocomp_pft_label_pa(io_idx_co_1st)= cpatch%nocomp_pft_label rio_area_pa(io_idx_co_1st) = cpatch%area - ! Patch level running means - call this%SetRMeanRestartVar(cpatch%tveg24, ir_tveg24_pa, io_idx_co_1st) - call this%SetRMeanRestartVar(cpatch%tveg_lpa, ir_tveglpa_pa, io_idx_co_1st) - call this%SetRMeanRestartVar(cpatch%tveg_longterm, ir_tveglongterm_pa, io_idx_co_1st) + ! Patch level running summaries + call this%SetRSummRestartVar(cpatch%tveg24, ir_tveg24_pa, io_idx_co_1st) + call this%SetRSummRestartVar(cpatch%tveg_lpa, ir_tveglpa_pa, io_idx_co_1st) + call this%SetRSummRestartVar(cpatch%tveg_longterm, ir_tveglongterm_pa, io_idx_co_1st) + + io_idx_pa_pft = io_idx_co_1st + do i_pft = 1, numpft + call this%SetRSummRestartVar(cpatch%btran24_ft(i_pft)%p, ir_btran24_pa_pft,io_idx_pa_pft) + io_idx_pa_pft = io_idx_pa_pft + 1 + end do if ( hlm_regeneration_model == TRS_regeneration ) then - call this%SetRMeanRestartVar(cpatch%seedling_layer_par24, ir_seedling_layer_par24_pa, io_idx_co_1st) - call this%SetRMeanRestartVar(cpatch%sdlng_mort_par, ir_sdlng_mort_par_pa,io_idx_co_1st) - call this%SetRMeanRestartVar(cpatch%sdlng2sap_par, ir_sdlng2sap_par_pa,io_idx_co_1st) + call this%SetRSummRestartVar(cpatch%seedling_layer_par24, ir_seedling_layer_par24_pa, io_idx_co_1st) + call this%SetRSummRestartVar(cpatch%sdlng_mort_par, ir_sdlng_mort_par_pa,io_idx_co_1st) + call this%SetRSummRestartVar(cpatch%sdlng2sap_par, ir_sdlng2sap_par_pa,io_idx_co_1st) io_idx_pa_pft = io_idx_co_1st do i_pft = 1, numpft - call this%SetRMeanRestartVar(cpatch%sdlng_mdd(i_pft)%p, ir_sdlng_mdd_pa,io_idx_pa_pft) - call this%SetRMeanRestartVar(cpatch%sdlng_emerg_smp(i_pft)%p, ir_sdlng_emerg_smp_pa,io_idx_pa_pft) + call this%SetRSummRestartVar(cpatch%sdlng_mdd(i_pft)%p, ir_sdlng_mdd_pa,io_idx_pa_pft) + call this%SetRSummRestartVar(cpatch%sdlng_emerg_smp(i_pft)%p, ir_sdlng_emerg_smp_pa,io_idx_pa_pft) io_idx_pa_pft = io_idx_pa_pft + 1 enddo end if @@ -2960,7 +3008,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in, bc_out) ! (Keeping as an example) ! Allocate running mean functions !allocate(new_cohort%tveg_lpa) - !call new_cohort%tveg_lpa%InitRMean(ema_lpa) + !call new_cohort%tveg_lpa%InitRSumm(ema_lpa) ! Update the previous @@ -3502,7 +3550,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end if ! (Keeping as an example) - !call this%GetRMeanRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) + !call this%GetRSummRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) ccohort%c_area = this%rvars(ir_c_area_co)%r81d(io_idx_co) ccohort%treelai = this%rvars(ir_treelai_co)%r81d(io_idx_co) @@ -3537,18 +3585,25 @@ subroutine get_restart_vectors(this, nc, nsites, sites) cpatch%ncl_p = this%rvars(ir_nclp_pa)%int1d(io_idx_co_1st) cpatch%zstar = this%rvars(ir_zstar_pa)%r81d(io_idx_co_1st) - call this%GetRMeanRestartVar(cpatch%tveg24, ir_tveg24_pa, io_idx_co_1st) - call this%GetRMeanRestartVar(cpatch%tveg_lpa, ir_tveglpa_pa, io_idx_co_1st) - call this%GetRMeanRestartVar(cpatch%tveg_longterm, ir_tveglongterm_pa, io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%tveg24, ir_tveg24_pa, io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%tveg_lpa, ir_tveglpa_pa, io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%tveg_longterm, ir_tveglongterm_pa, io_idx_co_1st) + + ! Daily summary for btran + io_idx_pa_pft = io_idx_co_1st + do pft = 1, numpft + call this%GetRSummRestartVar(cpatch%btran24_ft(pft)%p, ir_btran24_pa_pft, io_idx_pa_pft) + io_idx_pa_pft = io_idx_pa_pft + 1 + enddo if ( hlm_regeneration_model == TRS_regeneration ) then - call this%GetRMeanRestartVar(cpatch%seedling_layer_par24, ir_seedling_layer_par24_pa, io_idx_co_1st) - call this%GetRMeanRestartVar(cpatch%sdlng_mort_par, ir_sdlng_mort_par_pa,io_idx_co_1st) - call this%GetRMeanRestartVar(cpatch%sdlng2sap_par, ir_sdlng2sap_par_pa,io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%seedling_layer_par24, ir_seedling_layer_par24_pa, io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%sdlng_mort_par, ir_sdlng_mort_par_pa,io_idx_co_1st) + call this%GetRSummRestartVar(cpatch%sdlng2sap_par, ir_sdlng2sap_par_pa,io_idx_co_1st) io_idx_pa_pft = io_idx_co_1st do pft = 1, numpft - call this%GetRMeanRestartVar(cpatch%sdlng_mdd(pft)%p, ir_sdlng_mdd_pa,io_idx_pa_pft) - call this%GetRMeanRestartVar(cpatch%sdlng_emerg_smp(pft)%p, ir_sdlng_emerg_smp_pa,io_idx_pa_pft) + call this%GetRSummRestartVar(cpatch%sdlng_mdd(pft)%p, ir_sdlng_mdd_pa,io_idx_pa_pft) + call this%GetRSummRestartVar(cpatch%sdlng_emerg_smp(pft)%p, ir_sdlng_emerg_smp_pa,io_idx_pa_pft) io_idx_pa_pft = io_idx_pa_pft + 1 enddo end if diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningSummMod.F90 similarity index 55% rename from main/FatesRunningMeanMod.F90 rename to main/FatesRunningSummMod.F90 index 721030b974..cd75e6f7b2 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningSummMod.F90 @@ -1,8 +1,9 @@ -module FatesRunningMeanMod +module FatesRunningSummMod use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : fates_huge use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) use shr_log_mod, only : errMsg => shr_log_errMsg use FatesGlobals, only : endrun => fates_endrun @@ -29,7 +30,7 @@ module FatesRunningMeanMod ! memory period is, and if it should be zero'd ! These are globally defined on the proc. - type, public :: rmean_def_type + type, public :: rsumm_def_type real(r8) :: mem_period ! The total integration period (s) real(r8) :: up_period ! The period between updates (s) @@ -40,13 +41,13 @@ module FatesRunningMeanMod procedure :: define - end type rmean_def_type + end type rsumm_def_type ! This holds the time varying information for the mean ! which is instantiated on sites, patches, and cohorts - type, public :: rmean_type + type, public :: rsumm_type real(r8) :: c_mean ! The current mean value, if this ! is a moving window, it is the mean. @@ -59,6 +60,22 @@ module FatesRunningMeanMod ! as c_mean for moving windows, and for fixed windows ! it is the mean value when the time integration window ! last completed. + + real(r8) :: c_minimum ! The current minimum value. Depending on the time, + ! this will be a partial minimum as it needs to + ! compare with all to fully determine the minimum + + real(r8) :: l_minimum ! The latest reportable minimum value. This is the + ! minimum value for the last window period that was + ! fully assessed. + + real(r8) :: c_maximum ! The current maximum value. Depending on the time, + ! this will be a partial minimum as it needs to + ! compare with all to fully determine the maximum + + real(r8) :: l_maximum ! The latest reportable maximum value. This is the + ! maximum value for the last window period that was + ! fully assessed. integer :: c_index ! The number of values that have ! been added to the mean so far @@ -67,17 +84,19 @@ module FatesRunningMeanMod ! This points to the global structure that ! defines the nature of this mean/avg - type(rmean_def_type), pointer :: def_type + type(rsumm_def_type), pointer :: def_type contains procedure :: GetMean - procedure :: InitRMean - procedure :: UpdateRMean - procedure :: FuseRMean + procedure :: GetMin + procedure :: GetMax + procedure :: InitRSumm + procedure :: UpdateRSumm + procedure :: FuseRSumm procedure :: CopyFromDonor - end type rmean_type + end type rsumm_type logical, parameter :: debug = .true. @@ -87,37 +106,38 @@ module FatesRunningMeanMod ! Define the time methods that we want to have available to us - class(rmean_def_type), public, pointer :: ema_24hr ! Exponential moving average - 24hr window - class(rmean_def_type), public, pointer :: fixed_24hr ! Fixed, 24-hour window - class(rmean_def_type), public, pointer :: ema_lpa ! Exponential moving average - leaf photo acclimation - class(rmean_def_type), public, pointer :: ema_longterm ! Exponential moving average - long-term leaf photo acclimation - class(rmean_def_type), public, pointer :: ema_60day ! Exponential moving average, 60 day + class(rsumm_def_type), public, pointer :: ema_24hr ! Exponential moving average - 24hr window + class(rsumm_def_type), public, pointer :: fixed_24hr ! Fixed, 24-hour window + class(rsumm_def_type), public, pointer :: ema_lpa ! Exponential moving average - leaf photo acclimation + class(rsumm_def_type), public, pointer :: ema_longterm ! Exponential moving average - long-term leaf photo acclimation + class(rsumm_def_type), public, pointer :: ema_60day ! Exponential moving average, 60 day ! Updated daily - class(rmean_def_type), public, pointer :: ema_storemem ! EMA used for smoothing N/C and P/C storage - class(rmean_def_type), public, pointer :: ema_sdlng_mort_par ! EMA for seedling mort from light stress - class(rmean_def_type), public, pointer :: ema_sdlng2sap_par ! EMA for seedling to sapling transition rates + class(rsumm_def_type), public, pointer :: ema_storemem ! EMA used for smoothing N/C and P/C storage + class(rsumm_def_type), public, pointer :: ema_sdlng_mort_par ! EMA for seedling mort from light stress + class(rsumm_def_type), public, pointer :: ema_sdlng2sap_par ! EMA for seedling to sapling transition rates ! based in par - class(rmean_def_type), public, pointer :: ema_sdlng_emerg_h2o ! EMA for moisture-based seedling emergence - class(rmean_def_type), public, pointer :: ema_sdlng_mdd ! EMA for seedling moisture deficit days + class(rsumm_def_type), public, pointer :: ema_sdlng_emerg_h2o ! EMA for moisture-based seedling emergence + class(rsumm_def_type), public, pointer :: ema_sdlng_mdd ! EMA for seedling moisture deficit days ! If we want to have different running mean specs based on ! pft or other types of constants - type, public :: rmean_arr_type - class(rmean_type), pointer :: p - end type rmean_arr_type + type, public :: rsumm_arr_type + class(rsumm_type), pointer :: p + end type rsumm_arr_type contains subroutine define(this,mem_period,up_period,method) + ! This sub-routine sets parameters and method for this running summary structure - class(rmean_def_type) :: this + class(rsumm_def_type) :: this - real(r8),intent(in) :: mem_period - real(r8),intent(in) :: up_period - integer,intent(in) :: method + real(r8),intent(in) :: mem_period ! The total integration period (s) + real(r8),intent(in) :: up_period ! The period between updates (s) + integer,intent(in) :: method ! Is this a fixed or moving window? ! Check the memory and update periods if(debug) then @@ -143,42 +163,84 @@ end subroutine define ! ===================================================================================== function GetMean(this) + ! Function to retrieve the running average for this variable. - class(rmean_type) :: this + class(rsumm_type) :: this real(r8) :: GetMean - if(this%def_type%method .eq. moving_ema_window) then + select case(this%def_type%method) + case (moving_ema_window) if(this%c_index == 0 .and. debug) then write(fates_log(), *) 'attempting to get a running mean from a variable' write(fates_log(), *) 'that has not been given a value yet' call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end select GetMean = this%l_mean end function GetMean + + ! ===================================================================================== + function GetMin(this) + ! Function to retrieve the minimum value in this time window for this variable. + + class(rsumm_type) :: this + real(r8) :: GetMin + + select case(this%def_type%method) + case (moving_ema_window) + write(fates_log(), *) ' Attempting to get the minimum value from a variable' + write(fates_log(), *) 'that is integrated using exponential moving averages.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + case (fixed_window) + GetMin = this%l_minimum + end select + + + end function GetMin + ! ===================================================================================== - subroutine InitRMean(this,rmean_def,init_value,init_offset) + function GetMax(this) + ! Function to retrieve the maximum value in this time window for this variable. - class(rmean_type) :: this - type(rmean_def_type), target :: rmean_def - real(r8),intent(in),optional :: init_value - real(r8),intent(in),optional :: init_offset + class(rsumm_type) :: this + real(r8) :: GetMax + + select case(this%def_type%method) + case (moving_ema_window) + write(fates_log(), *) ' Attempting to get the maximum value from a variable' + write(fates_log(), *) 'that is integrated using exponential moving averages.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + case (fixed_window) + GetMax = this%l_maximum + end select + + end function GetMax + + ! ===================================================================================== + + subroutine InitRSumm(this,rsumm_def,init_value,init_offset) + ! Sub-routine to initialize the summary structure variables. ! If the initialization happens part-way through a fixed averaging window ! we need to account for this. The current method moves the position ! index to match the offset, and then assumes that the init_value provided ! was a constant over the offset period. - ! If the first value is offset, such that the we are a portion of the - ! way through the window, we need to account for this. + class(rsumm_type) :: this + type(rsumm_def_type), target :: rsumm_def + real(r8),intent(in),optional :: init_value + real(r8),intent(in),optional :: init_offset ! Point to the averaging type - this%def_type => rmean_def + this%def_type => rsumm_def - if(this%def_type%method .eq. fixed_window) then + select case(this%def_type%method) + case (fixed_window) if(debug) then if(.not.(present(init_offset).and.present(init_value)) )then @@ -189,11 +251,11 @@ subroutine InitRMean(this,rmean_def,init_value,init_offset) end if ! Check to see if the offset is an even increment of the update frequency - if( abs(real(nint(init_offset/rmean_def%up_period),r8)-(init_offset/rmean_def%up_period)) > nearzero ) then + if( abs(real(nint(init_offset/rsumm_def%up_period),r8)-(init_offset/rsumm_def%up_period)) > nearzero ) then write(fates_log(), *) 'when initializing a temporal mean on a fixed window' - write(fates_log(), *) 'the time offset must be an inrement of the update frequency' + write(fates_log(), *) 'the time offset must be an increment of the update frequency' write(fates_log(), *) 'offset: ',init_offset - write(fates_log(), *) 'up freq: ',rmean_def%up_period + write(fates_log(), *) 'up freq: ',rsumm_def%up_period call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -204,17 +266,26 @@ subroutine InitRMean(this,rmean_def,init_value,init_offset) end if - this%c_index = modulo(nint(init_offset/rmean_def%up_period),rmean_def%n_mem) - this%c_mean = real(this%c_index,r8)/real(rmean_def%n_mem,r8)*init_value + + ! If the first value is offset, such that the we are a portion of the + ! way through the window, we need to account for this. + this%c_index = modulo(nint(init_offset/rsumm_def%up_period),rsumm_def%n_mem) + this%c_mean = real(this%c_index,r8)/real(rsumm_def%n_mem,r8)*init_value this%l_mean = init_value - elseif(this%def_type%method .eq. moving_ema_window) then - - if(present(init_value))then + ! Minima and maxima are not weighted so they can both assigned the + ! initial value directly. + this%c_minimum = init_value + this%l_minimum = init_value + this%c_maximum = init_value + this%l_maximum = init_value + + case (moving_ema_window) + if (present(init_value)) then this%c_mean = init_value this%l_mean = init_value if(present(init_offset))then - this%c_index = min(nint(init_offset/rmean_def%up_period),rmean_def%n_mem) + this%c_index = min(nint(init_offset/rsumm_def%up_period),rsumm_def%n_mem) else this%c_index = 1 end if @@ -223,20 +294,26 @@ subroutine InitRMean(this,rmean_def,init_value,init_offset) this%l_mean = nan this%c_index = 0 end if - - end if - + + ! We do not track minima and maxima when handling exponential moving averages. + this%c_minimum = nan + this%l_minimum = nan + this%c_maximum = nan + this%l_maximum = nan + end select + return - end subroutine InitRMean + end subroutine InitRSumm ! ===================================================================================== subroutine CopyFromDonor(this, donor) + ! Sub-routine to transfer running summary information from one structure to another - class(rmean_type) :: this - class(rmean_type),intent(in) :: donor + class(rsumm_type) :: this ! Receptor structure + class(rsumm_type),intent(in) :: donor ! Donor structure if( .not.associated(this%def_type)) then write(fates_log(), *) 'Attempted to copy over running mean' @@ -246,11 +323,15 @@ subroutine CopyFromDonor(this, donor) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - this%c_mean = donor%c_mean - this%l_mean = donor%l_mean - this%c_index = donor%c_index + this%c_mean = donor%c_mean + this%l_mean = donor%l_mean + this%c_minimum = donor%c_minimum + this%l_minimum = donor%l_minimum + this%c_maximum = donor%c_maximum + this%l_maximum = donor%l_maximum + this%c_index = donor%c_index + - return end subroutine CopyFromDonor @@ -258,26 +339,28 @@ end subroutine CopyFromDonor ! ===================================================================================== - subroutine UpdateRMean(this, new_value) + subroutine UpdateRSumm(this, new_value) + ! Sub-routine that updates the running summary metrics. - class(rmean_type) :: this - real(r8) :: new_value ! The newest value added to the running mean - real(r8) :: wgt + class(rsumm_type) :: this + real(r8), intent(in) :: new_value ! The newest value added to the running mean + real(r8) :: wgt - if(this%def_type%method.eq.moving_ema_window) then + select case (this%def_type%method) + case (moving_ema_window) this%c_index = min(this%def_type%n_mem,this%c_index + 1) - + if(this%c_index==1) then this%c_mean = new_value else wgt = 1._r8/real(this%c_index,r8) this%c_mean = this%c_mean*(1._r8-wgt) + wgt*new_value end if - + this%l_mean = this%c_mean - - else + + case (fixed_window) ! If the last time we updated we had hit the ! end of the averaging memory period, and @@ -287,23 +370,29 @@ subroutine UpdateRMean(this, new_value) this%c_index = this%c_index + 1 wgt = this%def_type%up_period/this%def_type%mem_period this%c_mean = this%c_mean + new_value*wgt + this%c_minimum = min(this%c_minimum,new_value) + this%c_maximum = max(this%c_maximum,new_value) if(this%c_index == this%def_type%n_mem) then - this%l_mean = this%c_mean - this%c_mean = 0._r8 + this%l_mean = this%c_mean + this%c_mean = 0._r8 + this%l_minimum = this%c_minimum + this%c_minimum = fates_huge + this%l_maximum = this%c_maximum + this%c_maximum = -fates_huge this%c_index = 0 end if - - end if + end select return - end subroutine UpdateRmean + end subroutine UpdateRSumm ! ===================================================================================== - subroutine FuseRMean(this,donor,recip_wgt) - + subroutine FuseRSumm(this,donor,recip_wgt) + ! Sub-routine that fuses data from two running summary structures. + ! ! Rules for fusion: ! If both entities have valid means already, then you simply use the ! weight provided to combine them. @@ -313,18 +402,14 @@ subroutine FuseRMean(this,donor,recip_wgt) ! both. - class(rmean_type) :: this - class(rmean_type), pointer :: donor + class(rsumm_type) :: this + class(rsumm_type), pointer :: donor real(r8),intent(in) :: recip_wgt ! Weighting factor for recipient (0-1) - ! Unecessary - !if(this%def_type%n_mem .ne. donor%def_type%n_mem) then - ! write(fates_log(), *) 'memory size is somehow different during fusion' - ! write(fates_log(), *) 'of two running mean variables: '!,this%name,donor%name - ! call endrun(msg=errMsg(sourcefile, __LINE__)) - !end if - if(this%def_type%method .eq. fixed_window ) then + + select case (this%def_type%method) + case (fixed_window) if (this%c_index .ne. donor%c_index) then write(fates_log(), *) 'trying to fuse two fixed-window averages' write(fates_log(), *) 'that are at different points in the window?' @@ -333,21 +418,36 @@ subroutine FuseRMean(this,donor,recip_wgt) write(fates_log(), *) 'c_index', this%c_index, donor%c_index call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end select ! This last logic clause just simply prevents us from doing math ! on uninitialized values. If both are unintiailized, then ! leave the result as uninitialized - if( .not. (donor%c_index==0) ) then - - if(this%c_index==0) then - this%c_mean = donor%c_mean - this%l_mean = donor%l_mean - this%c_index = donor%c_index + if ( donor%c_index /= 0 ) then + + if ( this%c_index == 0 ) then + this%c_mean = donor%c_mean + this%l_mean = donor%l_mean + this%c_minimum = donor%c_minimum + this%l_minimum = donor%l_minimum + this%c_maximum = donor%c_maximum + this%l_maximum = donor%l_maximum + this%c_index = donor%c_index else ! Take the weighted mean between the two this%c_mean = this%c_mean*recip_wgt + donor%c_mean*(1._r8-recip_wgt) this%l_mean = this%l_mean*recip_wgt + donor%l_mean*(1._r8-recip_wgt) + + ! Minima and maxima are fused only for fixed window, otherwise they are + ! not computed and must remain nan + select case (this%def_type%method) + case (fixed_window) + this%c_minimum = this%c_minimum*recip_wgt + donor%c_minimum*(1._r8-recip_wgt) + this%l_minimum = this%l_minimum*recip_wgt + donor%l_minimum*(1._r8-recip_wgt) + this%c_maximum = this%c_maximum*recip_wgt + donor%c_maximum*(1._r8-recip_wgt) + this%l_maximum = this%l_maximum*recip_wgt + donor%l_maximum*(1._r8-recip_wgt) + end select + ! Update the index to the larger of the two this%c_index = max(this%c_index,donor%c_index) end if @@ -355,7 +455,7 @@ subroutine FuseRMean(this,donor,recip_wgt) end if return - end subroutine FuseRMean + end subroutine FuseRSumm -end module FatesRunningMeanMod +end module FatesRunningSummMod diff --git a/testing/testing_shr/FatesFactoryMod.F90 b/testing/testing_shr/FatesFactoryMod.F90 index aede3dbab9..d23d4d72fd 100644 --- a/testing/testing_shr/FatesFactoryMod.F90 +++ b/testing/testing_shr/FatesFactoryMod.F90 @@ -21,8 +21,8 @@ module FatesFactoryMod use EDParamsMod, only : nclmax use EDParamsMod, only : photo_temp_acclim_timescale use EDParamsMod, only : photo_temp_acclim_thome_time - use FatesRunningMeanMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm - use FatesRunningMeanMod, only : moving_ema_window, fixed_window + use FatesRunningSummMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm + use FatesRunningSummMod, only : moving_ema_window, fixed_window use EDCohortDynamicsMod, only : InitPRTObject use PRTParametersMod, only : prt_params use PRTGenericMod, only : element_pos