From 1d1b87fd93088ac074a5fa5ff9e59d90e1433faf Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 3 Mar 2021 17:12:44 -0500 Subject: [PATCH 01/21] First pass at adding classes and procedures for the running-mean type --- main/FatesRunningMeanMod.F90 | 232 +++++++++++++++++++++++++++++++++++ 1 file changed, 232 insertions(+) create mode 100644 main/FatesRunningMeanMod.F90 diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningMeanMod.F90 new file mode 100644 index 0000000000..f7b6ca611d --- /dev/null +++ b/main/FatesRunningMeanMod.F90 @@ -0,0 +1,232 @@ +module FatesRunningMeanMod + + + use FatesConstantsMod, only : nearzero + 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 + + + integer, parameter :: maxlen_varname = 8 + + type, public :: rmean_type + + real(r8), allocatable(:) :: mem_rmean ! Array storing the memory of the mean + real(r8) :: c_rmean ! The current mean value from the + ! available memory, as of the last update + integer :: c_index ! The current memory index as per the + ! last update + integer :: n_mem ! Total number of memory indices + logical :: filled ! Has enough time elapsed so all memory filled? + character(len=maxlen_varname) :: var_name ! A short name for this variable + ! for diagnostic purposes + + contains + + procedure :: InitRMean + procedure :: UpdateRMean + procedure :: FuseRMean + + end type rmean_type + +contains + + + subroutine InitRMean(this,name,mem_period,up_period) + + class(rmean_type) :: this + character(len=maxlen_varname) :: name ! The name of the new variable + real(r8) :: mem_period ! The period length in seconds that must be remembered + real(r8) :: up_period ! The period length in seconds that memory is updated + ! (i.e. the resolution of the memory) + + this%name = name + this%n_mem = nint(mem_period/up_period) + + if( abs(real(this%n_mem,r8)-mem_period/up_period) > nearzero ) then + write(fates_log(), *) 'While defining a running mean variable: ',this%var_name + write(fates_log(), *) 'an update and total memory period was specified' + write(fates_log(), *) 'where the update period is not an exact fraction of the period' + write(fates_log(), *) 'mem_period: ',mem_period + write(fates_log(), *) 'up_period: ',up_period + write(fates_log(), *) 'exiting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Otherwise we allocate + allocate(this%mem_rmean(this%n_mem)) + + ! Initialize with nonsense numbers + this%mem_rmean(:) = nan + + ! There are no memory spots filled with valid datapoints + this%filled = .false + + ! The current index of the memory is zero + this%c_index = 0 + + return + end subroutine InitRMean + + ! ===================================================================================== + + subroutine UpdateRMean(this, new_value) + + class(rmean_type) :: this + real(r8) :: new_value ! The newest value added to the running mean + + this%c_index = this%c_index + 1 + + ! Update the index of the memory array + ! and, determine if we have filled the memory yet + ! If this index is greater than our memory slots, + ! go back to the first + if(this%c_index==this%n_mem) then + this%filled = .true. + end if + + if(this%c_index>this%n_mem) this%c_index = 1 + + this%mem_rmean(this%c_index) = new_value + + ! Update the running mean value. It will return a value + ! if we have not filled in all the memory slots. To do this + ! it will take a mean over what is available + + if(this%filled) then + this%c_rmean = sum(this%mem_rmean)/real(this%n_mem,r8) + else + this%c_rmean = sum(this%mem_rmean(1:this%c_index))/real(this%c_index,r8) + end if + + + return + end subroutine UpdateRmean + + ! ===================================================================================== + + subroutine FuseRMean(this,donor,recip_wgt) + + ! When fusing the running mean of two entities, it is possible that they + ! may have a different amount of memory spaces filled (at least in FATES). This + ! is typical for newly created patches or cohorts, that litteraly just spawned + ! So what generally happens is we walk backwards from the current memory index + ! of both and take means where we can. In places where one entity has more + ! memory than the other, than we just use the value from the one that is there + + class(rmean_type) :: this + class(rmean_type), pointer :: donor + real(r8),intent(in) :: recip_wgt ! Weighting factor for recipient (0-1) + + integer :: r_id ! Loop index counter for the recipient (this) + integer :: d_id ! Loop index counter for the donor + + + if(this%n_mem .ne. donor%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%filled .and. donor%filled) then + + ! Both are filled, take averages of each and be sure to use relative positions + + d_id = donor%c_index + do r_id = 1,this%n_mem + this%mem_rmean(r_id) = this%mem_rmean(r_id)*(recip_wgt) + & + donor%mem_rmean(d_id)*(1._r8-recip_wgt) + d_id = d_id+ipos + if(d_id>donor%n_mem) d_id = 1 + end do + + elseif(this%filled .and. .not.donor%filled) then + + ! Here, we have only partial memory from the donor + ! we we iterate through the donor's memory + ! and average between the two in those spaces. Then + ! we leave the rest untouched because we accept + ! the values from the recipient. Also, keep the + ! recipient's index. + + r_id = this%c_index + do d_id = donor%c_index,1,-1 + this%mem_rmean(r_id) = this%mem_rmean(r_id)*(recip_wgt) + donor%mem_rmean(d_id)*(1._r8-recip_wgt) + r_id = r_id - 1 + if(r_id==0) r_id = this%n_mem + end do + + + elseif(.not.this%filled .and. donor%filled) then + + ! Here we only have partial memory from the recipient + ! so we iterate through the recipient's memory + ! and average between the two. Then we copy + ! over the values from the donor for the indices that we + ! didn't average because the donor has valid values + + d_id = donor%c_index + do r_id = this%c_index,1,-1 + this%mem_rmean(r_id) = this%mem_rmean(r_id)*(recip_wgt) + donor%mem_rmean(d_id)*(1._r8-recip_wgt) + d_id = d_id - 1 + if(d_id==0) d_id = donor%n_mem + end do + + do r_id = this%n_mem,this%c_index+1,-1 + this%mem_rmean(r_id) = donor%mem_rmean(d_id + d_id = d_id - 1 + end do + ! Pass the current index of the donor since that was filled + ! And also update the status to filled since it now + ! has all memory filled from the donor + this%c_index = donor%c_index + this%filled = .true. + + elseif(.not.this%filled .and. .not.donor%filled) then + + ! Here, neither is completely filled + + if( this%c_index>donor%c_index ) then + + ! In this case, leave all data as the recipient + ! except for where there is donor. Keep the recipient's + ! index and status, since it is larger and should remain unchanged + + r_id = this%c_index + do d_id = donor%c_index,1,-1 + this%mem_rmean(r_id) = this%mem_rmean(r_id)*(recip_wgt) + donor%mem_rmean(d_id)*(1._r8-recip_wgt) + end do + + else + + ! In this case, we do the same thing as the previous + ! clause, but just switch the roles, then + ! copy the donor to the recipient + ! Also transfer the index from the donor, since that was + ! higher and now reflects the filled memory + + d_id = donor%c_index + do r_id = this%c_index,1,-1 + donor%mem_rmean(r_id) = this%mem_rmean(r_id)*(recip_wgt) + donor%mem_rmean(d_id)*(1._r8-recip_wgt) + end do + this%mem_rmean(:) = donor%mem_reman + this%c_index = donor%c_index + end if + + + + end if + + + + ! take the mean of each + + + + + + return + end subroutine FuseRMean + + +end module FatesRunningMeanMod From 915440da268c9c553fe8b1bad161b9878a362144 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 12 Mar 2021 13:36:00 -0500 Subject: [PATCH 02/21] Adding running mean functions --- biogeochem/EDMortalityFunctionsMod.F90 | 5 +- biogeochem/EDPatchDynamicsMod.F90 | 14 +++- biogeochem/EDPhysiologyMod.F90 | 16 ++-- fire/SFMainMod.F90 | 2 +- main/EDInitMod.F90 | 1 - main/EDTypesMod.F90 | 8 +- main/FatesInterfaceMod.F90 | 64 +++++++++++--- main/FatesInterfaceTypesMod.F90 | 17 ++-- main/FatesRunningMeanMod.F90 | 111 +++++++++++++++---------- 9 files changed, 156 insertions(+), 82 deletions(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index fa0b933fc5..f188963767 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -62,7 +62,6 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor real(r8),intent(out) :: smort ! size dependent senescence term real(r8),intent(out) :: asmort ! age dependent senescence term - integer :: ifp real(r8) :: frac ! relativised stored carbohydrate real(r8) :: leaf_c_target ! target leaf biomass kgC real(r8) :: store_c @@ -173,8 +172,8 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor ! Eastern US carbon sink. Glob. Change Biol., 12, 2370-2390, ! doi: 10.1111/j.1365-2486.2006.01254.x - ifp = cohort_in%patchptr%patchno - temp_in_C = bc_in%t_veg24_pa(ifp) - tfrz + temp_in_C = cohort_in%patchptr%tveg24%get_mean() - tfrz + temp_dep_fraction = max(0.0_r8, min(1.0_r8, 1.0_r8 - (temp_in_C - & EDPftvarcon_inst%freezetol(cohort_in%pft))/frost_mort_buffer) ) frmort = EDPftvarcon_inst%mort_scalar_coldstress(cohort_in%pft) * temp_dep_fraction diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 5a49695d15..c7bc37ede7 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -45,6 +45,7 @@ module EDPatchDynamicsMod use FatesInterfaceTypesMod , only : bc_in_type use FatesInterfaceTypesMod , only : hlm_days_per_year use FatesInterfaceTypesMod , only : numpft + use FatesInterfaceTypesMod , only : hlm_stepsize use FatesGlobals , only : endrun => fates_endrun use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue, ifalse @@ -560,7 +561,7 @@ subroutine spawn_patches( currentSite, bc_in) allocate(new_patch_primary) call create_patch(currentSite, new_patch_primary, age, & - site_areadis_primary, primaryforest) + site_areadis_primary, primaryforest ) ! Initialize the litter pools to zero, these ! pools will be populated by looping over the existing patches @@ -1977,6 +1978,8 @@ subroutine create_patch(currentSite, new_patch, age, areap, label) allocate(new_patch%sabs_dif(hlm_numSWb)) allocate(new_patch%fragmentation_scaler(currentSite%nlevsoil)) + allocate(new_patch%tveg24) + call new_patch%tveg24%InitRMean(mem_period=86400._r8,up_period=hlm_stepsize) ! Litter ! Allocate, Zero Fluxes, and Initialize to "unset" values @@ -2478,6 +2481,9 @@ subroutine fuse_2_patches(csite, dp, rp) write(fates_log(),*) 'trying to fuse patches with different anthro_disturbance_label values' call endrun(msg=errMsg(sourcefile, __LINE__)) endif + + ! Weighted mean of the running means + call rp%tveg24%FuseRMean(dp%tveg24,rp%area*inv_sum_area) rp%fuel_eff_moist = (dp%fuel_eff_moist*dp%area + rp%fuel_eff_moist*rp%area) * inv_sum_area rp%livegrass = (dp%livegrass*dp%area + rp%livegrass*rp%area) * inv_sum_area @@ -2814,9 +2820,13 @@ subroutine dealloc_patch(cpatch) deallocate(cpatch%sabs_dir) deallocate(cpatch%sabs_dif) deallocate(cpatch%fragmentation_scaler) - end if + + ! Deallocate any running means + deallocate(cpatch%tveg24%mem) + deallocate(cpatch%tveg24) + return end subroutine dealloc_patch diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 32671e082e..0bf7f32df2 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -748,6 +748,7 @@ subroutine phenology( currentSite, bc_in ) ! Use the following layer index to calculate drought conditions ilayer_swater = minloc(abs(bc_in%z_sisl(:)-dphen_soil_depth),dim=1) + ! Parameter of drought decid leaf loss in mm in top layer...FIX(RF,032414) ! - this is arbitrary and poorly understood. Needs work. ED_ @@ -756,8 +757,8 @@ subroutine phenology( currentSite, bc_in ) temp_in_C = 0._r8 cpatch => CurrentSite%oldest_patch - do while(associated(cpatch)) - temp_in_C = temp_in_C + bc_in%t_veg24_pa(cpatch%patchno)*cpatch%area + do while(associated(cpatch)) + temp_in_C = temp_in_C + cpatch%tveg24%get_mean()*cpatch%area cpatch => cpatch%younger end do temp_in_C = temp_in_C * area_inv - tfrz @@ -2215,7 +2216,6 @@ subroutine fragmentation_scaler( currentPatch, bc_in) logical :: use_century_tfunc = .false. logical :: use_hlm_soil_scalar = .true. ! Use hlm input decomp fraction scalars integer :: j - integer :: ifp ! Index of a FATES Patch "ifp" real(r8) :: t_scalar ! temperature scalar real(r8) :: w_scalar ! moisture scalar real(r8) :: catanf ! hyperbolic temperature function from CENTURY @@ -2226,8 +2226,6 @@ subroutine fragmentation_scaler( currentPatch, bc_in) catanf(t1) = 11.75_r8 +(29.7_r8 / pi) * atan( pi * 0.031_r8 * ( t1 - 15.4_r8 )) catanf_30 = catanf(30._r8) - ifp = currentPatch%patchno - ! Use the hlm temp and moisture decomp fractions by default if ( use_hlm_soil_scalar ) then @@ -2239,17 +2237,17 @@ subroutine fragmentation_scaler( currentPatch, bc_in) if ( .not. use_century_tfunc ) then !calculate rate constant scalar for soil temperature,assuming that the base rate constants !are assigned for non-moisture limiting conditions at 25C. - if (bc_in%t_veg24_pa(ifp) >= tfrz) then - t_scalar = q10_mr**((bc_in%t_veg24_pa(ifp)-(tfrz+25._r8))/10._r8) + if (currentPatch%tveg24%get_mean() >= tfrz) then + t_scalar = q10_mr**((currentPatch%tveg24%get_mean()-(tfrz+25._r8))/10._r8) ! Q10**((t_soisno(c,j)-(tfrz+25._r8))/10._r8) else - t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((bc_in%t_veg24_pa(ifp)-tfrz)/10._r8)) + t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((currentPatch%tveg24%get_mean()-tfrz)/10._r8)) !Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-tfrz)/10._r8) endif else ! original century uses an arctangent function to calculate the ! temperature dependence of decomposition - t_scalar = max(catanf(bc_in%t_veg24_pa(ifp)-tfrz)/catanf_30,0.01_r8) + t_scalar = max(catanf(currentPatch%tveg24%get_mean()-tfrz)/catanf_30,0.01_r8) endif !Moisture Limitations diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index 127dfa43f9..155f6dc88d 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -142,7 +142,7 @@ subroutine fire_danger_index ( currentSite, bc_in) iofp = currentSite%oldest_patch%patchno - temp_in_C = bc_in%t_veg24_pa(iofp) - tfrz + temp_in_C = currentSite%oldest_patch%tveg24%get_mean() - tfrz rainfall = bc_in%precip24_pa(iofp)*sec_per_day rh = bc_in%relhumid24_pa(iofp) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 19a9d2cb00..951c041d60 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -462,7 +462,6 @@ subroutine init_patches( nsites, sites, bc_in) currentPatch%scorch_ht(:) = 0._r8 currentPatch%frac_burnt = 0._r8 currentPatch%burnt_frac_litter(:) = 0._r8 - currentPatch => currentPatch%older enddo enddo diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 110673f94a..94c4d77fcf 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -18,6 +18,7 @@ module EDTypesMod use FatesLitterMod, only : ncwd use FatesConstantsMod, only : n_anthro_disturbance_categories use FatesConstantsMod, only : days_per_year + use FatesRunningMeanMod, only : rmean_type implicit none private ! By default everything is private @@ -411,6 +412,12 @@ module EDTypesMod integer :: anthro_disturbance_label ! patch label for anthropogenic disturbance classification real(r8) :: age_since_anthro_disturbance ! average age for secondary forest since last anthropogenic disturbance + + ! 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) + + ! LEAF ORGANIZATION real(r8) :: pft_agb_profile(maxpft,n_dbh_bins) ! binned above ground biomass, for patch fusion: KgC/m2 real(r8) :: canopy_layer_tlai(nclmax) ! total leaf area index of each canopy layer @@ -684,7 +691,6 @@ module EDTypesMod type (ed_resources_management_type) :: resources_management ! resources_management at the site - ! INDICES real(r8) :: lat ! latitude: degrees real(r8) :: lon ! longitude: degrees diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 48aef9deed..244555616d 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -21,6 +21,8 @@ module FatesInterfaceMod use EDTypesMod , only : do_fates_salinity use EDTypesMod , only : numWaterMem use EDTypesMod , only : numlevsoil_max + use EDTypesMod , only : ed_site_type + use EDTypesMod , only : ed_patch_type use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue,ifalse use FatesConstantsMod , only : nearzero @@ -96,7 +98,8 @@ module FatesInterfaceMod public :: set_bcpconst public :: zero_bcs public :: set_bcs - + public :: UpdateFatesRMeansTStep + contains ! ==================================================================================== @@ -196,13 +199,13 @@ subroutine zero_bcs(fates,s) ! Input boundaries - fates%bc_in(s)%t_veg24_pa(:) = 0.0_r8 - fates%bc_in(s)%precip24_pa(:) = 0.0_r8 - fates%bc_in(s)%relhumid24_pa(:) = 0.0_r8 - fates%bc_in(s)%wind24_pa(:) = 0.0_r8 - fates%bc_in(s)%lightning24(:) = 0.0_r8 fates%bc_in(s)%pop_density(:) = 0.0_r8 + fates%bc_in(s)%precip24_pa(:) = 0.0_r8 + fates%bc_in(s)%relhumid24_pa(:) = 0.0_r8 + fates%bc_in(s)%wind24_pa(:) = 0.0_r8 + fates%bc_in(s)%tveg_pa(:) = 0.0_r8 + fates%bc_in(s)%solad_parb(:,:) = 0.0_r8 fates%bc_in(s)%solai_parb(:,:) = 0.0_r8 fates%bc_in(s)%smp_sl(:) = 0.0_r8 @@ -403,15 +406,13 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats) allocate(bc_in%t_scalar_sisl(nlevsoil_in)) ! Lightning (or successful ignitions) and population density + ! Fire related variables allocate(bc_in%lightning24(maxPatchesPerSite)) allocate(bc_in%pop_density(maxPatchesPerSite)) - - ! Vegetation Dynamics - allocate(bc_in%t_veg24_pa(maxPatchesPerSite)) - allocate(bc_in%wind24_pa(maxPatchesPerSite)) allocate(bc_in%relhumid24_pa(maxPatchesPerSite)) allocate(bc_in%precip24_pa(maxPatchesPerSite)) + allocate(bc_in%tveg_pa(maxPatchesPerSite)) ! Radiation allocate(bc_in%solad_parb(maxPatchesPerSite,hlm_numSWb)) @@ -429,6 +430,8 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats) allocate(bc_in%salinity_sl(nlevsoil_in)) endif + + ! Photosynthesis allocate(bc_in%filter_photo_pa(maxPatchesPerSite)) allocate(bc_in%dayl_factor_pa(maxPatchesPerSite)) @@ -1153,6 +1156,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_numlevgrnd = unset_int hlm_name = 'unset' hlm_hio_ignore_val = unset_double + hlm_stepsize = unset_double hlm_masterproc = unset_int hlm_ipedof = unset_int hlm_nu_com = 'unset' @@ -1360,6 +1364,14 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if + if( abs(hlm_stepsize-unset_double) sites(s)%oldest_patch + do while(associated(cpatch)) + ifp=ifp+1 + call cpatch%tveg24%UpdateRMean(bc_in(s)%tveg_pa(ifp)) + cpatch => cpatch%younger + enddo + enddo + + return + end subroutine UpdateFatesRMeansTStep + end module FatesInterfaceMod diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index 5b069709cc..67b36d7e6a 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -69,7 +69,9 @@ module FatesInterfaceTypesMod ! 0: none ! 1: p is on - + real(r8), public :: hlm_stepsize ! The step-size of the host land model (s) + ! moreover, this is the shortest main-model timestep + ! at which fates will be called on the main model integration loop real(r8), public :: hlm_hio_ignore_val ! This value can be flushed to history ! diagnostics, such that the @@ -335,23 +337,19 @@ module FatesInterfaceTypesMod real(r8),allocatable :: w_scalar_sisl(:) ! fraction by which decomposition is limited by moisture availability real(r8),allocatable :: t_scalar_sisl(:) ! fraction by which decomposition is limited by temperature - - ! Vegetation Dynamics - ! --------------------------------------------------------------------------------- + ! Fire Model ! 24-hour lightning or ignitions [#/km2/day] real(r8),allocatable :: lightning24(:) ! Population density [#/km2] real(r8),allocatable :: pop_density(:) - - ! Patch 24 hour vegetation temperature [K] - real(r8),allocatable :: t_veg24_pa(:) - ! Fire Model - ! Average precipitation over the last 24 hours [mm/s] real(r8), allocatable :: precip24_pa(:) + + ! Patch Vegetation temperature (K) + real(r8),allocatable :: tveg_pa(:) ! Average relative humidity over past 24 hours [-] real(r8), allocatable :: relhumid24_pa(:) @@ -359,7 +357,6 @@ module FatesInterfaceTypesMod ! Patch 24-hour running mean of wind (m/s ?) real(r8), allocatable :: wind24_pa(:) - ! Radiation variables for calculating sun/shade fractions ! --------------------------------------------------------------------------------- diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningMeanMod.F90 index f7b6ca611d..27e27ff573 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningMeanMod.F90 @@ -1,50 +1,74 @@ module FatesRunningMeanMod - use FatesConstantsMod, only : nearzero - 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 + use FatesConstantsMod, only : nearzero + use FatesConstantsMod, only : r8 => fates_r8 + 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 + use FatesGlobals, only : fates_log + implicit none + private integer, parameter :: maxlen_varname = 8 type, public :: rmean_type - real(r8), allocatable(:) :: mem_rmean ! Array storing the memory of the mean - real(r8) :: c_rmean ! The current mean value from the + real(r8), allocatable :: mem(:) ! Array storing the memory of the mean + real(r8) :: c_mean ! The current mean value from the ! available memory, as of the last update integer :: c_index ! The current memory index as per the ! last update integer :: n_mem ! Total number of memory indices logical :: filled ! Has enough time elapsed so all memory filled? - character(len=maxlen_varname) :: var_name ! A short name for this variable - ! for diagnostic purposes + !character(len=maxlen_varname) :: var_name ! A short name for this variable + ! ! for diagnostic purposes contains - + + procedure :: get_mean procedure :: InitRMean procedure :: UpdateRMean procedure :: FuseRMean end type rmean_type + + character(len=*), parameter, private :: sourcefile = & + __FILE__ contains - subroutine InitRMean(this,name,mem_period,up_period) + function get_mean(this) class(rmean_type) :: this - character(len=maxlen_varname) :: name ! The name of the new variable + real(r8) :: get_mean + + if(this%c_index == 0) then + write(fates_log(), *) 'attempting to get a running mean from a variable' + write(fates_log(), *) 'that has not experienced any time to accumluate memory' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + get_mean = this%c_mean + + end function get_mean + + + subroutine InitRMean(this,mem_period,up_period) !,init_value) + + class(rmean_type) :: this + !character(len=maxlen_varname) :: name ! The name of the new variable real(r8) :: mem_period ! The period length in seconds that must be remembered real(r8) :: up_period ! The period length in seconds that memory is updated - ! (i.e. the resolution of the memory) - - this%name = name + ! (i.e. the resolution of the memory) + !real(r8) :: init_value ! The first value to put in memory + + !this%name = name this%n_mem = nint(mem_period/up_period) if( abs(real(this%n_mem,r8)-mem_period/up_period) > nearzero ) then - write(fates_log(), *) 'While defining a running mean variable: ',this%var_name + write(fates_log(), *) 'While defining a running mean variable: '!,this%var_name write(fates_log(), *) 'an update and total memory period was specified' write(fates_log(), *) 'where the update period is not an exact fraction of the period' write(fates_log(), *) 'mem_period: ',mem_period @@ -54,17 +78,20 @@ subroutine InitRMean(this,name,mem_period,up_period) end if ! Otherwise we allocate - allocate(this%mem_rmean(this%n_mem)) + allocate(this%mem(this%n_mem)) ! Initialize with nonsense numbers - this%mem_rmean(:) = nan + this%mem(:) = nan ! There are no memory spots filled with valid datapoints - this%filled = .false + this%filled = .false. - ! The current index of the memory is zero + ! The current index of the memory is one this%c_index = 0 + !this%mem(1) = init_value + + return end subroutine InitRMean @@ -76,7 +103,7 @@ subroutine UpdateRMean(this, new_value) real(r8) :: new_value ! The newest value added to the running mean this%c_index = this%c_index + 1 - + ! Update the index of the memory array ! and, determine if we have filled the memory yet ! If this index is greater than our memory slots, @@ -87,19 +114,18 @@ subroutine UpdateRMean(this, new_value) if(this%c_index>this%n_mem) this%c_index = 1 - this%mem_rmean(this%c_index) = new_value + this%mem(this%c_index) = new_value ! Update the running mean value. It will return a value ! if we have not filled in all the memory slots. To do this ! it will take a mean over what is available if(this%filled) then - this%c_rmean = sum(this%mem_rmean)/real(this%n_mem,r8) + this%c_mean = sum(this%mem)/real(this%n_mem,r8) else - this%c_rmean = sum(this%mem_rmean(1:this%c_index))/real(this%c_index,r8) + this%c_mean = sum(this%mem(1:this%c_index))/real(this%c_index,r8) end if - - + return end subroutine UpdateRmean @@ -124,7 +150,7 @@ subroutine FuseRMean(this,donor,recip_wgt) if(this%n_mem .ne. donor%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 + write(fates_log(), *) 'of two running mean variables: '!,this%name,donor%name call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -134,9 +160,9 @@ subroutine FuseRMean(this,donor,recip_wgt) d_id = donor%c_index do r_id = 1,this%n_mem - this%mem_rmean(r_id) = this%mem_rmean(r_id)*(recip_wgt) + & - donor%mem_rmean(d_id)*(1._r8-recip_wgt) - d_id = d_id+ipos + this%mem(r_id) = this%mem(r_id)*(recip_wgt) + & + donor%mem(d_id)*(1._r8-recip_wgt) + d_id = d_id+1 if(d_id>donor%n_mem) d_id = 1 end do @@ -151,7 +177,7 @@ subroutine FuseRMean(this,donor,recip_wgt) r_id = this%c_index do d_id = donor%c_index,1,-1 - this%mem_rmean(r_id) = this%mem_rmean(r_id)*(recip_wgt) + donor%mem_rmean(d_id)*(1._r8-recip_wgt) + this%mem(r_id) = this%mem(r_id)*(recip_wgt) + donor%mem(d_id)*(1._r8-recip_wgt) r_id = r_id - 1 if(r_id==0) r_id = this%n_mem end do @@ -167,13 +193,13 @@ subroutine FuseRMean(this,donor,recip_wgt) d_id = donor%c_index do r_id = this%c_index,1,-1 - this%mem_rmean(r_id) = this%mem_rmean(r_id)*(recip_wgt) + donor%mem_rmean(d_id)*(1._r8-recip_wgt) + this%mem(r_id) = this%mem(r_id)*(recip_wgt) + donor%mem(d_id)*(1._r8-recip_wgt) d_id = d_id - 1 if(d_id==0) d_id = donor%n_mem end do do r_id = this%n_mem,this%c_index+1,-1 - this%mem_rmean(r_id) = donor%mem_rmean(d_id + this%mem(r_id) = donor%mem(d_id) d_id = d_id - 1 end do ! Pass the current index of the donor since that was filled @@ -194,7 +220,7 @@ subroutine FuseRMean(this,donor,recip_wgt) r_id = this%c_index do d_id = donor%c_index,1,-1 - this%mem_rmean(r_id) = this%mem_rmean(r_id)*(recip_wgt) + donor%mem_rmean(d_id)*(1._r8-recip_wgt) + this%mem(r_id) = this%mem(r_id)*(recip_wgt) + donor%mem(d_id)*(1._r8-recip_wgt) end do else @@ -207,22 +233,21 @@ subroutine FuseRMean(this,donor,recip_wgt) d_id = donor%c_index do r_id = this%c_index,1,-1 - donor%mem_rmean(r_id) = this%mem_rmean(r_id)*(recip_wgt) + donor%mem_rmean(d_id)*(1._r8-recip_wgt) + donor%mem(r_id) = this%mem(r_id)*(recip_wgt) + donor%mem(d_id)*(1._r8-recip_wgt) end do - this%mem_rmean(:) = donor%mem_reman + this%mem(:) = donor%mem this%c_index = donor%c_index end if - - end if - - - ! take the mean of each - - - + ! Update the mean based on the fusion + if(this%filled) then + this%c_mean = sum(this%mem)/real(this%n_mem,r8) + else + if(this%c_index>0) this%c_mean = & + sum(this%mem(1:this%c_index))/real(this%c_index,r8) + end if return From e6ea9b24f0651b22a560c1b40583c0acd586390a Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 19 Mar 2021 15:12:00 -0400 Subject: [PATCH 03/21] Re-worked the running mean functions: 1) abandoned maintaining all memory (ie simple moving average) in a window in favor of the exponential moving window, 2) added fixed window capabilities. --- biogeochem/EDPatchDynamicsMod.F90 | 6 +- main/FatesInterfaceMod.F90 | 16 +- main/FatesRunningMeanMod.F90 | 371 +++++++++++++++++------------- 3 files changed, 223 insertions(+), 170 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index c7bc37ede7..0b45d1b12c 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -85,7 +85,8 @@ module EDPatchDynamicsMod use SFParamsMod, only : SF_VAL_CWD_FRAC use EDParamsMod, only : logging_event_code use EDParamsMod, only : logging_export_frac - + use FatesRunningMeanMod, only : ema_24hr + ! CIME globals use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) use shr_log_mod , only : errMsg => shr_log_errMsg @@ -1979,7 +1980,7 @@ subroutine create_patch(currentSite, new_patch, age, areap, label) allocate(new_patch%fragmentation_scaler(currentSite%nlevsoil)) allocate(new_patch%tveg24) - call new_patch%tveg24%InitRMean(mem_period=86400._r8,up_period=hlm_stepsize) + call new_patch%tveg24%InitRMean(ema_24hr) ! No initial value ! Litter ! Allocate, Zero Fluxes, and Initialize to "unset" values @@ -2824,7 +2825,6 @@ subroutine dealloc_patch(cpatch) ! Deallocate any running means - deallocate(cpatch%tveg24%mem) deallocate(cpatch%tveg24) return diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 244555616d..d2656ef5ce 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -26,6 +26,7 @@ module FatesInterfaceMod use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue,ifalse use FatesConstantsMod , only : nearzero + use FatesConstantsMod , only : sec_per_day use FatesGlobals , only : fates_global_verbose use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun @@ -68,8 +69,11 @@ module FatesInterfaceMod use PRTInitParamsFatesMod , only : PRTCheckParams use PRTAllometricCarbonMod , only : InitPRTGlobalAllometricCarbon use PRTAllometricCNPMod , only : InitPRTGlobalAllometricCNP - - + use FatesRunningMeanMod , only : ema_24hr + use FatesRunningMeanMod , only : fixed_24hr + use FatesRunningMeanMod , only : moving_ema_window + use FatesRunningMeanMod , only : fixed_window + ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) @@ -797,6 +801,14 @@ subroutine SetFatesGlobalElements(use_fates) ! These will not be used if use_ed or use_fates is false call fates_history_maps() + + ! Instantiate the time-averaging method globals + allocate(ema_24hr) + call ema_24hr%define(sec_per_day, hlm_stepsize, moving_ema_window) + allocate(fixed_24hr) + call fixed_24hr%define(sec_per_day, hlm_stepsize, fixed_window) + + else ! If we are not using FATES, the cohort dimension is still diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningMeanMod.F90 index 27e27ff573..60fa1b0850 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningMeanMod.F90 @@ -12,18 +12,64 @@ module FatesRunningMeanMod private integer, parameter :: maxlen_varname = 8 + + ! These are flags that specify how the averaging window works. + ! Moving windows (default) can have an arbitrary size and update frequency) + ! and it is technically never reset, it just averages indefinitely. + ! But hourly, six-hourly, daily, monthly and yearly windows have pre-set + ! window sizes associated with their namesake, and more importantly, they + ! are zero'd at the beginning of the interval, and get equal average weighting + ! over their construction period. + + + integer, public, parameter :: moving_ema_window = 0 + integer, public, parameter :: fixed_window = 1 + + ! This type defines a type of mean. It does not + ! define the variable, but it defines how + ! often it is updated, how long its + ! memory period is, and if it should be zero'd + ! These are globally defined on the proc. + + type, public :: rmean_def_type + + real(r8) :: mem_period ! The total integration period (s) + real(r8) :: up_period ! The period between updates (s) + integer :: n_mem ! How many updates per integration period? + integer :: method ! Is this a fixed or moving window? + + contains + + procedure :: define + + end type rmean_def_type + + + ! This holds the time varying information for the mean + ! which is instantiated on sites, patches, and cohorts type, public :: rmean_type - real(r8), allocatable :: mem(:) ! Array storing the memory of the mean - real(r8) :: c_mean ! The current mean value from the - ! available memory, as of the last update - integer :: c_index ! The current memory index as per the - ! last update - integer :: n_mem ! Total number of memory indices - logical :: filled ! Has enough time elapsed so all memory filled? - !character(len=maxlen_varname) :: var_name ! A short name for this variable - ! ! for diagnostic purposes + real(r8) :: c_mean ! The current mean value, if this + ! is a moving window, its is the mean. + ! If this is a fixed window, it is only a partial mean + ! as the value uses equal update weights and is not + ! necessarily fully constructed. + + real(r8) :: l_mean ! The latest reportable mean value + ! this value is actually the same + ! as c_mean for moving windows, and for fixed windows + ! it is the mean value when the time integratino window + ! last completed. + + integer :: c_index ! The number of values that have + ! been added to the mean so far + ! once this is >= n_mem then + ! the ema weight hits its cap + + ! This points to the global structure that + ! defines the nature of this mean/avg + type(rmean_def_type), pointer :: def_type contains @@ -34,42 +80,31 @@ module FatesRunningMeanMod end type rmean_type + character(len=*), parameter, private :: sourcefile = & __FILE__ - -contains - function get_mean(this) + ! Define the time methods that we want to have available to us + + class(rmean_def_type), public, pointer :: ema_24hr + class(rmean_def_type), public, pointer :: fixed_24hr + +contains - class(rmean_type) :: this - real(r8) :: get_mean - if(this%c_index == 0) then - write(fates_log(), *) 'attempting to get a running mean from a variable' - write(fates_log(), *) 'that has not experienced any time to accumluate memory' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - get_mean = this%c_mean + subroutine define(this,mem_period,up_period,method) - end function get_mean - + class(rmean_def_type) :: this - subroutine InitRMean(this,mem_period,up_period) !,init_value) + real(r8),intent(in) :: mem_period + real(r8),intent(in) :: up_period + integer,intent(in) :: method - class(rmean_type) :: this - !character(len=maxlen_varname) :: name ! The name of the new variable - real(r8) :: mem_period ! The period length in seconds that must be remembered - real(r8) :: up_period ! The period length in seconds that memory is updated - ! (i.e. the resolution of the memory) - !real(r8) :: init_value ! The first value to put in memory - - !this%name = name - this%n_mem = nint(mem_period/up_period) - - if( abs(real(this%n_mem,r8)-mem_period/up_period) > nearzero ) then - write(fates_log(), *) 'While defining a running mean variable: '!,this%var_name - write(fates_log(), *) 'an update and total memory period was specified' + ! Check the memory and update periods + if( abs(nint(mem_period/up_period)-mem_period/up_period) > nearzero ) then + write(fates_log(), *) 'While defining a running mean definition' + write(fates_log(), *) 'an update and memory period was specified' write(fates_log(), *) 'where the update period is not an exact fraction of the period' write(fates_log(), *) 'mem_period: ',mem_period write(fates_log(), *) 'up_period: ',up_period @@ -77,20 +112,85 @@ subroutine InitRMean(this,mem_period,up_period) !,init_value) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - ! Otherwise we allocate - allocate(this%mem(this%n_mem)) + this%mem_period = mem_period + this%up_period = up_period + this%method = method + this%n_mem = nint(mem_period/up_period) + + return + end subroutine define + + ! ===================================================================================== + + function get_mean(this) + + class(rmean_type) :: this + real(r8) :: get_mean - ! Initialize with nonsense numbers - this%mem(:) = nan + if(this%def_type%method .eq. moving_ema_window) then + if(this%c_index == 0) then + write(fates_log(), *) 'attempting to get a running mean from a variable' + write(fates_log(), *) 'that has been given a value yet' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + else + if(this%c_index .ne. this%def_type%n_mem)then + write(fates_log(), *) 'attempting to get a mean over a fixed window' + write(fates_log(), *) 'at a time where the window has not completed' + write(fates_log(), *) 'its cycle yet' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + get_mean = this%l_mean + + end function get_mean + + ! ===================================================================================== + + subroutine InitRMean(this,rmean_def,init_value,init_offset) + + 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 - ! There are no memory spots filled with valid datapoints - this%filled = .false. + ! 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. - ! The current index of the memory is one - this%c_index = 0 + ! 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%mem(1) = init_value + ! Point to the averaging type + this%def_type => rmean_def + if(this%def_type%method .eq. fixed_window) then + + if(.not.(present(init_offset).and.present(init_value)) )then + write(fates_log(), *) 'when initializing a temporal mean on a fixed window' + write(fates_log(), *) 'there must be an initial value and a time offset' + write(fates_log(), *) 'specified.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + this%c_index = modulo(nint(init_offset/rmean_def%up_period)+1,rmean_def%n_mem) + this%c_mean = real(this%c_index,r8)/real(rmean_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 + this%c_mean = init_value + this%l_mean = init_value + this%c_index = 1 + else + this%c_mean = nan + this%l_mean = nan + this%c_index = 0 + end if + + end if return end subroutine InitRMean @@ -101,29 +201,41 @@ subroutine UpdateRMean(this, new_value) class(rmean_type) :: this real(r8) :: new_value ! The newest value added to the running mean + real(r8) :: wgt - this%c_index = this%c_index + 1 - - ! Update the index of the memory array - ! and, determine if we have filled the memory yet - ! If this index is greater than our memory slots, - ! go back to the first - if(this%c_index==this%n_mem) then - this%filled = .true. - end if + if(this%def_type%method.eq.moving_ema_window) then - if(this%c_index>this%n_mem) this%c_index = 1 - - this%mem(this%c_index) = new_value + 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 - ! Update the running mean value. It will return a value - ! if we have not filled in all the memory slots. To do this - ! it will take a mean over what is available + ! If the last time we updated we had hit the + ! end of the averaging memory period, and + ! we are not using an indefinite running + ! average, then zero things out - if(this%filled) then - this%c_mean = sum(this%mem)/real(this%n_mem,r8) - else - this%c_mean = sum(this%mem(1:this%c_index))/real(this%c_index,r8) + if(this%c_index == this%def_type%n_mem) then + this%c_mean = 0._r8 + this%c_index = 0 + end if + + 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 + + if(this%c_index == this%def_type%n_mem) then + this%l_mean = this%c_mean + end if + end if return @@ -133,123 +245,52 @@ end subroutine UpdateRmean subroutine FuseRMean(this,donor,recip_wgt) - ! When fusing the running mean of two entities, it is possible that they - ! may have a different amount of memory spaces filled (at least in FATES). This - ! is typical for newly created patches or cohorts, that litteraly just spawned - ! So what generally happens is we walk backwards from the current memory index - ! of both and take means where we can. In places where one entity has more - ! memory than the other, than we just use the value from the one that is there + ! Rules for fusion: + ! If both entities have valid means already, then you simply use the + ! weight provided to combine them. + ! If this is a moving average, then update the index to be the larger of + ! the two. + ! if this is a fixed window, check that the index is the same between + ! both. + class(rmean_type) :: this class(rmean_type), pointer :: donor real(r8),intent(in) :: recip_wgt ! Weighting factor for recipient (0-1) - integer :: r_id ! Loop index counter for the recipient (this) - integer :: d_id ! Loop index counter for the donor - - - if(this%n_mem .ne. donor%n_mem) then + 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%filled .and. donor%filled) then - - ! Both are filled, take averages of each and be sure to use relative positions - - d_id = donor%c_index - do r_id = 1,this%n_mem - this%mem(r_id) = this%mem(r_id)*(recip_wgt) + & - donor%mem(d_id)*(1._r8-recip_wgt) - d_id = d_id+1 - if(d_id>donor%n_mem) d_id = 1 - end do - - elseif(this%filled .and. .not.donor%filled) then - - ! Here, we have only partial memory from the donor - ! we we iterate through the donor's memory - ! and average between the two in those spaces. Then - ! we leave the rest untouched because we accept - ! the values from the recipient. Also, keep the - ! recipient's index. - - r_id = this%c_index - do d_id = donor%c_index,1,-1 - this%mem(r_id) = this%mem(r_id)*(recip_wgt) + donor%mem(d_id)*(1._r8-recip_wgt) - r_id = r_id - 1 - if(r_id==0) r_id = this%n_mem - end do - - - elseif(.not.this%filled .and. donor%filled) then - - ! Here we only have partial memory from the recipient - ! so we iterate through the recipient's memory - ! and average between the two. Then we copy - ! over the values from the donor for the indices that we - ! didn't average because the donor has valid values - - d_id = donor%c_index - do r_id = this%c_index,1,-1 - this%mem(r_id) = this%mem(r_id)*(recip_wgt) + donor%mem(d_id)*(1._r8-recip_wgt) - d_id = d_id - 1 - if(d_id==0) d_id = donor%n_mem - end do - - do r_id = this%n_mem,this%c_index+1,-1 - this%mem(r_id) = donor%mem(d_id) - d_id = d_id - 1 - end do - ! Pass the current index of the donor since that was filled - ! And also update the status to filled since it now - ! has all memory filled from the donor - this%c_index = donor%c_index - this%filled = .true. - - elseif(.not.this%filled .and. .not.donor%filled) then + if(this%def_type%method .eq. fixed_window ) then + 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?' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if - ! Here, neither is completely filled - - if( this%c_index>donor%c_index ) then - - ! In this case, leave all data as the recipient - ! except for where there is donor. Keep the recipient's - ! index and status, since it is larger and should remain unchanged - - r_id = this%c_index - do d_id = donor%c_index,1,-1 - this%mem(r_id) = this%mem(r_id)*(recip_wgt) + donor%mem(d_id)*(1._r8-recip_wgt) - end do - - else + ! 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 - ! In this case, we do the same thing as the previous - ! clause, but just switch the roles, then - ! copy the donor to the recipient - ! Also transfer the index from the donor, since that was - ! higher and now reflects the filled memory - - d_id = donor%c_index - do r_id = this%c_index,1,-1 - donor%mem(r_id) = this%mem(r_id)*(recip_wgt) + donor%mem(d_id)*(1._r8-recip_wgt) - end do - this%mem(:) = donor%mem + if(this%c_index==0) then + this%c_mean = donor%c_mean + this%l_mean = donor%l_mean 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) + ! Update the index to the larger of the two + this%c_index = max(this%c_index,donor%c_index) end if - - end if - ! Update the mean based on the fusion - if(this%filled) then - this%c_mean = sum(this%mem)/real(this%n_mem,r8) - else - if(this%c_index>0) this%c_mean = & - sum(this%mem(1:this%c_index))/real(this%c_index,r8) end if - - + return end subroutine FuseRMean From 7e7b075c446205796785dea8960064fdf6ca1b7b Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 30 Mar 2021 11:23:00 -0400 Subject: [PATCH 04/21] Added a copy function, a history variable and restart capabilities to running mean functions --- biogeochem/EDMortalityFunctionsMod.F90 | 2 +- biogeochem/EDPatchDynamicsMod.F90 | 8 ++ biogeochem/EDPhysiologyMod.F90 | 10 +-- fire/SFMainMod.F90 | 2 +- main/FatesHistoryInterfaceMod.F90 | 27 ++++++- main/FatesRestartInterfaceMod.F90 | 108 ++++++++++++++++++++++++- main/FatesRunningMeanMod.F90 | 38 +++++++-- 7 files changed, 178 insertions(+), 17 deletions(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index f188963767..ff788d8ef2 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -172,7 +172,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor ! Eastern US carbon sink. Glob. Change Biol., 12, 2370-2390, ! doi: 10.1111/j.1365-2486.2006.01254.x - temp_in_C = cohort_in%patchptr%tveg24%get_mean() - tfrz + temp_in_C = cohort_in%patchptr%tveg24%GetMean() - tfrz temp_dep_fraction = max(0.0_r8, min(1.0_r8, 1.0_r8 - (temp_in_C - & EDPftvarcon_inst%freezetol(cohort_in%pft))/frost_mort_buffer) ) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 0b45d1b12c..25c5decc79 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -668,6 +668,14 @@ subroutine spawn_patches( currentSite, bc_in) call mortality_litter_fluxes(currentSite, currentPatch, new_patch, patch_site_areadis) endif + + ! Copy any means or timers from the original patch to the new patch + ! These values will inherit all info from the original patch + ! -------------------------------------------------------------------------- + call new_patch%tveg24%CopyFromDonor(currentPatch%tveg24) + + + ! -------------------------------------------------------------------------- ! The newly formed patch from disturbance (new_patch), has now been given ! some litter from dead plants and pre-existing litter from the donor patches. diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 0bf7f32df2..10a2e98149 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -758,7 +758,7 @@ subroutine phenology( currentSite, bc_in ) temp_in_C = 0._r8 cpatch => CurrentSite%oldest_patch do while(associated(cpatch)) - temp_in_C = temp_in_C + cpatch%tveg24%get_mean()*cpatch%area + temp_in_C = temp_in_C + cpatch%tveg24%GetMean()*cpatch%area cpatch => cpatch%younger end do temp_in_C = temp_in_C * area_inv - tfrz @@ -2237,17 +2237,17 @@ subroutine fragmentation_scaler( currentPatch, bc_in) if ( .not. use_century_tfunc ) then !calculate rate constant scalar for soil temperature,assuming that the base rate constants !are assigned for non-moisture limiting conditions at 25C. - if (currentPatch%tveg24%get_mean() >= tfrz) then - t_scalar = q10_mr**((currentPatch%tveg24%get_mean()-(tfrz+25._r8))/10._r8) + if (currentPatch%tveg24%GetMean() >= tfrz) then + t_scalar = q10_mr**((currentPatch%tveg24%GetMean()-(tfrz+25._r8))/10._r8) ! Q10**((t_soisno(c,j)-(tfrz+25._r8))/10._r8) else - t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((currentPatch%tveg24%get_mean()-tfrz)/10._r8)) + t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((currentPatch%tveg24%GetMean()-tfrz)/10._r8)) !Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-tfrz)/10._r8) endif else ! original century uses an arctangent function to calculate the ! temperature dependence of decomposition - t_scalar = max(catanf(currentPatch%tveg24%get_mean()-tfrz)/catanf_30,0.01_r8) + t_scalar = max(catanf(currentPatch%tveg24%GetMean()-tfrz)/catanf_30,0.01_r8) endif !Moisture Limitations diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index 155f6dc88d..9e894df08f 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -142,7 +142,7 @@ subroutine fire_danger_index ( currentSite, bc_in) iofp = currentSite%oldest_patch%patchno - temp_in_C = currentSite%oldest_patch%tveg24%get_mean() - tfrz + temp_in_C = currentSite%oldest_patch%tveg24%GetMean() - tfrz rainfall = bc_in%precip24_pa(iofp)*sec_per_day rh = bc_in%relhumid24_pa(iofp) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index ed98301ee4..fa90a93080 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -527,7 +527,8 @@ module FatesHistoryInterfaceMod ! integer :: ih_fire_rate_of_spread_front_si_age integer :: ih_fire_intensity_si_age integer :: ih_fire_sum_fuel_si_age - + integer :: ih_tveg24_si_age + ! indices to (site x height) variables integer :: ih_canopy_height_dist_si_height integer :: ih_leaf_height_dist_si_height @@ -1763,7 +1764,6 @@ subroutine update_history_dyn(this,nc,nsites,sites) real(r8) :: npp_partition_error ! a check that the NPP partitions sum to carbon allocation real(r8) :: frac_canopy_in_bin ! fraction of a leaf's canopy that is within a given height bin real(r8) :: binbottom,bintop ! edges of height bins - real(r8) :: gpp_cached ! variable used to cache gpp value in previous time step; for C13 discrimination ! The following are all carbon states, turnover and net allocation flux variables @@ -2008,6 +2008,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! hio_fire_rate_of_spread_front_si_age => this%hvars(ih_fire_rate_of_spread_front_si_age)%r82d, & hio_fire_intensity_si_age => this%hvars(ih_fire_intensity_si_age)%r82d, & hio_fire_sum_fuel_si_age => this%hvars(ih_fire_sum_fuel_si_age)%r82d, & + hio_tveg24_si_age => this%hvars(ih_tveg24_si_age)%r82d, & hio_burnt_frac_litter_si_fuel => this%hvars(ih_burnt_frac_litter_si_fuel)%r82d, & hio_fuel_amount_si_fuel => this%hvars(ih_fuel_amount_si_fuel)%r82d, & hio_fuel_amount_age_fuel => this%hvars(ih_fuel_amount_age_fuel)%r82d, & @@ -2214,7 +2215,13 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_fire_sum_fuel_si_age(io_si, cpatch%age_class) = hio_fire_sum_fuel_si_age(io_si, cpatch%age_class) + & cpatch%sum_fuel * g_per_kg * cpatch%area * AREA_INV - + + if(cpatch%tveg24%c_index>0) then + hio_tveg24_si_age(io_si, cpatch%age_class) = & + hio_tveg24_si_age(io_si, cpatch%age_class) + & + cpatch%tveg24%GetMean()*cpatch%area + end if + if(associated(cpatch%tallest))then hio_trimming_si(io_si) = hio_trimming_si(io_si) + cpatch%tallest%canopy_trim * cpatch%area * AREA_INV endif @@ -2874,6 +2881,9 @@ subroutine update_history_dyn(this,nc,nsites,sites) if (hio_area_si_age(io_si, ipa2) .gt. tiny) then hio_lai_si_age(io_si, ipa2) = hio_lai_si_age(io_si, ipa2) / (hio_area_si_age(io_si, ipa2)*AREA) hio_ncl_si_age(io_si, ipa2) = hio_ncl_si_age(io_si, ipa2) / (hio_area_si_age(io_si, ipa2)*AREA) + + hio_tveg24_si_age(io_si, ipa2) = hio_tveg24_si_age(io_si, ipa2) / (hio_area_si_age(io_si, ipa2)*AREA) + do i_pft = 1, numpft iagepft = ipa2 + (i_pft-1) * nlevage hio_scorch_height_si_agepft(io_si, iagepft) = & @@ -3498,6 +3508,8 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) hio_c_lblayer_si(io_si) = hio_c_lblayer_si(io_si) + & cpatch%c_lblayer * cpatch%total_canopy_area + + ccohort => cpatch%shortest do while(associated(ccohort)) @@ -4549,6 +4561,15 @@ subroutine define_history_vars(this, initialize_variables) ivar=ivar, initialize=initialize_variables, index = ih_burnt_frac_litter_si_fuel ) + ! Running means + + call this%set_history_var(vname='TVEG24_AGE', units='Kelvin', & + long='24-hr running mean vegetation temperature by patch age', & + use_default='active', & + avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', flushval=hlm_hio_ignore_val, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_tveg24_si_age ) + + ! Litter Variables call this%set_history_var(vname='LITTER_IN', units='gC m-2 s-1', & diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 0b621dec0a..8a608d70d3 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -37,7 +37,7 @@ module FatesRestartInterfaceMod use FatesLitterMod, only : ndcmpy use PRTGenericMod, only : prt_global use PRTGenericMod, only : num_elements - + use FatesRunningMeanMod, only : rmean_type ! CIME GLOBALS use shr_log_mod , only : errMsg => shr_log_errMsg @@ -135,6 +135,8 @@ module FatesRestartInterfaceMod integer :: ir_gnd_alb_dif_pasb integer :: ir_gnd_alb_dir_pasb + ! Running Means + integer :: ir_tveg24patch_pa integer :: ir_ddbhdt_co integer :: ir_resp_tstep_co @@ -291,7 +293,9 @@ module FatesRestartInterfaceMod procedure, private :: GetCohortRealVector procedure, private :: SetCohortRealVector procedure, private :: RegisterCohortVector - + procedure, private :: DefineRMeanRestartVar + procedure, private :: GetRMeanRestartVar + procedure, private :: SetRMeanRestartVar end type fates_restart_interface_type @@ -1205,6 +1209,10 @@ subroutine define_restart_vars(this, initialize_variables) hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_promcflux_si ) + call this%DefineRMeanRestartVar(vname='fates_tveg24patch',vtype=cohort_r8, & + long_name='24-hour patch veg temp', & + units='K', initialize=initialize_variables,ivar=ivar, index = ir_tveg24patch_pa) + ! Register all of the PRT states and fluxes @@ -1218,7 +1226,90 @@ subroutine define_restart_vars(this, initialize_variables) this%num_restart_vars_ = ivar end subroutine define_restart_vars + + ! ===================================================================================== + + subroutine DefineRMeanRestartVar(this,vname,vtype,long_name,units,initialize,ivar,index) + + class(fates_restart_interface_type) :: this + character(len=*),intent(in) :: vname + character(len=*),intent(in) :: vtype + character(len=*),intent(in) :: long_name + character(len=*),intent(in) :: units + logical, intent(in) :: initialize + integer,intent(inout) :: ivar + integer,intent(inout) :: index + + integer :: dummy_index + + call this%set_restart_var(vname= trim(vname)//'_cmean', vtype=vtype, & + long_name=long_name//' current mean', & + units=units, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = index ) + + call this%set_restart_var(vname= trim(vname)//'_lmean', vtype=vtype, & + long_name=long_name//' latest mean', & + 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, & + hlms='CLM:ALM', initialize=initialize, ivar=ivar, index = dummy_index ) + + + return + end subroutine DefineRMeanRestartVar + + + ! ===================================================================================== + + subroutine GetRMeanRestartVar(this, rmean_var, ir_var_index, position_index) + + class(fates_restart_interface_type) , intent(inout) :: this + class(rmean_type), intent(inout) :: rmean_var + + integer,intent(in) :: ir_var_index + integer,intent(in) :: position_index + + integer :: i_pos ! vector position loop index + integer :: ir_pos_var ! global variable index + + + rmean_var%c_mean = this%rvars(ir_var_index)%r81d(position_index) + + rmean_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)) + + return + end subroutine GetRMeanRestartVar + + ! ======================================================================================= + + subroutine SetRMeanRestartVar(this, rmean_var, ir_var_index, position_index) + + class(fates_restart_interface_type) , intent(inout) :: this + class(rmean_type), intent(inout) :: rmean_var + + integer,intent(in) :: ir_var_index + integer,intent(in) :: 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) + + return + end subroutine SetRMeanRestartVar + + + ! ===================================================================================== subroutine DefinePRTRestartVars(this,initialize_variables,ivar) @@ -1410,6 +1501,12 @@ subroutine RegisterCohortVector(this,symbol_base, vtype, long_name_base, & end subroutine RegisterCohortVector + + + + + + ! ===================================================================================== subroutine GetCohortRealVector(this, state_vector, len_state_vector, & @@ -1918,6 +2015,9 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_patchdistturbcat_pa(io_idx_co_1st) = cpatch%anthro_disturbance_label rio_agesinceanthrodist_pa(io_idx_co_1st) = cpatch%age_since_anthro_disturbance rio_area_pa(io_idx_co_1st) = cpatch%area + + ! Patch level running means + call this%SetRMeanRestartVar(cpatch%tveg24, ir_tveg24patch_pa, io_idx_co_1st) ! set cohorts per patch for IO rio_ncohort_pa( io_idx_co_1st ) = cohortsperpatch @@ -2696,6 +2796,10 @@ subroutine get_restart_vectors(this, nc, nsites, sites) cpatch%solar_zenith_flag = ( rio_solar_zenith_flag_pa(io_idx_co_1st) .eq. itrue ) cpatch%solar_zenith_angle = rio_solar_zenith_angle_pa(io_idx_co_1st) + + call this%GetRMeanRestartVar(cpatch%tveg24, ir_tveg24patch_pa, io_idx_co_1st) + + ! set cohorts per patch for IO if ( debug ) then diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningMeanMod.F90 index 60fa1b0850..b535b6b678 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningMeanMod.F90 @@ -73,10 +73,11 @@ module FatesRunningMeanMod contains - procedure :: get_mean + procedure :: GetMean procedure :: InitRMean procedure :: UpdateRMean procedure :: FuseRMean + procedure :: CopyFromDonor end type rmean_type @@ -122,10 +123,10 @@ end subroutine define ! ===================================================================================== - function get_mean(this) + function GetMean(this) class(rmean_type) :: this - real(r8) :: get_mean + real(r8) :: GetMean if(this%def_type%method .eq. moving_ema_window) then if(this%c_index == 0) then @@ -141,9 +142,9 @@ function get_mean(this) call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if - get_mean = this%l_mean + GetMean = this%l_mean - end function get_mean + end function GetMean ! ===================================================================================== @@ -195,6 +196,33 @@ subroutine InitRMean(this,rmean_def,init_value,init_offset) return end subroutine InitRMean + + ! ===================================================================================== + + + subroutine CopyFromDonor(this, donor) + + class(rmean_type) :: this + class(rmean_type),intent(in) :: donor + + if( .not.associated(this%def_type)) then + write(fates_log(), *) 'Attempted to copy over running mean' + write(fates_log(), *) 'info from a donor into a new structure' + write(fates_log(), *) 'but the new structure did not have its' + write(fates_log(), *) 'def_type pointer associated' + 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 + + + return + end subroutine CopyFromDonor + + + ! ===================================================================================== subroutine UpdateRMean(this, new_value) From 8c9e4a47be6d81f393bb87dedaf75bacd7b0df24 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 2 Apr 2021 12:14:36 -0400 Subject: [PATCH 05/21] Added site-level running mean 24-veg temp diagnostic --- main/FatesHistoryInterfaceMod.F90 | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index fa90a93080..79bff307c5 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -528,6 +528,7 @@ module FatesHistoryInterfaceMod integer :: ih_fire_intensity_si_age integer :: ih_fire_sum_fuel_si_age integer :: ih_tveg24_si_age + integer :: ih_tveg24_si ! indices to (site x height) variables integer :: ih_canopy_height_dist_si_height @@ -2009,6 +2010,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_fire_intensity_si_age => this%hvars(ih_fire_intensity_si_age)%r82d, & hio_fire_sum_fuel_si_age => this%hvars(ih_fire_sum_fuel_si_age)%r82d, & hio_tveg24_si_age => this%hvars(ih_tveg24_si_age)%r82d, & + hio_tveg24_si => this%hvars(ih_tveg24_si)%r81d, & hio_burnt_frac_litter_si_fuel => this%hvars(ih_burnt_frac_litter_si_fuel)%r82d, & hio_fuel_amount_si_fuel => this%hvars(ih_fuel_amount_si_fuel)%r82d, & hio_fuel_amount_age_fuel => this%hvars(ih_fuel_amount_age_fuel)%r82d, & @@ -2152,6 +2154,9 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_harvest_carbonflux_si(io_si) = sites(s)%harvest_carbon_flux + hio_tveg24_si_age(io_si, :) = 0._r8 + hio_tveg24_si(io_si) = 0._r8 + ipa = 0 cpatch => sites(s)%oldest_patch do while(associated(cpatch)) @@ -2220,6 +2225,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_tveg24_si_age(io_si, cpatch%age_class) = & hio_tveg24_si_age(io_si, cpatch%age_class) + & cpatch%tveg24%GetMean()*cpatch%area + hio_tveg24_si(io_si) = hio_tveg24_si(io_si) + & + cpatch%tveg24%GetMean()*cpatch%area*area_inv end if if(associated(cpatch%tallest))then @@ -4564,11 +4571,16 @@ subroutine define_history_vars(this, initialize_variables) ! Running means call this%set_history_var(vname='TVEG24_AGE', units='Kelvin', & - long='24-hr running mean vegetation temperature by patch age', & + long='fates 24-hr running mean vegetation temperature by patch age', & use_default='active', & avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', flushval=hlm_hio_ignore_val, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_tveg24_si_age ) + call this%set_history_var(vname='TVEG24_SI', units='Kelvin', & + long='fates 24-hr running mean vegetation temperature by site', & + use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=hlm_hio_ignore_val, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_tveg24_si ) ! Litter Variables From b0f482f814008bc10b222acab240d911d78cc6f9 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 12 May 2021 10:14:30 -0400 Subject: [PATCH 06/21] Update main/FatesRunningMeanMod.F90 Co-authored-by: Charlie Koven --- main/FatesRunningMeanMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningMeanMod.F90 index b535b6b678..f1a30d8fa8 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningMeanMod.F90 @@ -131,7 +131,7 @@ function GetMean(this) if(this%def_type%method .eq. moving_ema_window) then if(this%c_index == 0) then write(fates_log(), *) 'attempting to get a running mean from a variable' - write(fates_log(), *) 'that has been given a value yet' + write(fates_log(), *) 'that has not been given a value yet' call endrun(msg=errMsg(sourcefile, __LINE__)) end if else From 7d9142353cc0da691d05e85f12f1e2953b6bb5e6 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 12 May 2021 10:15:03 -0400 Subject: [PATCH 07/21] Update main/FatesRunningMeanMod.F90 Co-authored-by: Charlie Koven --- main/FatesRunningMeanMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningMeanMod.F90 index f1a30d8fa8..43e5793b69 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningMeanMod.F90 @@ -51,7 +51,7 @@ module FatesRunningMeanMod type, public :: rmean_type real(r8) :: c_mean ! The current mean value, if this - ! is a moving window, its is the mean. + ! is a moving window, it is the mean. ! If this is a fixed window, it is only a partial mean ! as the value uses equal update weights and is not ! necessarily fully constructed. From d3b02c65ba4a69f2780805e3d2b1f6890d8a6b11 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 12 May 2021 10:15:15 -0400 Subject: [PATCH 08/21] Update main/FatesRunningMeanMod.F90 Co-authored-by: Charlie Koven --- main/FatesRunningMeanMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningMeanMod.F90 index 43e5793b69..98e2a07880 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningMeanMod.F90 @@ -59,7 +59,7 @@ module FatesRunningMeanMod real(r8) :: l_mean ! The latest reportable mean value ! this value is actually the same ! as c_mean for moving windows, and for fixed windows - ! it is the mean value when the time integratino window + ! it is the mean value when the time integration window ! last completed. integer :: c_index ! The number of values that have From 656a8ad926c62c8948361af13d23e3f7e8e4fa50 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 12 May 2021 10:47:13 -0400 Subject: [PATCH 09/21] Moving tveg24 to a fixed window average, part 1 --- biogeochem/EDPatchDynamicsMod.F90 | 11 +++++++++-- main/EDTypesMod.F90 | 2 +- main/FatesRunningMeanMod.F90 | 8 +++----- 3 files changed, 13 insertions(+), 8 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 25c5decc79..9d342cc6a6 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -85,7 +85,7 @@ module EDPatchDynamicsMod use SFParamsMod, only : SF_VAL_CWD_FRAC use EDParamsMod, only : logging_event_code use EDParamsMod, only : logging_export_frac - use FatesRunningMeanMod, only : ema_24hr + use FatesRunningMeanMod, only : ema_24hr, fixed_24hr ! CIME globals use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) @@ -1959,6 +1959,8 @@ end subroutine mortality_litter_fluxes subroutine create_patch(currentSite, new_patch, age, areap, label) + use FatesInterfaceTypesMod, only : hlm_current_tod,hlm_current_date,hlm_reference_date + ! ! !DESCRIPTION: ! Set default values for creating a new patch @@ -1988,8 +1990,13 @@ subroutine create_patch(currentSite, new_patch, age, areap, label) allocate(new_patch%fragmentation_scaler(currentSite%nlevsoil)) allocate(new_patch%tveg24) - call new_patch%tveg24%InitRMean(ema_24hr) ! No initial value + !call new_patch%tveg24%InitRMean(ema_24hr) ! No initial value + print*,hlm_current_tod,hlm_current_date,hlm_reference_date + stop + + call new_patch%tveg24%InitRMean(fixed_24hr ) + ! Litter ! Allocate, Zero Fluxes, and Initialize to "unset" values diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 94c4d77fcf..e4b274e5bc 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -415,7 +415,7 @@ module EDTypesMod ! 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 :: tveg24 ! 24-hour mean vegetation temperature (K) ! LEAF ORGANIZATION diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningMeanMod.F90 index 98e2a07880..0dc9fe3355 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningMeanMod.F90 @@ -11,18 +11,16 @@ module FatesRunningMeanMod implicit none private - integer, parameter :: maxlen_varname = 8 - ! These are flags that specify how the averaging window works. - ! Moving windows (default) can have an arbitrary size and update frequency) + ! Exponential moving average (EMA) windows have an arbitrary size and update frequency) ! and it is technically never reset, it just averages indefinitely. - ! But hourly, six-hourly, daily, monthly and yearly windows have pre-set + ! But hourly, six-hourly, daily, monthly and yearly fixed windows have pre-set ! window sizes associated with their namesake, and more importantly, they ! are zero'd at the beginning of the interval, and get equal average weighting ! over their construction period. - integer, public, parameter :: moving_ema_window = 0 + integer, public, parameter :: moving_ema_window = 0 ! (exponential moving average) integer, public, parameter :: fixed_window = 1 ! This type defines a type of mean. It does not From 7136c10f60633d61c6b253e02bb2739084b27d8e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 14 May 2021 15:27:19 -0400 Subject: [PATCH 10/21] Cleaned up some diagnostic print statements, finalizing vegtemp 24 rmean as using fixed window --- biogeochem/EDPatchDynamicsMod.F90 | 11 ++--- main/FatesRunningMeanMod.F90 | 78 ++++++++++++++++++------------- 2 files changed, 49 insertions(+), 40 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 9d342cc6a6..d3f5165674 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -580,7 +580,6 @@ subroutine spawn_patches( currentSite, bc_in) endif - ! next create patch to receive secondary forest area if ( site_areadis_secondary .gt. nearzero) then allocate(new_patch_secondary) @@ -1974,6 +1973,9 @@ subroutine create_patch(currentSite, new_patch, age, areap, label) real(r8), intent(in) :: areap ! initial area of this patch in m2. integer, intent(in) :: label ! anthropogenic disturbance label + real(r8), parameter :: temp_init_vegc = 15._r8 ! Until bc's are pointed to by sites + ! give veg temp a default temp. + ! !LOCAL VARIABLES: !--------------------------------------------------------------------- integer :: el ! element loop index @@ -1990,12 +1992,7 @@ subroutine create_patch(currentSite, new_patch, age, areap, label) allocate(new_patch%fragmentation_scaler(currentSite%nlevsoil)) allocate(new_patch%tveg24) - !call new_patch%tveg24%InitRMean(ema_24hr) ! No initial value - - print*,hlm_current_tod,hlm_current_date,hlm_reference_date - stop - - call new_patch%tveg24%InitRMean(fixed_24hr ) + call new_patch%tveg24%InitRMean(fixed_24hr,init_value=temp_init_vegc,init_offset=real(hlm_current_tod,r8) ) ! Litter ! Allocate, Zero Fluxes, and Initialize to "unset" values diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningMeanMod.F90 index 0dc9fe3355..5075417bea 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningMeanMod.F90 @@ -79,6 +79,8 @@ module FatesRunningMeanMod end type rmean_type + + logical, parameter :: debug = .true. character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -86,8 +88,8 @@ module FatesRunningMeanMod ! Define the time methods that we want to have available to us - class(rmean_def_type), public, pointer :: ema_24hr - class(rmean_def_type), public, pointer :: fixed_24hr + 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 contains @@ -101,16 +103,18 @@ subroutine define(this,mem_period,up_period,method) integer,intent(in) :: method ! Check the memory and update periods - if( abs(nint(mem_period/up_period)-mem_period/up_period) > nearzero ) then - write(fates_log(), *) 'While defining a running mean definition' - write(fates_log(), *) 'an update and memory period was specified' - write(fates_log(), *) 'where the update period is not an exact fraction of the period' - write(fates_log(), *) 'mem_period: ',mem_period - write(fates_log(), *) 'up_period: ',up_period - write(fates_log(), *) 'exiting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + if(debug) then + if( abs(nint(mem_period/up_period)-mem_period/up_period) > nearzero ) then + write(fates_log(), *) 'While defining a running mean definition' + write(fates_log(), *) 'an update and memory period was specified' + write(fates_log(), *) 'where the update period is not an exact fraction of the period' + write(fates_log(), *) 'mem_period: ',mem_period + write(fates_log(), *) 'up_period: ',up_period + write(fates_log(), *) 'exiting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if end if - + this%mem_period = mem_period this%up_period = up_period this%method = method @@ -127,18 +131,11 @@ function GetMean(this) real(r8) :: GetMean if(this%def_type%method .eq. moving_ema_window) then - if(this%c_index == 0) then + 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 - else - if(this%c_index .ne. this%def_type%n_mem)then - write(fates_log(), *) 'attempting to get a mean over a fixed window' - write(fates_log(), *) 'at a time where the window has not completed' - write(fates_log(), *) 'its cycle yet' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if end if GetMean = this%l_mean @@ -165,18 +162,35 @@ subroutine InitRMean(this,rmean_def,init_value,init_offset) this%def_type => rmean_def if(this%def_type%method .eq. fixed_window) then - - if(.not.(present(init_offset).and.present(init_value)) )then - write(fates_log(), *) 'when initializing a temporal mean on a fixed window' - write(fates_log(), *) 'there must be an initial value and a time offset' - write(fates_log(), *) 'specified.' - call endrun(msg=errMsg(sourcefile, __LINE__)) + + if(debug) then + if(.not.(present(init_offset).and.present(init_value)) )then + write(fates_log(), *) 'when initializing a temporal mean on a fixed window' + write(fates_log(), *) 'there must be an initial value and a time offset' + write(fates_log(), *) 'specified.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + 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 + 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(), *) 'offset: ',init_offset + write(fates_log(), *) 'up freq: ',rmean_def%up_period + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if(init_offset<-nearzero) then + write(fates_log(), *) 'offset must be positive: ',init_offset + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if - - this%c_index = modulo(nint(init_offset/rmean_def%up_period)+1,rmean_def%n_mem) + + 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 this%l_mean = init_value - + elseif(this%def_type%method .eq. moving_ema_window) then if(present(init_value))then @@ -250,18 +264,16 @@ subroutine UpdateRMean(this, new_value) ! average, then zero things out if(this%c_index == this%def_type%n_mem) then - this%c_mean = 0._r8 + this%l_mean = this%c_mean + this%c_mean = 0._r8 this%c_index = 0 + end if 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 - if(this%c_index == this%def_type%n_mem) then - this%l_mean = this%c_mean - end if - end if return From b887b062b9f0489eab2d4778d59ce038386bf349 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 24 Jun 2021 10:31:36 -0400 Subject: [PATCH 11/21] Adding both patch and cohort level photosynthesis acclimation vegetation temperature running mean. --- biogeochem/EDCohortDynamicsMod.F90 | 17 +++++++++++++++-- biogeochem/EDPatchDynamicsMod.F90 | 8 ++++++-- main/EDInitMod.F90 | 14 ++++++++------ main/EDTypesMod.F90 | 12 ++++++++++++ main/FatesInterfaceMod.F90 | 14 +++++++++++++- main/FatesRunningMeanMod.F90 | 1 + 6 files changed, 55 insertions(+), 11 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index b6714ee3e9..e50d168533 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -300,8 +300,10 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & call InitPRTBoundaryConditions(new_cohort) - - + ! Allocate running mean functions + allocate(new_cohort%tveg_lpa) + call new_cohort%tveg_lpa%InitRMean(ema_lpa) + ! Recuits do not have mortality rates, nor have they moved any ! carbon when they are created. They will bias our statistics ! until they have experienced a full day. We need a newly recruited flag. @@ -991,6 +993,9 @@ subroutine DeallocateCohort(currentCohort) ! ---------------------------------------------------------------------------------- type(ed_cohort_type),intent(inout) :: currentCohort + + ! Remove the running mean structure + deallocate(new_cohort%tveg_lpa) ! At this point, nothing should be pointing to current Cohort if (hlm_use_planthydro.eq.itrue) call DeallocateHydrCohort(currentCohort) @@ -1155,6 +1160,10 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) end do end if + ! Running mean fuses based on number density fraction just + ! like other variables + call currentCohort%tveg_lpa%FuseRMean(nextc%tveg_lpa,currentCohort%n/newn) + ! new cohort age is weighted mean of two cohorts currentCohort%coage = & (currentCohort%coage * (currentCohort%n/(currentCohort%n + nextc%n))) + & @@ -1786,6 +1795,7 @@ subroutine copy_cohort( currentCohort,copyc ) n%size_by_pft_class = o%size_by_pft_class n%coage_class = o%coage_class n%coage_by_pft_class = o%coage_by_pft_class + ! This transfers the PRT objects over. call n%prt%CopyPRTVartypes(o%prt) @@ -1795,6 +1805,9 @@ subroutine copy_cohort( currentCohort,copyc ) n%tpu25top = o%tpu25top n%kp25top = o%kp25top + ! Copy over running means + n%tveg_lpa%CopyFromDonor(o%tveg_lpa) + ! CARBON FLUXES n%gpp_acc_hold = o%gpp_acc_hold n%gpp_acc = o%gpp_acc diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index c9f5b2d16f..12da82bf23 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -85,7 +85,7 @@ module EDPatchDynamicsMod use SFParamsMod, only : SF_VAL_CWD_FRAC use EDParamsMod, only : logging_event_code use EDParamsMod, only : logging_export_frac - use FatesRunningMeanMod, only : ema_24hr, fixed_24hr + use FatesRunningMeanMod, only : ema_24hr, fixed_24hr, ema_lpa ! CIME globals use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) @@ -676,7 +676,7 @@ subroutine spawn_patches( currentSite, bc_in) ! These values will inherit all info from the original patch ! -------------------------------------------------------------------------- call new_patch%tveg24%CopyFromDonor(currentPatch%tveg24) - + call new_patch%tveg_lpa%CopyFromDonor(currentPatch%tveg_lpa) ! -------------------------------------------------------------------------- @@ -2013,6 +2013,8 @@ subroutine create_patch(currentSite, new_patch, age, areap, label) allocate(new_patch%tveg24) call new_patch%tveg24%InitRMean(fixed_24hr,init_value=temp_init_vegc,init_offset=real(hlm_current_tod,r8) ) + allocate(new_patch%tveg_lpa) + call new_patch%tveg_lpa%InitRmean(ema_lpa,init_value=temp_init_vegc) ! Litter ! Allocate, Zero Fluxes, and Initialize to "unset" values @@ -2517,6 +2519,7 @@ subroutine fuse_2_patches(csite, dp, rp) ! 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) rp%fuel_eff_moist = (dp%fuel_eff_moist*dp%area + rp%fuel_eff_moist*rp%area) * inv_sum_area rp%livegrass = (dp%livegrass*dp%area + rp%livegrass*rp%area) * inv_sum_area @@ -2858,6 +2861,7 @@ subroutine dealloc_patch(cpatch) ! Deallocate any running means deallocate(cpatch%tveg24) + deallocate(cpatch%tveg_lpa) return end subroutine dealloc_patch diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 6b3d861a28..8862e7237a 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -672,12 +672,14 @@ subroutine init_cohorts( site_in, patch_in, bc_in) endif !use_this_pft enddo !numpft - ! Zero the mass flux pools of the new cohorts -! temp_cohort => patch_in%tallest -! do while(associated(temp_cohort)) -! call temp_cohort%prt%ZeroRates() -! temp_cohort => temp_cohort%shorter -! end do + ! Pass patch level temperature to the new cohorts (this is a nominal 15C right now) + temp_cohort => patch_in%tallest + do while(associated(temp_cohort)) + call temp_cohort%tveg_lpa%UpdateRmean(patch_in%tveg_lpa%GetMean()) + temp_cohort => temp_cohort%shorter + end do + + call fuse_cohorts(site_in, patch_in,bc_in) call sort_cohorts(patch_in) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 9904e360e3..55b42e41d9 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -192,6 +192,12 @@ module EDTypesMod integer, public :: n_uptake_mode integer, public :: p_uptake_mode + + + ! Leaf photosynthetic temperature acclimation timescale [days] + ! (This variable may be moved to the FATES parameter file soon) + + real(r8), parameter, public :: leaf_photo_temp_acclim_days = 30._r8 !************************************ @@ -388,6 +394,12 @@ module EDTypesMod ! Hydraulics type(ed_cohort_hydr_type), pointer :: co_hydr ! All cohort hydraulics data, see FatesHydraulicsMemMod.F90 + + ! Running means + class(rmean_type), pointer :: tveg_lpa ! exponential moving average of leaf temperature at the + ! leaf photosynthetic acclimation time-scale + + end type ed_cohort_type !************************************ diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index ca0f95a03a..2b86c4f236 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -23,6 +23,7 @@ module FatesInterfaceMod use EDTypesMod , only : numlevsoil_max use EDTypesMod , only : ed_site_type use EDTypesMod , only : ed_patch_type + use EDTypesMod , only : leaf_photo_temp_acclim_days use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue,ifalse use FatesConstantsMod , only : nearzero @@ -72,6 +73,7 @@ module FatesInterfaceMod use PRTAllometricCNPMod , only : InitPRTGlobalAllometricCNP use FatesRunningMeanMod , only : ema_24hr use FatesRunningMeanMod , only : fixed_24hr + use FatesRunningMeanMod , only : ema_lpa use FatesRunningMeanMod , only : moving_ema_window use FatesRunningMeanMod , only : fixed_window @@ -865,7 +867,8 @@ subroutine SetFatesGlobalElements(use_fates) call ema_24hr%define(sec_per_day, hlm_stepsize, moving_ema_window) allocate(fixed_24hr) call fixed_24hr%define(sec_per_day, hlm_stepsize, fixed_window) - + allocate(ema_lpa) + call ema_lpa%define(leaf_photo_temp_acclim_days,hlm_stepsize,moving_ema_window) else @@ -1839,6 +1842,7 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in) type(bc_in_type), intent(in) :: bc_in(:) type(ed_patch_type), pointer :: cpatch + type(ed_cohort_type), pointer :: ccohort integer :: s, ifp @@ -1848,6 +1852,14 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in) do while(associated(cpatch)) ifp=ifp+1 call cpatch%tveg24%UpdateRMean(bc_in(s)%tveg_pa(ifp)) + call cpatch%tveg_lpa%UpdateRMean(bc_in(s)%tveg_pa(ifp)) + + ccohort => cpatch%tallest + do while (associated(ccohort)) + call ccohort%tveg_lpa%UpdateRMean(bc_in(s)%tveg_pa(ifp)) + ccohort => ccohort%shorter + end do + cpatch => cpatch%younger enddo enddo diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningMeanMod.F90 index 5075417bea..df1c12c7be 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningMeanMod.F90 @@ -90,6 +90,7 @@ module FatesRunningMeanMod 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 contains From 953c5eb666032b36d7f2604019988ed0f5f591d5 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 6 Jul 2021 11:20:19 -0400 Subject: [PATCH 12/21] Updating running mean functions to include parameter conrolled acclimation to temperature. Updated recruit initialization --- biogeochem/EDCohortDynamicsMod.F90 | 2 +- biogeophys/FatesPlantRespPhotosynthMod.F90 | 5 +++++ main/EDTypesMod.F90 | 10 ++-------- main/FatesInterfaceMod.F90 | 5 +++-- 4 files changed, 11 insertions(+), 11 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index e50d168533..8440e0b079 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -302,7 +302,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & ! Allocate running mean functions allocate(new_cohort%tveg_lpa) - call new_cohort%tveg_lpa%InitRMean(ema_lpa) + call new_cohort%tveg_lpa%InitRMean(ema_lpa,init_value=patchptr%tveg_lpa%GetMean()) ! Recuits do not have mortality rates, nor have they moved any ! carbon when they are created. They will bias our statistics diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 349f6473bf..644ad66d0e 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -355,8 +355,13 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if(currentPatch%countcohorts > 0.0)then ! Ignore empty patches + print*,"Ts: ",bc_in(s)%tveg_pa(ifp),cpatch%tveg_lpa%GetMean() + currentCohort => currentPatch%tallest do while (associated(currentCohort)) ! Cohort loop + + print*," c ",currentCohort%tveg_lpa%GetMean() + ! Identify the canopy layer (cl), functional type (ft) ! and the leaf layer (IV) for this cohort diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 55b42e41d9..7f053c832b 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -193,13 +193,6 @@ module EDTypesMod integer, public :: n_uptake_mode integer, public :: p_uptake_mode - - ! Leaf photosynthetic temperature acclimation timescale [days] - ! (This variable may be moved to the FATES parameter file soon) - - real(r8), parameter, public :: leaf_photo_temp_acclim_days = 30._r8 - - !************************************ !** COHORT type structure ** !************************************ @@ -430,7 +423,8 @@ module EDTypesMod ! 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 + ! leaf photosynthesis acclimation timescale ! LEAF ORGANIZATION real(r8) :: pft_agb_profile(maxpft,n_dbh_bins) ! binned above ground biomass, for patch fusion: KgC/m2 diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 2b86c4f236..7bb20b7f41 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -23,7 +23,6 @@ module FatesInterfaceMod use EDTypesMod , only : numlevsoil_max use EDTypesMod , only : ed_site_type use EDTypesMod , only : ed_patch_type - use EDTypesMod , only : leaf_photo_temp_acclim_days use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue,ifalse use FatesConstantsMod , only : nearzero @@ -40,6 +39,7 @@ module FatesInterfaceMod use EDParamsMod , only : FatesReportParams use EDParamsMod , only : bgc_soil_salinity use FatesPlantHydraulicsMod , only : InitHydroGlobals + use EDParamsMod , only : photo_temp_acclim_timescale use EDParamsMod , only : ED_val_history_sizeclass_bin_edges use EDParamsMod , only : ED_val_history_ageclass_bin_edges use EDParamsMod , only : ED_val_history_height_bin_edges @@ -868,7 +868,8 @@ subroutine SetFatesGlobalElements(use_fates) allocate(fixed_24hr) call fixed_24hr%define(sec_per_day, hlm_stepsize, fixed_window) allocate(ema_lpa) - call ema_lpa%define(leaf_photo_temp_acclim_days,hlm_stepsize,moving_ema_window) + call ema_lpa%define(photo_temp_acclim_timescale*sec_per_day, & + hlm_stepsize,moving_ema_window) else From 8fa8241d1a95699d74c257608397d8ddd9b98efa Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 7 Jul 2021 12:01:46 -0400 Subject: [PATCH 13/21] Various fixes to initialization and restarting cohort/patch level running means. --- biogeochem/EDCanopyStructureMod.F90 | 12 ++++++++++ biogeochem/EDCohortDynamicsMod.F90 | 5 ++-- biogeochem/EDPatchDynamicsMod.F90 | 14 +++++++---- main/EDTypesMod.F90 | 4 ++-- main/FatesInterfaceMod.F90 | 1 + main/FatesInventoryInitMod.F90 | 7 +++++- main/FatesRestartInterfaceMod.F90 | 36 ++++++++++++++++++++++------- 7 files changed, 61 insertions(+), 18 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 6aee6de8de..3e10966bb4 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -41,6 +41,7 @@ module EDCanopyStructureMod use PRTGenericMod, only : repro_organ use PRTGenericMod, only : struct_organ use PRTGenericMod, only : SetState + use FatesRunningMeanMod, only : ema_lpa ! CIME Globals @@ -671,6 +672,11 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr,bc_in) call InitHydrCohort(currentSite,copyc) endif + ! Initialize running means + allocate(copyc%tveg_lpa) + call copyc%tveg_lpa%InitRMean(ema_lpa,& + init_value=currentPatch%tveg_lpa%GetMean()) + call copy_cohort(currentCohort, copyc) newarea = currentCohort%c_area - cc_loss @@ -1123,6 +1129,12 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) if( hlm_use_planthydro.eq.itrue ) then call InitHydrCohort(CurrentSite,copyc) endif + + ! Initialize running means + allocate(copyc%tveg_lpa) + call copyc%tveg_lpa%InitRMean(ema_lpa,& + init_value=currentPatch%tveg_lpa%GetMean()) + call copy_cohort(currentCohort, copyc) !makes an identical copy... newarea = currentCohort%c_area - cc_gain !new area of existing cohort diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 8440e0b079..5f5999fbf1 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -16,6 +16,7 @@ module EDCohortDynamicsMod use FatesConstantsMod , only : fates_unset_r8 use FatesConstantsMod , only : nearzero use FatesConstantsMod , only : calloc_abs_error + use FatesRunningMeanMod , only : ema_lpa use FatesInterfaceTypesMod , only : hlm_days_per_year use FatesInterfaceTypesMod , only : nleafage use SFParamsMod , only : SF_val_CWD_frac @@ -995,7 +996,7 @@ subroutine DeallocateCohort(currentCohort) type(ed_cohort_type),intent(inout) :: currentCohort ! Remove the running mean structure - deallocate(new_cohort%tveg_lpa) + deallocate(currentCohort%tveg_lpa) ! At this point, nothing should be pointing to current Cohort if (hlm_use_planthydro.eq.itrue) call DeallocateHydrCohort(currentCohort) @@ -1806,7 +1807,7 @@ subroutine copy_cohort( currentCohort,copyc ) n%kp25top = o%kp25top ! Copy over running means - n%tveg_lpa%CopyFromDonor(o%tveg_lpa) + call n%tveg_lpa%CopyFromDonor(o%tveg_lpa) ! CARBON FLUXES n%gpp_acc_hold = o%gpp_acc_hold diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 25dd9e915f..209302b2e8 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -49,6 +49,7 @@ module EDPatchDynamicsMod use FatesGlobals , only : endrun => fates_endrun use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue, ifalse + use FatesConstantsMod , only : t_water_freeze_k_1atm use FatesPlantHydraulicsMod, only : InitHydrCohort use FatesPlantHydraulicsMod, only : AccumulateMortalityWaterStorage use FatesPlantHydraulicsMod, only : DeallocateHydrCohort @@ -698,7 +699,10 @@ subroutine spawn_patches( currentSite, bc_in) nc%prt => null() call InitPRTObject(nc%prt) call InitPRTBoundaryConditions(nc) - + ! Allocate running mean functions + allocate(nc%tveg_lpa) + call nc%tveg_lpa%InitRMean(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) + call zero_cohort(nc) ! nc is the new cohort that goes in the disturbed patch (new_patch)... currentCohort @@ -1992,8 +1996,8 @@ subroutine create_patch(currentSite, new_patch, age, areap, label) real(r8), intent(in) :: areap ! initial area of this patch in m2. integer, intent(in) :: label ! anthropogenic disturbance label - real(r8), parameter :: temp_init_vegc = 15._r8 ! Until bc's are pointed to by sites - ! give veg temp a default temp. + ! Until bc's are pointed to by sites give veg temp a default temp [K] + real(r8), parameter :: temp_init_veg = 15._r8+t_water_freeze_k_1atm ! !LOCAL VARIABLES: !--------------------------------------------------------------------- @@ -2011,9 +2015,9 @@ subroutine create_patch(currentSite, new_patch, age, areap, label) allocate(new_patch%fragmentation_scaler(currentSite%nlevsoil)) allocate(new_patch%tveg24) - call new_patch%tveg24%InitRMean(fixed_24hr,init_value=temp_init_vegc,init_offset=real(hlm_current_tod,r8) ) + call new_patch%tveg24%InitRMean(fixed_24hr,init_value=temp_init_veg,init_offset=real(hlm_current_tod,r8) ) allocate(new_patch%tveg_lpa) - call new_patch%tveg_lpa%InitRmean(ema_lpa,init_value=temp_init_vegc) + call new_patch%tveg_lpa%InitRmean(ema_lpa,init_value=temp_init_veg) ! Litter ! Allocate, Zero Fluxes, and Initialize to "unset" values diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 7f053c832b..77c9392a4b 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -390,7 +390,7 @@ module EDTypesMod ! Running means class(rmean_type), pointer :: tveg_lpa ! exponential moving average of leaf temperature at the - ! leaf photosynthetic acclimation time-scale + ! leaf photosynthetic acclimation time-scale [K] end type ed_cohort_type @@ -424,7 +424,7 @@ module EDTypesMod !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 - ! leaf photosynthesis acclimation timescale + ! leaf photosynthesis acclimation timescale [K] ! LEAF ORGANIZATION real(r8) :: pft_agb_profile(maxpft,n_dbh_bins) ! binned above ground biomass, for patch fusion: KgC/m2 diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 7bb20b7f41..21e9ca38cf 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -23,6 +23,7 @@ module FatesInterfaceMod use EDTypesMod , only : numlevsoil_max use EDTypesMod , only : ed_site_type use EDTypesMod , only : ed_patch_type + use EDTypesMod , only : ed_cohort_type use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue,ifalse use FatesConstantsMod , only : nearzero diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index efdebb8708..887310c105 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -62,6 +62,7 @@ module FatesInventoryInitMod use PRTGenericMod, only : phosphorus_element use PRTGenericMod, only : SetState use FatesConstantsMod, only : primaryforest + use FatesRunningMeanMod , only : ema_lpa use PRTGenericMod, only : StorageNutrientTarget implicit none @@ -1042,7 +1043,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & temp_cohort%structmemory = 0._r8 cstatus = leaves_on - stem_drop_fraction = EDPftvarcon_inst%phen_stem_drop_fraction(temp_cohort%pft) + stem_drop_fraction = EDPftvarcon_inst%phen_stem_drop_fraction(temp_cohort%pft) if( prt_params%season_decid(temp_cohort%pft) == itrue .and. & any(csite%cstatus == [phen_cstat_nevercold,phen_cstat_iscold])) then @@ -1069,6 +1070,10 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & prt_obj => null() call InitPRTObject(prt_obj) + ! Allocate running mean functions + allocate(temp_cohort%tveg_lpa) + call temp_cohort%tveg_lpa%InitRMean(ema_lpa,init_value=cpatch%tveg_lpa%GetMean()) + do el = 1,num_elements element_id = element_list(el) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index a0bc260066..04a17e1a1c 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -38,6 +38,7 @@ module FatesRestartInterfaceMod use PRTGenericMod, only : prt_global use PRTGenericMod, only : num_elements use FatesRunningMeanMod, only : rmean_type + use FatesRunningMeanMod, only : ema_lpa ! CIME GLOBALS use shr_log_mod , only : errMsg => shr_log_errMsg @@ -138,8 +139,10 @@ module FatesRestartInterfaceMod integer :: ir_gnd_alb_dir_pasb ! Running Means - integer :: ir_tveg24patch_pa - + integer :: ir_tveg24_pa + integer :: ir_tveglpa_pa + integer :: ir_tveglpa_co + integer :: ir_ddbhdt_co integer :: ir_resp_tstep_co integer :: ir_pft_co @@ -1217,9 +1220,16 @@ subroutine define_restart_vars(this, initialize_variables) call this%DefineRMeanRestartVar(vname='fates_tveg24patch',vtype=cohort_r8, & long_name='24-hour patch veg temp', & - units='K', initialize=initialize_variables,ivar=ivar, index = ir_tveg24patch_pa) - + units='K', initialize=initialize_variables,ivar=ivar, index = ir_tveg24_pa) + call this%DefineRMeanRestartVar(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_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) + ! Register all of the PRT states and fluxes @@ -1999,7 +2009,9 @@ subroutine set_restart_vectors(this,nc,nsites,sites) write(fates_log(),*) 'CLTV offsetNumCohorts II ',io_idx_co, & cohortsperpatch endif - + + call this%SetRMeanRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) + io_idx_co = io_idx_co + 1 ccohort => ccohort%taller @@ -2016,7 +2028,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_area_pa(io_idx_co_1st) = cpatch%area ! Patch level running means - call this%SetRMeanRestartVar(cpatch%tveg24, ir_tveg24patch_pa, io_idx_co_1st) + 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) ! set cohorts per patch for IO rio_ncohort_pa( io_idx_co_1st ) = cohortsperpatch @@ -2356,6 +2369,11 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in, bc_out) call InitHydrCohort(sites(s),new_cohort) end if + ! Allocate running mean functions + allocate(new_cohort%tveg_lpa) + call new_cohort%tveg_lpa%InitRMean(ema_lpa) + + ! Update the previous prev_cohort => new_cohort @@ -2766,6 +2784,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) call UpdatePlantPsiFTCFromTheta(ccohort,sites(s)%si_hydr) end if + + call this%GetRMeanRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) io_idx_co = io_idx_co + 1 @@ -2794,8 +2814,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) cpatch%solar_zenith_angle = rio_solar_zenith_angle_pa(io_idx_co_1st) - call this%GetRMeanRestartVar(cpatch%tveg24, ir_tveg24patch_pa, 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) ! set cohorts per patch for IO From da2d7e36e52c205a9883f2c3a65bb358d994bdd2 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 8 Jul 2021 13:03:37 -0400 Subject: [PATCH 14/21] Updates to leaf acclimation temperature running mean: history variables, and removing unnecessary (duplicative) boundary condition. --- biogeophys/FatesPlantRespPhotosynthMod.F90 | 5 --- main/EDInitMod.F90 | 17 +++++--- main/EDMainMod.F90 | 7 +++- main/EDTypesMod.F90 | 6 ++- main/FatesHistoryInterfaceMod.F90 | 16 ++++++++ main/FatesInterfaceMod.F90 | 47 ++++++++++++++++------ main/FatesInterfaceTypesMod.F90 | 3 -- 7 files changed, 74 insertions(+), 27 deletions(-) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 644ad66d0e..dbdfd0b3d4 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -355,14 +355,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if(currentPatch%countcohorts > 0.0)then ! Ignore empty patches - print*,"Ts: ",bc_in(s)%tveg_pa(ifp),cpatch%tveg_lpa%GetMean() - currentCohort => currentPatch%tallest do while (associated(currentCohort)) ! Cohort loop - print*," c ",currentCohort%tveg_lpa%GetMean() - - ! Identify the canopy layer (cl), functional type (ft) ! and the leaf layer (IV) for this cohort ft = currentCohort%pft diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 8862e7237a..b166c967de 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -68,6 +68,7 @@ module EDInitMod use PRTGenericMod, only : nitrogen_element use PRTGenericMod, only : phosphorus_element use PRTGenericMod, only : SetState + use FatesSizeAgeTypeIndicesMod,only : get_age_class_index ! CIME GLOBALS use shr_log_mod , only : errMsg => shr_log_errMsg @@ -130,6 +131,8 @@ subroutine init_site_vars( site_in, bc_in, bc_out ) allocate(site_in%area_pft(1:numpft)) allocate(site_in%use_this_pft(1:numpft)) + allocate(site_in%area_by_age(1:nlevage)) + do el=1,num_elements allocate(site_in%flux_diags(el)%leaf_litter_input(1:numpft)) allocate(site_in%flux_diags(el)%root_litter_input(1:numpft)) @@ -228,6 +231,8 @@ subroutine zero_site( site_in ) site_in%area_pft(:) = 0._r8 site_in%use_this_pft(:) = fates_unset_int + + site_in%area_by_age(:) = 0._r8 end subroutine zero_site ! ============================================================================ @@ -347,6 +352,7 @@ subroutine init_patches( nsites, sites, bc_in) ! !LOCAL VARIABLES: integer :: s integer :: el + integer :: ageclass real(r8) :: age !notional age of this patch ! dummy locals @@ -358,9 +364,7 @@ subroutine init_patches( nsites, sites, bc_in) type(ed_patch_type), pointer :: newp type(ed_patch_type), pointer :: currentPatch - ! List out some nominal patch values that are used for Near Bear Ground initializations - ! as well as initializing inventory - age = 0.0_r8 + ! --------------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------------- @@ -409,7 +413,9 @@ subroutine init_patches( nsites, sites, bc_in) sites(s)%oldest_patch => newp ! make new patch... - + ! List out some nominal patch values that are used for Near Bear Ground initializations + ! as well as initializing inventory + age = 0.0_r8 call create_patch(sites(s), newp, age, area, primaryforest) ! Initialize the litter pools to zero, these @@ -437,11 +443,11 @@ subroutine init_patches( nsites, sites, bc_in) end if - ! zero all the patch fire variables for the first timestep do s = 1, nsites currentPatch => sites(s)%youngest_patch do while(associated(currentPatch)) + ! zero all the patch fire variables for the first timestep currentPatch%litter_moisture(:) = 0._r8 currentPatch%fuel_eff_moist = 0._r8 currentPatch%livegrass = 0._r8 @@ -461,6 +467,7 @@ subroutine init_patches( nsites, sites, bc_in) currentPatch%scorch_ht(:) = 0._r8 currentPatch%frac_burnt = 0._r8 currentPatch%burnt_frac_litter(:) = 0._r8 + currentPatch => currentPatch%older enddo enddo diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 5deb2c5084..0518c1b34e 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -655,6 +655,8 @@ subroutine ed_update_site( currentSite, bc_in, bc_out ) call TotalBalanceCheck(currentSite,final_check_id) + currentSite%area_by_age(:) = 0._r8 + currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) @@ -666,7 +668,10 @@ subroutine ed_update_site( currentSite, bc_in, bc_out ) ! This cohort count is used in the photosynthesis loop call count_cohorts(currentPatch) - + ! Update the total area of by patch age class array + currentSite%area_by_age(currentPatch%age_class) = & + currentSite%area_by_age(currentPatch%age_class) + currentPatch%area + currentPatch => currentPatch%younger enddo diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 77c9392a4b..05a4178da9 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -707,7 +707,11 @@ module EDTypesMod ! Fixed Biogeography mode inputs real(r8), allocatable :: area_PFT(:) ! Area allocated to individual PFTs integer, allocatable :: use_this_pft(:) ! Is area_PFT > 0 ? (1=yes, 0=no) - + + ! Total area of patches in each age bin [m2] + real(r8), allocatable :: area_by_age(:) + + ! Mass Balance (allocation for each element) type(site_massbal_type), pointer :: mass_balance(:) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 480aefb180..cd9269849d 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -534,6 +534,8 @@ module FatesHistoryInterfaceMod integer :: ih_fire_sum_fuel_si_age integer :: ih_tveg24_si_age integer :: ih_tveg24_si + integer,public :: ih_tveglpa_si_age + integer,public :: ih_tveglpa_si ! indices to (site x height) variables integer :: ih_canopy_height_dist_si_height @@ -4624,6 +4626,20 @@ subroutine define_history_vars(this, initialize_variables) use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=hlm_hio_ignore_val, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_tveg24_si ) + + call this%set_history_var(vname='TVEGLPA_AGE', units='Kelvin', & + long='fates leaf photo-acclim running mean vegetation temperature by patch age', & + use_default='active', & + avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', flushval=hlm_hio_ignore_val, upfreq=2, & + ivar=ivar, initialize=initialize_variables, index = ih_tveglpa_si_age ) + + call this%set_history_var(vname='TVEGLPA_SI', units='Kelvin', & + long='fates leaf photo-acclim running mean vegetation temperature by site', & + use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=hlm_hio_ignore_val, upfreq=2, & + ivar=ivar, initialize=initialize_variables, index = ih_tveglpa_si ) + + ! Litter Variables diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 21e9ca38cf..d5273a9ca7 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -24,6 +24,7 @@ module FatesInterfaceMod use EDTypesMod , only : ed_site_type use EDTypesMod , only : ed_patch_type use EDTypesMod , only : ed_cohort_type + use EDTypesMod , only : area_inv use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue,ifalse use FatesConstantsMod , only : nearzero @@ -77,6 +78,8 @@ module FatesInterfaceMod use FatesRunningMeanMod , only : ema_lpa use FatesRunningMeanMod , only : moving_ema_window use FatesRunningMeanMod , only : fixed_window + use FatesHistoryInterfaceMod , only : fates_hist + use FatesHistoryInterfaceMod , only : ih_tveglpa_si_age,ih_tveglpa_si ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -249,7 +252,6 @@ subroutine zero_bcs(fates,s) fates%bc_in(s)%precip24_pa(:) = 0.0_r8 fates%bc_in(s)%relhumid24_pa(:) = 0.0_r8 fates%bc_in(s)%wind24_pa(:) = 0.0_r8 - fates%bc_in(s)%tveg_pa(:) = 0.0_r8 fates%bc_in(s)%solad_parb(:,:) = 0.0_r8 fates%bc_in(s)%solai_parb(:,:) = 0.0_r8 @@ -461,7 +463,6 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats) allocate(bc_in%wind24_pa(maxPatchesPerSite)) allocate(bc_in%relhumid24_pa(maxPatchesPerSite)) allocate(bc_in%precip24_pa(maxPatchesPerSite)) - allocate(bc_in%tveg_pa(maxPatchesPerSite)) ! Radiation allocate(bc_in%solad_parb(maxPatchesPerSite,hlm_numSWb)) @@ -871,7 +872,6 @@ subroutine SetFatesGlobalElements(use_fates) allocate(ema_lpa) call ema_lpa%define(photo_temp_acclim_timescale*sec_per_day, & hlm_stepsize,moving_ema_window) - else ! If we are not using FATES, the cohort dimension is still @@ -1845,28 +1845,51 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in) type(ed_patch_type), pointer :: cpatch type(ed_cohort_type), pointer :: ccohort - integer :: s, ifp - + integer :: s, ifp, io_si do s = 1,size(sites,dim=1) + ifp=0 cpatch => sites(s)%oldest_patch do while(associated(cpatch)) ifp=ifp+1 - call cpatch%tveg24%UpdateRMean(bc_in(s)%tveg_pa(ifp)) - call cpatch%tveg_lpa%UpdateRMean(bc_in(s)%tveg_pa(ifp)) + call cpatch%tveg24%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) + call cpatch%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) ccohort => cpatch%tallest do while (associated(ccohort)) - call ccohort%tveg_lpa%UpdateRMean(bc_in(s)%tveg_pa(ifp)) + call ccohort%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) ccohort => ccohort%shorter end do cpatch => cpatch%younger enddo - enddo - + end do + + ! Update running mean history variables + ! ------------------------------------------------------------------------------- + associate(hio_tveglpa_si_age => fates_hist%hvars(ih_tveglpa_si_age)%r82d, & + hio_tveglpa_si => fates_hist%hvars(ih_tveglpa_si)%r81d) + + do s = 1,size(sites,dim=1) + + io_si = sites(s)%h_gid + hio_tveglpa_si_age(io_si,:) = 0._r8 + hio_tveglpa_si(io_si) = 0._r8 + + cpatch => sites(s)%oldest_patch + do while(associated(cpatch)) + hio_tveglpa_si_age(io_si,cpatch%age_class) = & + hio_tveglpa_si_age(io_si,cpatch%age_class) + & + cpatch%tveg_lpa%GetMean()*cpatch%area/sites(s)%area_by_age(cpatch%age_class) + hio_tveglpa_si(io_si) = hio_tveglpa_si(io_si) + & + cpatch%tveg_lpa%GetMean()*cpatch%area*area_inv + cpatch => cpatch%younger + enddo + end do + end associate + return end subroutine UpdateFatesRMeansTStep - -end module FatesInterfaceMod + + end module FatesInterfaceMod diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index 88a2fbde5f..df04f8e4d3 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -351,9 +351,6 @@ module FatesInterfaceTypesMod ! Average precipitation over the last 24 hours [mm/s] real(r8), allocatable :: precip24_pa(:) - ! Patch Vegetation temperature (K) - real(r8),allocatable :: tveg_pa(:) - ! Average relative humidity over past 24 hours [-] real(r8), allocatable :: relhumid24_pa(:) From 6023f7218823e223a6b4b4f8f7cd698f8808bfe2 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 22 Oct 2021 13:29:41 -0400 Subject: [PATCH 15/21] Commented out cohort-level vegetation temp averages, as they are redundant with patch-level given how the HLMs handle temperature physics --- biogeochem/EDCanopyStructureMod.F90 | 14 ++++--- biogeochem/EDCohortDynamicsMod.F90 | 15 +++++--- biogeochem/EDPatchDynamicsMod.F90 | 6 ++- main/EDInitMod.F90 | 11 +++--- main/EDPftvarcon.F90 | 27 +++++++------ main/EDTypesMod.F90 | 6 ++- main/FatesInterfaceMod.F90 | 59 +++++++++++++++-------------- main/FatesInventoryInitMod.F90 | 5 ++- main/FatesRestartInterfaceMod.F90 | 28 +++++++++----- 9 files changed, 99 insertions(+), 72 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index c48d702000..5f8261ae26 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -669,10 +669,11 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr,bc_in) call InitHydrCohort(currentSite,copyc) endif + ! (keep as an example) ! Initialize running means - allocate(copyc%tveg_lpa) - call copyc%tveg_lpa%InitRMean(ema_lpa, & - init_value=currentPatch%tveg_lpa%GetMean()) + !allocate(copyc%tveg_lpa) + !call copyc%tveg_lpa%InitRMean(ema_lpa, & + ! init_value=currentPatch%tveg_lpa%GetMean()) call copy_cohort(currentCohort, copyc) @@ -1127,10 +1128,11 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) call InitHydrCohort(CurrentSite,copyc) endif + ! (keep as an example) ! Initialize running means - allocate(copyc%tveg_lpa) - call copyc%tveg_lpa%InitRMean(ema_lpa,& - init_value=currentPatch%tveg_lpa%GetMean()) + !allocate(copyc%tveg_lpa) + !call copyc%tveg_lpa%InitRMean(ema_lpa,& + ! init_value=currentPatch%tveg_lpa%GetMean()) call copy_cohort(currentCohort, copyc) !makes an identical copy... diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 4d9258b276..893c0e66db 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -310,8 +310,10 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & ! Allocate running mean functions - allocate(new_cohort%tveg_lpa) - call new_cohort%tveg_lpa%InitRMean(ema_lpa,init_value=patchptr%tveg_lpa%GetMean()) + + ! (Keeping as an example) + !! allocate(new_cohort%tveg_lpa) + !! call new_cohort%tveg_lpa%InitRMean(ema_lpa,init_value=patchptr%tveg_lpa%GetMean()) ! Recuits do not have mortality rates, nor have they moved any @@ -1004,8 +1006,9 @@ subroutine DeallocateCohort(currentCohort) type(ed_cohort_type),intent(inout) :: currentCohort + ! (Keeping as an example) ! Remove the running mean structure - deallocate(currentCohort%tveg_lpa) + ! deallocate(currentCohort%tveg_lpa) ! At this point, nothing should be pointing to current Cohort if (hlm_use_planthydro.eq.itrue) call DeallocateHydrCohort(currentCohort) @@ -1170,9 +1173,10 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) end do end if + ! (Keeping as an example) ! Running mean fuses based on number density fraction just ! like other variables - call currentCohort%tveg_lpa%FuseRMean(nextc%tveg_lpa,currentCohort%n/newn) + !!call currentCohort%tveg_lpa%FuseRMean(nextc%tveg_lpa,currentCohort%n/newn) ! new cohort age is weighted mean of two cohorts currentCohort%coage = & @@ -1816,8 +1820,9 @@ subroutine copy_cohort( currentCohort,copyc ) n%tpu25top = o%tpu25top n%kp25top = o%kp25top + ! (Keeping as an example) ! Copy over running means - call n%tveg_lpa%CopyFromDonor(o%tveg_lpa) + ! call n%tveg_lpa%CopyFromDonor(o%tveg_lpa) ! CARBON FLUXES n%gpp_acc_hold = o%gpp_acc_hold diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index fe53d69f0b..54dd928bd8 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -703,9 +703,11 @@ subroutine spawn_patches( currentSite, bc_in) nc%prt => null() call InitPRTObject(nc%prt) call InitPRTBoundaryConditions(nc) + + ! (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()) + !allocate(nc%tveg_lpa) + !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) call zero_cohort(nc) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 87fc570a72..6e35c38fe2 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -896,12 +896,13 @@ subroutine init_cohorts( site_in, patch_in, bc_in) endif !use_this_pft enddo !numpft + ! (Keeping as an example) ! Pass patch level temperature to the new cohorts (this is a nominal 15C right now) - temp_cohort => patch_in%tallest - do while(associated(temp_cohort)) - call temp_cohort%tveg_lpa%UpdateRmean(patch_in%tveg_lpa%GetMean()) - temp_cohort => temp_cohort%shorter - end do + !temp_cohort => patch_in%tallest + !do while(associated(temp_cohort)) + !call temp_cohort%tveg_lpa%UpdateRmean(patch_in%tveg_lpa%GetMean()) + !temp_cohort => temp_cohort%shorter + !end do call fuse_cohorts(site_in, patch_in,bc_in) call sort_cohorts(patch_in) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 89c0cf2ea5..46366d7275 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -1475,7 +1475,7 @@ subroutine FatesCheckParams(is_master) use FatesConstantsMod , only : fates_check_param_set use FatesConstantsMod , only : itrue, ifalse use EDParamsMod , only : logging_mechanical_frac, logging_collateral_frac, logging_direct_frac - use FatesInterfaceTypesMod , only : hlm_use_fixed_biogeog + use FatesInterfaceTypesMod , only : hlm_use_fixed_biogeog,hlm_use_sp ! Argument logical, intent(in) :: is_master ! Only log if this is the master proc @@ -1740,17 +1740,20 @@ subroutine FatesCheckParams(is_master) ! check that the host-fates PFT map adds to one along HLM dimension so that all the HLM area ! goes to a FATES PFT. Each FATES PFT can get < or > 1 of an HLM PFT. - do hlm_pft = 1,size( EDPftvarcon_inst%hlm_pft_map,2) - sumarea = sum(EDPftvarcon_inst%hlm_pft_map(1:npft,hlm_pft)) - if(abs(sumarea-1.0_r8).gt.nearzero)then - write(fates_log(),*) 'The distribution of this host land model PFT :',hlm_pft - write(fates_log(),*) 'into FATES PFTs, does not add up to 1.0.' - write(fates_log(),*) 'Error is:',sumarea-1.0_r8 - write(fates_log(),*) 'and the hlm_pft_map is:', EDPftvarcon_inst%hlm_pft_map(1:npft,hlm_pft) - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - end do !hlm_pft + if((hlm_use_fixed_biogeog.eq.itrue) .or. (hlm_use_sp.eq.itrue)) then + do hlm_pft = 1,size( EDPftvarcon_inst%hlm_pft_map,2) + sumarea = sum(EDPftvarcon_inst%hlm_pft_map(1:npft,hlm_pft)) + if(abs(sumarea-1.0_r8).gt.nearzero)then + write(fates_log(),*) 'The distribution of this host land model PFT :',hlm_pft + write(fates_log(),*) 'into FATES PFTs, does not add up to 1.0.' + write(fates_log(),*) 'Error is:',sumarea-1.0_r8 + write(fates_log(),*) 'and the hlm_pft_map is:', EDPftvarcon_inst%hlm_pft_map(1:npft,hlm_pft) + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do !hlm_pft + end if + end do !ipft diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 7ddc717986..bcbbd2a180 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -389,8 +389,10 @@ module EDTypesMod ! Running means - class(rmean_type), pointer :: tveg_lpa ! exponential moving average of leaf temperature at the - ! leaf photosynthetic acclimation time-scale [K] + + ! (keeping this in-code as an example) + !class(rmean_type), pointer :: tveg_lpa ! exponential moving average of leaf temperature at the + ! leaf photosynthetic acclimation time-scale [K] end type ed_cohort_type diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index b606cbdb9e..b288781c4d 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -147,6 +147,7 @@ module FatesInterfaceMod public :: zero_bcs public :: set_bcs public :: UpdateFatesRMeansTStep + public :: InitTimeAveragingGlobals contains @@ -730,7 +731,7 @@ subroutine SetFatesGlobalElements(use_fates) if (use_fates) then - ! first read the non-PFT parameters + ! Self explanatory, read the fates parameter file call FatesReadParameters() ! Identify the number of PFTs by evaluating a pft array @@ -873,15 +874,7 @@ subroutine SetFatesGlobalElements(use_fates) ! These will not be used if use_ed or use_fates is false call fates_history_maps() - - ! Instantiate the time-averaging method globals - allocate(ema_24hr) - call ema_24hr%define(sec_per_day, hlm_stepsize, moving_ema_window) - allocate(fixed_24hr) - call fixed_24hr%define(sec_per_day, hlm_stepsize, fixed_window) - allocate(ema_lpa) - call ema_lpa%define(photo_temp_acclim_timescale*sec_per_day, & - hlm_stepsize,moving_ema_window) + else ! If we are not using FATES, the cohort dimension is still @@ -898,6 +891,27 @@ subroutine SetFatesGlobalElements(use_fates) end subroutine SetFatesGlobalElements + ! ====================================================================== + + subroutine InitTimeAveragingGlobals() + + ! Instantiate the time-averaging method globals + ! NOTE: It may be possible in the future that the HLM model timesteps + ! are dynamic in time or space, in that case, these would no longer + ! be global constants. + + allocate(ema_24hr) + call ema_24hr%define(sec_per_day, hlm_stepsize, moving_ema_window) + allocate(fixed_24hr) + call fixed_24hr%define(sec_per_day, hlm_stepsize, fixed_window) + allocate(ema_lpa) + call ema_lpa%define(photo_temp_acclim_timescale*sec_per_day, & + hlm_stepsize,moving_ema_window) + + return + end subroutine InitTimeAveragingGlobals + + ! ====================================================================== subroutine InitPARTEHGlobals() @@ -1241,7 +1255,6 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_numlevgrnd = unset_int hlm_name = 'unset' hlm_hio_ignore_val = unset_double - hlm_stepsize = unset_double hlm_masterproc = unset_int hlm_ipedof = unset_int hlm_nu_com = 'unset' @@ -1451,14 +1464,6 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if( abs(hlm_stepsize-unset_double) cpatch%tallest - do while (associated(ccohort)) - call ccohort%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) - ccohort => ccohort%shorter - end do + ! (Keeping as an example) + !ccohort => cpatch%tallest + !do while (associated(ccohort)) + ! call ccohort%tveg_lpa%UpdateRMean(bc_in(s)%t_veg_pa(ifp)) + ! ccohort => ccohort%shorter + !end do cpatch => cpatch%younger enddo diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 3bf021d004..507f01dbee 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -1071,9 +1071,10 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & prt_obj => null() call InitPRTObject(prt_obj) + ! (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()) + !allocate(temp_cohort%tveg_lpa) + !call temp_cohort%tveg_lpa%InitRMean(ema_lpa,init_value=cpatch%tveg_lpa%GetMean()) do el = 1,num_elements diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 6436802bea..c641f5096f 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -147,7 +147,9 @@ module FatesRestartInterfaceMod ! Running Means integer :: ir_tveg24_pa integer :: ir_tveglpa_pa - integer :: ir_tveglpa_co + + ! (Keeping as an example) + !!integer :: ir_tveglpa_co integer :: ir_ddbhdt_co integer :: ir_resp_tstep_co @@ -591,6 +593,10 @@ subroutine define_restart_vars(this, initialize_variables) ivar=0 + print*,"ABOUT TO DEFINE RESTARTS" + stop + + ! ----------------------------------------------------------------------------------- ! Site level variables ! ----------------------------------------------------------------------------------- @@ -1266,10 +1272,11 @@ subroutine define_restart_vars(this, initialize_variables) call this%DefineRMeanRestartVar(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_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) + + ! (Keeping as an example) + !call this%DefineRMeanRestartVar(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) ! Register all of the PRT states and fluxes @@ -2050,7 +2057,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) cohortsperpatch endif - call this%SetRMeanRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) + ! (Keeping as an example) + ! call this%SetRMeanRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) io_idx_co = io_idx_co + 1 @@ -2431,9 +2439,10 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in, bc_out) call InitHydrCohort(sites(s),new_cohort) end if + ! (Keeping as an example) ! Allocate running mean functions - allocate(new_cohort%tveg_lpa) - call new_cohort%tveg_lpa%InitRMean(ema_lpa) + !allocate(new_cohort%tveg_lpa) + !call new_cohort%tveg_lpa%InitRMean(ema_lpa) ! Update the previous @@ -2851,7 +2860,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end if - call this%GetRMeanRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) + ! (Keeping as an example) + !call this%GetRMeanRestartVar(ccohort%tveg_lpa, ir_tveglpa_co, io_idx_co) if (hlm_use_sp .eq. itrue) then ccohort%c_area = this%rvars(ir_c_area_co)%r81d(io_idx_co) From 657b135bae79f264df88d1f7e5e6e988eec18984 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 22 Oct 2021 15:18:00 -0600 Subject: [PATCH 16/21] Removing debug print --- main/FatesRestartInterfaceMod.F90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index c641f5096f..f51a231ba3 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -593,10 +593,6 @@ subroutine define_restart_vars(this, initialize_variables) ivar=0 - print*,"ABOUT TO DEFINE RESTARTS" - stop - - ! ----------------------------------------------------------------------------------- ! Site level variables ! ----------------------------------------------------------------------------------- From 010f8d34005aa41fd167a8f42d4f0ab806c878a1 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 23 Nov 2021 21:29:22 -0500 Subject: [PATCH 17/21] Reorder fixed window update --- main/FatesRunningMeanMod.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/main/FatesRunningMeanMod.F90 b/main/FatesRunningMeanMod.F90 index df1c12c7be..7fa3bfd7cc 100644 --- a/main/FatesRunningMeanMod.F90 +++ b/main/FatesRunningMeanMod.F90 @@ -263,17 +263,17 @@ subroutine UpdateRMean(this, new_value) ! end of the averaging memory period, and ! we are not using an indefinite running ! average, then zero things out + + 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 if(this%c_index == this%def_type%n_mem) then this%l_mean = this%c_mean this%c_mean = 0._r8 this%c_index = 0 - end if - - 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 + end if From 7042adad8ba79a2936afefd0830d79e0f1c195ae Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 30 Nov 2021 13:42:32 -0500 Subject: [PATCH 18/21] Removing unnecessary history updates during running mean update --- main/FatesInterfaceMod.F90 | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index c7c20d6bc8..daef389587 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1901,29 +1901,6 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in) enddo end do - ! Update running mean history variables - ! ------------------------------------------------------------------------------- - associate(hio_tveglpa_si_age => fates_hist%hvars(ih_tveglpa_si_age)%r82d, & - hio_tveglpa_si => fates_hist%hvars(ih_tveglpa_si)%r81d) - - do s = 1,size(sites,dim=1) - - io_si = sites(s)%h_gid - hio_tveglpa_si_age(io_si,:) = 0._r8 - hio_tveglpa_si(io_si) = 0._r8 - - cpatch => sites(s)%oldest_patch - do while(associated(cpatch)) - hio_tveglpa_si_age(io_si,cpatch%age_class) = & - hio_tveglpa_si_age(io_si,cpatch%age_class) + & - cpatch%tveg_lpa%GetMean()*cpatch%area/sites(s)%area_by_age(cpatch%age_class) - hio_tveglpa_si(io_si) = hio_tveglpa_si(io_si) + & - cpatch%tveg_lpa%GetMean()*cpatch%area*area_inv - cpatch => cpatch%younger - enddo - end do - end associate - return end subroutine UpdateFatesRMeansTStep From faa6313a4bbb2431627ffaa03c8bdb642a4550e3 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 30 Nov 2021 14:11:26 -0500 Subject: [PATCH 19/21] Adding history diagnostics to fates vegetation temperature --- main/FatesHistoryInterfaceMod.F90 | 39 ++++++++++++++++++++++--------- 1 file changed, 28 insertions(+), 11 deletions(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index e36792cc97..44116774fe 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -9,6 +9,7 @@ module FatesHistoryInterfaceMod use FatesConstantsMod , only : mg_per_kg use FatesConstantsMod , only : pi_const use FatesConstantsMod , only : nearzero + use FatesConstantsMod , only : t_water_freeze_k_1atm use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun use EDTypesMod , only : nclmax @@ -254,10 +255,6 @@ module FatesHistoryInterfaceMod integer :: ih_pefflux_scpf integer :: ih_pneed_scpf - integer :: ih_daily_temp - integer :: ih_daily_rh - integer :: ih_daily_prec - integer :: ih_bdead_si integer :: ih_balive_si integer :: ih_agb_si @@ -301,7 +298,8 @@ module FatesHistoryInterfaceMod integer :: ih_scorch_height_si_agepft ! Indices to (site) variables - + integer :: ih_tveg24_si + integer :: ih_tveg_si integer :: ih_nep_si integer :: ih_hr_si @@ -3594,10 +3592,10 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) hio_fabi_sun_top_si_can => this%hvars(ih_fabi_sun_top_si_can)%r82d, & hio_fabi_sha_top_si_can => this%hvars(ih_fabi_sha_top_si_can)%r82d, & hio_parsun_top_si_can => this%hvars(ih_parsun_top_si_can)%r82d, & - hio_parsha_top_si_can => this%hvars(ih_parsha_top_si_can)%r82d & - ) - - + hio_parsha_top_si_can => this%hvars(ih_parsha_top_si_can)%r82d, & + hio_tveg24 => this%hvars(ih_tveg24_si)%r81d, & + hio_tveg => this%hvars(ih_tveg_si)%r81d) + ! Flush the relevant history variables call this%flush_hvars(nc,upfreq_in=2) @@ -3641,9 +3639,14 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) hio_c_lblayer_si(io_si) = hio_c_lblayer_si(io_si) + & cpatch%c_lblayer * cpatch%total_canopy_area / umol_per_mol - hio_rad_error_si(io_si) = hio_rad_error_si(io_si) + & + hio_rad_error_si(io_si) = hio_rad_error_si(io_si) + & cpatch%radiation_error * cpatch%area * AREA_INV - + + hio_tveg24(io_si) = hio_tveg24(io_si) + & + (bc_in(s)%t_veg24_pa(cpatch%patchno)- t_water_freeze_k_1atm)*cpatch%area*area_inv + hio_tveg(io_si) = hio_tveg(io_si) + & + (bc_in(s)%t_veg_pa(cpatch%patchno) - t_water_freeze_k_1atm)*cpatch%area*area_inv + ccohort => cpatch%shortest do while(associated(ccohort)) @@ -5106,6 +5109,20 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=2, & ivar=ivar, initialize=initialize_variables, index = ih_c_lblayer_si) + ! Temperature + + call this%set_history_var(vname='FATES_TVEG24', units='degree_Celsius', & + long='fates 24-hr running mean vegetation temperature by site', & + use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=2, & + ivar=ivar, initialize=initialize_variables, index = ih_tveg24_si ) + + call this%set_history_var(vname='FATES_TVEG', units='degree_Celsius', & + long='fates instantaneous mean vegetation temperature by site', & + use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=2, & + ivar=ivar, initialize=initialize_variables, index = ih_tveg_si ) + ! radiation error call this%set_history_var(vname='FATES_RAD_ERROR', units='W m-2 ', & From d4b753ddbef0070b280b12128ce0a42daee29ca6 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 30 Nov 2021 15:11:29 -0500 Subject: [PATCH 20/21] Adding in history variables for tveg and tveg24 --- main/FatesHistoryInterfaceMod.F90 | 2 +- main/FatesInterfaceMod.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 44116774fe..eca0bb829f 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -3643,7 +3643,7 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) cpatch%radiation_error * cpatch%area * AREA_INV hio_tveg24(io_si) = hio_tveg24(io_si) + & - (bc_in(s)%t_veg24_pa(cpatch%patchno)- t_water_freeze_k_1atm)*cpatch%area*area_inv + (cpatch%tveg24%GetMean()- t_water_freeze_k_1atm)*cpatch%area*area_inv hio_tveg(io_si) = hio_tveg(io_si) + & (bc_in(s)%t_veg_pa(cpatch%patchno) - t_water_freeze_k_1atm)*cpatch%area*area_inv diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index daef389587..0753c89829 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -79,7 +79,7 @@ module FatesInterfaceMod use FatesRunningMeanMod , only : moving_ema_window use FatesRunningMeanMod , only : fixed_window use FatesHistoryInterfaceMod , only : fates_hist - use FatesHistoryInterfaceMod , only : ih_tveglpa_si_age,ih_tveglpa_si + ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg From 0dfab4189bb0aa86df13ea3aedba48e804bfbf67 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 1 Dec 2021 14:27:46 -0500 Subject: [PATCH 21/21] Minor updates to the tveg24 history diagnostic --- main/FatesHistoryInterfaceMod.F90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index eca0bb829f..20918251d2 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -2063,6 +2063,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_cleafon_si => this%hvars(ih_cleafon_si)%r81d, & hio_dleafoff_si => this%hvars(ih_dleafoff_si)%r81d, & hio_dleafon_si => this%hvars(ih_dleafoff_si)%r81d, & + hio_tveg24 => this%hvars(ih_tveg24_si)%r81d, & hio_meanliqvol_si => this%hvars(ih_meanliqvol_si)%r81d, & hio_cbal_err_fates_si => this%hvars(ih_cbal_err_fates_si)%r81d, & hio_err_fates_si => this%hvars(ih_err_fates_si)%r82d ) @@ -2195,6 +2196,10 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_area_si_age(io_si,cpatch%age_class) = hio_area_si_age(io_si,cpatch%age_class) & + cpatch%area * AREA_INV + ! 24hr veg temperature + hio_tveg24(io_si) = hio_tveg24(io_si) + & + (cpatch%tveg24%GetMean()- t_water_freeze_k_1atm)*cpatch%area*AREA_INV + ! Increment some patch-age-resolved diagnostics hio_lai_si_age(io_si,cpatch%age_class) = hio_lai_si_age(io_si,cpatch%age_class) & @@ -3593,7 +3598,6 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) hio_fabi_sha_top_si_can => this%hvars(ih_fabi_sha_top_si_can)%r82d, & hio_parsun_top_si_can => this%hvars(ih_parsun_top_si_can)%r82d, & hio_parsha_top_si_can => this%hvars(ih_parsha_top_si_can)%r82d, & - hio_tveg24 => this%hvars(ih_tveg24_si)%r81d, & hio_tveg => this%hvars(ih_tveg_si)%r81d) ! Flush the relevant history variables @@ -3642,8 +3646,6 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) hio_rad_error_si(io_si) = hio_rad_error_si(io_si) + & cpatch%radiation_error * cpatch%area * AREA_INV - hio_tveg24(io_si) = hio_tveg24(io_si) + & - (cpatch%tveg24%GetMean()- t_water_freeze_k_1atm)*cpatch%area*area_inv hio_tveg(io_si) = hio_tveg(io_si) + & (bc_in(s)%t_veg_pa(cpatch%patchno) - t_water_freeze_k_1atm)*cpatch%area*area_inv @@ -5114,7 +5116,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_TVEG24', units='degree_Celsius', & long='fates 24-hr running mean vegetation temperature by site', & use_default='active', & - avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=2, & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_tveg24_si ) call this%set_history_var(vname='FATES_TVEG', units='degree_Celsius', &