NGEET · glemieux · Mar 21, 2023 · Mar 4, 2021 · Feb 3, 2023 · Feb 3, 2023
diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90
@@ -95,7 +95,7 @@ module EDPatchDynamicsMod
   use EDParamsMod,            only : maxpatch_primary
   use EDParamsMod,            only : maxpatch_secondary
   use EDParamsMod,            only : maxpatch_total
-  use FatesRunningMeanMod,    only : ema_24hr, fixed_24hr, ema_lpa
+  use FatesRunningMeanMod,    only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm
 
   ! CIME globals
   use shr_infnan_mod       , only : nan => shr_infnan_nan, assignment(=)
@@ -657,6 +657,7 @@ subroutine spawn_patches( currentSite, bc_in)
                       ! --------------------------------------------------------------------------
                       call new_patch%tveg24%CopyFromDonor(currentPatch%tveg24)
                       call new_patch%tveg_lpa%CopyFromDonor(currentPatch%tveg_lpa)
+                      call new_patch%tveg_longterm%CopyFromDonor(currentPatch%tveg_longterm)
 
 
                       ! --------------------------------------------------------------------------
@@ -2108,6 +2109,8 @@ subroutine create_patch(currentSite, new_patch, age, areap, label,nocomp_pft)
     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_veg)
+    allocate(new_patch%tveg_longterm)
+    call new_patch%tveg_longterm%InitRmean(ema_longterm,init_value=temp_init_veg)
 
     ! Litter
     ! Allocate, Zero Fluxes, and Initialize to "unset" values
@@ -2679,6 +2682,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)
+    call rp%tveg_longterm%FuseRMean(dp%tveg_longterm,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
@@ -3071,6 +3075,11 @@ subroutine dealloc_patch(cpatch)
        write(fates_log(),*) 'dealloc011: fail on deallocate(cpatch%tveg_lpa):'//trim(smsg)
        call endrun(msg=errMsg(sourcefile, __LINE__))
     endif
+    deallocate(cpatch%tveg_longterm, stat=istat, errmsg=smsg)
+    if (istat/=0) then
+       write(fates_log(),*) 'dealloc012: fail on deallocate(cpatch%tveg_longterm):'//trim(smsg)
+       call endrun(msg=errMsg(sourcefile, __LINE__))
+    endif
 
     return
   end subroutine dealloc_patch

diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90
@@ -32,6 +32,8 @@ module FATESPlantRespPhotosynthMod
   use FatesConstantsMod, only : fates_unset_r8
   use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm
   use FatesConstantsMod, only : nocomp_bareground
+  use FatesConstantsMod, only : photosynth_acclim_model_none
+  use FatesConstantsMod, only : photosynth_acclim_model_kumarathunge_etal_2019
   use FatesInterfaceTypesMod, only : hlm_use_planthydro
   use FatesInterfaceTypesMod, only : hlm_parteh_mode
   use FatesInterfaceTypesMod, only : numpft
@@ -55,10 +57,14 @@ module FATESPlantRespPhotosynthMod
   use PRTGenericMod,     only : store_organ
   use PRTGenericMod,     only : repro_organ
   use PRTGenericMod,     only : struct_organ
-  use EDParamsMod,       only : maintresp_nonleaf_baserate, stomatal_model, stomatal_assim_model
+  use EDParamsMod,       only : maintresp_nonleaf_baserate
+  use EDParamsMod,       only : stomatal_model
+  use EDParamsMod,       only : stomatal_assim_model
+  use EDParamsMod,       only : photo_tempsens_model
   use PRTParametersMod,  only : prt_params
-  use EDPftvarcon         , only : EDPftvarcon_inst
-
+  use EDPftvarcon      , only : EDPftvarcon_inst
+  use TemperatureType,   only : temperature_type
+
   ! CIME Globals
   use shr_log_mod , only      : errMsg => shr_log_errMsg
 
@@ -570,6 +576,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                                     currentCohort%kp25top,              &  ! in
                                     nscaler,                            &  ! in
                                     bc_in(s)%t_veg_pa(ifp),             &  ! in
+                                    currentPatch%tveg_lpa%GetMean(),    &  ! in
+                                    currentPatch%tveg_longterm%GetMean(),&  ! in
                                     btran_eff,                          &  ! in
                                     vcmax_z,                            &  ! out
                                     jmax_z,                             &  ! out
@@ -2186,6 +2194,8 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, &
    co2_rcurve_islope25top_ft, &
    nscaler,    &
    veg_tempk,      &
+   t_growth,   &
+   t_home,     &
    btran, &
    vcmax, &
    jmax, &
@@ -2219,6 +2229,8 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, &
    real(r8), intent(in) :: co2_rcurve_islope25top_ft ! initial slope of CO2 response curve
                                              ! (C4 plants) at 25C, canopy top, this pft
    real(r8), intent(in) :: veg_tempk           ! vegetation temperature
+   real(r8), intent(in) :: t_growth            ! T_growth (short-term running mean temperature) (K)
+   real(r8), intent(in) :: t_home              ! T_home (long-term running mean temperature) (K)
    real(r8), intent(in) :: btran           ! transpiration wetness factor (0 to 1)
 
    real(r8), intent(out) :: vcmax             ! maximum rate of carboxylation (umol co2/m**2/s)
@@ -2244,17 +2256,34 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, &
    real(r8) :: jmaxhd         ! deactivation energy for jmax (J/mol)
    real(r8) :: vcmaxse        ! entropy term for vcmax (J/mol/K)
    real(r8) :: jmaxse         ! entropy term for jmax (J/mol/K)
+   real(r8) :: t_growth_celsius ! average growing temperature
+   real(r8) :: t_home_celsius   ! average home temperature
+   real(r8) :: jvr            ! ratio of Jmax25 / Vcmax25
    real(r8) :: vcmaxc         ! scaling factor for high temperature inhibition (25 C = 1.0)
    real(r8) :: jmaxc          ! scaling factor for high temperature inhibition (25 C = 1.0)
 
-   vcmaxha = EDPftvarcon_inst%vcmaxha(FT)
-   jmaxha  = EDPftvarcon_inst%jmaxha(FT)
-
-   vcmaxhd = EDPftvarcon_inst%vcmaxhd(FT)
-   jmaxhd  = EDPftvarcon_inst%jmaxhd(FT)
-
-   vcmaxse = EDPftvarcon_inst%vcmaxse(FT)
-   jmaxse  = EDPftvarcon_inst%jmaxse(FT)
+   select case(photo_tempsens_model)
+   case (photosynth_acclim_model_none) !No temperature acclimation
+      vcmaxha = EDPftvarcon_inst%vcmaxha(FT)
+      jmaxha  = EDPftvarcon_inst%jmaxha(FT)
+      vcmaxhd = EDPftvarcon_inst%vcmaxhd(FT)
+      jmaxhd  = EDPftvarcon_inst%jmaxhd(FT)
+      vcmaxse = EDPftvarcon_inst%vcmaxse(FT)
+      jmaxse  = EDPftvarcon_inst%jmaxse(FT)
+   case (photosynth_acclim_model_kumarathunge_etal_2019) !Kumarathunge et al. temperature acclimation, Thome=30-year running mean
+      t_growth_celsius = t_growth-tfrz
+      t_home_celsius = t_home-tfrz
+      vcmaxha = (42.6_r8 + (1.14_r8*t_growth_celsius))*1e3_r8 !J/mol
+      jmaxha = 40.71_r8*1e3_r8 !J/mol
+      vcmaxhd = 200._r8*1e3_r8 !J/mol
+      jmaxhd = 200._r8*1e3_r8 !J/mol
+      vcmaxse = (645.13_r8 - (0.38_r8*t_growth_celsius))
+      jmaxse = 658.77_r8 - (0.84_r8*t_home_celsius) - 0.52_r8*(t_growth_celsius-t_home_celsius)
+      jvr = 2.56_r8 - (0.0375_r8*t_home_celsius)-(0.0202_r8*(t_growth_celsius-t_home_celsius))
+   case default
+      write (fates_log(),*)'error, incorrect leaf photosynthesis temperature acclimation model specified'
+      call endrun(msg=errMsg(sourcefile, __LINE__))
+   end select
 
    vcmaxc = fth25_f(vcmaxhd, vcmaxse)
    jmaxc  = fth25_f(jmaxhd, jmaxse)
@@ -2267,7 +2296,16 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, &
 
       ! Vcmax25top was already calculated to derive the nscaler function
       vcmax25 = vcmax25top_ft * nscaler
-      jmax25  = jmax25top_ft * nscaler
+      select case(photo_tempsens_model)
+      case (photosynth_acclim_model_none)
+         jmax25  = jmax25top_ft * nscaler
+      case (photosynth_acclim_model_kumarathunge_etal_2019) 
+         jmax25 = vcmax25*jvr
+      case default
+         write (fates_log(),*)'error, incorrect leaf photosynthesis temperature acclimation model specified'
+         call endrun(msg=errMsg(sourcefile, __LINE__))
+      end select
+
       co2_rcurve_islope25 = co2_rcurve_islope25top_ft * nscaler
 
       ! Adjust for temperature
@@ -2287,6 +2325,7 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, &
    vcmax = vcmax * btran
 
    return
+
  end subroutine LeafLayerBiophysicalRates
 
 subroutine lowstorage_maintresp_reduction(frac, pft, maintresp_reduction_factor)

diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90
@@ -23,19 +23,19 @@ module EDParamsMod
 
    real(r8),protected, public :: vai_top_bin_width           ! width in VAI units of uppermost leaf+stem
                                                              ! layer scattering element in each canopy layer [m2/m2]
-                                                             ! (NOT YET IMPLEMENTED)
    real(r8),protected, public :: vai_width_increase_factor   ! factor by which each leaf+stem scattering element
                                                              ! increases in VAI width (1 = uniform spacing)
-                                                             ! (NOT YET IMPLEMENTED)
    real(r8),protected, public :: photo_temp_acclim_timescale ! Length of the window for the exponential moving average (ema)
-                                                             ! of vegetation temperature used in photosynthesis
-                                                             ! temperature acclimation (NOT YET IMPLEMENTED)
-
+                                                             ! of vegetation temperature used in photosynthesis and respiration
+                                                             ! temperature acclimation [days]
+   real(r8),protected, public :: photo_temp_acclim_thome_time ! Length of the window for the long-term exponential moving average (ema)
+                                                              ! of vegetation temperature used in photosynthesis 
+                                                              ! T_home term in Kumarathunge parameterization [years]
    integer,protected, public :: maintresp_leaf_model  ! switch for choosing between leaf maintenance
                                                       ! respiration model. 1=Ryan (1991), 2=Atkin et al (2017)
    integer,protected, public :: photo_tempsens_model  ! switch for choosing the model that defines the temperature
                                                       ! sensitivity of photosynthetic parameters (vcmax, jmax).
-                                                      ! 1=non-acclimating (NOT YET IMPLEMENTED)
+                                                      ! 1=non-acclimating, 2=Kumarathunge et al., 2019
 
    real(r8),protected, public :: fates_mortality_disturbance_fraction ! the fraction of canopy mortality that results in disturbance
    real(r8),protected, public :: ED_val_comp_excln
@@ -94,6 +94,7 @@ module EDParamsMod
    integer, protected,allocatable,public :: hydr_htftype_node(:)
 
    character(len=param_string_length),parameter,public :: ED_name_photo_temp_acclim_timescale = "fates_leaf_photo_temp_acclim_timescale"
+   character(len=param_string_length),parameter,public :: ED_name_photo_temp_acclim_thome_time = "fates_leaf_photo_temp_acclim_thome_time"
    character(len=param_string_length),parameter,public :: name_photo_tempsens_model = "fates_leaf_photo_tempsens_model"
    character(len=param_string_length),parameter,public :: name_maintresp_model = "fates_maintresp_leaf_model"
    character(len=param_string_length),parameter,public :: ED_name_hydr_htftype_node = "fates_hydro_htftype_node"
@@ -264,6 +265,7 @@ subroutine FatesParamsInit()
     vai_top_bin_width                     = nan
     vai_width_increase_factor             = nan
     photo_temp_acclim_timescale           = nan
+    photo_temp_acclim_thome_time          = nan
     photo_tempsens_model                  = -9
     maintresp_leaf_model                  = -9
     fates_mortality_disturbance_fraction  = nan
@@ -348,6 +350,9 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=ED_name_photo_temp_acclim_timescale, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
+    call fates_params%RegisterParameter(name=ED_name_photo_temp_acclim_thome_time, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
     call fates_params%RegisterParameter(name=name_photo_tempsens_model,dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
@@ -443,7 +448,7 @@ subroutine FatesRegisterParams(fates_params)
 
     call fates_params%RegisterParameter(name=hydr_name_solver, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
-    
+
     call fates_params%RegisterParameter(name=hydr_name_kmax_rsurf1, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
@@ -551,6 +556,9 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetrieveParameter(name=ED_name_photo_temp_acclim_timescale, &
          data=photo_temp_acclim_timescale)
 
+    call fates_params%RetrieveParameter(name=ED_name_photo_temp_acclim_thome_time, &
+         data=photo_temp_acclim_thome_time)
+
     call fates_params%RetrieveParameter(name=name_photo_tempsens_model, &
          data=tmpreal)
     photo_tempsens_model = nint(tmpreal)
@@ -768,7 +776,8 @@ subroutine FatesReportParams(is_master)
         write(fates_log(),*) '-----------  FATES Scalar Parameters -----------------'
         write(fates_log(),fmt0) 'vai_top_bin_width = ',vai_top_bin_width
         write(fates_log(),fmt0) 'vai_width_increase_factor = ',vai_width_increase_factor
-        write(fates_log(),fmt0) 'photo_temp_acclim_timescale = ',photo_temp_acclim_timescale
+        write(fates_log(),fmt0) 'photo_temp_acclim_timescale (days) = ',photo_temp_acclim_timescale
+        write(fates_log(),fmt0) 'photo_temp_acclim_thome_time (years) = ',photo_temp_acclim_thome_time
         write(fates_log(),fmti) 'hydr_htftype_node = ',hydr_htftype_node
         write(fates_log(),fmt0) 'fates_mortality_disturbance_fraction = ',fates_mortality_disturbance_fraction
         write(fates_log(),fmt0) 'ED_val_comp_excln = ',ED_val_comp_excln

diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90
@@ -440,6 +440,8 @@ module EDTypesMod
      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 [K]
+     class(rmean_type), pointer :: tveg_longterm                ! Long-Term Running mean of vegetation temperature at the
+                                                                 ! leaf photosynthesis acclimation timescale [K] (i.e T_home)
 
      integer  ::  nocomp_pft_label                               ! Where nocomp is active, use this label for patch ID.
                                                                  ! Each patch ID corresponds to a pft number since each

diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90
@@ -65,6 +65,10 @@ module FatesConstantsMod
                                                             ! based on average age of global
                                                             ! secondary 1900s land in hurtt-2011
 
+  ! integer labels for specifying harvest units
+  integer, parameter, public :: photosynth_acclim_model_none = 1
+  integer, parameter, public :: photosynth_acclim_model_kumarathunge_etal_2019 = 2
+
   ! integer labels for specifying harvest units
   integer, parameter, public :: hlm_harvest_area_fraction = 1 ! Code for harvesting by area
   integer, parameter, public :: hlm_harvest_carbon = 2 ! Code for harvesting based on carbon extracted.