From 2d93f72802d1f412164f86d31d8fffa887535914 Mon Sep 17 00:00:00 2001 From: marius Date: Mon, 10 Oct 2022 14:05:32 +0200 Subject: [PATCH] implement hardening scheme, frost mortality impacts and hydraulic impacts --- biogeochem/EDCohortDynamicsMod.F90 | 8 ++ biogeochem/EDMortalityFunctionsMod.F90 | 159 +++++++++++++++++++-- biogeochem/EDPatchDynamicsMod.F90 | 2 +- biogeophys/FatesHydroWTFMod.F90 | 46 +++++- biogeophys/FatesPlantHydraulicsMod.F90 | 87 +++++++++-- biogeophys/FatesPlantRespPhotosynthMod.F90 | 63 ++++++-- main/EDInitMod.F90 | 6 + main/EDMainMod.F90 | 46 +++++- main/EDTypesMod.F90 | 9 +- main/FatesHistoryInterfaceMod.F90 | 62 ++++++++ main/FatesInterfaceMod.F90 | 46 ++++++ main/FatesInterfaceTypesMod.F90 | 19 ++- main/FatesRestartInterfaceMod.F90 | 27 ++++ 13 files changed, 538 insertions(+), 42 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 2d796a9048..d3d6c103f3 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -672,6 +672,8 @@ subroutine zero_cohort(cc_p) currentcohort%cambial_mort = 0._r8 currentCohort%c13disc_clm = 0._r8 currentCohort%c13disc_acc = 0._r8 + currentCohort%hard_level = -2._r8 + currentCohort%hard_level_prev = -2._r8 ! Daily nutrient fluxes are INTEGRATED over the course of the ! day. This variable MUST be zerod upon creation AND @@ -1381,6 +1383,10 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! their initization values, which should be the same for each if ( .not.currentCohort%isnew) then + currentCohort%hard_level = (currentCohort%n*currentCohort%hard_level + & + nextc%n*nextc%hard_level)/newn + currentCohort%hard_level_prev = (currentCohort%n*currentCohort%hard_level_prev + & + nextc%n*nextc%hard_level_prev)/newn currentCohort%seed_prod = (currentCohort%n*currentCohort%seed_prod + & nextc%n*nextc%seed_prod)/newn currentCohort%gpp_acc = (currentCohort%n*currentCohort%gpp_acc + & @@ -1868,6 +1874,8 @@ subroutine copy_cohort( currentCohort,copyc ) n%smort = o%smort n%asmort = o%asmort n%frmort = o%frmort + n%hard_level = o%hard_level + n%hard_level_prev = o%hard_level_prev ! logging mortalities, Yi Xu n%lmort_direct =o%lmort_direct diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 22e34e1abc..8bfa0bbe93 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -12,6 +12,7 @@ module EDMortalityFunctionsMod use EDTypesMod , only : ed_cohort_type use EDTypesMod , only : ed_site_type use EDTypesMod , only : ed_patch_type + use EDTypesMod , only : leaves_on use FatesConstantsMod , only : itrue,ifalse use FatesAllometryMod , only : bleaf use FatesAllometryMod , only : storage_fraction_of_target @@ -19,6 +20,12 @@ module EDMortalityFunctionsMod use FatesInterfaceTypesMod , only : hlm_use_ed_prescribed_phys use FatesInterfaceTypesMod , only : hlm_freq_day use FatesInterfaceTypesMod , only : hlm_use_planthydro + use FatesInterfaceTypesMod , only : hlm_model_day + use FatesInterfaceTypesMod , only : hlm_use_frosthard + use FatesInterfaceTypesMod , only : hlm_current_year + use FatesInterfaceTypesMod , only : hlm_current_month + use FatesInterfaceTypesMod , only : hlm_current_day + use PRTParametersMod , only : prt_params use EDLoggingMortalityMod , only : LoggingMortality_frac use EDParamsMod , only : fates_mortality_disturbance_fraction @@ -37,6 +44,7 @@ module EDMortalityFunctionsMod public :: mortality_rates public :: Mortality_Derivative public :: ExemptTreefallDist + public :: Hardening_scheme ! ============================================================================ @@ -49,7 +57,7 @@ module EDMortalityFunctionsMod - subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmort ) + subroutine mortality_rates( currentSite,cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmort ) ! ============================================================================ ! Calculate mortality rates from carbon storage, hydraulic cavitation, @@ -61,6 +69,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor use FatesConstantsMod, only : fates_check_param_set type (ed_cohort_type), intent(in) :: cohort_in + type (ed_site_type), intent(inout), target :: currentSite type (bc_in_type), intent(in) :: bc_in real(r8),intent(out) :: bmort ! background mortality : Fraction per year real(r8),intent(out) :: cmort ! carbon starvation mortality @@ -85,6 +94,10 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor real(r8) :: min_fmc_ar ! minimum fraction of maximum conductivity for absorbing root real(r8) :: min_fmc ! minimum fraction of maximum conductivity for whole plant real(r8) :: flc ! fractional loss of conductivity + real(r8) :: Tmin ! minimum daily average temperature + real(r8) :: hard_hydr ! reduced hydro mortality if hardened + real(r8) :: max_h !maximum hardiness level + real(r8), parameter :: min_h = -2.0_r8 !minimum hardiness level real(r8), parameter :: frost_mort_buffer = 5.0_r8 ! 5deg buffer for freezing mortality logical, parameter :: test_zero_mortality = .false. ! Developer test which ! may help to debug carbon imbalances @@ -140,9 +153,14 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor min_fmc = min(min_fmc_ag, min_fmc_tr) min_fmc = min(min_fmc, min_fmc_ar) flc = 1.0_r8-min_fmc - if(flc >= hf_flc_threshold .and. hf_flc_threshold < 1.0_r8 )then - hmort = (flc-hf_flc_threshold)/(1.0_r8-hf_flc_threshold) * & - EDPftvarcon_inst%mort_scalar_hydrfailure(cohort_in%pft) + if(flc >= hf_flc_threshold .and. hf_flc_threshold < 1.0_r8 )then + if (hlm_use_hydrohard .eq. itrue .and. currentSite%hard_level2(cohort_in%pft) < -3._r8) then + max_h=min(max(EDPftvarcon_inst%freezetol(cohort_in%pft),max(currentSite%hardtemp,-60._r8)-10._r8),min_h) + hard_hydr=EDPftvarcon_inst%mort_scalar_hydrfailure(cohort_in%pft)*( ((currentSite%hard_level2(cohort_in%pft)-max_h)/(min_h-max_h) ) *0.5_r8+0.5_r8) + else + hard_hydr=EDPftvarcon_inst%mort_scalar_hydrfailure(cohort_in%pft) + end if + hmort = (flc-hf_flc_threshold)/(1.0_r8-hf_flc_threshold) * hard_hydr else hmort = 0.0_r8 endif @@ -180,11 +198,20 @@ 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%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) ) - frmort = EDPftvarcon_inst%mort_scalar_coldstress(cohort_in%pft) * temp_dep_fraction + if (hlm_use_frosthard .eq. itrue) then !Hardiness level dependent frost mortality + Tmin=bc_in%tmin24_si-273.15_r8 + temp_dep_fraction = max(0.0_r8, min(1.0_r8,(-Tmin + & + max(currentSite%hard_level2(cohort_in%pft),EDPftvarcon_inst%mort_scalar_coldstress(cohort_in%pft)))/frost_mort_buffer)) + if (nint(hlm_model_day)<185) then + temp_dep_fraction=0._r8 + endif + frmort = EDPftvarcon_inst%mort_scalar_coldstress(cohort_in%pft)*temp_dep_fraction + else + temp_in_C = bc_in%t_veg24_pa(ifp) - 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 + endif !mortality_rates = bmort + hmort + cmort @@ -248,7 +275,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr ! Mortality for trees in the understorey. !if trees are in the canopy, then their death is 'disturbance'. This probably needs a different terminology - call mortality_rates(currentCohort,bc_in,cmort,hmort,bmort,frmort,smort, asmort) + call mortality_rates(currentSite,currentCohort,bc_in,cmort,hmort,bmort,frmort,smort, asmort) call LoggingMortality_frac(ipft, currentCohort%dbh, currentCohort%canopy_layer, & currentCohort%lmort_direct, & currentCohort%lmort_collateral, & @@ -293,7 +320,119 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr return end subroutine Mortality_Derivative +! ============================================================================ + + subroutine Hardening_scheme(currentSite,currentPatch,cohort_in,bc_in ) + ! + ! !DESCRIPTION: + ! Hardening module from Rammig et al. 2010. + ! Controls the unrealistic root water release by reducing root confuctivity, controls freezing mortality + + ! !ARGUMENTS + type (ed_site_type), intent(inout), target :: currentSite + type (ed_patch_type), intent(inout), target :: currentPatch + type (ed_cohort_type), intent(inout) :: cohort_in + type (bc_in_type), intent(in) :: bc_in + + + ! !LOCAL VARIABLES: + real(r8) :: Tmean ! Daily average temperature °C + real(r8) :: Tmin ! Daily minimum temperature °C + real(r8) :: max_h ! Maximum hardiness level + real(r8) :: max_h_dehard ! Maximum hardiness level for dehardening function + real(r8), parameter :: min_h = -2.0_r8 ! Minimum hardiness level from Bigras for Picea abies (°C) + real(r8) :: target_h ! Target hardiness + real(r8) :: rate_h ! Hardening rate + real(r8) :: rate_dh ! Dehardening rate + real(r8) :: hard_level_prev ! Temporary variable for the previous time-step hardiness level + real(r8) :: dayl_thresh ! Threshold for the onset of the hardening only period + integer :: ipft ! pft index + + Tmean=bc_in%temp24_si-273.15_r8 + if (Tmean<-200.0_r8) then + Tmean=0.0_r8 + endif + Tmin=bc_in%tmin24_si-273.15_r8 + + max_h=min(max(EDPftvarcon_inst%freezetol(cohort_in%pft),max(currentSite%hardtemp,-60._r8)-10._r8),min_h) + + !Calculation of the target hardiness + if (Tmean <= max_h/1.5_r8) then ! + target_h=max_h + else if (Tmean>= (6.0_r8-max_h/6._r8)) then + target_h=min_h + else + target_h = -sin((3.14159_r8*(0.5_r8+(Tmean-max_h/1.5_r8)/(-max_h/1.5_r8+(6.0_r8-max_h/6._r8)))))*(min_h-max_h)/2._r8-(min_h-max_h)/2._r8+min_h + end if + !Calculation of the hardening rate + if (Tmean <= max_h/2) then + rate_h=(max_h-min_h)/-31.11_r8+0.1_r8 + else if (Tmean >= 20.0_r8) then + rate_h=0.1_r8 + else + rate_h = sin((3.14159_r8*(0.5+(Tmean-max_h/2._r8)/(-max_h/2._r8+20._r8))))*((max_h-min_h)/-62.22_r8)+(((max_h-min_h)/-62.22_r8)+0.1_r8) + end if + !Calculation of the dehardening rate + if (max_h>-5.111_r8) then !this if steatement so that there is always possibility for 0.5 dehardening rate at 12 degreesC. + max_h_dehard=-5.111_r8 + else + max_h_dehard=max_h + end if + ipft = cohort_in%pft + if (prt_params%season_decid(ipft) == itrue) then + if (Tmean <= 5.0_r8) then + rate_dh=0.0_r8 + else if (Tmean >= 15.0_r8) then + rate_dh=5.0_r8*(max_h_dehard-min_h)/-31.11_r8 + else + rate_dh=(Tmean-5.0_r8)*((max_h_dehard-min_h)/-62.22_r8) + end if + else + if (Tmean <= 2.5_r8) then + rate_dh=0.0_r8 + else if (Tmean >= 12.5_r8) then + rate_dh=5.0_r8*(max_h_dehard-min_h)/-31.11_r8 + else + rate_dh=(Tmean-2.5_r8)*((max_h_dehard-min_h)/-62.22_r8) + end if + end if + + dayl_thresh= 42000.0_r8 + ( (-30.0_r8 - max(-60.0_r8,min(0.0_r8,currentSite%hardtemp)) )/15.0_r8) * 4500.0_r8 + !Hardening calculation + cohort_in%hard_level_prev = cohort_in%hard_level + + if (cohort_in%hard_level_prev + rate_dh > min_h) then + cohort_in%hard_level = min_h + else if (cohort_in%hard_level_prev >= target_h) then + cohort_in%hard_level = cohort_in%hard_level_prev - rate_h + else if (cohort_in%hard_level_prev < target_h) then + cohort_in%hard_level = cohort_in%hard_level_prev + rate_dh + end if + if (bc_in%dayl_si <= dayl_thresh .and. bc_in%dayl_si < bc_in%prev_dayl_si) then + if (cohort_in%hard_level_prev >= target_h) then + cohort_in%hard_level = cohort_in%hard_level_prev - rate_h + else + cohort_in%hard_level = cohort_in%hard_level_prev ! now dehardening is not possible in autumn but the hardiness level is also allowed to remain steady + end if + else if (bc_in%dayl_si > dayl_thresh .and. bc_in%dayl_si < bc_in%prev_dayl_si) then + cohort_in%hard_level = min_h + end if + ipft = cohort_in%pft + if (prt_params%season_decid(ipft) == itrue .and. cohort_in%status_coh == leaves_on .and. bc_in%dayl_si > bc_in%prev_dayl_si .and. & + (nint(hlm_model_day) >= currentSite%cleafondate .or. nint(hlm_model_day) >= currentSite%dleafondate)) then + cohort_in%hard_level = min_h + end if + if (cohort_in%hard_level > min_h) then + cohort_in%hard_level = min_h + end if + if (cohort_in%hard_level < max_h) then + cohort_in%hard_level = max_h + end if + + return + end subroutine Hardening_scheme + ! ============================================================================ function ExemptTreefallDist(ccohort) result(is_exempt) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 141aaad52e..e84f0640bf 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -207,7 +207,7 @@ subroutine disturbance_rates( site_in, bc_in) ! Mortality for trees in the understorey. currentCohort%patchptr => currentPatch - call mortality_rates(currentCohort,bc_in,cmort,hmort,bmort,frmort,smort,asmort) + call mortality_rates(site_in,currentCohort,bc_in,cmort,hmort,bmort,frmort,smort,asmort) currentCohort%dmort = cmort+hmort+bmort+frmort+smort+asmort call carea_allom(currentCohort%dbh,currentCohort%n,site_in%spread,currentCohort%pft, & currentCohort%c_area) diff --git a/biogeophys/FatesHydroWTFMod.F90 b/biogeophys/FatesHydroWTFMod.F90 index 8356b2b165..1b331019be 100644 --- a/biogeophys/FatesHydroWTFMod.F90 +++ b/biogeophys/FatesHydroWTFMod.F90 @@ -72,6 +72,7 @@ module FatesHydroWTFMod procedure :: dpsidth_from_th => dpsidth_from_th_base procedure :: set_wrf_param => set_wrf_param_base procedure :: get_thsat => get_thsat_base + procedure :: set_wrf_hard => set_wrf_hard_base ! All brands of WRFs have access to these tools to operate ! above and below sat and residual, should they want to @@ -226,6 +227,7 @@ module FatesHydroWTFMod procedure :: set_wrf_param => set_wrf_param_tfs procedure :: get_thsat => get_thsat_tfs procedure :: bisect_pv + procedure :: set_wrf_hard => set_wrf_cohort_hardening end type wrf_type_tfs ! Water Conductivity Function @@ -417,6 +419,14 @@ function dftcdpsi_from_psi_base(this,psi) result(dftcdpsi) write(fates_log(),*) 'check how the class pointer was setup' call endrun(msg=errMsg(sourcefile, __LINE__)) end function dftcdpsi_from_psi_base + subroutine set_wrf_hard_base(this,params_in) + class(wrf_type) :: this + real(r8),intent(in) :: params_in(:) + write(fates_log(),*) 'The base water retention function' + write(fates_log(),*) 'should never be actualized' + write(fates_log(),*) 'check how the class pointer was setup' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end subroutine set_wrf_hard_base ! ===================================================================================== ! Van Genuchten Functions are defined here @@ -747,10 +757,14 @@ function th_from_psi_cch(this,psi) result(th) class(wrf_type_cch) :: this real(r8), intent(in) :: psi real(r8) :: th - + real(r8) :: thmin + thmin=this%th_sat*(-20._r8/this%psi_sat)**(-1.0_r8/this%beta) !We want to avoid Soil water potentials below -25 + if(psi>this%psi_max) then ! Linear range for extreme values th = this%th_max + (psi-this%psi_max)/this%dpsidth_max + else if (psi<-20._r8) then + th = thmin + thmin/5._r8*(psi+20._r8) else th = this%th_sat*(psi/this%psi_sat)**(-1.0_r8/this%beta) end if @@ -764,9 +778,13 @@ function psi_from_th_cch(this,th) result(psi) class(wrf_type_cch) :: this real(r8),intent(in) :: th real(r8) :: psi - + real(r8) :: thmin + thmin=this%th_sat*(-20._r8/this%psi_sat)**(-1.0_r8/this%beta) + if(th>this%th_max) then psi = this%psi_max + this%dpsidth_max*(th-max_sf_interp*this%th_sat) + else if (th < thmin) then + psi = -20._r8 - 5._r8 * (thmin-th)/thmin else psi = this%psi_sat*(th/this%th_sat)**(-this%beta) end if @@ -778,12 +796,16 @@ end function psi_from_th_cch function dpsidth_from_th_cch(this,th) result(dpsidth) class(wrf_type_cch) :: this - real(r8),intent(in) :: th - real(r8) :: dpsidth - + real(r8),intent(in) :: th + real(r8) :: dpsidth + real(r8) :: thmin + thmin=this%th_sat*(-20._r8/this%psi_sat)**(-1.0_r8/this%beta) + ! Differentiate: if(th>this%th_max) then dpsidth = this%dpsidth_max + else if (th < thmin ) then + dpsidth = 5._r8/thmin else dpsidth = -this%beta*this%psi_sat/this%th_sat * (th/this%th_sat)**(-this%beta-1._r8) end if @@ -1494,6 +1516,20 @@ subroutine set_wrf_param_tfs(this,params_in) return end subroutine set_wrf_param_tfs + ! ===================================================================================== + ! Set the hardening changed variables in wrf + ! PV curve changes if plants harden or deharden + + subroutine set_wrf_cohort_hardening(this,params_in) + + class(wrf_type_tfs) :: this + real(r8), intent(in) :: params_in(:) + + this%pinot = params_in(1) + this%epsil = params_in(2) + + return + end subroutine set_wrf_cohort_hardening ! ===================================================================================== function get_thsat_tfs(this) result(th_sat) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 2225022e9f..b3eb801c69 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -27,6 +27,7 @@ module FatesPlantHydraulicsMod use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : fates_huge use FatesConstantsMod, only : denh2o => dens_fresh_liquid_water + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm use FatesConstantsMod, only : grav_earth use FatesConstantsMod, only : ifalse, itrue use FatesConstantsMod, only : pi_const @@ -60,6 +61,7 @@ module FatesPlantHydraulicsMod use FatesInterfaceTypesMod , only : bc_in_type use FatesInterfaceTypesMod , only : bc_out_type use FatesInterfaceTypesMod , only : hlm_use_planthydro + use FatesInterfaceTypesMod , only : hlm_use_hydrohard use FatesInterfaceTypesMod , only : hlm_ipedof use FatesInterfaceTypesMod , only : numpft use FatesInterfaceTypesMod , only : nlevsclass @@ -1300,7 +1302,9 @@ subroutine FuseCohortHydraulics(currentSite,currentCohort, nextCohort, bc_in, ne do k=1,n_hypool_ag vol_c1 = currentCohort%n*ccohort_hydr%th_ag(k)*ccohort_hydr%v_ag_init(k) vol_c2 = nextCohort%n*ncohort_hydr%th_ag(k)*ncohort_hydr%v_ag(k) - ccohort_hydr%th_ag(k) = (vol_c1+vol_c2)/(ccohort_hydr%v_ag(k)*newn) + if( (currentCohort%status_coh == leaves_on) .or. (k > n_hypool_leaf) ) then !Bug fix + ccohort_hydr%th_ag(k) = (vol_c1+vol_c2)/(ccohort_hydr%v_ag(k)*newn) + end if end do vol_c1 = currentCohort%n*ccohort_hydr%th_troot*ccohort_hydr%v_troot_init @@ -2425,7 +2429,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) type (ed_cohort_type), pointer :: ccohort ! current cohort pointer type(ed_site_hydr_type), pointer :: csite_hydr ! site hydraulics pointer type(ed_cohort_hydr_type), pointer :: ccohort_hydr ! cohort hydraulics pointer - + class(wrf_type_tfs), pointer :: wrf_tfs ! Local arrays ! accumulated water content change over all cohorts in a column [m3 m-3] @@ -2462,6 +2466,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) real(r8) :: sumcheck ! used to debug mass balance in soil horizon diagnostics integer :: nlevrhiz ! local for number of rhizosphere levels integer :: sc ! size class index + integer :: pm ! plant media index real(r8) :: lat,lon ! latitude and longitude of site real(r8) :: eff_por ! effective porosity real(r8) :: h2osoi_liqvol ! liquid water content [m3/m3] @@ -2472,6 +2477,8 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) real(r8) :: sum_l_aroot ! sum of root length of cohort, for disaggregation real(r8) :: rootfr ! fraction of root mass in soil layer, for disaggregation real(r8) :: z_fr ! Maximum fine root depth, used in disaggregation + real(r8) :: pinot_hard ! modified pinot parameter to make PV curve vary + real(r8) :: epsil_hard ! modified epsil parameter to make PV curve vary integer, parameter :: soilz_disagg = 0 ! disaggregate rhizosphere layers based on depth integer, parameter :: soilk_disagg = 1 ! disaggregate rhizosphere layers based on conductance @@ -2579,7 +2586,19 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) ccohort_hydr => ccohort%co_hydr ft = ccohort%pft - + + !update hardening for each cohort in BC hydraulics loop. + if (hlm_use_hydrohard .eq. itrue .and. ccohort%hard_level<-3._r8) then + if (ccohort%hard_level .ne. ccohort%hard_level_prev) then + do pm = 1,n_hypool_plant + pinot_hard=EDPftvarcon_inst%hydr_pinot_node(ft,pm)-(1._r8-(ccohort%hard_level+70._r8)/67._r8)*0.5_r8 ! solute pinot max change if hardening is -1 + epsil_hard=EDPftvarcon_inst%hydr_epsil_node(ft,pm)+(1._r8-(ccohort%hard_level+70._r8)/67._r8)*10._r8! pressure epsil max change if hardening is +15 + call wrf_plant(pm,ft)%p%set_wrf_hard([pinot_hard,epsil_hard]) + end do + end if + ccohort%hard_level_prev=ccohort%hard_level + end if + ! Relative transpiration of this cohort from the whole patch ! Note that g_sb_laweight / gscan_patch is the weighting that gives cohort contribution per area ! [mm H2O/plant/s] = [mm H2O/ m2 / s] * [m2 / patch] * [cohort/plant] * [patch/cohort] @@ -2632,14 +2651,14 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) call MatSolve2D(csite_hydr,ccohort,ccohort_hydr, & dtime,qflx_tran_veg_indiv, & sapflow,rootuptake(1:nlevrhiz),wb_err_plant,dwat_plant, & - dth_layershell_col) + dth_layershell_col,sites(s)%hard_level2(ft)) elseif(hydr_solver == hydr_solver_2DPicard) then call PicardSolve2D(csite_hydr,ccohort,ccohort_hydr, & dtime,qflx_tran_veg_indiv, & sapflow,rootuptake(1:nlevrhiz),wb_err_plant,dwat_plant, & - dth_layershell_col,csite_hydr%num_nodes) + dth_layershell_col,csite_hydr%num_nodes,sites(s)%hard_level2(ft)) elseif(hydr_solver == hydr_solver_1DTaylor ) then @@ -4750,7 +4769,7 @@ end subroutine Hydraulics_Tridiagonal subroutine MatSolve2D(csite_hydr,cohort,cohort_hydr, & tmx,qtop, & sapflow,rootuptake,wb_err_plant , dwat_plant, & - dth_layershell_site) + dth_layershell_site,hard_level) ! --------------------------------------------------------------------------------- @@ -4790,6 +4809,7 @@ subroutine MatSolve2D(csite_hydr,cohort,cohort_hydr, & type(ed_cohort_type) , intent(inout), target :: cohort real(r8),intent(in) :: tmx ! time interval to integrate over [s] real(r8),intent(in) :: qtop + real(r8),intent(in) :: hard_level real(r8),intent(out) :: sapflow ! time integrated mass flux between transp-root and stem [kg] real(r8),intent(out) :: rootuptake(:) ! time integrated mass flux between rhizosphere and aroot [kg] @@ -5125,7 +5145,7 @@ subroutine MatSolve2D(csite_hydr,cohort,cohort_hydr, & ! of each connection. This IS dependant on total potential h_node ! because of the root-soil radial conductance. - call SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_up) + call SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_up, hard_level) ! calculate boundary fluxes do icnx=1,csite_hydr%num_connections @@ -5526,7 +5546,7 @@ end function SumBetweenDepths subroutine PicardSolve2D(csite_hydr,cohort,cohort_hydr, & tmx,qtop, & sapflow,rootuptake,wb_err_plant , dwat_plant, & - dth_layershell_site,nnode) + dth_layershell_site,nnode,hard_level) ! --------------------------------------------------------------------------------- @@ -5565,6 +5585,7 @@ subroutine PicardSolve2D(csite_hydr,cohort,cohort_hydr, & type(ed_cohort_type) , intent(inout), target :: cohort real(r8),intent(in) :: tmx ! time interval to integrate over [s] real(r8),intent(in) :: qtop + real(r8),intent(in) :: hard_level integer :: nnode !total number of nodes real(r8),intent(out) :: sapflow ! time integrated mass flux between transp-root and stem [kg] real(r8),intent(out) :: rootuptake(:) ! time integrated mass flux between rhizosphere and aroot [kg] @@ -5848,7 +5869,7 @@ subroutine PicardSolve2D(csite_hydr,cohort,cohort_hydr, & ! of each connection. This IS dependant on total potential h_node ! because of the root-soil radial conductance. - call SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_up) + call SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_up, hard_level) ! calculate boundary fluxes do icnx=1,csite_hydr%num_connections @@ -5869,7 +5890,9 @@ subroutine PicardSolve2D(csite_hydr,cohort,cohort_hydr, & cfl_max = max(cfl_max,abs(k_eff*(h_node(id_dn) -h_node(id_up)))*dtime/volx/denh2o) enddo !Top node - cfl_max = max(cfl_max, abs(qtop * dtime/v_node(1)/denh2o)) + if(cohort%status_coh == leaves_on) then !bug fix + cfl_max = max(cfl_max, abs(qtop * dtime/v_node(1)/denh2o)) + end if ! To avoid extreme large clf_max due to large qtop from small gw weight cfl_max = min(20._r8,cfl_max) @@ -5956,7 +5979,7 @@ subroutine PicardSolve2D(csite_hydr,cohort,cohort_hydr, & ! of each connection. This IS dependant on total potential h_node ! because of the root-soil radial conductance. - call SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_up) + call SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_up, hard_level) ! calculate boundary fluxes do icnx=1,csite_hydr%num_connections @@ -6180,7 +6203,7 @@ end subroutine PicardSolve2D ! ===================================================================================== -subroutine SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_up) +subroutine SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_up, hard_level) ! ------------------------------------------------------------------------------- ! This subroutine sets the maximum conductances @@ -6195,11 +6218,16 @@ subroutine SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_ type(ed_site_hydr_type), intent(in),target :: csite_hydr type(ed_cohort_hydr_type), intent(in),target :: cohort_hydr real(r8),intent(in) :: h_node(:) ! Total (matric+height) potential at each node (Mpa) + real(r8),intent(in) :: hard_level ! Hardiness level of plants real(r8),intent(out) :: kmax_dn(:) ! Max conductance of downstream sides of connections (kg s-1 MPa-1) real(r8),intent(out) :: kmax_up(:) ! Max conductance of upstream sides of connections (kg s-1 MPa-1) real(r8):: aroot_frac_plant ! Fraction of the cohort's fine-roots ! out of the total in the current layer + real(r8):: k_factor + real(r8):: hard_rate_temporary + real(r8):: soil_rate_temporary + real(r8):: hard_max integer :: icnx ! connection index integer :: inode ! node index integer :: istem ! stem index @@ -6209,6 +6237,7 @@ subroutine SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_ kmax_dn(:) = fates_unset_r8 kmax_up(:) = fates_unset_r8 + k_factor=3._r8 ! Factor by which we reduce conductivity when plants are cold acclimated ! Set leaf to stem connections (only 1 leaf layer ! this will break if we have multiple, as there would ! need to be assumptions about which compartment @@ -6217,18 +6246,37 @@ subroutine SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_ kmax_dn(icnx) = cohort_hydr%kmax_petiole_to_leaf kmax_up(icnx) = cohort_hydr%kmax_stem_upper(1) + if (hlm_use_hydrohard .eq. itrue) then + if (hard_level<-3_r8) then + hard_rate_temporary=((hard_level+3._r8)/k_factor) + kmax_dn(icnx)=kmax_dn(icnx)*10**hard_rate_temporary + kmax_up(icnx)=kmax_up(icnx)*10**hard_rate_temporary + end if + end if ! Stem to stem connections do istem = 1,n_hypool_stem-1 icnx = icnx + 1 kmax_dn(icnx) = cohort_hydr%kmax_stem_lower(istem) kmax_up(icnx) = cohort_hydr%kmax_stem_upper(istem+1) + if (hlm_use_hydrohard .eq. itrue) then + if (hard_level<-3_r8) then + kmax_dn(icnx)=kmax_dn(icnx)*10**hard_rate_temporary + kmax_up(icnx)=kmax_up(icnx)*10**hard_rate_temporary + end if + end if enddo ! Path is between lowest stem and transporting root icnx = icnx + 1 kmax_dn(icnx) = cohort_hydr%kmax_stem_lower(n_hypool_stem) kmax_up(icnx) = cohort_hydr%kmax_troot_upper - + if (hlm_use_hydrohard .eq. itrue) then + if (hard_level<-3_r8) then + kmax_dn(icnx)=kmax_dn(icnx)*10**hard_rate_temporary + kmax_up(icnx)=kmax_up(icnx)*10**hard_rate_temporary + end if + end if + ! Path is between the transporting root and the absorbing roots inode = n_hypool_ag do j = 1,csite_hydr%nlevrhiz @@ -6242,6 +6290,12 @@ subroutine SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_ kmax_dn(icnx) = cohort_hydr%kmax_troot_lower(j) kmax_up(icnx) = cohort_hydr%kmax_aroot_upper(j) + if (hlm_use_hydrohard .eq. itrue) then + if (hard_level<-3_r8) then + kmax_dn(icnx)=kmax_dn(icnx)*10**hard_rate_temporary + kmax_up(icnx)=kmax_up(icnx)*10**hard_rate_temporary + end if + end if elseif( k == 2) then ! aroot-soil ! Special case. Maximum conductance depends on the @@ -6255,7 +6309,12 @@ subroutine SetMaxCondConnections(csite_hydr, cohort_hydr, h_node, kmax_dn, kmax_ 1._r8/cohort_hydr%kmax_aroot_radial_out(j)) end if kmax_up(icnx) = csite_hydr%kmax_upper_shell(j,1)*aroot_frac_plant - + if (hlm_use_hydrohard .eq. itrue) then + if (hard_level<-3_r8) then + kmax_dn(icnx)=kmax_dn(icnx)*10**hard_rate_temporary + kmax_up(icnx)=kmax_up(icnx)*10**hard_rate_temporary + end if + endif else ! soil - soil kmax_dn(icnx) = csite_hydr%kmax_lower_shell(j,k-2)*aroot_frac_plant kmax_up(icnx) = csite_hydr%kmax_upper_shell(j,k-1)*aroot_frac_plant diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 23722bab33..c87b22b1f8 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -34,9 +34,13 @@ module FATESPlantRespPhotosynthMod use FatesInterfaceTypesMod, only : hlm_parteh_mode use FatesInterfaceTypesMod, only : numpft use FatesInterfaceTypesMod, only : nleafage + use FatesInterfaceTypesMod, only : hlm_use_hydrohard + use FatesInterfaceTypesMod, only : hlm_use_planthydro + use FatesInterfaceTypesMod, only : hlm_model_day use EDTypesMod, only : maxpft use EDTypesMod, only : nlevleaf use EDTypesMod, only : nclmax + use EDTypesMod, only : leaves_on use PRTGenericMod, only : max_nleafage use EDTypesMod, only : do_fates_salinity use EDParamsMod, only : q10_mr @@ -239,6 +243,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) real(r8) :: lai_current ! the LAI in the current leaf layer real(r8) :: cumulative_lai ! the cumulative LAI, top down, to the leaf layer of interest real(r8) :: leaf_psi ! leaf xylem matric potential [MPa] (only meaningful/used w/ hydro) + real(r8) :: stomatal_intercept_hard ! stomatal intercept influenced by the hardiness + real(r8) :: hard_rate_temporary! temporary hardiness level + real(r8) :: k_factor ! rate at which we reduce stomatal parameters real(r8), allocatable :: rootfr_ft(:,:) ! Root fractions per depth and PFT ! ----------------------------------------------------------------------------------- @@ -423,7 +430,18 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if (hlm_use_planthydro.eq.itrue ) then - stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentCohort%co_hydr%btran ) + stomatal_intercept_hard = stomatal_intercept(ft) + if (hlm_use_hydrohard.eq.itrue .and. sites(s)%hard_level2(ft) < -3._r8 ) then + if (prt_params%season_decid(ft) == itrue .and. currentCohort%status_coh == leaves_on .and. & + (nint(hlm_model_day) >= sites(s)%cleafondate .or. nint(hlm_model_day) >= sites(s)%dleafondate)) then + else + k_factor=20._r8 + hard_rate_temporary=((sites(s)%hard_level2(ft)+3._r8)/k_factor) + stomatal_intercept_hard = stomatal_intercept(ft)*10**hard_rate_temporary + endif + end if + + stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept_hard*currentCohort%co_hydr%btran ) btran_eff = currentCohort%co_hydr%btran ! dinc_vai(:) is the total vegetation area index of each "leaf" layer @@ -541,7 +559,11 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! Part IX: This call calculates the actual photosynthesis for the ! leaf layer, as well as the stomatal resistance and the net assimilated carbon. - call LeafLayerPhotosynthesis(currentPatch%f_sun(cl,ft,iv), & ! in + call LeafLayerPhotosynthesis(sites(s)%hard_level2(ft), & ! in + currentCohort%status_coh , & ! in + sites(s)%cleafondate , & ! in + sites(s)%dleafondate , & ! in + currentPatch%f_sun(cl,ft,iv) , & ! in currentPatch%ed_parsun_z(cl,ft,iv), & ! in currentPatch%ed_parsha_z(cl,ft,iv), & ! in currentPatch%ed_laisun_z(cl,ft,iv), & ! in @@ -858,7 +880,11 @@ end subroutine FatesPlantRespPhotosynthDrive ! ======================================================================================= -subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in +subroutine LeafLayerPhotosynthesis(hard_level, & ! in + status_coh, & ! in + cleafondate, & ! in + dleafondate, & ! in + f_sun_lsl, & ! in parsun_lsl, & ! in parsha_lsl, & ! in laisun_lsl, & ! in @@ -982,6 +1008,14 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in real(r8) :: init_co2_inter_c ! First guess intercellular co2 specific to C path real(r8) :: term ! intermediate variable in Medlyn stomatal conductance model real(r8) :: vpd ! water vapor deficit in Medlyn stomatal model (KPa) + real(r8) :: hard_level ! hardiness level of plants + real(r8) :: stomatal_intercept_hard + real(r8) :: hard_rate_temporary + real(r8) :: k_factor + real(r8) :: bb_slope_hard + integer :: status_coh + integer :: dleafondate + integer :: cleafondate ! Parameters @@ -1010,7 +1044,20 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in associate( bb_slope => EDPftvarcon_inst%bb_slope ,& ! slope of BB relationship, unitless medlyn_slope=> EDPftvarcon_inst%medlyn_slope , & ! Slope for Medlyn stomatal conductance model method, the unit is KPa^0.5 stomatal_intercept=> EDPftvarcon_inst%stomatal_intercept ) !Unstressed minimum stomatal conductance, the unit is umol/m**2/s - + + stomatal_intercept_hard = stomatal_intercept(ft) + bb_slope_hard= bb_slope(ft) + if (hlm_use_hydrohard.eq.itrue .and. hard_level < -3._r8 ) then ! reduce stomatal conductance with hardiness + if (prt_params%season_decid(ft) == itrue .and. status_coh == leaves_on .and. & + (nint(hlm_model_day) >= cleafondate .or. nint(hlm_model_day) >= dleafondate)) then + else + k_factor=20._r8 + hard_rate_temporary=((hard_level+3._r8)/k_factor) + stomatal_intercept_hard = stomatal_intercept(ft)*10**hard_rate_temporary + bb_slope_hard= bb_slope(ft)*10**hard_rate_temporary + endif + end if + ! photosynthetic pathway: 0. = c4, 1. = c3 c3c4_path_index = nint(EDPftvarcon_inst%c3psn(ft)) @@ -1185,9 +1232,9 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in else if ( stomatal_model == ballberry_model ) then !stomatal conductance calculated from Ball et al. (1987) aquad = leaf_co2_ppress - bquad = leaf_co2_ppress*(gb_mol - stomatal_intercept_btran) - bb_slope(ft) * a_gs * can_press + bquad = leaf_co2_ppress*(gb_mol - stomatal_intercept_btran) - bb_slope_hard * a_gs * can_press cquad = -gb_mol*(leaf_co2_ppress*stomatal_intercept_btran + & - bb_slope(ft)*anet*can_press * ceair/ veg_esat ) + bb_slope_hard*anet*can_press * ceair/ veg_esat ) call quadratic_f (aquad, bquad, cquad, r1, r2) gs_mol = max(r1,r2) @@ -1261,7 +1308,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! Compare with Ball-Berry model: gs_mol = m * an * hs/leaf_co2_ppress*p + b else if ( stomatal_model == 1 ) then hs = (gb_mol*ceair + gs_mol* veg_esat ) / ((gb_mol+gs_mol)*veg_esat ) - gs_mol_err = bb_slope(ft)*max(anet, 0._r8)*hs/leaf_co2_ppress*can_press + stomatal_intercept_btran + gs_mol_err = bb_slope_hard*max(anet, 0._r8)*hs/leaf_co2_ppress*can_press + stomatal_intercept_btran end if if (abs(gs_mol-gs_mol_err) > 1.e-01_r8) then @@ -1290,7 +1337,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in psn_out = 0._r8 anet_av_out = 0._r8 - rstoma_out = min(rsmax0,cf/(stem_cuticle_loss_frac*stomatal_intercept(ft))) + rstoma_out = min(rsmax0,cf/(stem_cuticle_loss_frac*stomatal_intercept_hard)) c13disc_z = 0.0_r8 end if if_leafarea !is there leaf area? diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index d3177ca72a..9cda1325b6 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -198,6 +198,9 @@ subroutine zero_site( site_in ) site_in%cstatus = fates_unset_int ! are leaves in this pixel on or off? site_in%dstatus = fates_unset_int site_in%grow_deg_days = nan ! growing degree days + site_in%hardtemp = -2._r8 + site_in%hard_level2(:) = -2._r8 + site_in%Tmin_24_fates = 0.0_r8 site_in%snow_depth = nan site_in%nchilldays = fates_unset_int site_in%ncolddays = fates_unset_int @@ -327,6 +330,9 @@ subroutine set_site_properties( nsites, sites,bc_in ) sites(s)%dleafoffdate = dleafoff - hlm_day_of_year sites(s)%dleafondate = dleafon - hlm_day_of_year sites(s)%grow_deg_days = GDD + sites(s)%hardtemp = -2._r8 + sites(s)%hard_level2(1:numpft) = -2._r8 + sites(s)%Tmin_24_fates = 0.0_r8 sites(s)%water_memory(1:numWaterMem) = watermem sites(s)%vegtemp_memory(1:num_vegtemp_mem) = 0._r8 diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 7b9d623413..237e0c8201 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -26,6 +26,10 @@ module EDMainMod use FatesInterfaceTypesMod , only : hlm_masterproc use FatesInterfaceTypesMod , only : numpft use FatesInterfaceTypesMod , only : hlm_use_nocomp + use FatesInterfaceTypesMod , only : hlm_use_hydrohard + use FatesInterfaceTypesMod , only : hlm_use_frosthard + use FatesInterfaceTypesMod , only : hlm_model_day + use FatesHydroWTFMod , only : wrf_type, wrf_type_tfs use PRTGenericMod , only : prt_carbon_allom_hyp use PRTGenericMod , only : prt_cnp_flex_allom_hyp use PRTGenericMod , only : nitrogen_element @@ -83,7 +87,7 @@ module EDMainMod use EDPatchDynamicsMod , only : get_frac_site_primary use FatesGlobals , only : endrun => fates_endrun use ChecksBalancesMod , only : SiteMassStock - use EDMortalityFunctionsMod , only : Mortality_Derivative + use EDMortalityFunctionsMod , only : Mortality_Derivative,Hardening_scheme use EDTypesMod , only : AREA_INV use PRTGenericMod, only : carbon12_element use PRTGenericMod, only : all_carbon_elements @@ -304,13 +308,15 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) ! ! !USES: use FatesInterfaceTypesMod, only : hlm_use_cohort_age_tracking - use FatesConstantsMod, only : itrue + use FatesConstantsMod , only : itrue + use EDTypesMod , only : maxpft ! !ARGUMENTS: type(ed_site_type) , intent(inout) :: currentSite type(bc_in_type) , intent(in) :: bc_in type(bc_out_type) , intent(inout) :: bc_out + class(wrf_type_tfs), pointer :: wrf_tfs ! ! !LOCAL VARIABLES: type(site_massbal_type), pointer :: site_cmass @@ -328,12 +334,48 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) logical :: is_drought ! logical for if the plant (site) is in a drought state real(r8) :: delta_dbh ! correction for dbh real(r8) :: delta_hite ! correction for hite + real(r8) :: ncohort_pft(maxpft) + real(r8) :: number_fraction_pft real(r8) :: current_npp ! place holder for calculating npp each year in prescribed physiology mode !----------------------------------------------------------------------- real(r8) :: frac_site_primary + if (hlm_use_hydrohard.eq.itrue .or. hlm_use_frosthard.eq.itrue) then + currentSite%Tmin_24_fates=bc_in%tmin24_si-273.15_r8 + if (nint(hlm_model_day)>=366) then + currentSite%hardtemp=bc_in%t_mean_5yr_si-273.15_r8 + else if (nint(hlm_model_day)<366) then + currentSite%hardtemp=bc_in%t_min_yr_inst_si-273.15_r8 + end if + currentPatch => currentSite%youngest_patch + do while(associated(currentPatch)) + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + if ( (.not. currentCohort%isnew) .or. currentCohort%hard_level>-1.0_r8 ) then + ft = currentCohort%pft + !hard_level will be updated, ED_ecosystem_dynamics is called once a day at beginning of day + call Hardening_scheme( currentSite, currentPatch, currentCohort, bc_in ) + currentSite%hard_level2(ft) = currentCohort%hard_level + endif + currentCohort => currentCohort%taller + enddo ! cohort loop + currentPatch => currentPatch%older + end do !patch loop + + currentPatch => currentSite%youngest_patch + do while(associated(currentPatch)) + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + ft = currentCohort%pft + currentCohort%hard_level = currentSite%hard_level2(ft) + currentCohort => currentCohort%taller + enddo ! cohort loop + currentPatch => currentPatch%older + end do !patch loop + end if + call get_frac_site_primary(currentSite, frac_site_primary) ! Set a pointer to this sites carbon12 mass balance diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 86bb36e31e..ee98c30142 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -353,6 +353,8 @@ module EDTypesMod ! forest patch). fraction /per logging activity real(r8) :: seed_prod ! diagnostic seed production rate [kgC/plant/day] + real(r8) :: hard_level ! Hardiness state (K) + real(r8) :: hard_level_prev ! Hardiness state previous (K) ! NITROGEN POOLS ! ---------------------------------------------------------------------------------- @@ -750,7 +752,10 @@ module EDTypesMod real(r8) :: fdi ! daily probability an ignition event will start a fire real(r8) :: NF ! daily ignitions in km2 real(r8) :: NF_successful ! daily ignitions in km2 that actually lead to fire - + real(r8) :: hardtemp ! Temperature index for hardiness calculation + real(r8) :: hard_level2(1:maxpft) ! Temporary hardiness variable + real(r8) :: Tmin_24_fates + ! PLANT HYDRAULICS type(ed_site_hydr_type), pointer :: si_hydr @@ -1094,6 +1099,8 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%cambial_mort = ', ccohort%cambial_mort write(fates_log(),*) 'co%size_class = ', ccohort%size_class write(fates_log(),*) 'co%size_by_pft_class = ', ccohort%size_by_pft_class + write(fates_log(),*) 'co%hard_level = ', ccohort%hard_level + write(fates_log(),*) 'co%hard_level_prev = ', ccohort%hard_level_prev if (associated(ccohort%co_hydr) ) then call dump_cohort_hydr(ccohort) endif diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 391d5f87d9..e5b63ce62c 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -36,6 +36,8 @@ module FatesHistoryInterfaceMod use FatesHistoryVariableType , only : fates_history_variable_type use FatesInterfaceTypesMod , only : hlm_hio_ignore_val use FatesInterfaceTypesMod , only : hlm_use_planthydro + use FatesInterfaceTypesMod , only : hlm_use_hydrohard + use FatesInterfaceTypesMod , only : hlm_use_frosthard use FatesInterfaceTypesMod , only : hlm_use_ed_st3 use FatesInterfaceTypesMod , only : hlm_use_cohort_age_tracking use FatesInterfaceTypesMod , only : numpft @@ -311,6 +313,9 @@ module FatesHistoryInterfaceMod integer :: ih_tveg_si integer :: ih_nep_si integer :: ih_hr_si + integer :: ih_hardtemp_si + integer :: ih_Tmin_24_fates_si + integer :: ih_hardlevel_si_pft integer :: ih_c_stomata_si integer :: ih_c_lblayer_si @@ -1882,6 +1887,9 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_npp_froot_si => this%hvars(ih_npp_froot_si)%r81d, & hio_npp_croot_si => this%hvars(ih_npp_croot_si)%r81d, & hio_npp_stor_si => this%hvars(ih_npp_stor_si)%r81d, & + hio_hardlevel_si_pft => this%hvars(ih_hardlevel_si_pft)%r82d, & + hio_hardtemp_si => this%hvars(ih_hardtemp_si)%r81d, & + hio_Tmin_24_fates_si => this%hvars(ih_Tmin_24_fates_si)%r81d, & hio_bstor_canopy_si_scpf => this%hvars(ih_bstor_canopy_si_scpf)%r82d, & hio_bstor_understory_si_scpf => this%hvars(ih_bstor_understory_si_scpf)%r82d, & hio_bleaf_canopy_si_scpf => this%hvars(ih_bleaf_canopy_si_scpf)%r82d, & @@ -2132,6 +2140,33 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! Fire danger index (FDI) (0-1) hio_fire_fdi_si(io_si) = sites(s)%FDI + if (hlm_use_hydrohard .eq. itrue .or. hlm_use_frosthard .eq. itrue) then + hio_hardtemp_si(io_si) = sites(s)%hardtemp + hio_Tmin_24_fates_si(io_si) = sites(s)%Tmin_24_fates + ncohort_pft(:) = 0.0_r8 + ! Normalization counters + cpatch => sites(s)%oldest_patch + do while(associated(cpatch)) + ccohort => cpatch%shortest + do while(associated(ccohort)) + if ( .not. ccohort%isnew ) then + ft = ccohort%pft + ncohort_pft(ft) = ncohort_pft(ft) + ccohort%n + end if + ccohort => ccohort%taller + enddo ! cohort loop + cpatch => cpatch%younger + end do !patch loop + do ft = 1, numpft + if (ncohort_pft(ft)>0._r8)then + hio_hardlevel_si_pft(io_si,ft)=0._r8 + endif + enddo + else + hio_Tmin_24_fates_si(io_si) = 0.0_r8 + hio_hardtemp_si(io_si) = -2._r8 + end if + ! If hydraulics are turned on, track the error terms associated with ! dynamics [kg/m2] if(hlm_use_planthydro.eq.itrue)then @@ -2465,6 +2500,16 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! have any meaning, otherwise they are just inialization values notnew: if( .not.(ccohort%isnew) ) then + if (hlm_use_hydrohard .eq. itrue .or. hlm_use_frosthard .eq. itrue) then + ft = ccohort%pft + number_fraction_pft = (ccohort%n / ncohort_pft(ft)) + + hio_hardlevel_si_pft(io_si,ft) = hio_hardlevel_si_pft(io_si,ft) + & + ccohort%hard_level * number_fraction_pft + else + hio_hardlevel_si_pft(io_si,ft) = -2._r8 + end if + ! Turnover pools [kgC/day] * [day/yr] = [kgC/yr] sapw_m_turnover = ccohort%prt%GetTurnover(sapw_organ, carbon12_element) * days_per_year store_m_turnover = ccohort%prt%GetTurnover(store_organ, carbon12_element) * days_per_year @@ -4439,6 +4484,23 @@ subroutine define_history_vars(this, initialize_variables) ! Site level counting variables + + + call this%set_history_var(vname='PFThardiness', units='degC', & + long='Hardiness level of vegetation', use_default='active', & + avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_hardlevel_si_pft) + + call this%set_history_var(vname='hardtemp', units='DegC', & + long='5yr running mean temperature for hardiness', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_hardtemp_si ) + + call this%set_history_var(vname='Tmin_24_fates', units='DegC', & + long='Tmin_24_fates', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_Tmin_24_fates_si ) + call this%set_history_var(vname='FATES_NPATCHES', units='', & long='total number of patches per site', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', & diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 0ca5aafb63..0065e7fac6 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -257,6 +257,12 @@ subroutine zero_bcs(fates,s) ! Input boundaries + fates%bc_in(s)%temp24_si = 0.0_r8 + fates%bc_in(s)%t_mean_5yr_si = 0.0_r8 + fates%bc_in(s)%t_min_yr_inst_si = 0.0_r8 + fates%bc_in(s)%tmin24_si = 0.0_r8 + fates%bc_in(s)%dayl_si = 0.0_r8 + fates%bc_in(s)%prev_dayl_si = 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 @@ -1372,6 +1378,8 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_sf_successful_ignitions_def = unset_int hlm_sf_anthro_ignitions_def = unset_int hlm_use_planthydro = unset_int + hlm_use_hydrohard = unset_int + hlm_use_frosthard = unset_int hlm_use_lu_harvest = unset_int hlm_num_lu_harvest_cats = unset_int hlm_use_cohort_age_tracking = unset_int @@ -1419,6 +1427,32 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(), *) '' write(fates_log(), *) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' end if + + if ( .not.((hlm_use_hydrohard.eq.1).or.(hlm_use_hydrohard.eq.0)) ) then + if (fates_global_verbose()) then + write(fates_log(), *) 'The FATES namelist hydrohard flag must be 0 or 1, exiting' + end if + call endrun(msg=errMsg(sourcefile, __LINE__)) + elseif (hlm_use_hydrohard.eq.1 ) then + write(fates_log(), *) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(fates_log(), *) '' + write(fates_log(), *) ' use_fates_hydrohard is an EXPERIMENTAL FEATURE ' + write(fates_log(), *) '' + write(fates_log(), *) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + end if + + if ( .not.((hlm_use_frosthard.eq.1).or.(hlm_use_frosthard.eq.0)) ) then + if (fates_global_verbose()) then + write(fates_log(), *) 'The FATES namelist frosthard flag must be 0 or 1, exiting' + end if + call endrun(msg=errMsg(sourcefile, __LINE__)) + elseif (hlm_use_frosthard.eq.1 ) then + write(fates_log(), *) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(fates_log(), *) '' + write(fates_log(), *) ' use_fates_frosthard is an EXPERIMENTAL FEATURE ' + write(fates_log(), *) '' + write(fates_log(), *) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + end if if ( (hlm_use_lu_harvest .lt. 0).or.(hlm_use_lu_harvest .gt. 1) ) then write(fates_log(), *) 'The FATES lu_harvest flag must be 0 or 1, exiting' @@ -1782,6 +1816,18 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) if (fates_global_verbose()) then write(fates_log(),*) 'Transfering hlm_use_planthydro= ',ival,' to FATES' end if + + case('use_hydrohard') + hlm_use_hydrohard = ival + if (fates_global_verbose()) then + write(fates_log(),*) 'Transfering hlm_use_hydrohard= ',ival,' to FATES' + end if + + case('use_frosthard') + hlm_use_frosthard = ival + if (fates_global_verbose()) then + write(fates_log(),*) 'Transfering hlm_use_frosthard= ',ival,' to FATES' + end if case('use_lu_harvest') hlm_use_lu_harvest = ival diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index f08fba4ce8..4b0146abf0 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -140,7 +140,17 @@ module FatesInterfaceTypesMod integer, public :: hlm_use_planthydro ! This flag signals whether or not to use ! plant hydraulics (bchristo/xu methods) ! 1 = TRUE, 0 = FALSE - ! THIS IS CURRENTLY NOT SUPPORTED + ! THIS IS CURRENTLY NOT SUPPORTED + + integer, public :: hlm_use_hydrohard ! This flag signals whether or not to use + ! hydraulic hardening + ! 1 = TRUE, 0 = FALSE + ! THIS IS CURRENTLY NOT SUPPORTED + + integer, public :: hlm_use_frosthard ! This flag signals whether or not to use + ! frost hardening mortality + ! 1 = TRUE, 0 = FALSE + ! THIS IS CURRENTLY NOT SUPPORTED integer, public :: hlm_use_cohort_age_tracking ! This flag signals whether or not to use ! cohort age tracking. 1 = TRUE, 0 = FALSE @@ -482,6 +492,13 @@ module FatesInterfaceTypesMod real(r8) :: snow_depth_si ! Depth of snow in snowy areas of site (m) real(r8) :: frac_sno_eff_si ! Fraction of ground covered by snow (0-1) + real(r8) :: temp24_si + real(r8) :: t_mean_5yr_si + real(r8) :: t_min_yr_inst_si + real(r8) :: tmin24_si + real(r8) :: dayl_si + real(r8) :: prev_dayl_si + ! Hydrology variables for BTRAN ! --------------------------------------------------------------------------------- diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index a720d629f3..e298c4f4bd 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -20,6 +20,8 @@ module FatesRestartInterfaceMod use FatesInterfaceTypesMod, only : hlm_use_planthydro use FatesInterfaceTypesMod, only : hlm_use_sp use FatesInterfaceTypesMod, only : fates_maxElementsPerSite + use FatesInterfaceTypesMod, only : hlm_use_hydrohard + use FatesInterfaceTypesMod, only : hlm_use_frosthard use EDCohortDynamicsMod, only : UpdateCohortBioPhysRates use FatesHydraulicsMemMod, only : nshell use FatesHydraulicsMemMod, only : n_hypool_ag @@ -90,6 +92,8 @@ module FatesRestartInterfaceMod integer :: ir_phenmodeldate_si integer :: ir_acc_ni_si integer :: ir_gdd_si + integer :: ir_hard_level_co + integer :: ir_hard_level_prev_co integer :: ir_snow_depth_si integer :: ir_trunk_product_si integer :: ir_ncohort_pa @@ -683,6 +687,15 @@ subroutine define_restart_vars(this, initialize_variables) ! 1D cohort Variables ! ----------------------------------------------------------------------------------- + if (hlm_use_hydrohard .eq. itrue .or. hlm_use_frosthard .eq. itrue) then + call this%set_restart_var(vname='fates_hard_level', vtype=cohort_r8, & + long_name='hard_level', units='DegC', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hard_level_co ) + call this%set_restart_var(vname='fates_hard_level_prev', vtype=cohort_r8, & + long_name='hard_level', units='DegC', flushval = flushinvalid, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hard_level_prev_co ) + end if + call this%set_restart_var(vname='fates_seed_prod', vtype=cohort_r8, & long_name='fates cohort - seed production', units='kgC/plant', flushval = flushinvalid, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_seed_prod_co ) @@ -1746,6 +1759,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_phenmodeldate_si => this%rvars(ir_phenmodeldate_si)%int1d, & rio_acc_ni_si => this%rvars(ir_acc_ni_si)%r81d, & rio_gdd_si => this%rvars(ir_gdd_si)%r81d, & + rio_hard_level_co => this%rvars(ir_hard_level_co)%r81d, & + rio_hard_level_prev_co => this%rvars(ir_hard_level_prev_co)%r81d, & rio_snow_depth_si => this%rvars(ir_snow_depth_si)%r81d, & rio_trunk_product_si => this%rvars(ir_trunk_product_si)%r81d, & rio_ncohort_pa => this%rvars(ir_ncohort_pa)%int1d, & @@ -2005,6 +2020,11 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_asmort_co(io_idx_co) = ccohort%asmort rio_frmort_co(io_idx_co) = ccohort%frmort + if (hlm_use_hydrohard .eq. itrue .or. hlm_use_frosthard .eq. itrue) then + rio_hard_level_co(io_idx_co) = ccohort%hard_level + rio_hard_level_prev_co(io_idx_co) = ccohort%hard_level_prev + end if + ! Nutrient uptake/efflux rio_daily_no3_uptake_co(io_idx_co) = ccohort%daily_no3_uptake rio_daily_nh4_uptake_co(io_idx_co) = ccohort%daily_nh4_uptake @@ -2579,6 +2599,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_phenmodeldate_si => this%rvars(ir_phenmodeldate_si)%int1d, & rio_acc_ni_si => this%rvars(ir_acc_ni_si)%r81d, & rio_gdd_si => this%rvars(ir_gdd_si)%r81d, & + rio_hard_level_co => this%rvars(ir_hard_level_co)%r81d, & + rio_hard_level_prev_co => this%rvars(ir_hard_level_prev_co)%r81d, & rio_snow_depth_si => this%rvars(ir_snow_depth_si)%r81d, & rio_trunk_product_si => this%rvars(ir_trunk_product_si)%r81d, & rio_ncohort_pa => this%rvars(ir_ncohort_pa)%int1d, & @@ -2809,6 +2831,11 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%asmort = rio_asmort_co(io_idx_co) ccohort%frmort = rio_frmort_co(io_idx_co) + if (hlm_use_hydrohard .eq. itrue .or. hlm_use_frosthard .eq. itrue) then + ccohort%hard_level = rio_hard_level_co(io_idx_co) + ccohort%hard_level_prev = rio_hard_level_prev_co(io_idx_co) + end if + ! Nutrient uptake / efflux ccohort%daily_nh4_uptake = rio_daily_nh4_uptake_co(io_idx_co) ccohort%daily_no3_uptake = rio_daily_no3_uptake_co(io_idx_co)