diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 12cab3db61..2baf1b322a 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -198,10 +198,6 @@ module EDPhysiologyMod ! a residual amount of leaves, which may create ! computational problems. The current threshold ! is the same used in ED-2.2. - - real(r8), parameter :: smp_lwr_bound = -1000000._r8 ! Imposed soil matric potential lower bound for - ! frozen or excessively dry soils, used when - ! computing water stress. ! ============================================================================ contains @@ -911,13 +907,11 @@ subroutine phenology( currentSite, bc_in ) ! Phenology. ! ! !USES: - use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm use EDParamsMod, only : ED_val_phen_a, ED_val_phen_b, ED_val_phen_c use EDParamsMod, only : ED_val_phen_chiltemp use EDParamsMod, only : ED_val_phen_mindayson use EDParamsMod, only : ED_val_phen_ncolddayslim use EDParamsMod, only : ED_val_phen_coldtemp - use EDBtranMod, only : check_layer_water ! ! !ARGUMENTS: type(ed_site_type), intent(inout), target :: currentSite @@ -926,25 +920,15 @@ subroutine phenology( currentSite, bc_in ) ! ! !LOCAL VARIABLES: - type(fates_patch_type),pointer :: cpatch integer :: model_day_int ! integer model day 1 - inf integer :: ncolddays ! no days underneath the threshold for leaf drop - integer :: i_wmem ! Loop counter for water mem days integer :: i_tmem ! Loop counter for veg temp mem days integer :: ipft ! plant functional type index - integer :: j ! Soil layer index real(r8) :: mean_10day_liqvol ! mean soil liquid volume over last 10 days [m3/m3] real(r8) :: mean_10day_smp ! mean soil matric potential over last 10 days [mm] - real(r8) :: leaf_c ! leaf carbon [kg] - real(r8) :: fnrt_c ! fineroot carbon [kg] - real(r8) :: sapw_c ! sapwood carbon [kg] - real(r8) :: store_c ! storage carbon [kg] - real(r8) :: struct_c ! structure carbon [kg] real(r8) :: gdd_threshold ! GDD accumulation function, - real(r8) :: rootfrac_notop ! Total rooting fraction excluding the top soil layer integer :: ncdstart ! beginning of counting period for chilling degree days. integer :: gddstart ! beginning of counting period for growing degree days. - integer :: nlevroot ! Number of rooting levels to consider real(r8) :: temp_in_C ! daily averaged temperature in celsius real(r8) :: temp_wgt ! canopy area weighting factor for daily average ! vegetation temperature calculation @@ -986,38 +970,16 @@ subroutine phenology( currentSite, bc_in ) logical :: last_flush_long_ago ! Has it been a very long time since last flushing? - ! This is the integer model day. The first day of the simulation is 1, and it - ! continues monotonically, indefinitely - ! Advance it. (this should be a global, no reason - ! for site level, but we don't have global scalars in the - ! restart file) - currentSite%phen_model_date = currentSite%phen_model_date + 1 + ! This is the elapsed number of days since the very beginning of the simulation time model_day_int = currentSite%phen_model_date - ! Parameter of drought decid leaf loss in mm in top layer...FIX(RF,032414) - ! - this is arbitrary and poorly understood. Needs work. ED_ !Parameters: defaults from Botta et al. 2000 GCB,6 709-725 !Parameters, default from from SDGVM model of senesence - temp_in_C = 0._r8 - temp_wgt = 0._r8 - cpatch => CurrentSite%oldest_patch - do while(associated(cpatch)) - temp_in_C = temp_in_C + cpatch%tveg24%GetMean()*cpatch%total_canopy_area - temp_wgt = temp_wgt + cpatch%total_canopy_area - cpatch => cpatch%younger - end do - if(temp_wgt>nearzero)then - temp_in_C = temp_in_C/temp_wgt - tfrz - else - ! If there is no canopy area, we use the veg temperature - ! of the first patch, which is the forcing air temperature - ! as defined in CLM/ELM. The forcing air temperature - ! should be the same among all patches. (Although - ! it is unlikely there are more than 1 in this scenario) - temp_in_C = CurrentSite%oldest_patch%tveg24%GetMean() - tfrz - end if + ! Retrieve average canopy temperature of the current day + temp_in_C = currentSite%vegtemp_memory(1) + !-----------------Cold Phenology--------------------! @@ -1047,10 +1009,6 @@ subroutine phenology( currentSite, bc_in ) ! -68 + 638 * (-0.001 * ncd) gdd_threshold = ED_val_phen_a + ED_val_phen_b*exp(ED_val_phen_c*real(currentSite%nchilldays,r8)) - !Accumulate temperature of last 10 days. - currentSite%vegtemp_memory(2:num_vegtemp_mem) = currentSite%vegtemp_memory(1:num_vegtemp_mem-1) - currentSite%vegtemp_memory(1) = temp_in_C - !count number of days for leaves off ncolddays = 0 do i_tmem = 1,num_vegtemp_mem @@ -1188,53 +1146,6 @@ subroutine phenology( currentSite, bc_in ) phen_moist_threshold = prt_params%phen_moist_threshold (ipft) phen_doff_time = prt_params%phen_doff_time (ipft) - - ! Update soil moisture information memory (we always track the last 10 days) - do i_wmem = numWaterMem,2,-1 !shift memory to previous day, to make room for current day - currentSite%liqvol_memory(i_wmem,ipft) = currentSite%liqvol_memory(i_wmem-1,ipft) - currentSite%smp_memory (i_wmem,ipft) = currentSite%smp_memory (i_wmem-1,ipft) - end do - - ! Find the rooting depth distribution for PFT - call set_root_fraction( currentSite%rootfrac_scr, ipft, currentSite%zi_soil, & - bc_in%max_rooting_depth_index_col ) - nlevroot = max(2,min(ubound(currentSite%zi_soil,1),bc_in%max_rooting_depth_index_col)) - - ! The top most layer is typically very thin (~ 2cm) and dries rather quickly. Despite - ! being thin, it can have a non-negligible rooting fraction (e.g., using - ! exponential_2p_root_profile with default parameters make the top layer to contain - ! about 7% of the total fine root density). To avoid overestimating dryness, we - ! ignore the top layer when calculating the memory. - rootfrac_notop = sum(currentSite%rootfrac_scr(2:nlevroot)) - if ( rootfrac_notop <= nearzero ) then - ! Unlikely, but just in case all roots are in the first layer, we use the second - ! layer the second layer (to avoid FPE issues). - currentSite%rootfrac_scr(2) = 1.0_r8 - rootfrac_notop = 1.0_r8 - end if - - ! Set the memory to be the weighted average of the soil properties, using the - ! root fraction of each layer (except the topmost one) as the weighting factor. - - currentSite%liqvol_memory(1,ipft) = sum( bc_in%h2o_liqvol_sl (2:nlevroot) * & - currentSite%rootfrac_scr(2:nlevroot) ) / & - rootfrac_notop - currentSite%smp_memory (1,ipft) = 0._r8 - do j = 2,nlevroot - if(check_layer_water(bc_in%h2o_liqvol_sl(j),bc_in%tempk_sl(j)) ) then - currentSite%smp_memory (1,ipft) = currentSite%smp_memory (1,ipft) + & - bc_in%smp_sl (j) * & - currentSite%rootfrac_scr(j) / & - rootfrac_notop - else - ! Nominal extreme suction for frozen or unreasonably dry soil - currentSite%smp_memory (1,ipft) = currentSite%smp_memory (1,ipft) + & - smp_lwr_bound * & - currentSite%rootfrac_scr(j) / & - rootfrac_notop - end if - end do - ! Calculate the mean soil moisture ( liquid volume (m3/m3) and matric potential (mm)) ! over the last 10 days mean_10day_liqvol = sum(currentSite%liqvol_memory(1:numWaterMem,ipft)) / & diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index bfba6d7d56..12f1cd79d8 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -106,6 +106,7 @@ module EDMainMod use EDPftvarcon, only : EDPftvarcon_inst use FatesHistoryInterfaceMod, only : fates_hist use FatesLandUseChangeMod, only: FatesGrazing + use FatesCumulativeMemoryMod, only : UpdateCumulativeMemoryVars ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -168,6 +169,10 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out) end do call currentSite%flux_diags%ZeroFluxDiags() + ! Call a routine that will compute cumulative variables and "memory" averages for + ! a suite of variables. These variables are mostly used for leaf phenology, but they + ! may be useful for other components (disturbances, mortality, management). + call UpdateCumulativeMemoryVars(currentSite,bc_in) ! Call a routine that simply identifies if logging should occur ! This is limited to a global event until more structured event handling is enabled diff --git a/main/FatesCumulativeMemoryMod.F90 b/main/FatesCumulativeMemoryMod.F90 new file mode 100644 index 0000000000..56b93e657d --- /dev/null +++ b/main/FatesCumulativeMemoryMod.F90 @@ -0,0 +1,228 @@ +module FatesCumulativeMemoryMod + !---~--- + ! This module contains procedures that update multiple cumulative/memory variables + ! that drive leaf phenology, and potentially mortality, disturbances, and management. + ! These sub-routines are related to FatesRunningMeanMod, however, they need a separate + ! module because many of these procedures require access to EDTypesMod, which in turn + ! depends on FatesRunningMeanMod. + !---~--- + + + + use EDBtranMod , only : check_layer_water + use EDTypesMod , only : area_inv + use EDTypesMod , only : ed_site_type + use EDTypesMod , only : num_vegtemp_mem + use EDTypesMod , only : numWaterMem + use EDTypesMod , only : phen_cstat_iscold + use FatesAllometryMod , only : set_root_fraction + use FatesConstantsMod , only : ndays_per_year + use FatesConstantsMod , only : nearzero + use FatesConstantsMod , only : r8 => fates_r8 + use FatesConstantsMod , only : tfrz => t_water_freeze_k_1atm + use FatesInterfaceTypesMod, only : bc_in_type + use FatesInterfaceTypesMod, only : hlm_day_of_year + use FatesInterfaceTypesMod, only : numpft + use FatesPatchMod , only : fates_patch_type + + implicit none + private + + !---~--- + ! The main sub-routine that updates all cumulative and memory-related variables. + ! This is the only sub-routine that needs to be public, the specific procedures that + ! update time, temperature and moisture memories are called by this sub-routine. + !---~--- + public :: UpdateCumulativeMemoryVars + + !---~--- + ! Imposed soil matric potential lower bound for frozen or excessively dry soils, + ! used when computing water stress. + !---~--- + real(r8), public, parameter :: smp_lwr_bound = -1000000._r8 + +contains + + subroutine UpdateCumulativeMemoryVars(currentSite,bc_in) + !---~--- + ! This sub-routine updates multiple variables that retain cumulative or "memory" + ! related information. These are mostly used for phenology, but some of the variables + ! can be useful for disturbances, mortality and management. + !---~--- + + ! Arguments + type(ed_site_type), intent(inout), target :: currentSite + type(bc_in_type), intent(in) :: bc_in + + + ! Update elapsed time + call UpdatePhenologyDate(currentSite) + + ! Update temperature-related variables. Currently this is just the temperature + ! memory, but it will eventually contain growing degree days and number of + ! chilling days. + call UpdateCumulativeThermal(currentSite,bc_in) + + ! Update moisture-related memory variables. + call UpdateMemoryMoisture(currentSite,bc_in) + + return + end subroutine UpdateCumulativeMemoryVars + + + + + + + subroutine UpdatePhenologyDate(currentSite) + !---~--- + ! This sub-routine simply updates the number of elapsed days in this simulation. + ! The first day of the simulation is 1, and it continues monotonically, + ! indefinitely. + !---~--- + + ! Arguments + type(ed_site_type), intent(inout), target :: currentSite + + !---~--- + ! Advance elapsed time. The only reason this is a site variable instead of a + ! global variable is that we need to save this information to the restart file, + ! and we do not have global scalars in the restart file. + !---~--- + currentSite%phen_model_date = currentSite%phen_model_date + 1 + + return + end subroutine UpdatePhenologyDate + + + + + + + subroutine UpdateCumulativeThermal(currentSite,bc_in) + !---~--- + ! Currently this sub-routine updates the canopy temperature memory variable only. + ! In the future, this will also host updates to growing degree days and number of + ! chilling days. + !---~--- + + ! Arguments + type(ed_site_type), intent(inout), target :: currentSite + type(bc_in_type), intent(in) :: bc_in + + ! Local variables + type(fates_patch_type), pointer :: cpatch ! Current patch + real(r8) :: temp_in_C ! Daily averaged temperature in degC + real(r8) :: temp_wgt ! Canopy area weighting factor for daily average + ! vegetation temperature calculation + + ! Find the average temperature for the previous 24 hours + temp_in_C = 0._r8 + temp_wgt = 0._r8 + cpatch => CurrentSite%oldest_patch + do while( associated(cpatch) ) + temp_in_C = temp_in_C + cpatch%tveg24%GetMean()*cpatch%total_canopy_area + temp_wgt = temp_wgt + cpatch%total_canopy_area + cpatch => cpatch%younger + end do + if ( temp_wgt > nearzero ) then + temp_in_C = temp_in_C/temp_wgt - tfrz + else + ! If there is no canopy area, we use the vegetation temperature + ! of the first patch, which is the forcing air temperature + ! as defined in CLM/ELM. The forcing air temperature + ! should be the same across all patches. (Although + ! it is unlikely there are more than 1 in this scenario) + temp_in_C = CurrentSite%oldest_patch%tveg24%GetMean() - tfrz + end if + + + ! Record temperature for the last 10 days to the memory variable. + currentSite%vegtemp_memory(2:num_vegtemp_mem) = currentSite%vegtemp_memory(1:num_vegtemp_mem-1) + currentSite%vegtemp_memory(1) = temp_in_C + + + return + end subroutine UpdateCumulativeThermal + + + + + + + subroutine UpdateMemoryMoisture(currentSite,bc_in) + !---~--- + ! This sub-routine updates the moisture-related memory variables. These variables + ! are used for leaf phenology, and may become useful for other modules too. + !---~--- + + ! Arguments + type(ed_site_type), intent(inout), target :: currentSite + type(bc_in_type), intent(in) :: bc_in + + ! Local variables + integer :: ipft ! Plant Functional Type index + integer :: i_wmem ! Loop counter for water memory days + integer :: j ! Soil layer index + integer :: nlevroot ! Number of rooting levels to consider + real(r8) :: rootfrac_notop ! Total rooting fraction excluding the top soil layer + + + + ! The soil memory variables are defined for each PFT + pft_memory_loop: do ipft=1,numpft + + ! Update soil moisture information memory (we always track the last 10 days). + ! We shift the memory from days to the previous day, and make room for current day + do i_wmem = numWaterMem,2,-1 + currentSite%liqvol_memory(i_wmem,ipft) = currentSite%liqvol_memory(i_wmem-1,ipft) + currentSite%smp_memory (i_wmem,ipft) = currentSite%smp_memory (i_wmem-1,ipft) + end do + + ! Find the rooting depth distribution for PFT + call set_root_fraction( currentSite%rootfrac_scr, ipft, currentSite%zi_soil, & + bc_in%max_rooting_depth_index_col ) + nlevroot = max(2,min(ubound(currentSite%zi_soil,1),bc_in%max_rooting_depth_index_col)) + + + ! The top most layer is typically very thin (~ 2cm) and dries rather quickly. Despite + ! being thin, it can have a non-negligible rooting fraction (e.g., using + ! exponential_2p_root_profile with default parameters make the top layer to contain + ! about 7% of the total fine root density). To avoid overestimating dryness, we + ! ignore the top layer when calculating the memory. + rootfrac_notop = sum(currentSite%rootfrac_scr(2:nlevroot)) + if ( rootfrac_notop <= nearzero ) then + ! Unlikely, but just in case all roots are in the first layer, we use the second + ! layer the second layer (to avoid FPE issues). + currentSite%rootfrac_scr(2) = 1.0_r8 + rootfrac_notop = 1.0_r8 + end if + + ! Set the memory to be the weighted average of the liquid volumetric soil water, + ! using the root fraction of each layer (except the topmost one) as the weighting + ! factor. + currentSite%liqvol_memory(1,ipft) = & + sum( bc_in%h2o_liqvol_sl(2:nlevroot) * currentSite%rootfrac_scr(2:nlevroot) ) / & + rootfrac_notop + + ! For the soil matric potential, we must check whether or not to include the layer + ! based on the layer temperature. + currentSite%smp_memory (1,ipft) = 0._r8 + root_loop: do j = 2,nlevroot + if ( check_layer_water(bc_in%h2o_liqvol_sl(j),bc_in%tempk_sl(j)) ) then + currentSite%smp_memory (1,ipft) = & + currentSite%smp_memory (1,ipft) + & + bc_in%smp_sl(j) * currentSite%rootfrac_scr(j) / rootfrac_notop + else + ! Nominal extreme suction for frozen or unreasonably dry soil + currentSite%smp_memory (1,ipft) = & + currentSite%smp_memory (1,ipft) + & + smp_lwr_bound * currentSite%rootfrac_scr(j) / rootfrac_notop + end if + end do root_loop + end do pft_memory_loop + + + return + end subroutine UpdateMemoryMoisture +end module FatesCumulativeMemoryMod diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 86d3bbeddb..7a2527e9ee 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -3079,8 +3079,6 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) real(r8) :: storep_understory_scpf(numpft*nlevsclass) real(r8) :: storec_canopy_scpf(numpft*nlevsclass) real(r8) :: storec_understory_scpf(numpft*nlevsclass) - real(r8) :: a_sapw ! sapwood area [m^2] - real(r8) :: c_sapw ! sapwood biomass [kgC] integer :: i_dist, j_dist @@ -8895,6 +8893,13 @@ subroutine define_history_vars(this, initialize_variables) ! over the short timestep. We turn off these variables when we want ! to save time (and some space) + call this%set_history_var(vname='FATES_NPP_AP', units='kg m-2 s-1', & + long='net primary productivity by age bin in kg carbon per m2 per second'// & + this%per_ageclass_norm_info('FATES_PATCHAREA/FATES_PATCHAREA_AP'), & + use_default='inactive', avgflag='A', vtype=site_age_r8, & + hlms='CLM:ALM', upfreq=group_hifr_complx, ivar=ivar, initialize=initialize_variables, & + index = ih_npp_si_age) + call this%set_history_var(vname='FATES_GPP_AP', units='kg m-2 s-1', & long='gross primary productivity by age bin in kg carbon per m2 per second'// & this%per_ageclass_norm_info('FATES_PATCHAREA/FATES_PATCHAREA_AP'), & diff --git a/testing/testing_shr/FatesFactoryMod.F90 b/testing/testing_shr/FatesFactoryMod.F90 index 81c0288a36..92a359a7b6 100644 --- a/testing/testing_shr/FatesFactoryMod.F90 +++ b/testing/testing_shr/FatesFactoryMod.F90 @@ -112,7 +112,7 @@ subroutine InitializeGlobals(step_size) do i = 2,nlevleaf dlower_vai(i) = dlower_vai(i-1) + dinc_vai(i-1) end do - + end subroutine InitializeGlobals !---------------------------------------------------------------------------------------