Track climatology of leaf temp in grain-filling crop phase

src/biogeochem/CMakeLists.txt

@@ -10,6 +10,7 @@ list(APPEND clm_sources
   DustEmisLeung2023.F90
   DustEmisZender2003.F90
   DustEmisFactory.F90
+  CropHeatStress.F90
   CropReprPoolsMod.F90
   CropType.F90
   CNVegStateType.F90

src/biogeochem/CNPhenologyMod.F90

@@ -2014,8 +2014,9 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
     use clm_time_manager , only : get_prev_calday, get_curr_days_per_year, is_beg_curr_year
     use clm_time_manager , only : get_average_days_per_year
     use clm_time_manager , only : get_prev_date
-    use clm_time_manager , only : is_doy_in_interval, is_end_curr_day
+    use clm_time_manager , only : is_doy_in_interval, is_end_curr_day, is_beg_curr_day
     use clm_time_manager , only : get_doy_tomorrow
+    use CropHeatStress   , only : crop_heatstress_reset, crop_heatstress_ndays, calc_HS_factor     !added by SdR
     use pftconMod        , only : ntmp_corn, nswheat, nwwheat, ntmp_soybean
     use pftconMod        , only : nirrig_tmp_corn, nirrig_swheat, nirrig_wwheat, nirrig_tmp_soybean
     use pftconMod        , only : ntrp_corn, nsugarcane, ntrp_soybean, ncotton, nrice
@@ -2100,7 +2101,8 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
          a5tmin            =>    temperature_inst%t_a5min_patch                , & ! Input:  [real(r8) (:) ]  5-day running mean of min 2-m temperature         
          a10tmin           =>    temperature_inst%t_a10min_patch               , & ! Input:  [real(r8) (:) ]  10-day running mean of min 2-m temperature        
          gdd020            =>    temperature_inst%gdd020_patch                 , & ! Input:  [real(r8) (:) ]  20 yr mean of gdd0                                
-         gdd820            =>    temperature_inst%gdd820_patch                 , & ! Input:  [real(r8) (:) ]  20 yr mean of gdd8                                
+         gdd820            =>    temperature_inst%gdd820_patch                 , & ! Input:  [real(r8) (:) ]  20 yr mean of gdd8
+         t_veg_day         =>    temperature_inst%t_veg_day_patch             , & ! Input:  [real(r8)  (:) ] daily vegetation temperature; added by SdR
 
          fertnitro         =>    crop_inst%fertnitro_patch                     , & ! Input:  [real(r8) (:) ]  fertilizer nitrogen
          hui               =>    crop_inst%hui_patch                           , & ! Input:  [real(r8) (:) ]  crop patch heat unit index (growing degree-days); set to 0 at sowing and accumulated until harvest
@@ -2110,6 +2112,9 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
          vf                =>    crop_inst%vf_patch                            , & ! Output: [real(r8) (:) ]  vernalization factor                              
          sowing_count      =>    crop_inst%sowing_count                        , & ! Inout:  [integer  (:) ]  number of sowing events this year for this patch
          harvest_count     =>    crop_inst%harvest_count                       , & ! Inout:  [integer  (:) ]  number of harvest events this year for this patch
+         HS_ndays          =>    crop_inst%HS_ndays_patch                      , & ! Input:  [real(r8) (:) ]  number of crop heat stressed days; added by SdR
+         HS_factor         =>    crop_inst%HS_factor_patch                     , & ! Input:  [real(r8) (:) ]  heat stress factor ; added by SdR
+         heatwave_crop     =>    crop_inst%heatwave_crop_patch                 , & ! Input:  [real(r8) (:) ]  check if heatwave condition is true; added by SdR
          peaklai           =>    cnveg_state_inst%peaklai_patch                , & ! Output: [integer  (:) ]  1: max allowed lai; 0: not at max
          tlai              =>    canopystate_inst%tlai_patch                   , & ! Input:  [real(r8) (:) ]  one-sided leaf area index, no burying by snow     
 
@@ -2157,6 +2162,12 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
          g = patch%gridcell(p)
          h = inhemi(p)
 
+         ! added by SdR
+         if (is_beg_curr_day()) then
+            call crop_heatstress_ndays(HS_ndays(p), heatwave_crop(p), t_veg_day(p), croplive(p))
+            call calc_HS_factor(HS_factor(p), HS_ndays(p), t_veg_day(p),  croplive(p))
+         end if
+
          ! background litterfall and transfer rates; long growing season factor
 
          bglfr(p) = 0._r8 ! this value changes later in a crop's life cycle
@@ -2560,6 +2571,9 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
                   crop_inst%harvest_reason_thisyr_patch(p, harvest_count(p)) = harvest_reason
                endif
 
+               ! Reset heat-stress variables
+               call crop_heatstress_reset(HS_ndays(p), heatwave_crop(p))
+
                croplive(p) = .false.     ! no re-entry in greater if-block
                cphase(p) = cphase_harvest
                if (tlai(p) > 0._r8) then ! plant had emerged before harvest
@@ -2593,7 +2607,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
 
             else if (hui(p) >= huigrain(p)) then
                cphase(p) = cphase_grainfill
-               bglfr(p) = 1._r8/(leaf_long(ivt(p))*avg_dayspyr*secspday)
+               bglfr(p) = (1._r8/(leaf_long(ivt(p))*avg_dayspyr*secspday)) * HS_factor(p)
             end if
 
             ! continue fertilizer application while in phase 2;
@@ -2676,6 +2690,8 @@ subroutine CropPhase(bounds, num_pcropp, filter_pcropp, &
     do fp = 1, num_pcropp
        p = filter_pcropp(fp)
 
+
+
        if (croplive(p)) then
           ! Start with cphase_planted, but this might get changed in the later
           ! conditional blocks.

src/biogeochem/CropHeatStress.F90

@@ -0,0 +1,136 @@
+module CropHeatStress
+
+#include "shr_assert.h"
+
+  !-----------------------------------------------------------------------
+  ! !MODULE: CropHeatStress
+  !
+  ! !DESCRIPTION:
+  ! Module for implementing the effect of heat stress on crop production
+  ! Added by SdR. 
+  !
+  ! !USES:
+  !------------------------------------------------------------------------
+  use shr_kind_mod      , only : r8 => shr_kind_r8
+  use shr_log_mod       , only : errMsg => shr_log_errMsg
+  use shr_sys_mod       , only : shr_sys_flush
+  use shr_infnan_mod    , only : isnan => shr_infnan_isnan, isinf => shr_infnan_isinf, nan => shr_infnan_nan, assignment(=)
+  use clm_varcon, only : spval
+  !
+  implicit none
+  save
+  private
+  !
+  ! !PUBLIC MEMBER FUNCTIONS:
+  public  :: crop_heatstress_ndays       ! SdR: checks for number of days above Tcrit for crop heat stress
+  public  :: calc_HS_factor              ! SdR: calculates heat stress magnitude affecting leaf area decline (grainfill phase)
+  public  :: crop_heatstress_reset 
+
+  !
+  ! !PUBLIC FOR UNIT TESTING
+  real(r8), public, parameter :: tcrit = 302.15_r8
+  real(r8), public, parameter :: tmax = 318.15_r8
+  real(r8), public, parameter :: HS_ndays_min = 3._r8
+
+  character(len=*), parameter, private :: sourcefile = &
+       __FILE__
+  !------------------------------------------------------------------------
+
+contains
+
+  subroutine crop_heatstress_reset(HS_ndays, heatwave_crop)
+
+    ! !DESCRIPTION
+    ! Unsets variables related to crop heat stress
+
+    ! !ARGUMENTS:
+    real(r8), intent(inout)    :: HS_ndays ! number of crop heat stress days (ndays) should be integer at final implementation
+    real(r8), intent(inout)    :: heatwave_crop ! keep track if heatwave is activated
+
+    HS_ndays = 0._r8
+    heatwave_crop = 0._r8
+
+  end subroutine crop_heatstress_reset
+
+  !------------------------------------------------------------------------
+  subroutine crop_heatstress_ndays(HS_ndays, heatwave_crop, t_veg_day, croplive)
+
+    ! !DESCRIPTION:
+    ! added by SdR for heat stress implementation
+    ! function to keep track of critical temperature for crops that is exceeded at the end of each day
+    ! needs a minimum of 3 consecutive days before heat stress is activated
+
+    ! !ARGUMENTS:
+    real(r8),        intent(inout)    :: HS_ndays              ! number of crop heat stress days (ndays) should be integer at final implementation
+    real(r8),        intent(inout)    :: heatwave_crop         ! keep track if heatwave is activated
+    real(r8),        intent(in)       :: t_veg_day
+    logical,         intent(in)       :: croplive              ! crop between sowing and harvest
+
+    !----------------------------------------------------------------------
+
+    ! No heat stress if crop isn't alive
+    if (.not. croplive) then
+       call crop_heatstress_reset(HS_ndays, heatwave_crop)
+       return
+    end if
+
+    ! Don't do anything if it's not a real temperature
+    if (t_veg_day > spval / 1000._r8) then
+       return
+    end if
+
+    ! check if tcrit is exceeded and count days
+    if (t_veg_day >= tcrit) then
+         HS_ndays = HS_ndays + 1.0_r8
+    else
+         HS_ndays = 0.0_r8
+    end if
+
+    ! check if a heatwave is occurring
+    if (HS_ndays >= (HS_ndays_min - 0.2_r8)) then
+         heatwave_crop = 1.0_r8 
+    else
+         heatwave_crop = 0.0_r8
+    end if
+
+  end subroutine crop_heatstress_ndays
+
+
+  subroutine calc_HS_factor(HS_factor, HS_ndays, t_veg_day, croplive)
+
+    ! !DESCRIPTION:
+    ! function to calculate heat stress instensity by applying a factor to bglfr (increasing LAI decline). function based on Apsim-Nwheat model Asseng et al. 2011:https://doi.org/10.1111/j.1365-2486.2010.02262.x
+    ! WIP: different tcrit, tmax values for different crops or climatologies
+
+    ! !ARGUMENTS:
+    real(r8),        intent(inout)     :: HS_factor         ! keep track if heatwave is activated
+    real(r8),        intent(in)        :: HS_ndays          ! number of crop heat stress days (ndays) should be integer at final implementation
+    real(r8),        intent(in)        :: t_veg_day         ! daily vegetation temperature (Kelvin)
+    logical,         intent(in)        :: croplive          ! crop between sowing and harvest
+
+    ! !LOCAL VARIABLES:
+    integer  :: day_min
+
+    !-----------------------------------------------------------------------
+
+    day_min = 3
+
+    !check  if stress occurs
+    if (HS_ndays == day_min .and. croplive .and. t_veg_day >= tcrit) then
+       ! onset heatwave
+       HS_factor = 1.5_r8
+    else if (HS_ndays > day_min .and. croplive .and. t_veg_day >= tcrit) then
+      if (t_veg_day <= tmax ) then
+          HS_factor = 4 - (1 - (t_veg_day - tcrit)/2)
+      else
+          HS_factor = 11._r8
+      end if
+    else
+       HS_factor = 1._r8
+    end if
+
+
+  end subroutine calc_HS_factor 
+
+
+end module CropHeatStress

src/biogeochem/CropType.F90

@@ -73,6 +73,10 @@ module CropType
      ! achieved before the GDD threshold for grain fill has been reached; see CropPhenology().
      real(r8), pointer :: hui_patch               (:)   ! crop patch heat unit index (ddays)
      real(r8), pointer :: gddaccum_patch          (:)   ! patch growing degree-days from planting (air) (ddays)
+     ! added by SdR as part of HS implementation (20-08-24)
+     real(r8) , pointer :: HS_ndays_patch           (:)   ! patch day count for heat stress
+     real(r8) , pointer :: heatwave_crop_patch      (:)   ! check if heatwave is activated
+     real(r8) , pointer :: HS_factor_patch          (:)   ! patch day count for heat stress
 
    contains
      ! Public routines
@@ -255,6 +259,11 @@ subroutine InitAllocate(this, bounds)
     allocate(this%sowing_count(begp:endp)) ; this%sowing_count(:) = 0
     allocate(this%harvest_count(begp:endp)) ; this%harvest_count(:) = 0
 
+    ! added by SdR as part of heat stress implementation (22-08-24)
+    allocate(this%HS_ndays_patch           (begp:endp))                      ; this%HS_ndays_patch           (:)   = 0.0_r8
+    allocate(this%heatwave_crop_patch      (begp:endp))                      ; this%heatwave_crop_patch      (:)   = 0.0_r8
+    allocate(this%HS_factor_patch          (begp:endp))                      ; this%HS_factor_patch          (:)   = 1.0_r8
+
   end subroutine InitAllocate
 
   !-----------------------------------------------------------------------
@@ -364,6 +373,22 @@ subroutine InitHistory(this, bounds)
          avgflag='I', long_name='Reason for each crop harvest; should only be output annually', &
          ptr_patch=this%harvest_reason_thisyr_patch, default='inactive')
 
+    ! added by SdR for heat stress implementation
+    this%HS_ndays_patch(begp:endp) = 0._r8
+    call hist_addfld1d (fname='HS_NDAYS', units='ndays', &
+         avgflag='X', long_name='number of days with daily crop temperature above critical', &
+         ptr_patch=this%HS_ndays_patch, default='inactive')
+
+    this%heatwave_crop_patch(begp:endp) = 0._r8
+    call hist_addfld1d (fname='HW', units='boolean', &
+         avgflag='I', long_name='crop heatwave activated', &
+         ptr_patch=this%heatwave_crop_patch, default='inactive')
+
+    this%HS_factor_patch(begp:endp) = spval
+    call hist_addfld1d (fname='HSF', units='unitless', &
+         avgflag='A', long_name='crop stress factor', &
+         ptr_patch=this%HS_factor_patch, default='inactive')
+
     this%gdd20_baseline_patch(begp:endp) = spval
     call hist_addfld1d (fname='GDD20_BASELINE', units='ddays', &
          avgflag='I', long_name='Baseline mean growing-degree days accumulated during accumulation period (from input)', &
@@ -635,6 +660,16 @@ subroutine Restart(this, bounds, ncid, cnveg_state_inst, flag)
        end if
        deallocate(temp1d)
 
+       ! added by SdR as part of heatstress implementation
+       call restartvar(ncid=ncid, flag=flag,  varname='HS_ndays_patch',xtype=ncd_double, &
+            dim1name='pft', long_name='number of heatstressed days crop', &
+            units='ndays', &
+            interpinic_flag='interp', readvar=readvar, data=this%HS_ndays_patch)
+       call restartvar(ncid=ncid, flag=flag,  varname='HS_factor_patch',xtype=ncd_double, &
+            dim1name='pft', long_name='heat stress factor', &
+            units='unitless', &
+            interpinic_flag='interp', readvar=readvar, data=this%HS_factor_patch)
+
        call restartvar(ncid=ncid, flag=flag,  varname='harvdate', xtype=ncd_int,  &
             dim1name='pft', long_name='harvest date', units='jday', nvalid_range=(/1,366/), &
             interpinic_flag='interp', readvar=readvar, data=this%harvdate_patch)

src/biogeochem/test/CMakeLists.txt

@@ -3,3 +3,4 @@ add_subdirectory(CNVegComputeSeed_test)
 add_subdirectory(CNPhenology_test)
 add_subdirectory(Latbaset_test)
 add_subdirectory(DustEmis_test)
+add_subdirectory(CropHeatStress_test)

src/biogeochem/test/CropHeatStress_test/CMakeLists.txt

@@ -0,0 +1,6 @@
+set (pfunit_sources
+  test_CropHeatStress.pf)
+
+add_pfunit_ctest(CropHeatStress
+  TEST_SOURCES "${pfunit_sources}"
+  LINK_LIBRARIES clm csm_share esmf_wrf_timemgr)