diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm
index c2fdfd6ff8..204a4de350 100755
--- a/bld/CLMBuildNamelist.pm
+++ b/bld/CLMBuildNamelist.pm
@@ -1006,6 +1006,13 @@ sub setup_cmdl_bgc {
if ( &value_is_true($nl->get_value($var)) && $nl_flags->{'soil_decomp_method'} ne "CENTURYKoven2013" ) {
$log->fatal_error("$var can only be on with CENTURYKoven2013 soil decomposition");
}
+
+ # Set use_nvmovement
+ $var = "use_nvmovement";
+ add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var);
+ if ( &value_is_true($nl->get_value($var)) && !&value_is_true($nl_flags->{'use_nitrif_denitrif'}) ) {
+ $log->fatal_error("$var cannot be on with use_nitrif_denitrif = .false.");
+ }
} # end bgc
diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml
index 1d9c85c33c..20ab70bb72 100644
--- a/bld/namelist_files/namelist_defaults_ctsm.xml
+++ b/bld/namelist_files/namelist_defaults_ctsm.xml
@@ -2431,6 +2431,7 @@ lnd/clm2/surfdata_esmf/NEON/ctsm5.3.0/surfdata_1x1_NEON_TOOL_hist_2000_78pfts_c2
.true.
.false.
.true.
+.false.
FvCB1980
diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml
index af3550ff36..69a243bd27 100644
--- a/bld/namelist_files/namelist_definition_ctsm.xml
+++ b/bld/namelist_files/namelist_definition_ctsm.xml
@@ -1002,6 +1002,12 @@ How LUNA and Photosynthesis (if needed) will get Leaf nitrogen content
lnc_opt = true get from leaf N from CN model
lnc_opt = false get based on LAI and fixed CN ratio from parameter file
+
+ use_nvmovement = true use soil nitrogen vertical movement
+ use_nvmovement = false do not use soil nitrogen vertical movement
+ use_nvmovement cannot be true while use_nitrif_denitrif is false
+
diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml
index ab53ea371d..1ea0dc67a2 100644
--- a/cime_config/testdefs/testlist_clm.xml
+++ b/cime_config/testdefs/testlist_clm.xml
@@ -595,6 +595,24 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/cime_config/testdefs/testmods_dirs/clm/nvmovement/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/nvmovement/include_user_mods
new file mode 100644
index 0000000000..fe0e18cf88
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/nvmovement/include_user_mods
@@ -0,0 +1 @@
+../default
diff --git a/cime_config/testdefs/testmods_dirs/clm/nvmovement/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/nvmovement/user_nl_clm
new file mode 100644
index 0000000000..258d3bbdf4
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/nvmovement/user_nl_clm
@@ -0,0 +1 @@
+ use_nvmovement = .true.
diff --git a/doc/ChangeLog b/doc/ChangeLog
index 428ab62539..69dc644440 100644
--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -1,4 +1,80 @@
===============================================================
+Tag name: ctsm5.3.072
+Originator(s): jinmuluo (Jinmu Luo, Cornell University)
+Date: Tue 26 Aug 2025 02:53:26 PM MDT
+One-line Summary: New vertical movement scheme for soil nitrate
+
+Purpose and description of changes
+----------------------------------
+
+ Introducing a vertical movement scheme for soil nitrate but setting its use to .false. by default.
+ Details about the new scheme in the PR #2992.
+
+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.]
+
+[ ] clm6_0
+
+[ ] clm5_0
+
+[ ] ctsm5_0-nwp
+
+[ ] clm4_5
+
+
+Bugs fixed
+----------
+List of CTSM issues fixed (include CTSM Issue # and description) [one per line]:
+ No corresponding github issue
+
+Notes of particular relevance for users
+---------------------------------------
+Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables):
+ use_nvmovement can be set by user
+
+Changes made to namelist defaults (e.g., changed parameter values):
+ use_nvmovement defaults to .false.
+
+Changes to documentation:
+ Planned for Jinmu's NSF-NCAR visit this Fall
+
+Notes of particular relevance for developers:
+---------------------------------------------
+Changes to tests or testing:
+ New tests with nvmovement testmods, one with matrixcnOn and one with it off
+ ERP_D_Ld5.f10_f10_mg37.I1850Clm60BgcCrop.derecho_intel.clm-nvmovement--clm-matrixcnOn
+ ERP_D.f10_f10_mg37.IHistClm60Bgc.derecho_intel.clm-decStart--clm-nvmovement
+
+Testing summary:
+----------------
+
+ [PASS means all tests PASS; OK means tests PASS other than expected fails.]
+
+ build-namelist tests (if CLMBuildNamelist.pm has changed):
+
+ derecho - PASS with all commits up to and including 1cf9859
+
+ regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing):
+
+ derecho ----- OK
+ izumi ------- OK
+
+Answer changes
+--------------
+Changes answers relative to baseline: No
+
+Other details
+-------------
+Pull Requests that document the changes (include PR ids):
+ https://github.com/ESCOMP/ctsm/pull/2992
+
+===============================================================
+===============================================================
Tag name: ctsm5.3.071
Originator(s): samrabin (Sam Rabin, UCAR/TSS)
Date: Fri Aug 22 13:49:25 MDT 2025
diff --git a/doc/ChangeSum b/doc/ChangeSum
index f01334b0e4..b8340937d7 100644
--- a/doc/ChangeSum
+++ b/doc/ChangeSum
@@ -1,5 +1,6 @@
Tag Who Date Summary
============================================================================================================================
+ ctsm5.3.072 jinmuluo 08/26/2025 New vertical movement scheme for soil nitrate
ctsm5.3.071 samrabin 08/22/2025 Merge b4b-dev to master
ctsm5.3.070 glemieux 08/22/2025 Update default FATES parameter file and add FATES managed fire namelist option
ctsm5.3.069 samrabin 08/12/2025 Add SystemTests to run subset_data and then CTSM
diff --git a/src/biogeochem/CNDriverMod.F90 b/src/biogeochem/CNDriverMod.F90
index 9a51d9f616..d14e619b03 100644
--- a/src/biogeochem/CNDriverMod.F90
+++ b/src/biogeochem/CNDriverMod.F90
@@ -44,6 +44,7 @@ module CNDriverMod
use SoilWaterRetentionCurveMod , only : soil_water_retention_curve_type
use CLMFatesInterfaceMod , only : hlm_fates_interface_type
use CropReprPoolsMod , only : nrepr
+ use SoilHydrologyType , only: soilhydrology_type
!
! !PUBLIC TYPES:
implicit none
@@ -1013,7 +1014,7 @@ subroutine CNDriverLeaching(bounds, &
c13_cnveg_carbonstate_inst,c14_cnveg_carbonstate_inst, &
c13_cnveg_carbonflux_inst,c14_cnveg_carbonflux_inst, &
c13_soilbiogeochem_carbonstate_inst,c14_soilbiogeochem_carbonstate_inst,&
- c13_soilbiogeochem_carbonflux_inst,c14_soilbiogeochem_carbonflux_inst)
+ c13_soilbiogeochem_carbonflux_inst,c14_soilbiogeochem_carbonflux_inst, soilhydrology_inst)
!
! !DESCRIPTION:
! Update the nitrogen leaching rate as a function of soluble mineral N and total soil water outflow.
@@ -1028,6 +1029,8 @@ subroutine CNDriverLeaching(bounds, &
use clm_time_manager , only: is_first_step_of_this_run_segment,is_beg_curr_year,is_end_curr_year,get_curr_date
use CNSharedParamsMod , only: use_matrixcn
use SoilBiogeochemDecompCascadeConType, only: use_soil_matrixcn
+ use SoilNitrogenMovementMod , only: SoilNitrogenMovement
+ use clm_varctl, only : use_nvmovement
!
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds
@@ -1045,6 +1048,7 @@ subroutine CNDriverLeaching(bounds, &
type(cnveg_carbonflux_type) , intent(inout) :: cnveg_carbonflux_inst
type(cnveg_carbonstate_type) , intent(inout) :: cnveg_carbonstate_inst
type(soilstate_type) , intent(inout) :: soilstate_inst
+ type(soilhydrology_type) , intent(in) :: soilhydrology_inst
type(soilbiogeochem_state_type) , intent(inout) :: soilbiogeochem_state_inst
type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst
type(soilbiogeochem_carbonstate_type) , intent(inout) :: soilbiogeochem_carbonstate_inst
@@ -1062,8 +1066,16 @@ subroutine CNDriverLeaching(bounds, &
type(soilbiogeochem_carbonflux_type) , intent(inout) :: c14_soilbiogeochem_carbonflux_inst
integer p,fp,yr,mon,day,sec
!-----------------------------------------------------------------------
-
- ! Mineral nitrogen dynamics (deposition, fixation, leaching)
+
+ ! soil nitrate fast aqueous movement, leaching will be evaluted here
+ if (use_nitrif_denitrif .and. use_nvmovement) then
+ call t_startf('SoilNitrogenMovementMod')
+ call SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterstatebulk_inst, &
+ soilstate_inst, soilhydrology_inst, soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst)
+ call t_stopf('SoilNitrogenMovementMod')
+ end if
+
+ ! Mineral nitrogen dynamics: deposition, fixation. If use_nvmoment false, also leaching.
call t_startf('SoilBiogeochemNLeaching')
call SoilBiogeochemNLeaching(bounds, num_bgc_soilc, filter_bgc_soilc, &
diff --git a/src/biogeochem/CNVegetationFacade.F90 b/src/biogeochem/CNVegetationFacade.F90
index b47237690d..9ff26dd56d 100644
--- a/src/biogeochem/CNVegetationFacade.F90
+++ b/src/biogeochem/CNVegetationFacade.F90
@@ -96,6 +96,7 @@ module CNVegetationFacade
use SoilBiogeochemPrecisionControlMod , only: SoilBiogeochemPrecisionControl
use SoilWaterRetentionCurveMod , only : soil_water_retention_curve_type
use CLMFatesInterfaceMod , only : hlm_fates_interface_type
+ use SoilHydrologyType , only : soilhydrology_type
!
implicit none
private
@@ -1074,7 +1075,7 @@ subroutine EcosystemDynamicsPostDrainage(this, bounds, num_allc, filter_allc, &
soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, &
c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, &
c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, &
- soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst)
+ soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst, soilhydrology_inst)
!
! !DESCRIPTION:
! Do the main science for CN vegetation that needs to be done after hydrology-drainage
@@ -1113,6 +1114,7 @@ subroutine EcosystemDynamicsPostDrainage(this, bounds, num_allc, filter_allc, &
type(soilbiogeochem_carbonstate_type) , intent(inout) :: c14_soilbiogeochem_carbonstate_inst
type(soilbiogeochem_nitrogenflux_type) , intent(inout) :: soilbiogeochem_nitrogenflux_inst
type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst
+ type(soilhydrology_type) , intent(in) :: soilhydrology_inst
!
! !LOCAL VARIABLES:
@@ -1135,7 +1137,7 @@ subroutine EcosystemDynamicsPostDrainage(this, bounds, num_allc, filter_allc, &
this%c13_cnveg_carbonstate_inst,this%c14_cnveg_carbonstate_inst, &
this%c13_cnveg_carbonflux_inst,this%c14_cnveg_carbonflux_inst, &
c13_soilbiogeochem_carbonstate_inst,c14_soilbiogeochem_carbonstate_inst,&
- c13_soilbiogeochem_carbonflux_inst,c14_soilbiogeochem_carbonflux_inst)
+ c13_soilbiogeochem_carbonflux_inst,c14_soilbiogeochem_carbonflux_inst, soilhydrology_inst)
! Set controls on very low values in critical state variables
diff --git a/src/biogeophys/SoilHydrologyType.F90 b/src/biogeophys/SoilHydrologyType.F90
index 07ad2ca45b..9be9073bc0 100644
--- a/src/biogeophys/SoilHydrologyType.F90
+++ b/src/biogeophys/SoilHydrologyType.F90
@@ -51,6 +51,8 @@ Module SoilHydrologyType
real(r8), pointer :: top_ice_col (:) ! col VIC ice len in top layers
real(r8), pointer :: top_moist_limited_col(:) ! col VIC soil moisture in top layers, limited to no greater than top_max_moist_col
real(r8), pointer :: ice_col (:,:) ! col VIC soil ice (kg/m2) for VIC soil layers
+ real(r8), pointer :: qout_col (:,:) ! flux of water out of soil layer [mm H2O/s]
+ real(r8), pointer :: qin_col (:,:) ! flux of water into soil layer [mm H2O/s]
contains
@@ -142,6 +144,8 @@ subroutine InitAllocate(this, bounds)
allocate(this%top_ice_col (begc:endc)) ; this%top_ice_col (:) = nan
allocate(this%top_moist_limited_col(begc:endc)) ; this%top_moist_limited_col(:) = nan
allocate(this%ice_col (begc:endc,nlayert)) ; this%ice_col (:,:) = nan
+ allocate(this%qout_col (begc:endc,nlevsoi)) ; this%qout_col (:,:) = nan
+ allocate(this%qin_col (begc:endc,nlevsoi)) ; this%qin_col (:,:) = nan
end subroutine InitAllocate
@@ -149,7 +153,7 @@ end subroutine InitAllocate
subroutine InitHistory(this, bounds, use_aquifer_layer)
!
! !USES:
- use histFileMod , only : hist_addfld1d
+ use histFileMod , only : hist_addfld1d, hist_addfld2d
!
! !ARGUMENTS:
class(soilhydrology_type) :: this
@@ -192,6 +196,16 @@ subroutine InitHistory(this, bounds, use_aquifer_layer)
avgflag='A', long_name='perched water table depth (natural vegetated and crop landunits only)', &
ptr_col=this%zwt_perched_col, l2g_scale_type='veg')
+ this%qout_col(begc:endc, :) = spval
+ call hist_addfld2d (fname="QOUT", units='mm H2O/s', type2d='levsoi', &
+ avgflag='A', long_name='flux of water out of soil layer', &
+ ptr_col=this%qout_col, default='inactive')
+
+ this%qin_col(begc:endc, :) = spval
+ call hist_addfld2d (fname="QIN", units='mm H2O/s', type2d='levsoi', &
+ avgflag='A', long_name='flux of water into soil layer', &
+ ptr_col=this%qin_col, default='inactive')
+
end subroutine InitHistory
!-----------------------------------------------------------------------
@@ -287,7 +301,6 @@ subroutine Restart(this, bounds, ncid, flag)
character(len=*) , intent(in) :: flag ! 'read' or 'write'
!
! !LOCAL VARIABLES:
- integer :: j,c ! indices
logical :: readvar ! determine if variable is on initial file
!-----------------------------------------------------------------------
diff --git a/src/biogeophys/SoilWaterMovementMod.F90 b/src/biogeophys/SoilWaterMovementMod.F90
index 85bcf42c5e..dd4f48090f 100644
--- a/src/biogeophys/SoilWaterMovementMod.F90
+++ b/src/biogeophys/SoilWaterMovementMod.F90
@@ -579,7 +579,8 @@ subroutine soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, &
zwt => soilhydrology_inst%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m)
icefrac => soilhydrology_inst%icefrac_col , & ! Input: [real(r8) (:,:) ] fraction of ice
hkdepth => soilhydrology_inst%hkdepth_col , & ! Input: [real(r8) (:) ] decay factor (m)
-
+ qout_col => soilhydrology_inst%qout_col , & ! Output: [real(r8) (:,:) ] soil water out of the bottom, mm h2o/s
+ qin_col => soilhydrology_inst%qin_col , & ! Output: [real(r8) (:,:) ] soil water into the bottom, mm h2o/s
smpmin => soilstate_inst%smpmin_col , & ! Input: [real(r8) (:) ] restriction for min of soil potential (mm)
watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity)
hksat => soilstate_inst%hksat_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity at saturation (mm H2O /s)
@@ -787,7 +788,8 @@ subroutine soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, &
amx(c,j) = 0._r8
bmx(c,j) = dzmm(c,j)*(sdamp+1._r8/dtime) + dqodw1(c,j)
cmx(c,j) = dqodw2(c,j)
-
+ qin_col(c,j) = qin(c,j)
+ qout_col(c,j) = qout(c,j)
end do
! Nodes j=2 to j=nlevsoi-1
@@ -811,7 +813,8 @@ subroutine soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, &
amx(c,j) = -dqidw0(c,j)
bmx(c,j) = dzmm(c,j)/dtime - dqidw1(c,j) + dqodw1(c,j)
cmx(c,j) = dqodw2(c,j)
-
+ qin_col(c,j) = qin(c,j)
+ qout_col(c,j) = qout(c,j)
end do
end do
@@ -833,6 +836,8 @@ subroutine soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, &
amx(c,j) = -dqidw0(c,j)
bmx(c,j) = dzmm(c,j)/dtime - dqidw1(c,j) + dqodw1(c,j)
cmx(c,j) = 0._r8
+ qin_col(c,j) = qin(c,j)
+ qout_col(c,j) = qout(c,j)
! next set up aquifer layer; hydrologically inactive
rmx(c,j+1) = 0._r8
@@ -883,6 +888,8 @@ subroutine soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, &
amx(c,j+1) = -dqidw0(c,j+1)
bmx(c,j+1) = dzmm(c,j+1)/dtime - dqidw1(c,j+1) + dqodw1(c,j+1)
cmx(c,j+1) = 0._r8
+ qin_col(c,j) = qin(c,j)
+ qout_col(c,j) = qout(c,j)
endif
end do
@@ -1154,7 +1161,8 @@ subroutine soilwater_moisture_form(bounds, num_hydrologyc, &
qcharge => soilhydrology_inst%qcharge_col , & ! Input: [real(r8) (:) ] aquifer recharge rate (mm/s)
zwt => soilhydrology_inst%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m)
-
+ qout_col => soilhydrology_inst%qout_col , & ! Output: [real(r8) (:,:) ] soil water out of the bottom, mm h2o/s
+ qin_col => soilhydrology_inst%qin_col , & ! Output: [real(r8) (:,:) ] soil water into the bottom, mm h2o/s
watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity)
smp_l => soilstate_inst%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matrix potential [mm]
hk_l => soilstate_inst%hk_l_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity (mm/s)
@@ -1402,7 +1410,8 @@ subroutine soilwater_moisture_form(bounds, num_hydrologyc, &
! call endrun(subname // ':: negative soil moisture values found!')
endif
end do
-
+ qin_col(c,1:nlayers) = qin(c,1:nlayers)
+ qout_col(c,1:nlayers) = qout(c,1:nlayers)
end do ! spatial loop
diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90
index 2cea242a67..a436b75a97 100644
--- a/src/main/clm_driver.F90
+++ b/src/main/clm_driver.F90
@@ -1157,7 +1157,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, &
c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, &
c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, &
- soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst)
+ soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst, soilhydrology_inst)
call t_stopf('EcosysDynPostDrainage')
end if
diff --git a/src/main/clm_varcon.F90 b/src/main/clm_varcon.F90
index e45f5c440e..985660142a 100644
--- a/src/main/clm_varcon.F90
+++ b/src/main/clm_varcon.F90
@@ -127,6 +127,7 @@ module clm_varcon
real(r8), public, parameter :: g_to_mg = 1.0e3_r8 ! coefficient to convert g to mg
real(r8), public, parameter :: cm3_to_m3 = 1.0e-6_r8 ! coefficient to convert cm3 to m3
real(r8), public, parameter :: pct_to_frac = 1.0e-2_r8 ! coefficient to convert % to fraction
+ real(r8), public, parameter :: mmh2o_to_m3h2o_per_m2 = 1.0e-3_r8 ! coefficient to convert mm H2O to m3 H2O/m2
!!! C13
real(r8), public, parameter :: preind_atm_del13c = -6.0_r8 ! preindustrial value for atmospheric del13C
diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90
index dc323893f2..9a0d2901b9 100644
--- a/src/main/clm_varctl.F90
+++ b/src/main/clm_varctl.F90
@@ -522,6 +522,7 @@ module clm_varctl
logical, public :: use_lch4 = .true.
logical, public :: use_nitrif_denitrif = .true.
+ logical, public :: use_nvmovement = .false.
logical, public :: use_extralakelayers = .false.
logical, public :: use_vichydro = .false.
logical, public :: use_cn = .false.
diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90
index 13f86edf78..089503dc8b 100644
--- a/src/main/controlMod.F90
+++ b/src/main/controlMod.F90
@@ -319,7 +319,7 @@ subroutine control_init(dtime)
use_lch4, use_nitrif_denitrif, use_extralakelayers, &
use_vichydro, use_cn, use_cndv, use_crop, use_fertilizer, &
use_grainproduct, use_snicar_frc, use_vancouver, use_mexicocity, use_noio, &
- use_nguardrail, crop_residue_removal_frac, flush_gdd20
+ use_nguardrail, crop_residue_removal_frac, flush_gdd20, use_nvmovement
! SNICAR
namelist /clm_inparm/ &
@@ -730,6 +730,7 @@ subroutine control_spmd()
call mpi_bcast (use_lch4, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (use_nitrif_denitrif, 1, MPI_LOGICAL, 0, mpicom, ier)
+ call mpi_bcast (use_nvmovement, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (use_extralakelayers, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (use_vichydro, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (use_cn, 1, MPI_LOGICAL, 0, mpicom, ier)
@@ -1028,6 +1029,7 @@ subroutine control_print ()
write(iulog,*) 'process control parameters:'
write(iulog,*) ' use_lch4 = ', use_lch4
write(iulog,*) ' use_nitrif_denitrif = ', use_nitrif_denitrif
+ write(iulog,*) ' use_nvmovement = ', use_nvmovement
write(iulog,*) ' use_extralakelayers = ', use_extralakelayers
write(iulog,*) ' use_vichydro = ', use_vichydro
write(iulog,*) ' use_excess_ice = ', use_excess_ice
diff --git a/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 b/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90
index a646feb1d7..f0a42b379c 100644
--- a/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90
+++ b/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90
@@ -10,7 +10,7 @@ module SoilBiogeochemNLeachingMod
use abortutils , only : endrun
use decompMod , only : bounds_type
use clm_varcon , only : dzsoi_decomp, zisoi
- use clm_varctl , only : use_nitrif_denitrif
+ use clm_varctl , only : use_nitrif_denitrif, use_nvmovement
use SoilBiogeochemNitrogenStateType , only : soilbiogeochem_nitrogenstate_type
use SoilBiogeochemNitrogenFluxType , only : soilbiogeochem_nitrogenflux_type
use WaterStateBulkType , only : waterstatebulk_type
@@ -212,16 +212,21 @@ subroutine SoilBiogeochemNLeaching(bounds, num_bgc_soilc, filter_bgc_soilc, &
if (h2osoi_liq(c,j) > 0._r8) then
disn_conc = (sf_no3 * smin_no3_vr(c,j) * col%dz(c,j) )/(h2osoi_liq(c,j) )
end if
- !
- ! calculate the N leaching flux as a function of the dissolved
- ! concentration and the sub-surface drainage flux
- smin_no3_leached_vr(c,j) = disn_conc * drain_tot(c) * h2osoi_liq(c,j) / ( tot_water(c) * col%dz(c,j) )
- !
- ! ensure that leaching rate isn't larger than soil N pool
- smin_no3_leached_vr(c,j) = min(smin_no3_leached_vr(c,j), smin_no3_vr(c,j) / dt )
- !
- ! limit the leaching flux to a positive value
- smin_no3_leached_vr(c,j) = max(smin_no3_leached_vr(c,j), 0._r8)
+ ! Evaluating the leaching flux in this module, if use_nvmovement is not true
+ ! and in SoilNitrogenMovementMod if use_nvmovement is true
+ if (.not. use_nvmovement) then
+ !
+ ! calculate the N leaching flux as a function of the dissolved
+ ! concentration and the sub-surface drainage flux
+ smin_no3_leached_vr(c,j) = disn_conc * drain_tot(c) * h2osoi_liq(c,j) / ( tot_water(c) * col%dz(c,j) )
+ !
+ ! ensure that leaching rate isn't larger than soil N pool
+ smin_no3_leached_vr(c,j) = min(smin_no3_leached_vr(c,j), smin_no3_vr(c,j) / dt )
+ !
+ ! limit the leaching flux to a positive value
+ smin_no3_leached_vr(c,j) = max(smin_no3_leached_vr(c,j), 0._r8)
+ end if
+
!
!
! calculate the N loss from surface runoff, assuming a shallow mixing of surface waters into soil and removal based on runoff
@@ -242,13 +247,17 @@ subroutine SoilBiogeochemNLeaching(bounds, num_bgc_soilc, filter_bgc_soilc, &
! limit the flux to a positive value
smin_no3_runoff_vr(c,j) = max(smin_no3_runoff_vr(c,j), 0._r8)
- ! limit the flux based on current smin_no3 state
- ! only let at most the assumed soluble fraction
- ! of smin_no3 be leached on any given timestep
- smin_no3_leached_vr(c,j) = min(smin_no3_leached_vr(c,j), (sf_no3 * smin_no3_vr(c,j))/dt)
+ ! Evaluating the leaching flux in this module, if use_nvmovement is not true
+ ! and in SoilNitrogenMovementMod if use_nvmovement is true
+ if (.not. use_nvmovement) then
+ ! limit the flux based on current smin_no3 state
+ ! only let at most the assumed soluble fraction
+ ! of smin_no3 be leached on any given timestep
+ smin_no3_leached_vr(c,j) = min(smin_no3_leached_vr(c,j), (sf_no3 * smin_no3_vr(c,j))/dt)
- ! limit the flux to a positive value
- smin_no3_leached_vr(c,j) = max(smin_no3_leached_vr(c,j), 0._r8)
+ ! limit the flux to a positive value
+ smin_no3_leached_vr(c,j) = max(smin_no3_leached_vr(c,j), 0._r8)
+ end if
end do
end do
diff --git a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 b/src/soilbiogeochem/SoilNitrogenMovementMod.F90
new file mode 100644
index 0000000000..cc352fd827
--- /dev/null
+++ b/src/soilbiogeochem/SoilNitrogenMovementMod.F90
@@ -0,0 +1,247 @@
+module SoilNitrogenMovementMod
+
+ !------------------------------------------------------------------------
+ ! DESCRIPTION
+ ! implementation of Soil-Water-Atmosphere-Plant (SWAP3.2)
+ ! and Pantakar 1980 algorithm to Community Land Model
+ ! SWAP3.2 website: https://edepot.wur.nl/39776
+ ! Author: Jinmu Luo, Cornell EAS, April 1 2024
+
+ use decompMod , only : bounds_type
+ use shr_kind_mod , only : r8 => shr_kind_r8
+ use shr_infnan_mod , only : isnan => shr_infnan_isnan
+ use shr_infnan_mod , only : isinf => shr_infnan_isinf
+ use clm_varctl , only : iulog
+ use spmdMod , only : masterproc
+ use abortutils , only : endrun
+ use clm_time_manager , only : get_step_size_real
+ use clm_time_manager , only : get_curr_date
+ use SoilBiogeochemNitrogenFluxType , only : soilbiogeochem_nitrogenflux_type
+ use SoilBiogeochemNitrogenStateType , only : soilbiogeochem_nitrogenstate_type
+ use WaterStateBulkType , only : waterstatebulk_type
+ use SoilStatetype , only : soilstate_type
+ use SoilHydrologyType , only : soilhydrology_type
+ use ColumnType , only : col
+
+ !
+ implicit none
+ private
+ !
+ ! PUBLIC MEMBER FUNCTIONS
+ public SoilNitrogenMovement
+
+ character(len=*), parameter, private :: sourcefile = &
+ __FILE__
+ !------------------------------------------------------------------------
+
+ contains
+
+ !------------------------------------------------------------------------
+ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterstatebulk_inst, &
+ soilstate_inst, soilhydrology_inst, soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst)
+ !
+ ! implementation of the advection-diffusion algorithm in Patankar1980
+ !
+ ! This part of code is only designed for fast aqueous transport and leaching of inorganic nitrate
+ ! no sources and other sinks are included in this module
+ ! leaching flux is taken out of soil pool
+ ! Author: Jinmu Luo, Cornell EAS, Apr 11 2024
+ !
+ !USES:
+ use decompMod , only : bounds_type
+ use clm_varpar , only : nlevdecomp, nlevgrnd
+ use clm_time_manager , only : get_step_size_real, get_curr_date
+ use clm_varcon , only : zsoi, zisoi, dzsoi_decomp, mmh2o_to_m3h2o_per_m2
+ use ColumnType , only : col
+ use clm_varctl , only : use_bedrock
+ use TridiagonalMod , only : Tridiagonal
+ !ARGUMENTS:
+ type(bounds_type) , intent(in) :: bounds ! bounds
+ integer , intent(in) :: num_bgc_soilc ! number of soil columns in filter
+ integer , intent(in) :: filter_bgc_soilc(:) ! filter for soil columns
+ type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst
+ type(soilstate_type) , intent(in) :: soilstate_inst
+ type(soilhydrology_type) , intent(in) :: soilhydrology_inst
+ type(soilbiogeochem_nitrogenflux_type) , intent(inout) :: soilbiogeochem_nitrogenflux_inst
+ type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst
+
+ !LOCAL VARIABLES:
+ integer :: c,fc,j ! do loop indices
+ integer :: year, mon, day, tod
+ integer :: jtop(bounds%begc:bounds%endc) ! top level at each column
+ real(r8) :: dtime ! land model time step (sec)
+ real(r8) :: wafc, wafc2 ! Fraction of water that is liquid by mass
+ real(r8) :: dispersion_length = 0.1_r8 ! dispersion length (m), Jury et al., 1991
+ real(r8) :: theta, thetasat ! soil water and soil water at the saturation level
+ real(r8) :: no3_diffusivity_in_water = 1.7e-9_r8 ! Molecular diffusivity of NO3- in water, m2/s
+ real(r8) :: dissolve_frac = 1.0_r8 ! dissolve fraction
+ real(r8) :: flux_component_gridpoint_ahead ! A function in Patankar 1980, figure 5.6
+ real(r8) :: peclet_num ! Peclet number in Patankar 1980, foumula 5.18
+ real(r8) :: peclet_num_in, peclet_num_out ! temporary Peclet numbers
+ real(r8) :: qflx_in, qflx_out ! water fluxes same as qin, qout but in m3 H2O/m2/s
+ real(r8) :: dz_node(1:nlevdecomp+1) ! difference between nodes
+ real(r8) :: mass_old(bounds%begc:bounds%endc) ! Temporal column mass, g/m2
+ real(r8) :: mass_new(bounds%begc:bounds%endc) ! Temporal column mass, g/m2
+ real(r8) :: swliq(bounds%begc:bounds%endc,1:nlevdecomp) ! volumetric liquid soil water [m3/m3], hardwired to 1 for non-transport layers and layers below bedrock
+ real(r8) :: total_diffusivity(bounds%begc:bounds%endc,1:nlevdecomp+1) ! Total diffusivity
+ real(r8) :: a_tri(bounds%begc:bounds%endc,0:nlevdecomp+1) ! "a" vector for tridiagonal matrix
+ real(r8) :: b_tri(bounds%begc:bounds%endc,0:nlevdecomp+1) ! "b" vector for tridiagonal matrix
+ real(r8) :: c_tri(bounds%begc:bounds%endc,0:nlevdecomp+1) ! "c" vector for tridiagonal matrix
+ real(r8) :: r_tri(bounds%begc:bounds%endc,0:nlevdecomp+1) ! "r" vector for tridiagonal solution
+ real(r8) :: conc_trcr(bounds%begc:bounds%endc,0:nlevdecomp+1) ! temporary for concentration, g/m3H2O
+
+ ! set up the A function, table 5.2 in Patankar 1980 has multiple A function.
+ ! Notes: According to Table 5.2, here we use the "Power Law" version of the function
+ ! A is a dimensionless coefficient described in equation 5.37 of Patankar (1980)
+ ! The same identical function appears in CLM's SoilBiogeochemLittVertTranspMod.F90 as aaa(pe)
+ ! Patankar (1980) is posted here: https://github.com/ESCOMP/CTSM/pull/2992#discussion_r2294809728
+ flux_component_gridpoint_ahead(peclet_num) = max(0._r8, ( 1._r8 - 0.1_r8 * abs(peclet_num) )**5 )
+
+ associate(&
+ h2osoi_vol => waterstatebulk_inst%h2osoi_vol_col , & ! Input: [real(r8) (:,:) ] volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3]
+ h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] col liquid water (kg/m2)
+ h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] col ice lens (kg/m2)
+ watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity)
+ qout => soilhydrology_inst%qout_col , & ! Input: [real(r8) (:,:) ] soil water out of the bottom, mm h2o/s
+ qin => soilhydrology_inst%qin_col , & ! Input: [real(r8) (:,:) ] soil water into the bottom, mm h2o/s
+ smin_no3_vr => soilbiogeochem_nitrogenstate_inst%smin_no3_vr_col , & ! Inout: [real(r8) (:,:) ] soil nitrate concentration, gN/m3
+ smin_no3_leached_vr => soilbiogeochem_nitrogenflux_inst%smin_no3_leached_vr_col & ! Output: [real(r8) (:,:) ] rate of mineral NO3 leaching (gN/m3/s)
+ )
+
+ !Get the size of model time step
+ dtime = get_step_size_real()
+ call get_curr_date(year, mon, day, tod)
+
+ ! Preparing for the necessary parameters, like Delta z, Jtop, and the total diffusivity coefficients.
+ do j = 1, nlevdecomp + 1
+ if (j < nlevdecomp) then
+ dz_node(j) = zsoi(j+1) - zsoi(j)
+ end if
+ do fc = 1, num_bgc_soilc
+ c = filter_bgc_soilc(fc)
+ jtop(c) = 0
+ if (j < nlevdecomp) then
+ wafc = h2osoi_liq(c,j)/(h2osoi_liq(c,j) + h2osoi_ice(c,j))
+ wafc2 = h2osoi_liq(c,j+1)/(h2osoi_liq(c,j+1) + h2osoi_ice(c,j+1))
+ swliq(c,j) = wafc * h2osoi_vol(c,j)
+ swliq(c,j+1) = wafc2 * h2osoi_vol(c,j+1)
+ if (swliq(c,j) == 0._r8 .or. swliq(c,j+1) == 0._r8) then
+ theta = 0._r8
+ thetasat = 1._r8
+ else
+ theta = swliq(c,j) + dzsoi_decomp(j)/2 * (swliq(c,j+1) - swliq(c,j))/dz_node(j)
+ thetasat = watsat(c,j) + dzsoi_decomp(j)/2 * (watsat(c,j+1) - watsat(c,j))/dz_node(j)
+ end if
+ ! here we refer the j as the interface of j + zj/2
+ total_diffusivity(c,j) = theta * (no3_diffusivity_in_water * theta**(7/3) * thetasat**(-2))
+ total_diffusivity(c,j) = total_diffusivity(c,j) + dispersion_length * mmh2o_to_m3h2o_per_m2 * abs(qout(c,j))
+ else
+ !no gradient for the last layer
+ total_diffusivity(c,j) = total_diffusivity(c,j-1)
+ end if
+ end do ! Loop for columns
+ end do ! Loop for depths
+ dz_node(nlevdecomp) = dz_node(nlevdecomp-1)
+ dz_node(nlevdecomp+1) = dz_node(nlevdecomp)
+
+ ! Calculate the tridiagonal matrix, using Crank-Nicholson soluiton
+ ! dc/dt = (1-alpha)dF(t+1)/dx + alpha*F(t)/dx,
+ ! alpha=0
+ do j = 0, nlevdecomp + 1
+ do fc = 1, num_bgc_soilc
+ c = filter_bgc_soilc(fc)
+ if ( j==0 .or. j==nlevdecomp+1) then
+ !atmosphere and bottom layer, no concentration gradient here
+ conc_trcr(c,j) = 0._r8
+ a_tri(c,j) = 0._r8
+ b_tri(c,j) = 1._r8
+ c_tri(c,j) = 0._r8
+ r_tri(c,j) = 0._r8
+ elseif (swliq(c,j) == 0._r8 .or. j > col%nbedrock(c)) then
+ ! extremely dry condition and layers beneath the bedrock, no aqueous transport of nitrate
+ conc_trcr(c,j) = dissolve_frac * smin_no3_vr(c,j)
+ a_tri(c,j) = 0._r8
+ b_tri(c,j) = 1._r8
+ c_tri(c,j) = 0._r8
+ r_tri(c,j) = conc_trcr(c,j)
+ swliq(c,j) = 1.0_r8 ! change swliq into 1 to be used in the update session below
+ elseif ( j == 1) then
+ ! topmost soil layer, flux only interacts with the layer below it
+ conc_trcr(c,j) = dissolve_frac * smin_no3_vr(c,j)/swliq(c,j)
+ qflx_out = qout(c,j) * mmh2o_to_m3h2o_per_m2
+ peclet_num_out = qflx_out * dz_node(j) / total_diffusivity(c,j)
+ a_tri(c,j) = 0._r8
+ c_tri(c,j) = -total_diffusivity(c,j) / dz_node(j) * flux_component_gridpoint_ahead(peclet_num_out) - max(-qflx_out, 0._r8)
+ b_tri(c,j) = total_diffusivity(c,j) / dz_node(j) * flux_component_gridpoint_ahead(peclet_num_out) + max(qflx_out, 0._r8) + swliq(c,j) / dtime * dzsoi_decomp(j)
+ r_tri(c,j) = conc_trcr(c,j)/dtime*swliq(c,j)*dzsoi_decomp(j)
+ elseif ( j == col%nbedrock(c)) then
+ ! Assume the bottom layer concentration is always zero
+ ! This method count the loss at this layer as the leaching flux at the bottom
+ conc_trcr(c,j) = 0._r8
+ a_tri(c,j) = 0._r8
+ b_tri(c,j) = 1._r8
+ c_tri(c,j) = 0._r8
+ r_tri(c,j) = conc_trcr(c,j)/dtime*swliq(c,j)*dzsoi_decomp(j)
+ else
+ ! Active layers from second one to bedrock-1, concentration should be in gN/m3Water
+ conc_trcr(c,j) = dissolve_frac * smin_no3_vr(c,j)/swliq(c,j)
+ qflx_in = qin(c,j) * mmh2o_to_m3h2o_per_m2 ! mm H2O/s to m3 H2O/m2/s
+ qflx_out = qout(c,j) * mmh2o_to_m3h2o_per_m2
+ peclet_num_in = qflx_in * dz_node(j-1) / total_diffusivity(c,j-1)
+ peclet_num_out = qflx_out * dz_node(j) / total_diffusivity(c,j)
+ a_tri(c,j) = -total_diffusivity(c,j-1) / dz_node(j-1) * flux_component_gridpoint_ahead(peclet_num_in) - max(qflx_in, 0._r8)
+ c_tri(c,j) = -total_diffusivity(c,j) / dz_node(j) * flux_component_gridpoint_ahead(peclet_num_out) - max(-qflx_out, 0._r8)
+ b_tri(c,j) = total_diffusivity(c,j-1) / dz_node(j-1) * flux_component_gridpoint_ahead(peclet_num_in) + max(-qflx_in, 0._r8) + &
+ total_diffusivity(c,j) / dz_node(j) * flux_component_gridpoint_ahead(peclet_num_out) + max(qflx_out, 0._r8) + swliq(c,j) / dtime * dzsoi_decomp(j)
+ r_tri(c,j) = conc_trcr(c,j)/dtime*swliq(c,j)*dzsoi_decomp(j)
+ end if
+ end do ! Loop for columns
+ end do ! Loop for depths
+
+ ! solve the tridiagonal matrix
+ call Tridiagonal(bounds, 0, nlevdecomp+1, &
+ jtop(bounds%begc:bounds%endc), &
+ num_bgc_soilc, filter_bgc_soilc, &
+ a_tri(bounds%begc:bounds%endc, :), &
+ b_tri(bounds%begc:bounds%endc, :), &
+ c_tri(bounds%begc:bounds%endc, :), &
+ r_tri(bounds%begc:bounds%endc, :), &
+ conc_trcr(bounds%begc:bounds%endc,0:nlevdecomp+1))
+
+ ! Calculate the leaching flux
+ mass_old(bounds%begc:bounds%endc) = 0._r8
+ mass_new(bounds%begc:bounds%endc) = 0._r8
+ do j = 1, nlevdecomp
+ do fc = 1, num_bgc_soilc
+ c = filter_bgc_soilc(fc)
+ smin_no3_leached_vr(c,j) = 0._r8
+ mass_old(c) = mass_old(c) + smin_no3_vr(c,j)*dzsoi_decomp(j)
+ mass_new(c) = mass_new(c) + (smin_no3_vr(c,j) * (1._r8 - dissolve_frac) + conc_trcr(c,j) * swliq(c,j)) * dzsoi_decomp(j)
+ end do
+ end do
+
+ do fc = 1, num_bgc_soilc
+ c = filter_bgc_soilc(fc)
+ ! g/m3/sec, leaching mass is at the layer above the bedrock
+ smin_no3_leached_vr(c, col%nbedrock(c)) = max(0._r8, (mass_old(c) - mass_new(c))/dzsoi_decomp(col%nbedrock(c))/dtime)
+ end do
+
+ ! Update the pools of interest
+ do j = 1, nlevdecomp
+ do fc = 1, num_bgc_soilc
+ c = filter_bgc_soilc(fc)
+ smin_no3_vr(c,j) = smin_no3_vr(c,j) - smin_no3_vr(c,j) * dissolve_frac + conc_trcr(c,j) * swliq(c,j)
+ ! Return this leaching flux back to smin_no3 pool, and update will be finished in CNNStateUpdate3Mod
+ if( j == col%nbedrock(c) ) then
+ smin_no3_vr(c,j) = smin_no3_vr(c,j) + smin_no3_leached_vr(c,j)*dtime
+ end if
+ end do ! loop for columns
+ end do !loop for depths
+
+ end associate
+
+
+ end subroutine SoilNitrogenMovement
+
+
+end module SoilNitrogenMovementMod