Merge branch 'fix_parameters_gregorian'

Use avg days per year when converting param units

When converting parameter units from per-year to per-second, use average
days per year instead of current number of days per year. This is
relevant when running with a Gregorian calendar.

See #1612 for details.

Resolves #1612 (Some uses of get_days_per_year should use the
average number of days in a year, not the number of days in the current
year)

doc/ChangeLog

@@ -1,4 +1,93 @@
 ===============================================================
+Tag name: ctsm5.1.dev080
+Originator(s): sacks (Bill Sacks)
+Date: Thu Feb 24 15:26:02 MST 2022
+One-line Summary: Use avg days per year when converting param units
+
+Purpose and description of changes
+----------------------------------
+
+When converting parameter units from per-year to per-second, use average
+days per year instead of current number of days per year. This is
+relevant when running with a Gregorian calendar.
+
+See https://github.com/ESCOMP/CTSM/issues/1612 for details.
+
+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.]
+
+[ ] clm5_1
+
+[ ] clm5_0
+
+[ ] ctsm5_0-nwp
+
+[ ] clm4_5
+
+
+Bugs fixed or introduced
+------------------------
+Issues fixed (include CTSM Issue #):
+- Resolves ESCOMP/CTSM#1612 (Some uses of get_days_per_year should use
+  the average number of days in a year, not the number of days in the
+  current year)
+
+Known bugs introduced in this tag (include issue #):
+- ESCOMP/CTSM#1624 (Change get_average_days_per_year to use
+  ESMF_CalendarGet)
+
+Notes of particular relevance for developers:
+---------------------------------------------
+Caveats for developers (e.g., code that is duplicated that requires double maintenance):
+- Once ESMF supports it, we should change get_average_days_per_year to
+  use ESMF_CalendarGet (see ESCOMP/CTSM#1624)
+
+Testing summary:
+----------------
+
+    cheyenne ---- OK
+    izumi ------- OK
+
+Test
+ERS_Ly3_P72x2_Vmct.f10_f10_mg37.IHistClm50BgcCropG.cheyenne_intel.clm-cropMonthOutput
+initially failed COMPARE_base_rest and BASELINE comparisons; rerunning
+solved the issue.
+
+Answer changes
+--------------
+
+Changes answers relative to baseline: YES
+
+  Summarize any changes to answers, i.e.,
+    - what code configurations: Gregorian cases with BGC
+    - what platforms/compilers: all
+    - nature of change (roundoff; larger than roundoff/same climate; new climate):
+      larger than roundoff / same climate
+
+      Changes a few BGC-related parameters by a small amount (< 1/365 in
+      a relative sense) for Gregorian cases
+
+      Changes answers for these tests:
+      - SMS_Ly5_Mmpi-serial.1x1_smallvilleIA.IHistClm50BgcCropQianRs.izumi_gnu.clm-gregorian_cropMonthOutput
+      - DAE_C2_D_Lh12.f10_f10_mg37.I2000Clm50BgcCrop.cheyenne_intel.clm-DA_multidrv
+      - DAE_N2_D_Lh12_Vmct.f10_f10_mg37.I2000Clm50BgcCrop.cheyenne_intel.clm-DA_multidrv      
+
+      In principle, might change answers by roundoff for NOLEAP BGC
+      tests on some machines / compilers, but that wasn't seen in any
+      aux_clm testing.
+
+Other details
+-------------
+Pull Requests that document the changes (include PR ids):
+https://github.com/ESCOMP/CTSM/pull/1625
+
+===============================================================
+===============================================================
 Tag name: ctsm5.1.dev079
 Originator(s): sacks (Bill Sacks)
 Date: Thu Feb 24 14:40:58 MST 2022

doc/ChangeSum

@@ -1,5 +1,6 @@
 Tag                   Who      Date  Summary
 ============================================================================================================================
+    ctsm5.1.dev080    sacks 02/24/2022 Use avg days per year when converting param units
     ctsm5.1.dev079    sacks 02/24/2022 Changes to CropPhenology timing
     ctsm5.1.dev078    sacks 02/24/2022 Rework single-point testing
     ctsm5.1.dev077   rgknox 02/22/2022 Updates to FATES API, including removal of patch dimensions from fates history and using soil instead of ground layers for fates history

src/biogeochem/CNC14DecayMod.F90

@@ -5,7 +5,7 @@ module CNC14DecayMod
   !
   ! !USES:
   use shr_kind_mod                       , only : r8 => shr_kind_r8
-  use clm_time_manager                   , only : get_step_size_real, get_curr_days_per_year
+  use clm_time_manager                   , only : get_step_size_real, get_average_days_per_year
   use clm_varpar                         , only : nlevdecomp, ndecomp_pools
   use clm_varcon                         , only : secspday
   use clm_varctl                         , only : spinup_state
@@ -87,7 +87,7 @@ subroutine C14Decay( bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, &
 
       ! set time steps
       dt = get_step_size_real()
-      days_per_year = get_curr_days_per_year()
+      days_per_year = get_average_days_per_year()
 
       half_life = 5730._r8 * secspday * days_per_year
       decay_const = - log(0.5_r8) / half_life

src/biogeochem/CNGapMortalityMod.F90

@@ -88,7 +88,7 @@ subroutine CNGapMortality (bounds, num_soilc, filter_soilc, num_soilp, filter_so
     ! Gap-phase mortality routine for coupled carbon-nitrogen code (CN)
     !
     ! !USES:
-    use clm_time_manager , only: get_curr_days_per_year
+    use clm_time_manager , only: get_average_days_per_year
     use clm_varpar       , only: nlevdecomp_full
     use clm_varcon       , only: secspday
     use clm_varctl       , only: use_cndv, spinup_state
@@ -179,7 +179,7 @@ subroutine CNGapMortality (bounds, num_soilc, filter_soilc, num_soilp, filter_so
 
          end if
 
-         m  = am/(get_curr_days_per_year() * secspday)
+         m  = am/(get_average_days_per_year() * secspday)
 
          !------------------------------------------------------
          ! patch-level gap mortality carbon fluxes

src/biogeochem/CNPhenologyMod.F90

@@ -605,7 +605,7 @@ subroutine CNEvergreenPhenology (num_soilp, filter_soilp , &
     !
     ! !USES:
     use clm_varcon       , only : secspday
-    use clm_time_manager , only : get_curr_days_per_year
+    use clm_time_manager , only : get_average_days_per_year
     use clm_varctl       , only : CN_evergreen_phenology_opt   
     !
     ! !ARGUMENTS:
@@ -618,7 +618,7 @@ subroutine CNEvergreenPhenology (num_soilp, filter_soilp , &
     type(cnveg_nitrogenflux_type)  , intent(inout) :: cnveg_nitrogenflux_inst   
     !
     ! !LOCAL VARIABLES:
-    real(r8):: dayspyr                    ! Days per year
+    real(r8):: avg_dayspyr                ! Average days per year
     integer :: p                          ! indices
     integer :: fp                         ! lake filter patch index
 
@@ -693,12 +693,12 @@ subroutine CNEvergreenPhenology (num_soilp, filter_soilp , &
          lgsf       => cnveg_state_inst%lgsf_patch    & ! Output: [real(r8) (:) ]  long growing season factor [0-1]                  
          )
 
-      dayspyr   = get_curr_days_per_year()
+      avg_dayspyr = get_average_days_per_year()
 
       do fp = 1,num_soilp
          p = filter_soilp(fp)
          if (evergreen(ivt(p)) == 1._r8) then
-            bglfr(p) = 1._r8/(leaf_long(ivt(p)) * dayspyr * secspday)
+            bglfr(p) = 1._r8/(leaf_long(ivt(p)) * avg_dayspyr * secspday)
             bgtr(p)  = 0._r8
             lgsf(p)  = 0._r8
          end if
@@ -1220,7 +1220,7 @@ subroutine CNStressDecidPhenology (num_soilp, filter_soilp , &
     ! per year.
     !
     ! !USES:
-    use clm_time_manager , only : get_curr_days_per_year
+    use clm_time_manager , only : get_average_days_per_year
     use CNSharedParamsMod, only : use_fun
     use clm_varcon       , only : secspday
     use shr_const_mod    , only : SHR_CONST_TKFRZ, SHR_CONST_PI
@@ -1243,7 +1243,7 @@ subroutine CNStressDecidPhenology (num_soilp, filter_soilp , &
     real(r8),parameter :: secspqtrday = secspday / 4  ! seconds per quarter day
     integer :: g,c,p           ! indices
     integer :: fp              ! lake filter patch index
-    real(r8):: dayspyr         ! days per year
+    real(r8):: avg_dayspyr     ! average days per year
     real(r8):: crit_onset_gdd  ! degree days for onset trigger
     real(r8):: soilt           ! temperature of top soil layer
     real(r8):: psi             ! water stress of top soil layer
@@ -1338,8 +1338,7 @@ subroutine CNStressDecidPhenology (num_soilp, filter_soilp , &
          deadcrootn_storage_to_xfer          =>    cnveg_nitrogenflux_inst%deadcrootn_storage_to_xfer_patch      & ! Output:  [real(r8) (:)   ]                                                    
          )
 
-      ! set time steps
-      dayspyr = get_curr_days_per_year()
+      avg_dayspyr = get_average_days_per_year()
 
       ! specify rain threshold for leaf onset
       rain_threshold = 20._r8
@@ -1588,7 +1587,7 @@ subroutine CNStressDecidPhenology (num_soilp, filter_soilp , &
             ! calculate long growing season factor (lgsf)
             ! only begin to calculate a lgsf greater than 0.0 once the number
             ! of days active exceeds days/year.
-            lgsf(p) = max(min(3.0_r8*(days_active(p)-leaf_long(ivt(p))*dayspyr )/dayspyr, 1._r8),0._r8)
+            lgsf(p) = max(min(3.0_r8*(days_active(p)-leaf_long(ivt(p))*avg_dayspyr )/avg_dayspyr, 1._r8),0._r8)
             ! RosieF. 5 Nov 2015.  Changed this such that the increase in leaf turnover is faster after
             ! trees enter the 'fake evergreen' state. Otherwise, they have a whole year of 
             ! cheating, with less litterfall than they should have, resulting in very high LAI. 
@@ -1603,7 +1602,7 @@ subroutine CNStressDecidPhenology (num_soilp, filter_soilp , &
                ! calculate the background litterfall rate (bglfr)
                ! in units 1/s, based on leaf longevity (yrs) and correction for long growing season
 
-               bglfr(p) = (1._r8/(leaf_long(ivt(p))*dayspyr*secspday))*lgsf(p)
+               bglfr(p) = (1._r8/(leaf_long(ivt(p))*avg_dayspyr*secspday))*lgsf(p)
             end if
 
             ! set background transfer rate when active but not in the phenological onset period
@@ -1614,7 +1613,7 @@ subroutine CNStressDecidPhenology (num_soilp, filter_soilp , &
                ! in complete turnover of the storage pools in one year at steady state,
                ! once lgsf has reached 1.0 (after 730 days active).
 
-               bgtr(p) = (1._r8/(dayspyr*secspday))*lgsf(p)
+               bgtr(p) = (1._r8/(avg_dayspyr*secspday))*lgsf(p)
 
                ! set carbon fluxes for shifting storage pools to transfer pools
 
@@ -1661,6 +1660,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
 
     ! !USES:
     use clm_time_manager , only : get_prev_date, get_prev_calday, get_curr_days_per_year
+    use clm_time_manager , only : get_average_days_per_year
     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
@@ -1699,7 +1699,8 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
     integer g         ! gridcell indices
     integer h         ! hemisphere indices
     integer idpp      ! number of days past planting
-    real(r8) dayspyr  ! days per year
+    real(r8) dayspyr  ! days per year in this year
+    real(r8) avg_dayspyr ! average number of days per year
     real(r8) crmcorn  ! comparitive relative maturity for corn
     real(r8) ndays_on ! number of days to fertilize
     logical do_plant_normal ! are the normal planting rules defined and satisfied?
@@ -1765,6 +1766,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
 
       ! get time info
       dayspyr = get_curr_days_per_year()
+      avg_dayspyr = get_average_days_per_year()
       jday    = get_prev_calday()
       call get_prev_date(kyr, kmo, kda, mcsec)
 
@@ -2119,7 +2121,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
 
             else if (hui(p) >= huigrain(p)) then
                cphase(p) = 3._r8
-               bglfr(p) = 1._r8/(leaf_long(ivt(p))*dayspyr*secspday)
+               bglfr(p) = 1._r8/(leaf_long(ivt(p))*avg_dayspyr*secspday)
             end if
 
             ! continue fertilizer application while in phase 2;

src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90

@@ -516,7 +516,7 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, &
     !  written by C. Koven based on original CLM4 decomposition cascade
     !
     ! !USES:
-    use clm_time_manager , only : get_curr_days_per_year
+    use clm_time_manager , only : get_average_days_per_year
     use shr_const_mod    , only : SHR_CONST_PI
     use clm_varcon       , only : secspday
     !
@@ -592,7 +592,7 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, &
               errMsg(sourcefile, __LINE__))
       endif
 
-      days_per_year = get_curr_days_per_year()
+      days_per_year = get_average_days_per_year()
 
       ! set "Q10" parameter
       Q10 = CNParamsShareInst%Q10

src/soilbiogeochem/SoilBiogeochemDecompCascadeMIMICSMod.F90

@@ -760,7 +760,7 @@ subroutine decomp_rates_mimics(bounds, num_soilc, filter_soilc, &
     ! decomposition cascade model
     !
     ! !USES:
-    use clm_time_manager , only : get_curr_days_per_year
+    use clm_time_manager , only : get_average_days_per_year
     use clm_varcon       , only : secspday, secsphr, tfrz
     use clm_varcon       , only : g_to_mg, cm3_to_m3
     !
@@ -873,7 +873,7 @@ subroutine decomp_rates_mimics(bounds, num_soilc, filter_soilc, &
 
       mino2lim = CNParamsShareInst%mino2lim
 
-      days_per_year = get_curr_days_per_year()
+      days_per_year = get_average_days_per_year()
 
 !     ! Set "decomp_depth_efolding" parameter
 !     decomp_depth_efolding = CNParamsShareInst%decomp_depth_efolding

src/utils/clm_time_manager.F90

@@ -38,6 +38,7 @@ module clm_time_manager
         get_prev_calday,          &! return calendar day at beginning of current timestep
         get_calday,               &! return calendar day from input date
         get_calendar,             &! return calendar
+        get_average_days_per_year,&! return the average number of days per year for the given calendar
         get_curr_days_per_year,   &! return the days per year for year as of the end of the current time step
         get_prev_days_per_year,   &! return the days per year for year as of the beginning of the current time step
         get_curr_yearfrac,        &! return the fractional position in the current year, as of the end of the current timestep
@@ -123,6 +124,9 @@ module clm_time_manager
    private :: TimeGetymd
    private :: check_timemgr_initialized
 
+   character(len=*), parameter, private :: sourcefile = &
+        __FILE__
+
    !=========================================================================================
 contains
   !=========================================================================================
@@ -1258,6 +1262,63 @@ end function get_calendar
 
   !=========================================================================================
 
+  real(r8) function get_average_days_per_year()
+
+    !---------------------------------------------------------------------------------
+    ! Get the average number of days per year for the given calendar.
+    !
+    ! This should be used, for example, when converting a parameter from units of
+    ! per-year to units of per-second (so that the parameter will have a fixed, constant
+    ! value rather than a slightly different value on leap years vs. non-leap years).
+
+    real(r8) :: avg_days_per_year
+    real(r8) :: curr_days_per_year
+
+    real(r8), parameter :: days_per_year_noleap = 365._r8
+
+    ! From the definition of ESMF_CALKIND_GREGORIAN in
+    ! https://earthsystemmodeling.org/docs/release/latest/ESMF_refdoc/node6.html: "In the
+    ! Gregorian calendar every fourth year is a leap year in which February has 29 and not
+    ! 28 days; however, years divisible by 100 are not leap years unless they are also
+    ! divisible by 400." This results in an average number of days per year of 365.2425.
+    real(r8), parameter :: days_per_year_gregorian = 365.2425_r8
+
+    character(len=*), parameter :: subname = 'get_average_days_per_year'
+    !---------------------------------------------------------------------------------
+
+    ! BUG(wjs, 2022-02-01, ESCOMP/CTSM#1624) Ideally we would use ESMF_CalendarGet here,
+    ! but that currently isn't possible (see notes in issue 1624 for details)
+    if (to_upper(calendar) == NO_LEAP_C) then
+       avg_days_per_year = days_per_year_noleap
+    else if (to_upper(calendar) == GREGORIAN_C) then
+       avg_days_per_year = days_per_year_gregorian
+    else
+       call shr_sys_abort(subname//' ERROR: unrecognized calendar specified= '//trim(calendar))
+    end if
+
+    ! Paranoia: Since we're using a hard-coded value, let's make sure that the user hasn't
+    ! done some customizations to the calendar that change the days per year from what we
+    ! expect: Compare the hard-coded value with the number of days per year in the
+    ! current year, which comes from the actual ESMF calendar; the two should be close.
+    ! (This check can be removed once we address issue 1624, making the results of this
+    ! function depend on the actual ESMF calendar instead of a hard-coded value.)
+    curr_days_per_year = get_curr_days_per_year()
+    if (abs(avg_days_per_year - curr_days_per_year) > 1._r8) then
+       write(iulog,*) 'ERROR: hard-coded average days per year differs by more than expected'
+       write(iulog,*) 'from current days per year. Are you using a non-standard calendar?'
+       write(iulog,*) 'avg_days_per_year (hard-coded)          = ', avg_days_per_year
+       write(iulog,*) 'curr_days_per_year (from ESMF calendar) = ', curr_days_per_year
+       write(iulog,*) 'You can fix this by changing the hard-coded parameters in '//subname
+       write(iulog,*) 'in file: '//sourcefile
+       call shr_sys_abort(subname//' ERROR: hard-coded average days per year differs by more than expected')
+    end if
+
+    get_average_days_per_year = avg_days_per_year
+
+  end function get_average_days_per_year
+
+  !=========================================================================================
+
   integer function get_curr_days_per_year( offset )
 
     !---------------------------------------------------------------------------------