Refactor 20-year running means of crop GDD accumulation

doc/ChangeLog

@@ -1,4 +1,67 @@
 ===============================================================
+Tag name: ctsm5.1.dev161
+Originator(s): samrabin (Sam Rabin, UCAR/TSS, [email protected])
+Date: Thu Jan  4 09:00:28 MST 2024
+One-line Summary: Refactor 20-year running means of crop GDD accumulation
+
+Purpose and description of changes
+----------------------------------
+
+Three variables track the 20-year running mean of GDD accumulation (base temperatures 0, 8, and 10°C) during the "growing season" (April through September in the Northern Hemisphere, October through March in the Southern Hemisphere). This PR refactors those to use accumulMod, resolving overly-strong weighting of the first few years after a crop becomes active.
+
+
+Significant changes to scientifically-supported configurations
+--------------------------------------------------------------
+
+Does this tag change answers significantly for any of the following physics configurations?
+(Details of any changes will be given in the "Answer changes" section below.)
+
+    [Put an [X] in the box for any configuration with significant answer changes.]
+
+[X] clm5_1
+
+[X] clm5_0
+
+[X] ctsm5_0-nwp
+
+[X] clm4_5
+
+
+Bugs fixed
+----------
+CTSM issues fixed (include CTSM Issue #):
+* Fixes #75.
+
+
+Testing summary:
+----------------
+
+  regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing):
+
+    derecho ----- DIFF
+    izumi ------- DIFF
+
+
+Answer changes
+--------------
+
+Changes answers relative to baseline: YES
+
+  Summarize any changes to answers, i.e.,
+    - what code configurations: All crop configurations
+    - what platforms/compilers: All
+    - nature of change (roundoff; larger than roundoff/same climate; new climate): larger than roundoff/same climate
+
+GDD020, GDD820, and GDD1020 will differ, most strongly in the first few years after a crop becomes active. This will have downstream effects on lots of variables, since those are used in determining sowing date and maturity requirements.
+
+Other details
+-------------
+
+Pull Requests that document the changes (include PR ids):
+* #2060: Refactor 20-year running means of crop GDD accumulation (https://github.com/ESCOMP/CTSM/pull/2060)
+
+===============================================================
+===============================================================
 Tag name: ctsm5.1.dev160
 Originator(s): glemieux (Gregory Lemieux, LBNL, [email protected])
 Date: Sat 30 Dec 2023 11:23:47 PM MST

doc/ChangeSum

@@ -1,5 +1,6 @@
 Tag                   Who      Date  Summary
 ============================================================================================================================
+    ctsm5.1.dev161 samrabin 01/04/2024 Refactor 20-year running means of crop GDD accumulation
     ctsm5.1.dev160 glemieux 12/30/2023 FATES landuse version 1
     ctsm5.1.dev159 multiple 12/12/2023 Various BFB fixes and updates
     ctsm5.1.dev158     erik 12/07/2023 First tag with testing moved to Derecho and working PE-layouts for Derecho

src/biogeochem/CNPhenologyMod.F90

@@ -369,7 +369,7 @@ subroutine CNPhenology (bounds, num_soilc, filter_soilc, num_soilp, &
 
 
     if ( phase == 1 ) then
-       call CNPhenologyClimate(num_soilp, filter_soilp, num_pcropp, filter_pcropp, &
+       call CNPhenologyClimate(num_soilp, filter_soilp, &
             temperature_inst, cnveg_state_inst, crop_inst)
 
        call CNEvergreenPhenology(num_soilp, filter_soilp, &
@@ -531,7 +531,7 @@ subroutine CNPhenologyInit(bounds)
   end subroutine CNPhenologyInit
 
   !-----------------------------------------------------------------------
-  subroutine CNPhenologyClimate (num_soilp, filter_soilp, num_pcropp, filter_pcropp, &
+  subroutine CNPhenologyClimate (num_soilp, filter_soilp, &
        temperature_inst, cnveg_state_inst, crop_inst)
     !
     ! !DESCRIPTION:
@@ -544,8 +544,6 @@ subroutine CNPhenologyClimate (num_soilp, filter_soilp, num_pcropp, filter_pcrop
     ! !ARGUMENTS:
     integer                , intent(in)    :: num_soilp       ! number of soil patches in filter
     integer                , intent(in)    :: filter_soilp(:) ! filter for soil patches
-    integer                , intent(in)    :: num_pcropp      ! number of prognostic crops in filter
-    integer                , intent(in)    :: filter_pcropp(:)! filter for prognostic crop patches
     type(temperature_type) , intent(inout) :: temperature_inst
     type(cnveg_state_type) , intent(inout) :: cnveg_state_inst
     type(crop_type)        , intent(inout) :: crop_inst
@@ -554,24 +552,12 @@ subroutine CNPhenologyClimate (num_soilp, filter_soilp, num_pcropp, filter_pcrop
     integer  :: p       ! indices
     integer  :: fp      ! filter patch index
     real(r8) :: dayspyr ! days per year (days)
-    integer  :: kyr     ! current year
-    integer  :: kmo     ! month of year  (1, ..., 12)
-    integer  :: kda     ! day of month   (1, ..., 31)
-    integer  :: mcsec   ! seconds of day (0, ..., seconds/day)
-    real(r8), parameter :: yravg   = 20.0_r8      ! length of years to average for gdd
-    real(r8), parameter :: yravgm1 = yravg-1.0_r8 ! minus 1 of above
     !-----------------------------------------------------------------------
 
     associate(                                                & 
          nyrs_crop_active => crop_inst%nyrs_crop_active_patch,   & ! InOut:  [integer (:)  ]  number of years this crop patch has been active
 
          t_ref2m        => temperature_inst%t_ref2m_patch ,   & ! Input:  [real(r8) (:) ]  2m air temperature (K)
-         gdd0           => temperature_inst%gdd0_patch    ,   & ! Output: [real(r8) (:) ]  growing deg. days base 0 deg C (ddays)            
-         gdd8           => temperature_inst%gdd8_patch    ,   & ! Output: [real(r8) (:) ]     "     "    "    "   8  "  "    "               
-         gdd10          => temperature_inst%gdd10_patch   ,   & ! Output: [real(r8) (:) ]     "     "    "    "  10  "  "    "               
-         gdd020         => temperature_inst%gdd020_patch  ,   & ! Output: [real(r8) (:) ]  20-yr mean of gdd0 (ddays)                        
-         gdd820         => temperature_inst%gdd820_patch  ,   & ! Output: [real(r8) (:) ]  20-yr mean of gdd8 (ddays)                        
-         gdd1020        => temperature_inst%gdd1020_patch ,   & ! Output: [real(r8) (:) ]  20-yr mean of gdd10 (ddays)                       
 
          tempavg_t2m    => cnveg_state_inst%tempavg_t2m_patch & ! Output: [real(r8) (:) ]  temp. avg 2m air temperature (K)                  
          )
@@ -585,41 +571,6 @@ subroutine CNPhenologyClimate (num_soilp, filter_soilp, num_pcropp, filter_pcrop
          tempavg_t2m(p) = tempavg_t2m(p) + t_ref2m(p) * (fracday/dayspyr)
       end do
 
-      !
-      ! The following crop related steps are done here rather than CropPhenology
-      ! so that they will be completed each time-step rather than with doalb.
-      !
-      ! NOTE(wjs, 2022-02-03) The above comment about doalb no longer applies, because
-      ! there is no longer a doalb conditional around the CropPhenology call. Therefore,
-      ! we could move these calculations into CropPhenology if it made sense to do so.
-      !
-      ! The following lines come from ibis's climate.f + stats.f
-      ! gdd SUMMATIONS ARE RELATIVE TO THE PLANTING DATE (see subr. updateAccFlds)
-
-      if (num_pcropp > 0) then
-         ! get time-related info
-         call get_curr_date(kyr, kmo, kda, mcsec)
-      end if
-
-      do fp = 1,num_pcropp
-         p = filter_pcropp(fp)
-         if (kmo == 1 .and. kda == 1 .and. nyrs_crop_active(p) == 0) then ! YR 1:
-            gdd020(p)  = 0._r8                             ! set gdd..20 variables to 0
-            gdd820(p)  = 0._r8                             ! and crops will not be planted
-            gdd1020(p) = 0._r8
-         end if
-         if (kmo == 1 .and. kda == 1 .and. mcsec == 0) then        ! <-- END of EVERY YR:
-            if (nyrs_crop_active(p) == 1) then                     ! <-- END of YR 1
-               gdd020(p)  = gdd0(p)                                ! <-- END of YR 1
-               gdd820(p)  = gdd8(p)                                ! <-- END of YR 1
-               gdd1020(p) = gdd10(p)                               ! <-- END of YR 1
-            end if                                                 ! <-- END of YR 1
-            gdd020(p)  = (yravgm1* gdd020(p)  + gdd0(p))  / yravg  ! gdd..20 must be long term avgs
-            gdd820(p)  = (yravgm1* gdd820(p)  + gdd8(p))  / yravg  ! so ignore results for yrs 1 & 2
-            gdd1020(p) = (yravgm1* gdd1020(p) + gdd10(p)) / yravg 
-         end if
-      end do
-
     end associate
 
   end subroutine CNPhenologyClimate

src/biogeophys/TemperatureType.F90

@@ -1215,6 +1215,17 @@ subroutine InitAccBuffer (this, bounds)
        call init_accum_field (name='GDD10', units='K', &
             desc='growing degree-days base 10C from planting', accum_type='runaccum', accum_period=not_used,  &
             subgrid_type='pft', numlev=1, init_value=0._r8)
+
+       ! 20-year running means (20*365 days)
+       call init_accum_field (name='GDD020', units='K', &
+            desc='20-year running mean of growing degree days base 0C from planting', accum_type='runmean', accum_period=-20*365, &
+            subgrid_type='pft', numlev=1, init_value=0._r8)
+       call init_accum_field (name='GDD820', units='K', &
+            desc='20-year running mean of growing degree days base 8C from planting', accum_type='runmean', accum_period=-20*365, &
+            subgrid_type='pft', numlev=1, init_value=0._r8)
+       call init_accum_field (name='GDD1020', units='K', &
+            desc='20-year running mean of growing degree days base 10C from planting', accum_type='runmean', accum_period=-20*365, &
+            subgrid_type='pft', numlev=1, init_value=0._r8)
 
     end if
 
@@ -1330,6 +1341,15 @@ subroutine InitAccVars(this, bounds)
        call extract_accum_field ('GDD10', rbufslp, nstep)
        this%gdd10_patch(begp:endp) = rbufslp(begp:endp)
 
+       call extract_accum_field ('GDD020', rbufslp, nstep)
+       this%gdd020_patch(begp:endp) = rbufslp(begp:endp)
+
+       call extract_accum_field ('GDD820', rbufslp, nstep)
+       this%gdd820_patch(begp:endp) = rbufslp(begp:endp)
+
+       call extract_accum_field ('GDD1020', rbufslp, nstep)
+       this%gdd1020_patch(begp:endp) = rbufslp(begp:endp)
+
     end if
 
     deallocate(rbufslp)
@@ -1342,7 +1362,7 @@ subroutine UpdateAccVars (this, bounds)
     !
     ! USES
     use shr_const_mod    , only : SHR_CONST_CDAY, SHR_CONST_TKFRZ
-    use clm_time_manager , only : get_step_size, get_nstep, is_end_curr_day, get_curr_date
+    use clm_time_manager , only : get_step_size, get_nstep, is_end_curr_day, get_curr_date, is_end_curr_year
     use accumulMod       , only : update_accum_field, extract_accum_field, accumResetVal
     use CNSharedParamsMod, only : upper_soil_layer
     !
@@ -1576,11 +1596,42 @@ subroutine UpdateAccVars (this, bounds)
        call update_accum_field  ('GDD10', rbufslp, nstep)
        call extract_accum_field ('GDD10', this%gdd10_patch, nstep)
 
+       ! Accumulate and extract running 20-year means
+       if (is_end_curr_year()) then
+          call update_accum_field  ('GDD020', this%gdd0_patch, nstep)
+          call extract_accum_field ('GDD020', this%gdd020_patch, nstep)
+          call update_accum_field  ('GDD820', this%gdd8_patch, nstep)
+          call extract_accum_field ('GDD820', this%gdd820_patch, nstep)
+          call update_accum_field  ('GDD1020', this%gdd10_patch, nstep)
+          call extract_accum_field ('GDD1020', this%gdd1020_patch, nstep)
+       end if
+
     end if
 
     deallocate(rbufslp)
     deallocate(rbufslc)
 
   end subroutine UpdateAccVars
 
+  subroutine Clean(this)
+     !
+     ! !DESCRIPTION:
+     ! Finalize this instance
+     !
+     ! !USES:
+     !
+     ! !ARGUMENTS:
+     class(temperature_type), intent(inout) :: this
+     !
+     ! !LOCAL VARIABLES:
+
+     character(len=*), parameter :: subname = 'Clean'
+     !-----------------------------------------------------------------------
+
+     deallocate(this%gdd020_patch)
+     deallocate(this%gdd820_patch)
+     deallocate(this%gdd1020_patch)
+
+   end subroutine Clean
+
 end module TemperatureType