Add soil tillage for crops

bld/CLMBuildNamelist.pm

@@ -1633,6 +1633,11 @@ sub process_namelist_inline_logic {
   ###############################
   setup_logic_crop_inparm($opts,  $nl_flags, $definition, $defaults, $nl);
 
+  ###############################
+  # namelist group: tillage     #
+  ###############################
+  setup_logic_tillage($opts,  $nl_flags, $definition, $defaults, $nl);
+
   ###############################
   # namelist group: ch4par_in   #
   ###############################
@@ -2241,6 +2246,17 @@ sub setup_logic_crop_inparm {
   }
 }
 
+#-------------------------------------------------------------------------------
+
+sub setup_logic_tillage {
+  my ($opts, $nl_flags, $definition, $defaults, $nl) = @_;
+
+  my $tillage_mode = remove_leading_and_trailing_quotes( $nl->get_value( "tillage_mode" ) );
+  if ( $tillage_mode ne "off" && $tillage_mode ne "" && not &value_is_true($nl->get_value('use_crop')) ) {
+      $log->fatal_error( "Tillage only works on crop columns, so use_crop must be true if tillage is enabled." );
+  }
+}
+
 #-------------------------------------------------------------------------------
 sub error_if_set {
    # do a fatal_error and exit if any of the input variable names are set
@@ -4538,7 +4554,7 @@ sub write_output_files {
                soil_resis_inparm  bgc_shared canopyfluxes_inparm aerosol
                clmu_inparm clm_soilstate_inparm clm_nitrogen clm_snowhydrology_inparm
                cnprecision_inparm clm_glacier_behavior crop_inparm irrigation_inparm
-               surfacealbedo_inparm water_tracers_inparm);
+               surfacealbedo_inparm water_tracers_inparm tillage_inparm);
 
   #@groups = qw(clm_inparm clm_canopyhydrology_inparm clm_soilhydrology_inparm
   #             finidat_consistency_checks dynpft_consistency_checks);

bld/namelist_files/namelist_defaults_ctsm.xml

@@ -484,9 +484,9 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
 <!-- The default filenames are given relative to the root directory
      for the CLM2 data in the CESM distribution -->
 <!-- Plant function types (relative to {csmdata}) -->
-<paramfile phys="clm5_1" >lnd/clm2/paramdata/ctsm51_params.c231117.nc</paramfile>
-<paramfile phys="clm5_0" >lnd/clm2/paramdata/clm50_params.c231117.nc</paramfile>
-<paramfile phys="clm4_5" >lnd/clm2/paramdata/clm45_params.c231117.nc</paramfile>
+<paramfile phys="clm5_1" >lnd/clm2/paramdata/ctsm51_params.c240105.nc</paramfile>
+<paramfile phys="clm5_0" >lnd/clm2/paramdata/clm50_params.c240105.nc</paramfile>
+<paramfile phys="clm4_5" >lnd/clm2/paramdata/clm45_params.c240105.nc</paramfile>
 
 <!-- ================================================================== -->
 <!-- FATES default parameter file                                       -->
@@ -2796,4 +2796,12 @@ use_crop=".true.">lnd/clm2/surfdata_map/ctsm5.1.dev052/landuse.timeseries_mpasa1
 <stream_meshfile_exice    use_excess_ice=".true.">lnd/clm2/paramdata/exice_init_0.125x0.125_ESMFmesh_cdf5_c20220802.nc</stream_meshfile_exice>
 <stream_mapalgo_exice     use_excess_ice=".true.">bilinear</stream_mapalgo_exice>
 
+<!-- =========================================  -->
+<!-- Defaults for tillage                       -->
+<!-- =========================================  -->
+
+<tillage_mode>off</tillage_mode>
+<use_original_tillage_phases>.false.</use_original_tillage_phases>
+<max_tillage_depth>0.26d00</max_tillage_depth>
+
 </namelist_defaults>

bld/namelist_files/namelist_definition_ctsm.xml

@@ -2967,5 +2967,23 @@ Mapping method from excess ice input stream data to the model resolution
     none     = no interpolation
 </entry>
 
+<!-- ========================================================================================  -->
+<!-- Namelist items controlling tillage                                                        -->
+<!-- ========================================================================================  -->
+
+<entry id="tillage_mode" type="char*8" category="physics"
+       group="tillage_inparm" valid_values="off,low,high">
+Whether to till crop soil, and if so, with what intensity.
+</entry>
+
+<entry id="use_original_tillage_phases" type="logical" category="physics"
+       group="tillage_inparm" valid_values="" value=".false.">
+Toggle to use original (Graham et al. 2021) tillage logic, with bug for seasons crossing into a new calendar year
+</entry>
+
+<entry id="max_tillage_depth" type="real" category="physics"
+       group="tillage_inparm" valid_values="" value="0.26d00">
+Maximum depth to till soil (m). Default 0.26; original (Graham et al., 2021) value was unintentionally 0.32.
+</entry>
 
 </namelist_definition>

cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/user_nl_clm

@@ -1,2 +1,2 @@
-paramfile = '$DIN_LOC_ROOT/lnd/clm2/paramdata/ctsm51_ciso_cwd_hr_params.c231117.nc'
+paramfile = '$DIN_LOC_ROOT/lnd/clm2/paramdata/ctsm51_ciso_cwd_hr_params.c240105.nc'
 hist_fincl1 = 'CWDC_HR','C13_CWDC_HR','C14_CWDC_HR','CWD_HR_L2','CWD_HR_L2_vr','CWD_HR_L3','CWD_HR_L3_vr'

cime_config/testdefs/testmods_dirs/clm/till/include_user_mods

@@ -0,0 +1 @@
+../crop

cime_config/testdefs/testmods_dirs/clm/till/user_nl_clm

@@ -0,0 +1,2 @@
+tillage_mode = 'high'
+

doc/source/tech_note/Crop_Irrigation/CLM50_Tech_Note_Crop_Irrigation.rst

@@ -39,7 +39,8 @@ These updates appear in detail in the sections below. Many also appear in :ref:`
 Available new features since the CLM5 release
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 - Addition of bioenergy crops
-- Ablity to customize crop calendars (sowing windows/dates, maturity requirements) using stream files
+- Ability to customize crop calendars (sowing windows/dates, maturity requirements) using stream files
+- Cropland soil tillage
 
 .. _The crop model:
 
@@ -708,6 +709,12 @@ Separate reproductive pool
 ''''''''''''''''''''''''''
 One notable difference between natural vegetation and crops is the presence of reproductive carbon and nitrogen pools. Accounting for the reproductive pools helps determine whether crops are performing reasonably through yield calculations. The reproductive pool is maintained similarly to the leaf, stem, and fine root pools, but allocation of carbon and nitrogen does not begin until the grain fill stage of crop development. Equation :eq:`25.5` describes the carbon and nitrogen allocation coefficients to the reproductive pool. In CLM5BGCCROP, as allocation declines in stem, leaf, and root pools (see section :numref:`Grain fill to harvest`) during the grain fill stage of growth, increasing amounts of carbon and nitrogen are available for grain development.
 
+.. _Tillage:
+
+Tillage
+'''''''
+Tillage is represented as an enhancement of the decomposition rate coefficient; see section :numref:`decomp_mgmt_modifiers`.
+
 .. _The irrigation model:
 
 The irrigation model

doc/source/tech_note/Decomposition/CLM50_Tech_Note_Decomposition.rst

@@ -209,6 +209,35 @@ The combined decomposition rate scalar (:math:`{r}_{total}`,unitless) is:
 
    r_{total} =r_{tsoil} r_{water} r_{oxygen} r_{depth} .
 
+.. _decomp_mgmt_modifiers:
+
+Management modifiers on decomposition rate
+--------------------------------------------------
+
+Tillage of cropland soil is represented as an additional rate scalar that depends on tillage intensity (default off), soil pool, and time since planting :ref:`(Graham et al., 2021) <Grahametal2021>`. The tillage enhancement is strongest in the first 14 days after planting (idpp < 15), weaker in the next 30 days (15 ≤ idpp < 45), weaker still in the next 30 days (45 ≤ idpp < 75), and nonexistent after that (idpp ≥ 75).
+
+.. list-table:: Tillage decomposition rate scalars. Values in each cell represent enhancement in different periods of days past planting: [0, 14], [15, 44], [45, 74].
+   :header-rows: 1
+
+   * - \
+     - low
+     - high
+   * - Litter 2 (cel_lit)
+     - 1.5, 1.5, 1.1
+     - 1.8, 1.5, 1.1
+   * - Litter 3 (lig_lit)
+     - 1.5, 1.5, 1.1
+     - 1.8, 1.5, 1.1
+   * - SOM 1 (act_som)
+     - 1.0, 1.0, 1.0
+     - 1.2, 1.0, 1.0
+   * - SOM 2 (slo_som)
+     - 3.0, 1.6, 1.3
+     - 4.8, 3.5, 2.5
+   * - SOM 3 (pas_som)
+     - 3.0, 1.6, 1.3
+     - 4.8, 3.5, 2.5
+
 N-limitation of Decomposition Fluxes
 -----------------------------------------
 

doc/source/tech_note/References/CLM50_Tech_Note_References.rst

@@ -463,6 +463,10 @@ Gotangco Castillo C., Levis S., and Thornton P. 2012. Evaluation of the new CNDV
 
 Graham, S.T., Famiglietti, J.S., and Maidment, D.R. 1999. Five-minute, 1/2°, and 1° data sets of continental watersheds and river networks for use in regional and global hydrologic and climate system modeling studies. Water Resour. Res. 35:583-587.
 
+.. _Grahametal2021:
+
+Graham, M. W., Thomas, R. Q., Lombardozzi, D. L., & O'Rourke, M. E. (2021). Modest capacity of no-till farming to offset emissions over 21st century. Environmental Research Letters, 16(5), 054055. doi: 10.1088/1748-9326/abe6c6
+
 .. _Gravenetal2017:
 
 Graven, H., C. E. Allison, D. M. Etheridge, S. Hammer, R. F. Keeling, I. Levin, H. A. J. Meijer, M. Rubino, P. P. Tans, C. M. Trudinger, B. H. Vaughn and J. W. C. White, 2017. Compiled records of carbon isotopes in atmospheric CO2 for historical simulations in CMIP6, Geoscientific Model Development, in review. doi: 10.5194/gmd-2017-166.

doc/source/users_guide/running-special-cases/Running-with-tillage.rst

@@ -0,0 +1,38 @@
+.. _running-with-tillage:
+
+.. include:: ../substitutions.rst
+
+=====================
+ Running with tillage
+=====================
+
+
+Cropland tillage (Sect. :numref:`decomp_mgmt_modifiers`) can be toggled by specifying a value of ``'low'`` (low intensity) or ``'high'`` (high intensity) for the  ``tillage_mode`` namelist option. By default this option is ``'off'``.
+
+Depth of tillage can be changed with the ``max_tillage_depth`` parameter (meters; default 0.26).
+
+Tillage multipliers for different soil pools and time since planting are defined on the parameter file, in variables ``bgc_till_decompk_multipliers`` (for CENTURY soil) and ``mimics_till_decompk_multipliers`` (for MIMICS soil). These variables were originally added with the script at ``tools/contrib/add_tillage_to_paramsfile.py``, which can be modified as needed to change tillage multipliers.
+
+
+Example: Crop simulation with tillage
+-------------------------------------
+::
+
+   > cime/scripts/create_newcase -case IHistClm51BgcCrop_till -res f19_g17_gl4 -compset IHistClm51BgcCrop
+
+
+   > cd IHistClm51BgcCrop_till
+   > ./case.setup
+
+   # turn on tillage ('low' or 'high'; default 'off')
+   > echo "tillage_mode = 'high'" >> user_nl_clm
+
+Reverting fixes relative to original tillage implementation
+-----------------------------------------------------------
+
+The current implementation of tillage in CTSM is based on work by :ref:`Graham et al. (2021) <Grahametal2021>`, but with fixes to how days after planting is calculated and to default tillage depth. To run without those changes:
+
+::
+
+   > echo "use_original_tillage_phases = .true." >> user_nl_clm
+   > echo "max_tillage_depth = 0.32" >> user_nl_clm

src/biogeochem/CNDriverMod.F90

@@ -343,12 +343,14 @@ subroutine CNDriverNoLeaching(bounds,
     call t_startf('DecompRate')
     if (decomp_method == century_decomp) then
        call decomp_rate_constants_bgc(bounds, num_bgc_soilc, filter_bgc_soilc, &
-            soilstate_inst, temperature_inst, ch4_inst, soilbiogeochem_carbonflux_inst)
+            soilstate_inst, temperature_inst, ch4_inst, soilbiogeochem_carbonflux_inst, &
+            cnveg_state_inst%idop_patch)
     else if (decomp_method == mimics_decomp) then
        call decomp_rates_mimics(bounds, num_bgc_soilc, filter_bgc_soilc, &
             num_bgc_vegp, filter_bgc_vegp, clm_fates, &
             soilstate_inst, temperature_inst, cnveg_carbonflux_inst, ch4_inst, &
-            soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst)
+            soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, &
+            cnveg_state_inst%idop_patch)
     end if
     call t_stopf('DecompRate')
 

src/biogeochem/CNPhenologyMod.F90

@@ -52,6 +52,7 @@ module CNPhenologyMod
   public :: CNPhenologyInit      ! Initialization
   public :: CNPhenology          ! Update
   public :: CropPhase            ! Get the current phase of each crop patch
+  public :: DaysPastPlanting     ! Get how many days it's been since crop was planted
 
   ! !PUBLIC for unit testing
   public :: CNPhenologySetNML         ! Set the namelist setttings explicitly for unit tests
@@ -2255,12 +2256,7 @@ subroutine CropPhenology(num_pcropp, filter_pcropp                     , &
             end if
 
             ! days past planting may determine harvest
-
-            if (jday >= idop(p)) then
-               idpp = jday - idop(p)
-            else
-               idpp = int(dayspyr) + jday - idop(p)
-            end if
+            idpp = DaysPastPlanting(idop(p), jday)
 
             ! onset_counter initialized to zero when .not. croplive
             ! offset_counter relevant only at time step of harvest
@@ -2772,6 +2768,38 @@ subroutine PlantCrop(p, leafcn_in, jday, kyr, do_plant_normal, &
 
   end subroutine PlantCrop
 
+  !-----------------------------------------------------------------------
+  function DaysPastPlanting(idop, jday_in)
+    ! !USES:
+    use clm_time_manager, only : get_prev_calday, get_curr_days_per_year
+    !
+    ! !ARGUMENTS:
+    integer,           intent(in) :: idop ! patch day of planting
+    integer, optional, intent(in) :: jday_in ! julian day of the year
+    !
+    ! !LOCAL VARIABLES
+    integer :: DaysPastPlanting
+    integer :: jday
+
+    ! Must use separate jday_in and jday because we can't redefine an intent(in)
+    ! variable, even if it wasn't provided in the function call.
+    if (present(jday_in)) then
+       jday = jday_in
+    else
+       ! Use prev instead of curr to avoid jday=1 in last timestep of year
+       jday = get_prev_calday()
+    end if
+
+    if (jday >= idop) then
+       DaysPastPlanting = jday - idop
+    else
+      ! As long as crops have at most a 365-day growing season, using get_curr_days_per_year()
+      ! should give the same result of this function as using get_prev_days_per_year().
+       DaysPastPlanting = jday - idop + get_curr_days_per_year()
+    end if
+
+  end function DaysPastPlanting
+
   !-----------------------------------------------------------------------
   subroutine vernalization(p, &
        canopystate_inst, temperature_inst, waterdiagnosticbulk_inst, cnveg_state_inst, crop_inst, &

src/main/readParamsMod.F90

@@ -39,6 +39,7 @@ subroutine readParameters (photosyns_inst)
     use SoilBiogeochemLittVertTranspMod   , only : readSoilBiogeochemLittVertTranspParams => readParams
     use SoilBiogeochemPotentialMod        , only : readSoilBiogeochemPotentialParams      => readParams
     use SoilBiogeochemDecompMod           , only : readSoilBiogeochemDecompParams         => readParams
+    use TillageMod                        , only : readTillageParams                      => readParams
     use SoilBiogeochemDecompCascadeMIMICSMod, only : readSoilBiogeochemDecompMimicsParams => readParams
     use SoilBiogeochemDecompCascadeBGCMod , only : readSoilBiogeochemDecompBgcParams      => readParams
     use ch4Mod                            , only : readCH4Params                          => readParams
@@ -105,6 +106,7 @@ subroutine readParameters (photosyns_inst)
           call readSoilBiogeochemDecompBgcParams(ncid)
        end if
        call readSoilBiogeochemDecompParams(ncid)
+       call readTillageParams(ncid, NLFilename_in)
        call readSoilBiogeochemLittVertTranspParams(ncid)
        call readSoilBiogeochemNitrifDenitrifParams(ncid)
        call readSoilBiogeochemNLeachingParams(ncid)

src/soilbiogeochem/CMakeLists.txt

@@ -7,6 +7,7 @@ list(APPEND clm_sources
   SoilBiogeochemStateType.F90
   SoilBiogeochemNitrogenStateType.F90
   SoilBiogeochemNitrogenFluxType.F90
+  TillageMod.F90
   )
 
 sourcelist_to_parent(clm_sources)

src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90

@@ -229,8 +229,6 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i
     !  initialize rate constants and decomposition pathways following the decomposition cascade of the BGC model.
     !  written by C. Koven 
     !
-    ! !USES:
-    !
     ! !ARGUMENTS:
     type(bounds_type)               , intent(in)    :: bounds  
     type(soilbiogeochem_state_type) , intent(inout) :: soilbiogeochem_state_inst
@@ -511,7 +509,8 @@ end subroutine init_decompcascade_bgc
 
   !-----------------------------------------------------------------------
   subroutine decomp_rate_constants_bgc(bounds, num_bgc_soilc, filter_bgc_soilc, &
-       soilstate_inst, temperature_inst, ch4_inst, soilbiogeochem_carbonflux_inst)
+       soilstate_inst, temperature_inst, ch4_inst, soilbiogeochem_carbonflux_inst, &
+       idop)
     !
     ! !DESCRIPTION:
     !  calculate rate constants and decomposition pathways for the CENTURY decomposition cascade model
@@ -521,6 +520,9 @@ subroutine decomp_rate_constants_bgc(bounds, num_bgc_soilc, filter_bgc_soilc, &
     use clm_time_manager , only : get_average_days_per_year, get_step_size
     use shr_const_mod    , only : SHR_CONST_PI
     use clm_varcon       , only : secspday
+    use TillageMod       , only : get_do_tillage
+    use TillageMod       , only : get_apply_tillage_multipliers
+    use landunit_varcon  , only : istcrop
     !
     ! !ARGUMENTS:
     type(bounds_type)                    , intent(in)    :: bounds          
@@ -530,6 +532,7 @@ subroutine decomp_rate_constants_bgc(bounds, num_bgc_soilc, filter_bgc_soilc, &
     type(temperature_type)               , intent(in)    :: temperature_inst
     type(ch4_type)                       , intent(in)    :: ch4_inst
     type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst
+    integer, optional                    , intent(in)    :: idop(:) ! patch day of planting
     !
     ! !LOCAL VARIABLES:
     real(r8), parameter :: eps = 1.e-6_r8
@@ -595,6 +598,10 @@ subroutine decomp_rate_constants_bgc(bounds, num_bgc_soilc, filter_bgc_soilc, &
               errMsg(sourcefile, __LINE__))
       endif
 
+      if (get_do_tillage() .and. .not. present(idop)) then
+         call endrun("Do not enable tillage without providing idop to decomp_rate_constants_bgc().")
+      end if
+
       days_per_year = get_average_days_per_year()
       dt = real( get_step_size(), r8 )
 
@@ -895,6 +902,12 @@ subroutine decomp_rate_constants_bgc(bounds, num_bgc_soilc, filter_bgc_soilc, &
                decomp_k(c,j,i_cwd) = k_frag * t_scalar(c,j) * w_scalar(c,j) * &
                   depth_scalar(c,j) * o_scalar(c,j) * spinup_geogterm_cwd(c)
             end if
+
+            ! Tillage
+            if (get_do_tillage()) then
+               call get_apply_tillage_multipliers(idop, c, j, decomp_k(c,j,:))
+            end if
+
             ! Above into soil matrix
             if(use_soil_matrixcn)then
                ! same for cwd but only if fates is not enabled; fates handles CWD
@@ -904,6 +917,7 @@ subroutine decomp_rate_constants_bgc(bounds, num_bgc_soilc, filter_bgc_soilc, &
             end if !use_soil_matrixcn
          end do
       end do
+
       pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l1s1) = 1.0_r8
       pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l2s1) = 1.0_r8
       pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l3s2) = 1.0_r8