diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 17d64cd59f..cdfceb7ccb 100755 --- a/bld/CLMBuildNamelist.pm +++ b/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); diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 3997066a0b..6206391794 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -484,9 +484,9 @@ attributes from the config_cache.xml file (with keys converted to upper-case). -lnd/clm2/paramdata/ctsm51_params.c231117.nc -lnd/clm2/paramdata/clm50_params.c231117.nc -lnd/clm2/paramdata/clm45_params.c231117.nc +lnd/clm2/paramdata/ctsm51_params.c240105.nc +lnd/clm2/paramdata/clm50_params.c240105.nc +lnd/clm2/paramdata/clm45_params.c240105.nc @@ -2796,4 +2796,12 @@ use_crop=".true.">lnd/clm2/surfdata_map/ctsm5.1.dev052/landuse.timeseries_mpasa1 lnd/clm2/paramdata/exice_init_0.125x0.125_ESMFmesh_cdf5_c20220802.nc bilinear + + + + +off +.false. +0.26d00 + diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 97310ebe80..f763dcb968 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/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 + + + + + +Whether to till crop soil, and if so, with what intensity. + + + +Toggle to use original (Graham et al. 2021) tillage logic, with bug for seasons crossing into a new calendar year + + + +Maximum depth to till soil (m). Default 0.26; original (Graham et al., 2021) value was unintentionally 0.32. + diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index 5e15eba958..3b51a596b9 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -1787,7 +1787,7 @@ - + @@ -2296,7 +2296,7 @@ - + @@ -2305,7 +2305,7 @@ - + diff --git a/cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/user_nl_clm index 534cbb84b8..4e073859be 100644 --- a/cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/user_nl_clm +++ b/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' diff --git a/cime_config/testdefs/testmods_dirs/clm/till/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/till/include_user_mods new file mode 100644 index 0000000000..23ea3745e6 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/till/include_user_mods @@ -0,0 +1 @@ +../crop diff --git a/cime_config/testdefs/testmods_dirs/clm/till/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/till/user_nl_clm new file mode 100644 index 0000000000..f94c6df309 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/till/user_nl_clm @@ -0,0 +1,2 @@ +tillage_mode = 'high' + diff --git a/doc/source/tech_note/Crop_Irrigation/CLM50_Tech_Note_Crop_Irrigation.rst b/doc/source/tech_note/Crop_Irrigation/CLM50_Tech_Note_Crop_Irrigation.rst index 2addbd2507..6ac4bc85c6 100755 --- a/doc/source/tech_note/Crop_Irrigation/CLM50_Tech_Note_Crop_Irrigation.rst +++ b/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 diff --git a/doc/source/tech_note/Decomposition/CLM50_Tech_Note_Decomposition.rst b/doc/source/tech_note/Decomposition/CLM50_Tech_Note_Decomposition.rst index c0970d2bd5..bf6d52ee45 100644 --- a/doc/source/tech_note/Decomposition/CLM50_Tech_Note_Decomposition.rst +++ b/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) `. 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 ----------------------------------------- diff --git a/doc/source/tech_note/References/CLM50_Tech_Note_References.rst b/doc/source/tech_note/References/CLM50_Tech_Note_References.rst index a1f84e055b..b824f705bd 100644 --- a/doc/source/tech_note/References/CLM50_Tech_Note_References.rst +++ b/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. diff --git a/doc/source/users_guide/running-special-cases/Running-with-tillage.rst b/doc/source/users_guide/running-special-cases/Running-with-tillage.rst new file mode 100644 index 0000000000..77309aea07 --- /dev/null +++ b/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) `, 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 diff --git a/src/biogeochem/CNDriverMod.F90 b/src/biogeochem/CNDriverMod.F90 index 4ab1e1f51e..bee506c8cb 100644 --- a/src/biogeochem/CNDriverMod.F90 +++ b/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') diff --git a/src/biogeochem/CNPhenologyMod.F90 b/src/biogeochem/CNPhenologyMod.F90 index 42a533d965..ac8e0e7b72 100644 --- a/src/biogeochem/CNPhenologyMod.F90 +++ b/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, & diff --git a/src/main/readParamsMod.F90 b/src/main/readParamsMod.F90 index 2d60b29be5..c11a741198 100644 --- a/src/main/readParamsMod.F90 +++ b/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) diff --git a/src/soilbiogeochem/CMakeLists.txt b/src/soilbiogeochem/CMakeLists.txt index f4545a6a76..e2baa2d1b2 100644 --- a/src/soilbiogeochem/CMakeLists.txt +++ b/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) diff --git a/src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90 b/src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90 index e1b6e68aad..f4785fc4d8 100644 --- a/src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90 +++ b/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 diff --git a/src/soilbiogeochem/SoilBiogeochemDecompCascadeMIMICSMod.F90 b/src/soilbiogeochem/SoilBiogeochemDecompCascadeMIMICSMod.F90 index 7185febca6..65091755f5 100644 --- a/src/soilbiogeochem/SoilBiogeochemDecompCascadeMIMICSMod.F90 +++ b/src/soilbiogeochem/SoilBiogeochemDecompCascadeMIMICSMod.F90 @@ -755,7 +755,8 @@ end subroutine init_decompcascade_mimics subroutine decomp_rates_mimics(bounds, num_bgc_soilc, filter_bgc_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_bgc_soilc, filter_bgc_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_bgc_soilc, filter_bgc_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 @@ -881,6 +885,10 @@ subroutine decomp_rates_mimics(bounds, num_bgc_soilc, filter_bgc_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() @@ -1298,6 +1306,11 @@ subroutine decomp_rates_mimics(bounds, num_bgc_soilc, filter_bgc_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 diff --git a/src/soilbiogeochem/TillageMod.F90 b/src/soilbiogeochem/TillageMod.F90 new file mode 100644 index 0000000000..4a24daf4c2 --- /dev/null +++ b/src/soilbiogeochem/TillageMod.F90 @@ -0,0 +1,341 @@ +module TillageMod + !----------------------------------------------------------------------- + ! !DESCRIPTION: + ! Module for soil tillage. + ! + ! !USES: + use shr_kind_mod , only : r8 => shr_kind_r8, CS => shr_kind_CS + use shr_log_mod , only : errMsg => shr_log_errMsg + use abortutils , only : endrun + use clm_varctl , only : iulog + use clm_varpar , only : ndecomp_pools + use ColumnType , only : col + use PatchType , only : patch + ! + implicit none + private + ! !PUBLIC MEMBER PROCEDURES + public :: readParams + public :: get_do_tillage + public :: get_apply_tillage_multipliers + public :: get_fraction_tilled + ! !PUBLIC DATA MEMBERS + character(len=CS), public :: tillage_mode ! off, low, high + ! + ! !PRIVATE DATA MEMBERS + integer :: tillage_intensity + integer, parameter :: tillage_off = 0 + integer, parameter :: tillage_low = 1 + integer, parameter :: tillage_high = 2 + logical :: use_original_tillage_phases ! Use buggy tillage phase determination? + real(r8), pointer :: tillage_mults_allphases(:,:) ! (ndecomp_pools, ntill_stages_max) + integer, parameter :: ntill_stages_max = 3 ! How many different tillage phases are there? (Not including all-1 phases.) + integer, parameter :: ntill_intensities_max = 2 ! How many different tillage intensities are allowed (other than "off")? + real(r8) :: max_tillage_depth ! Maximum depth to till (m) + +!============================================================================== +contains +!============================================================================== + + subroutine readParams_namelist(NLFilename) + ! + ! Read namelist parameters related to tillage. + ! + ! !USES: + use spmdMod , only : masterproc, mpicom + use clm_nlUtilsMod , only : find_nlgroup_name + use shr_mpi_mod , only : shr_mpi_bcast + ! !ARGUMENTS: + character(len=*), intent(in) :: NLFilename ! Namelist filename + ! + ! !LOCAL VARIABLES + integer :: nu_nml ! unit for namelist file + integer :: nml_error ! namelist i/o error flag + character(*), parameter :: subname = "('readParams_namelist')" + + namelist /tillage_inparm/ & + tillage_mode, & + use_original_tillage_phases, & + max_tillage_depth + + ! Default values + tillage_mode = 'off' + use_original_tillage_phases = .false. + max_tillage_depth = 0.26_r8 ! Graham et al. (2021) unintentionally used 0.32 + + ! Read tillage namelist + if (masterproc) then + open(newunit=nu_nml, file=trim(NLFilename), status='old', iostat=nml_error ) + call find_nlgroup_name(nu_nml, 'tillage_inparm', status=nml_error) + if (nml_error == 0) then + read(nu_nml, nml=tillage_inparm, iostat=nml_error) + if (nml_error /= 0) then + call endrun(subname // ':: ERROR reading tillage namelist') + end if + else + call endrun(subname // ':: ERROR finding tillage namelist') + end if + close(nu_nml) + endif + call shr_mpi_bcast(tillage_mode, mpicom) + call shr_mpi_bcast(use_original_tillage_phases , mpicom) + call shr_mpi_bcast(max_tillage_depth, mpicom) + + if (masterproc) then + write(iulog,*) ' ' + write(iulog,*) 'tillage settings:' + write(iulog,*) ' tillage_mode = ',tillage_mode + write(iulog,*) ' use_original_tillage_phases = ',use_original_tillage_phases + write(iulog,*) ' max_tillage_depth = ',max_tillage_depth + endif + + ! Assign these + if (tillage_mode == "off") then + tillage_intensity = tillage_off + else if (tillage_mode == "low") then + tillage_intensity = tillage_low + else if (tillage_mode == "high") then + tillage_intensity = tillage_high + else + call endrun(subname // ':: ERROR Unrecognized tillage_mode') + end if + + end subroutine readParams_namelist + + + subroutine readParams_netcdf(ncid) + ! !DESCRIPTION: + ! + ! Read paramfile parameters to be used in tillage. + ! + ! !USES + use ncdio_pio , only : file_desc_t, ncd_io + use clm_varpar, only : ndecomp_pools_max + use SoilBiogeochemDecompCascadeConType, only : no_soil_decomp, century_decomp, mimics_decomp, decomp_method + ! + ! !ARGUMENTS: + type(file_desc_t),intent(inout) :: ncid ! pio netCDF file id + ! + ! !LOCAL VARIABLES: + character(len=32) :: subname = 'readParams_netcdf' + character(len=100) :: errCode = 'Error reading tillage params ' + logical :: readv ! has variable been read in or not + real(r8), allocatable :: tempr(:,:,:) ! temporary to read in constant + character(len=100) :: tString ! temp. var for reading + character(len=3) :: decomp_method_str + + ! Initialize tillage multipliers as all 1, and exit if not tilling + allocate(tillage_mults_allphases(ndecomp_pools, ntill_stages_max)) + tillage_mults_allphases(:,:) = 1.0_r8 + if (.not. get_do_tillage()) then + return + end if + + ! Handle decomposition method + select case( decomp_method ) + case( no_soil_decomp ) + return + case( century_decomp ) + tString = 'bgc_till_decompk_multipliers' + case( mimics_decomp ) + tString = 'mimics_till_decompk_multipliers' + case default + write(decomp_method_str, '(I3)') decomp_method + call endrun('Bad decomp_method = '//decomp_method_str ) + end select + + ! Read off of netcdf file + allocate(tempr(ntill_intensities_max,ndecomp_pools_max,ntill_stages_max)) + call ncd_io(trim(tString), tempr, 'read', ncid, readvar = readv, posNOTonfile = .true.) + if (.not. readv) then + call endrun(msg=trim(errCode)//trim(tString)//errMsg(__FILE__, __LINE__)) + end if + + ! Save + tillage_mults_allphases = tempr(tillage_intensity,1:ndecomp_pools,:) + + end subroutine readParams_netcdf + + + subroutine readParams(ncid, NLFilename) + ! !USES + use ncdio_pio , only : file_desc_t + ! + ! !ARGUMENTS: + type(file_desc_t),intent(inout) :: ncid ! pio netCDF file id + character(len=*), intent(in) :: NLFilename ! Namelist filename + + call readParams_namelist(NLFilename) + call readParams_netcdf(ncid) + + end subroutine readParams + + + function get_do_tillage() + logical :: get_do_tillage + get_do_tillage = tillage_intensity > tillage_off + end function get_do_tillage + + + subroutine get_tillage_multipliers(tillage_mults, idop) + ! !DESCRIPTION: + ! + ! Get the tillage effective multiplier if prognostic crops are on and + ! tillage is turned on. Created by Sam Levis. Modified by Michael Graham + ! to use days past planting (idpp). + ! + ! Modified by Sam Rabin to fix a bug where idpp wasn't actually used, which + ! would affect growing seasons that crossed over into a new calendar year. + ! Previous behavior can be requested with namelist variable use_original_tillage_phases. + ! + ! Original code had two versions depending on cell's GDP, but this seems to + ! have been only an initial effort that was (a) never published and (b) not + ! completely developed. + ! + ! !USES: + use clm_time_manager, only : get_curr_calday, get_curr_days_per_year + use pftconMod , only : ntmp_corn, nirrig_tmp_corn, ntmp_soybean, nirrig_tmp_soybean + use CNPhenologyMod , only : DaysPastPlanting + ! !ARGUMENTS: + real(r8) , intent(inout) :: tillage_mults(:) ! tillage multipliers for this patch + integer , intent(in) :: idop ! patch day of planting + ! + ! !LOCAL VARIABLES: + ! + integer :: day ! julian day + integer :: idpp ! days past planting + integer :: phase ! which tillage phase are we in? + real(r8) dayspyr ! days per year + !----------------------------------------------------------------------- + + ! info from DAYCENT (Melannie Hartman CSU) + ! temp. cereals: P 30 d bef, C 15 d bef, D on day of planting + ! corn, soy : P C D & HW-7 30 d aftr + + phase = 0 + + if (use_original_tillage_phases) then + day = get_curr_calday() + if (day >= idop .and. day < idop+15) then ! based on Point Chisel Tandem Disk multipliers + phase = 1 + else if (day >= idop+15 .and. day < idop+45) then ! based on Field and Row Cultivator multipliers + phase = 2 + else if (day >= idop+45 .and. day < idop+75) then ! based on Rod Weed Row Planter + phase = 3 + end if + else + idpp = DaysPastPlanting(idop) + if (idpp < 15) then ! based on Point Chisel Tandem Disk multipliers + phase = 1 + else if (idpp < 45) then ! based on Field and Row Cultivator multipliers + phase = 2 + else if (idpp < 75) then ! based on Rod Weed Row Planter + phase = 3 + end if + end if + + tillage_mults(:) = 1._r8 + if (phase > 0) then + if (phase > ntill_stages_max) then + call endrun(msg='Tillage phase > ntill_stages_max') + end if + tillage_mults = tillage_mults_allphases(:, phase) + end if + + end subroutine get_tillage_multipliers + + + function get_fraction_tilled(layer_bottom, layer_thickness, max_tillage_depth_gft) result(fraction_tilled) + ! !ARGUMENTS + real(r8), intent(in) :: layer_bottom ! Soil layer interface (between j and j+1) depth (zisoi) + real(r8), intent(in) :: layer_thickness ! Soil layer thickness (dzsoi_decomp) + real(r8) :: max_tillage_depth_gft ! Maximum tillage depth + ! !LOCAL VARIABLES + real(r8) :: layer_top + ! !RESULT + real(r8) :: fraction_tilled ! Fraction of this layer that's within the tillage depth + + ! If the top of the layer is below the max tillage depth, do not till. + layer_top = layer_bottom - layer_thickness + if (layer_top > max_tillage_depth_gft) then + fraction_tilled = 0._r8 + return + end if + + ! Handle zero-thickness layers. This may not be necessary. + if (layer_thickness == 0._r8) then + if (layer_bottom <= max_tillage_depth_gft) then + fraction_tilled = 1._r8 + else + fraction_tilled = 0._r8 + end if + return + end if + + fraction_tilled = max(0._r8, min(1._r8, (max_tillage_depth_gft - layer_top) / layer_thickness)) + + end function get_fraction_tilled + + + subroutine get_apply_tillage_multipliers(idop, c, j, decomp_k) + ! !DESCRIPTION: + ! + ! Multiply decomposition rate constants by tillage coefficients. + ! Written by Sam Rabin, based on original code by Michael Graham. + ! + ! !USES + use pftconMod , only : npcropmin + use clm_varcon, only : zisoi, dzsoi_decomp + use landunit_varcon , only : istcrop + use PatchType , only : patch + ! + ! !ARGUMENTS: + integer , intent(in) :: idop(:) ! patch day of planting + integer , intent(in) :: c ! index of column this is being called for + integer , intent(in) :: j ! index of soil layer this is being called for + real(r8), dimension(:), intent(inout) :: decomp_k ! Output: [real(r8) (:) ] rate constant for decomposition (1./sec) + ! + ! !LOCAL VARIABLES + integer :: p + real :: sumwt ! sum of all patch weights, to check + real(r8), dimension(ndecomp_pools) :: tillage_mults + real(r8), dimension(ndecomp_pools) :: tillage_mults_1patch + real(r8) :: fraction_tilled ! Fraction of this layer that's within the tillage depth + + ! Skip tillage if column is inactive or this layer doesn't get tilled + fraction_tilled = get_fraction_tilled(zisoi(j), dzsoi_decomp(j), max_tillage_depth) + if (.not. col%active(c) .or. fraction_tilled == 0._r8 .or. col%lun_itype(c) /= istcrop) then + return + end if + + ! Initialize tillage multipliers to 0. We will loop through all patches in column, + ! adding patch-weighted multipliers to this. + tillage_mults(:) = 0.0_r8 + + sumwt = 0.0_r8 + do p = col%patchi(c),col%patchf(c) + if (patch%active(p) .and. patch%wtcol(p) /= 0._r8) then + if (patch%itype(p) < npcropmin) then + ! Do not till generic crops + tillage_mults_1patch(:) = 1._r8 + else + call get_tillage_multipliers(tillage_mults_1patch, idop(p)) + end if + tillage_mults = tillage_mults + tillage_mults_1patch * patch%wtcol(p) + sumwt = sumwt + patch%wtcol(p) + end if + end do + if (abs(1.0_r8 - sumwt) > 1.e-6_r8) then + call endrun('ERROR Active crop patch weights does not sum to 1') + end if + + ! Adjust tillage_mults to consider fraction of this layer that's within tillage depth. + tillage_mults = tillage_mults * fraction_tilled & + + 1._r8 * (1._r8 - fraction_tilled) + + ! Apply + decomp_k = decomp_k * tillage_mults(:) + + end subroutine get_apply_tillage_multipliers + +end module TillageMod diff --git a/src/soilbiogeochem/test/CMakeLists.txt b/src/soilbiogeochem/test/CMakeLists.txt index 988cb531b5..22d8fba60f 100644 --- a/src/soilbiogeochem/test/CMakeLists.txt +++ b/src/soilbiogeochem/test/CMakeLists.txt @@ -1 +1,2 @@ add_subdirectory(ACSpinup_test) +add_subdirectory(tillage_test) diff --git a/src/soilbiogeochem/test/tillage_test/CMakeLists.txt b/src/soilbiogeochem/test/tillage_test/CMakeLists.txt new file mode 100644 index 0000000000..fbc550dd03 --- /dev/null +++ b/src/soilbiogeochem/test/tillage_test/CMakeLists.txt @@ -0,0 +1,3 @@ +add_pfunit_ctest(tillage + TEST_SOURCES "test_tillage.pf" + LINK_LIBRARIES clm csm_share esmf_wrf_timemgr) diff --git a/src/soilbiogeochem/test/tillage_test/test_tillage.pf b/src/soilbiogeochem/test/tillage_test/test_tillage.pf new file mode 100644 index 0000000000..77721aedd4 --- /dev/null +++ b/src/soilbiogeochem/test/tillage_test/test_tillage.pf @@ -0,0 +1,54 @@ +module test_tillage + + ! Tests of the functions in TillageMod + + use funit + + use shr_kind_mod , only : r8 => shr_kind_r8 + use TillageMod, only : get_fraction_tilled + + implicit none + save + + real(r8), parameter :: tol = 1.e-8_r8 + +contains + + @Test + subroutine test_get_fraction_tilled() + + integer, parameter :: nlayers = 5 + real(r8) :: zisoi(nlayers) ! Depth of soil interfaces (bottom of each layer) + real(r8) :: dzsoi_decomp(nlayers) ! Thickness of soil layers + integer :: j ! Soil layer + real(r8), parameter :: max_tillage_depth = 0.21_r8 + + zisoi = [0.01_r8, 0.05_r8, 0.1_r8, 0.5_r8, 1._r8] + dzsoi_decomp = [zisoi(1) - 0._r8, & + zisoi(2) - zisoi(1), & + zisoi(3) - zisoi(2), & + zisoi(4) - zisoi(3), & + zisoi(5) - zisoi(4)] + + @assertEqual(1._r8, get_fraction_tilled(zisoi(1), dzsoi_decomp(1), max_tillage_depth)) + @assertEqual(1._r8, get_fraction_tilled(zisoi(2), dzsoi_decomp(2), max_tillage_depth)) + @assertEqual(1._r8, get_fraction_tilled(zisoi(3), dzsoi_decomp(3), max_tillage_depth)) + @assertEqual(0.275_r8, get_fraction_tilled(zisoi(4), dzsoi_decomp(4), max_tillage_depth), tolerance=tol) + @assertEqual(0._r8, get_fraction_tilled(zisoi(5), dzsoi_decomp(5), max_tillage_depth)) + + end subroutine test_get_fraction_tilled + + + @Test + + subroutine test_get_fraction_tilled_0thickness() + real(r8), parameter :: max_tillage_depth = 0.5_r8 + + @assertEqual(1._r8, get_fraction_tilled(0.4_r8, 0._r8, max_tillage_depth)) + @assertEqual(1._r8, get_fraction_tilled(0.5_r8, 0._r8, max_tillage_depth)) + @assertEqual(0._r8, get_fraction_tilled(0.6_r8, 0._r8, max_tillage_depth)) + end subroutine test_get_fraction_tilled_0thickness + + + +end module test_tillage diff --git a/tools/contrib/add_tillage_to_paramsfile.py b/tools/contrib/add_tillage_to_paramsfile.py new file mode 100644 index 0000000000..ebea329ef6 --- /dev/null +++ b/tools/contrib/add_tillage_to_paramsfile.py @@ -0,0 +1,171 @@ +import xarray as xr +import numpy as np +import os +import subprocess +import argparse +import sys + +def make_dataarray(np_array, decomp_model, ntill_intensities_max, ndecomp_pools_max, ntill_stages_max): + intensities_dim = xr.IndexVariable( + dims = "ntill_intensities_max", + data = np.arange(ntill_intensities_max), + ) + pools_dim = xr.IndexVariable( + dims = "ndecomp_pools_max", + data = np.arange(ndecomp_pools_max), + ) + stages_dim = xr.IndexVariable( + dims = "ntill_stages_max", + data = np.arange(ntill_stages_max), + ) + + # Name DataArray + if decomp_model.lower() == "mimics": + da_name = "mimics" + elif decomp_model.lower() == "century": + da_name = "bgc" + da_name += "_till_decompk_multipliers" + + da = xr.DataArray( + data = np_array, + dims = { + "ntill_intensities_max": intensities_dim, + "ndecomp_pools_max": pools_dim, + "ntill_stages_max": stages_dim + }, + name = da_name, + attrs = { + "long_name": f"Value by which decomp_k should be multiplied during tillage with {decomp_model} soil", + "units": "unitless", + } + ) + + # netCDF variable needs dimensions reversed from how they're specified in code + da = da.transpose() + + return da + +def main(file_in, + file_out): + # Get git info + thisDir = os.path.dirname(__file__) + git_status = subprocess.run( + ["git", "status"], + capture_output=True, + ) + git_status = git_status.stdout.decode() + repo_is_clean = "working tree clean" in git_status + if not repo_is_clean: + print("WARNING: Repo not clean; will not save params file.") + print(git_status) + git_log = subprocess.run( + ["git", "log", "-1"], + capture_output=True, + ) + git_log = git_log.stdout.decode() + + + # Set up dimensions + ds0 = xr.open_dataset(file_in) + ntill_intensities_max = 2 + ndecomp_pools_max = ds0.dims["ndecomp_pools_max"] + ntill_stages_max = 3 + tillage_shape_ips = (ntill_intensities_max, ndecomp_pools_max, ntill_stages_max) + tillage_shape_ps = (ndecomp_pools_max, ntill_stages_max) + + + # Fill CENTURY array + # Define pool indices + i_litr_min = 0 # 1 in FORTRAN, but Python is 0-indexed + i_met_lit = i_litr_min + i_cel_lit = i_met_lit + 1 + i_lig_lit = i_cel_lit + 1 + i_act_som = i_lig_lit + 1 + i_slo_som = i_act_som + 1 + i_pas_som = i_slo_som + 1 + tillage_century = np.full(tillage_shape_ips, 1.0) + tillage_century_lo = np.full(tillage_shape_ps, 1.0) + tillage_century_lo[i_act_som,:] = np.array([1.0, 1.0, 1.0]) + tillage_century_lo[i_slo_som,:] = np.array([3.0, 1.6, 1.3]) + tillage_century_lo[i_pas_som,:] = np.array([3.0, 1.6, 1.3]) + tillage_century_lo[i_cel_lit,:] = np.array([1.5, 1.5, 1.1]) + tillage_century_lo[i_lig_lit,:] = np.array([1.5, 1.5, 1.1]) + tillage_century[0,:,:] = tillage_century_lo + tillage_century_hi = np.full(tillage_shape_ps, 1.0) + tillage_century_hi[i_act_som,:] = np.array([1.2, 1.0, 1.0]) + tillage_century_hi[i_slo_som,:] = np.array([4.8, 3.5, 2.5]) + tillage_century_hi[i_pas_som,:] = np.array([4.8, 3.5, 2.5]) + tillage_century_hi[i_cel_lit,:] = np.array([1.8, 1.5, 1.1]) + tillage_century_hi[i_lig_lit,:] = np.array([1.8, 1.5, 1.1]) + tillage_century[1,:,:] = tillage_century_hi + + + # Fill MIMICS array + i_litr_min = 1 + i_met_lit = i_litr_min + i_str_lit = i_met_lit + 1 + i_avl_som = i_str_lit + 1 + i_chem_som = i_avl_som + 1 + i_phys_som = i_chem_som + 1 + tillage_mimics = np.full(tillage_shape_ips, 1.0) + tillage_mimics[:,i_avl_som,:] = tillage_century[:,i_act_som,:] + tillage_mimics[:,i_chem_som,:] = tillage_century[:,i_slo_som,:] + tillage_mimics[:,i_phys_som,:] = tillage_century[:,i_pas_som,:] + if not np.array_equal(tillage_century[:,i_cel_lit,:], tillage_century[:,i_lig_lit,:]): + raise RuntimeError("How to combine 2 CENTURY litter pools into 1 MIMICS litter pool?") + tillage_mimics[:,i_str_lit,:] = tillage_century[:,i_cel_lit,:] + + # Make DataArrays + tillage_century_da = make_dataarray( + tillage_century, "CENTURY", + ntill_intensities_max, ndecomp_pools_max, ntill_stages_max) + tillage_mimics_da = make_dataarray( + tillage_mimics, "MIMICS", + ntill_intensities_max, ndecomp_pools_max, ntill_stages_max) + + if not repo_is_clean: + raise RuntimeError("Clean up git repo before trying to save!") + + ds1 = ds0.copy() + ds0.close() + + ds1[tillage_century_da.name] = tillage_century_da + ds1[tillage_mimics_da.name] = tillage_mimics_da + ds1.attrs['latest_git_log'] = git_log + + ds1.to_netcdf(file_out, format="NETCDF3_CLASSIC") + + +if __name__ == "__main__": + ############################### + ### Process input arguments ### + ############################### + parser = argparse.ArgumentParser( + description="Adds tillage parameters to a CLM parameter file (netCDF)." + ) + + # Define arguments + parser.add_argument( + "-i", + "--input-file", + help="Parameter file (netCDF) to which you wish to add tillage parameters.", + type=str, + required=True, + ) + parser.add_argument( + "-o", + "--output-file", + help="Output parameter file.", + type=str, + required=True, + ) + + # Get arguments + args = parser.parse_args(sys.argv[1:]) + + ########### + ### Run ### + ########### + main(os.path.realpath(args.input_file), + os.path.realpath(args.output_file), + ) \ No newline at end of file