Add soil tillage for crops

bld/CLMBuildNamelist.pm

@@ -1629,6 +1629,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   #
   ###############################
@@ -2183,6 +2188,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( "It doesn't make sense to use tillage with use_crop false" );
+  }
+}
+
 #-------------------------------------------------------------------------------
 sub error_if_set {
    # do a fatal_error and exit if any of the input variable names are set
@@ -4412,7 +4428,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

@@ -2723,4 +2723,11 @@ 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>.false.</use_original_tillage>
+
 </namelist_defaults>

bld/namelist_files/namelist_definition_ctsm.xml

@@ -2881,5 +2881,18 @@ 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" type="logical" category="physics"
+       group="tillage_inparm" valid_values="" value=".false.">
+Toggle to use original tillage logic
+</entry>
 
 </namelist_definition>

src/biogeochem/CNDriverMod.F90

@@ -340,12 +340,14 @@ subroutine CNDriverNoLeaching(bounds,
     call t_startf('DecompRate')
     if (decomp_method == century_decomp) then
        call decomp_rate_constants_bgc(bounds, num_soilc, filter_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_soilc, filter_soilc, &
             num_soilp, filter_soilp, 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/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)
        call readSoilBiogeochemLittVertTranspParams(ncid)
        call readSoilBiogeochemNitrifDenitrifParams(ncid)
        call readSoilBiogeochemNLeachingParams(ncid)

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_soilc, filter_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_soilc, filter_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_soilc, filter_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_soilc, filter_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_soilc, filter_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_soilc, filter_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

src/soilbiogeochem/SoilBiogeochemDecompCascadeMIMICSMod.F90

@@ -755,7 +755,8 @@ end subroutine init_decompcascade_mimics
   subroutine decomp_rates_mimics(bounds, num_soilc, filter_soilc, &
        num_soilp, filter_soilp, clm_fates, &
        soilstate_inst, temperature_inst, cnveg_carbonflux_inst, &
-       ch4_inst, soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst)
+       ch4_inst, soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, &
+       idop)
     !
     ! !DESCRIPTION:
     ! Calculate rates and decomposition pathways for the MIMICS
@@ -768,6 +769,8 @@ subroutine decomp_rates_mimics(bounds, num_soilc, filter_soilc, &
     use subgridAveMod    , only : p2c
     use PatchType        , only : patch
     use pftconMod        , only : pftname
+    use TillageMod       , only : get_do_tillage
+    use TillageMod       , only : get_apply_tillage_multipliers
     !
     ! !ARGUMENTS:
     type(bounds_type)                    , intent(in)    :: bounds          
@@ -782,6 +785,7 @@ subroutine decomp_rates_mimics(bounds, num_soilc, filter_soilc, &
     type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst
     type(soilbiogeochem_carbonstate_type), intent(in)    :: soilbiogeochem_carbonstate_inst
     type(hlm_fates_interface_type)       , intent(inout) :: clm_fates
+    integer, optional                    , intent(in)    :: idop(:) ! patch day of planting
     !
     ! !LOCAL VARIABLES:
     real(r8), parameter :: eps = 1.e-6_r8
@@ -883,6 +887,10 @@ subroutine decomp_rates_mimics(bounds, num_soilc, filter_soilc, &
          spinup_factor  => decomp_cascade_con%spinup_factor              & ! Input:  [real(r8)          (:)     ]  factor for AD spinup associated with each pool           
          )
 
+      if (get_do_tillage() .and. .not. present(idop)) then
+         call endrun("Do not enable tillage without providing idop to decomp_rate_constants_mimics().")
+      end if
+
       mino2lim = CNParamsShareInst%mino2lim
 
       days_per_year = get_average_days_per_year()
@@ -1301,6 +1309,11 @@ subroutine decomp_rates_mimics(bounds, num_soilc, filter_soilc, &
             if (.not. use_fates) then
                decomp_k(c,j,i_cwd) = k_frag * w_d_o_scalars  ! * 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
          end do
       end do