From f0199d7f04f3e148f3c1e87715ffe97bf21d3371 Mon Sep 17 00:00:00 2001 From: Jinmu Luo Date: Thu, 6 Mar 2025 12:41:00 -0700 Subject: [PATCH 01/29] Merging vertical N leaching with CTSM5.3.027 --- src/biogeochem/CNDriverMod.F90 | 17 +- src/biogeochem/CNVegetationFacade.F90 | 6 +- src/biogeophys/SoilHydrologyType.F90 | 16 +- src/biogeophys/SoilWaterMovementMod.F90 | 19 +- src/main/clm_driver.F90 | 2 +- .../SoilBiogeochemNLeachingMod.F90 | 40 +-- .../SoilNitrogenMovementMod.F90 | 237 ++++++++++++++++++ 7 files changed, 309 insertions(+), 28 deletions(-) create mode 100644 src/soilbiogeochem/SoilNitrogenMovementMod.F90 diff --git a/src/biogeochem/CNDriverMod.F90 b/src/biogeochem/CNDriverMod.F90 index 0274fcc87f..b40ebff792 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 @@ -1016,7 +1017,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. @@ -1031,6 +1032,7 @@ 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_nvmovement ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds @@ -1048,6 +1050,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 @@ -1065,8 +1068,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, leaching[only if use_nvmoment=.false.]) 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 dd20ce50a2..817adba03e 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 @@ -1061,7 +1062,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 @@ -1100,6 +1101,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: @@ -1122,7 +1124,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..5f12c23d07 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) + + 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 of soil layer', & + ptr_col=this%qin_col) + end subroutine InitHistory !----------------------------------------------------------------------- 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 f93143d9e3..fde7f1b5a2 100644 --- a/src/main/clm_driver.F90 +++ b/src/main/clm_driver.F90 @@ -1116,7 +1116,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/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 b/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 index a646feb1d7..e5b85db612 100644 --- a/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 +++ b/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 @@ -83,6 +83,7 @@ subroutine SoilBiogeochemNLeaching(bounds, num_bgc_soilc, filter_bgc_soilc, & ! !USES: use clm_varpar , only : nlevdecomp, nlevsoi use clm_time_manager , only : get_step_size_real + use SoilNitrogenMovementMod , only: use_nvmovement ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds @@ -212,16 +213,20 @@ 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 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,16 @@ 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 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..a626c17c4f --- /dev/null +++ b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 @@ -0,0 +1,237 @@ +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 + + ! !!!!!!! Author note, I currently using a flag to control this leaching module !!!!!! + logical, parameter, public :: use_nvmovement = .TRUE. + + contains + + + !# main code here + 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 + 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 ! H2Oliq/(H2Oliq + H2Oice) + real(r8) :: Ldis = 0.1_r8 ! dispersion length (m), Jury et al., 1991 + real(r8) :: pcl1, pcl2, qflx1, qflx2, theta, thetasat, nerror + real(r8) :: Dw = 1.7e-9_r8 ! Molecular diffusivity of NO3- in water, m2/s + real(r8) :: disp = 1.0_r8 ! dissove percentage + real(r8) :: afunc ! A function in Patankar 1980, figure 5.6 + real(r8) :: pcl ! Peclet number in Patankar 1980, foumula 5.18 + 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], 1 for nonewater + real(r8) :: total_D(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. + afunc(pcl) = max(0._r8, ( 1._r8 - 0.1_r8 * abs(pcl) )**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_D(c,j) = theta*(Dw * theta**(7/3) * thetasat**(-2)) + total_D(c,j) = total_D(c,j) + Ldis*0.001*abs(qout(c,j)) + else + !no gradient for the last layer + total_D(c,j) = total_D(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) = disp * 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) = disp * smin_no3_vr(c,j)/swliq(c,j) + qflx1 = qin(c,j) * 0.001_r8 ! mmH2O/sec to m3H2O/m2/sec + qflx2 = qout(c,j) * 0.001_r8 + pcl1 = qflx1*dz_node(j)/total_D(c,j) + pcl2 = qflx2*dz_node(j)/total_D(c,j) + a_tri(c,j) = 0._r8 + c_tri(c,j) = -total_D(c,j)/dz_node(j) * afunc(pcl2) - max(-qflx2, 0._r8) + b_tri(c,j) = total_D(c,j)/dz_node(j) * afunc(pcl2) + max(qflx2, 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; c_tri(c,j) = 0._r8; b_tri(c,j) = 1._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) = disp * smin_no3_vr(c,j)/swliq(c,j) + qflx1 = qin(c,j) * 0.001_r8 ! mmH2O/sec to m3H2O/m2/sec + qflx2 = qout(c,j) * 0.001_r8 + pcl1 = qflx1*dz_node(j-1)/total_D(c,j-1) + pcl2 = qflx2*dz_node(j)/total_D(c,j) + a_tri(c,j) = -total_D(c,j-1)/dz_node(j-1) * afunc(pcl1) - max(qflx1, 0._r8) + c_tri(c,j) = -total_D(c,j)/dz_node(j) * afunc(pcl2) - max(-qflx2, 0._r8) + b_tri(c,j) = total_D(c,j-1)/dz_node(j-1) * afunc(pcl1) + max(-qflx1, 0._r8) + & + total_D(c,j)/dz_node(j) * afunc(pcl2) + max(qflx2, 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-disp) + 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) + nerror = mass_old(c) - mass_new(c) + ! 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)*disp + 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 From 0b2e352360c835c960cc809e9df81135fca05668 Mon Sep 17 00:00:00 2001 From: Jinmu Luo Date: Tue, 6 May 2025 14:42:08 -0600 Subject: [PATCH 02/29] change the address in spin-up check tool --- tools/contrib/SpinupStability_BGC_v10.ncl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tools/contrib/SpinupStability_BGC_v10.ncl b/tools/contrib/SpinupStability_BGC_v10.ncl index 5ed7516455..1f4be8505e 100644 --- a/tools/contrib/SpinupStability_BGC_v10.ncl +++ b/tools/contrib/SpinupStability_BGC_v10.ncl @@ -35,8 +35,8 @@ begin ;=======================================================================; ; GLOBAL EXAMPLE - caseid = "clm50_release-clm5.0.15_2deg_GSWP3V1_1850AD" - username = "oleson" + caseid = "ctsm53027_LEACHING_pAD" + username = "jinmuluo" annual_hist = True region = "Global" ; Global, Arctic, or SPT (single point) subper = 20 ; Subsampling period in years @@ -54,7 +54,7 @@ begin ; do_plot = False ;=======================================================================; - data_dir = "/glade/scratch/"+username+"/archive/"+caseid+"/lnd/hist/" + data_dir = "/glade/derecho/scratch/"+username+"/archive/"+caseid+"/lnd/hist/" if ( systemfunc("test -d "+data_dir+"; echo $?" ) .ne. 0 )then print( "Input directory does not exist or not found: "+data_dir ); print( "Make sure username and caseid and base directory is set correctly" ) From c417cbb524e2d82657f3b6e6bf4658db2168ec25 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Mon, 28 Jul 2025 11:14:14 -0600 Subject: [PATCH 03/29] Revert SpinupStability_BGC_v10.ncl changes as not pertinent to this PR --- tools/contrib/SpinupStability_BGC_v10.ncl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tools/contrib/SpinupStability_BGC_v10.ncl b/tools/contrib/SpinupStability_BGC_v10.ncl index 1f4be8505e..5ed7516455 100644 --- a/tools/contrib/SpinupStability_BGC_v10.ncl +++ b/tools/contrib/SpinupStability_BGC_v10.ncl @@ -35,8 +35,8 @@ begin ;=======================================================================; ; GLOBAL EXAMPLE - caseid = "ctsm53027_LEACHING_pAD" - username = "jinmuluo" + caseid = "clm50_release-clm5.0.15_2deg_GSWP3V1_1850AD" + username = "oleson" annual_hist = True region = "Global" ; Global, Arctic, or SPT (single point) subper = 20 ; Subsampling period in years @@ -54,7 +54,7 @@ begin ; do_plot = False ;=======================================================================; - data_dir = "/glade/derecho/scratch/"+username+"/archive/"+caseid+"/lnd/hist/" + data_dir = "/glade/scratch/"+username+"/archive/"+caseid+"/lnd/hist/" if ( systemfunc("test -d "+data_dir+"; echo $?" ) .ne. 0 )then print( "Input directory does not exist or not found: "+data_dir ); print( "Make sure username and caseid and base directory is set correctly" ) From 614487dcb12efa960e325fd792cc15495c476296 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Mon, 4 Aug 2025 18:14:57 -0600 Subject: [PATCH 04/29] Make use_nvmovement a namelist parameter and make it default .false. --- bld/CLMBuildNamelist.pm | 7 +++++++ bld/namelist_files/namelist_defaults_ctsm.xml | 1 + bld/namelist_files/namelist_definition_ctsm.xml | 6 ++++++ src/biogeochem/CNDriverMod.F90 | 3 ++- src/main/clm_varctl.F90 | 1 + src/main/controlMod.F90 | 4 +++- src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 | 3 +-- src/soilbiogeochem/SoilNitrogenMovementMod.F90 | 10 +++++----- 8 files changed, 26 insertions(+), 9 deletions(-) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 4cd02532f8..898b5e16ca 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -1005,6 +1005,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 and check that it is not true at the same time as use_soil_matrixcn + $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->get_value('use_soil_matrixcn')) ) { + $log->fatal_error("$var can only be on with use_soil_matrixcn = .false."); + } } # end bgc diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index a435bb3b4a..222f1a822d 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -2429,6 +2429,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 03c9ba420e..9316947a56 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -991,6 +991,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_soil_matrixcn is true + diff --git a/src/biogeochem/CNDriverMod.F90 b/src/biogeochem/CNDriverMod.F90 index 82b606e70d..e442708c7c 100644 --- a/src/biogeochem/CNDriverMod.F90 +++ b/src/biogeochem/CNDriverMod.F90 @@ -1029,7 +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_nvmovement + use SoilNitrogenMovementMod , only: SoilNitrogenMovement + use clm_varctl, only : use_nvmovement ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index 41978ae695..6a34714896 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -521,6 +521,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 6d363a9a6e..20688ca61f 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -318,7 +318,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/ & @@ -729,6 +729,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) @@ -1026,6 +1027,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 e5b85db612..604a63d3c0 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 @@ -83,7 +83,6 @@ subroutine SoilBiogeochemNLeaching(bounds, num_bgc_soilc, filter_bgc_soilc, & ! !USES: use clm_varpar , only : nlevdecomp, nlevsoi use clm_time_manager , only : get_step_size_real - use SoilNitrogenMovementMod , only: use_nvmovement ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds diff --git a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 index a626c17c4f..68649dedea 100644 --- a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 +++ b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 @@ -30,13 +30,13 @@ module SoilNitrogenMovementMod ! PUBLIC MEMBER FUNCTIONS public SoilNitrogenMovement - ! !!!!!!! Author note, I currently using a flag to control this leaching module !!!!!! - logical, parameter, public :: use_nvmovement = .TRUE. - - contains + character(len=*), parameter, private :: sourcefile = & + __FILE__ + !------------------------------------------------------------------------ + contains - !# main code here + !------------------------------------------------------------------------ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterstatebulk_inst, & soilstate_inst, soilhydrology_inst, soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst) ! From 8d34bc76af6e3554ccb9cb2cca6e98861406c6ec Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Tue, 5 Aug 2025 16:13:11 -0600 Subject: [PATCH 05/29] Change comment to make consistent with the code --- src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 b/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 index 604a63d3c0..57b919e447 100644 --- a/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 +++ b/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 @@ -212,7 +212,7 @@ 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 - ! Evaluating the leaching flux in SoilNitrogenMovementMod, if use_nvmovement is true + ! Evaluating the leaching flux in SoilNitrogenMovementMod, if use_nvmovement is not true if (.not. use_nvmovement) then ! ! calculate the N leaching flux as a function of the dissolved @@ -246,7 +246,7 @@ 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) - ! Evaluating the leaching flux in SoilNitrogenMovementMod, if use_nvmovement is true + ! Evaluating the leaching flux in SoilNitrogenMovementMod, if use_nvmovement is not true if (.not. use_nvmovement) then ! limit the flux based on current smin_no3 state ! only let at most the assumed soluble fraction From d175bb7d7a8c07266dfb3f8e276cd1ededf237f6 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Tue, 5 Aug 2025 17:07:52 -0600 Subject: [PATCH 06/29] QIN & QOUT needed for b4b restarts and not needed in history --- src/biogeophys/SoilHydrologyType.F90 | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/src/biogeophys/SoilHydrologyType.F90 b/src/biogeophys/SoilHydrologyType.F90 index 5f12c23d07..e03390fac9 100644 --- a/src/biogeophys/SoilHydrologyType.F90 +++ b/src/biogeophys/SoilHydrologyType.F90 @@ -144,8 +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 + allocate(this%qout_col (begc:endc,nlevgrnd)) ; this%qout_col (:,:) = nan + allocate(this%qin_col (begc:endc,nlevgrnd)) ; this%qin_col (:,:) = nan end subroutine InitAllocate @@ -199,12 +199,12 @@ subroutine InitHistory(this, bounds, use_aquifer_layer) 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) + 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 of soil layer', & - ptr_col=this%qin_col) + avgflag='A', long_name='flux of water into soil layer', & + ptr_col=this%qin_col, default='inactive') end subroutine InitHistory @@ -301,8 +301,8 @@ 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 + real(r8), pointer :: ptr2d(:,:) ! temp. pointers for slicing larger arrays !----------------------------------------------------------------------- call restartvar(ncid=ncid, flag=flag, varname='FROST_TABLE', xtype=ncd_double, & @@ -326,6 +326,20 @@ subroutine Restart(this, bounds, ncid, flag) this%zwt_perched_col(bounds%begc:bounds%endc) = col%zi(bounds%begc:bounds%endc,nlevsoi) end if + ptr2d => this%qout_col + call restartvar(ncid=ncid, flag=flag, varname='QOUT', xtype=ncd_double, & + dim1name='column', dim2name='levgrnd', switchdim=.true., & + long_name='flux of water out of soil layer', units='mm h2o/s', fill_value=spval, & + scale_by_thickness=.false., & + interpinic_flag='skip', readvar=readvar, data=ptr2d) + + ptr2d => this%qin_col + call restartvar(ncid=ncid, flag=flag, varname='QIN', xtype=ncd_double, & + dim1name='column', dim2name='levgrnd', switchdim=.true., & + long_name='flux of water into soil layer', units='mm h2o/s', fill_value=spval, & + scale_by_thickness=.false., & + interpinic_flag='skip', readvar=readvar, data=ptr2d) + end subroutine Restart !----------------------------------------------------------------------- From 9e62984e83e75b9b3eb67b20d2bcc16532487ad9 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Tue, 5 Aug 2025 17:21:10 -0600 Subject: [PATCH 07/29] Correct comment update committed earlier today --- src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 b/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 index 57b919e447..f0a42b379c 100644 --- a/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 +++ b/src/soilbiogeochem/SoilBiogeochemNLeachingMod.F90 @@ -212,7 +212,8 @@ 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 - ! Evaluating the leaching flux in SoilNitrogenMovementMod, if use_nvmovement is not true + ! 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 @@ -246,7 +247,8 @@ 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) - ! Evaluating the leaching flux in SoilNitrogenMovementMod, if use_nvmovement is not true + ! 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 From 5b9c0e25096306553165186fcfc08432834afe04 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Tue, 5 Aug 2025 17:35:37 -0600 Subject: [PATCH 08/29] Reverting mistake: QIN & QOUT do not help with bit-for-bit restart --- src/biogeophys/SoilHydrologyType.F90 | 19 ++----------------- 1 file changed, 2 insertions(+), 17 deletions(-) diff --git a/src/biogeophys/SoilHydrologyType.F90 b/src/biogeophys/SoilHydrologyType.F90 index e03390fac9..138044154f 100644 --- a/src/biogeophys/SoilHydrologyType.F90 +++ b/src/biogeophys/SoilHydrologyType.F90 @@ -144,8 +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,nlevgrnd)) ; this%qout_col (:,:) = nan - allocate(this%qin_col (begc:endc,nlevgrnd)) ; this%qin_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 @@ -302,7 +302,6 @@ subroutine Restart(this, bounds, ncid, flag) ! ! !LOCAL VARIABLES: logical :: readvar ! determine if variable is on initial file - real(r8), pointer :: ptr2d(:,:) ! temp. pointers for slicing larger arrays !----------------------------------------------------------------------- call restartvar(ncid=ncid, flag=flag, varname='FROST_TABLE', xtype=ncd_double, & @@ -326,20 +325,6 @@ subroutine Restart(this, bounds, ncid, flag) this%zwt_perched_col(bounds%begc:bounds%endc) = col%zi(bounds%begc:bounds%endc,nlevsoi) end if - ptr2d => this%qout_col - call restartvar(ncid=ncid, flag=flag, varname='QOUT', xtype=ncd_double, & - dim1name='column', dim2name='levgrnd', switchdim=.true., & - long_name='flux of water out of soil layer', units='mm h2o/s', fill_value=spval, & - scale_by_thickness=.false., & - interpinic_flag='skip', readvar=readvar, data=ptr2d) - - ptr2d => this%qin_col - call restartvar(ncid=ncid, flag=flag, varname='QIN', xtype=ncd_double, & - dim1name='column', dim2name='levgrnd', switchdim=.true., & - long_name='flux of water into soil layer', units='mm h2o/s', fill_value=spval, & - scale_by_thickness=.false., & - interpinic_flag='skip', readvar=readvar, data=ptr2d) - end subroutine Restart !----------------------------------------------------------------------- From 93ba4bb359fee24bb88cabd627cb4ed6458bbadd Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Wed, 6 Aug 2025 14:29:59 -0600 Subject: [PATCH 09/29] New testdefs for aux_clm test with use_nvmovement = .true. --- cime_config/testdefs/testlist_clm.xml | 9 +++++++++ .../testmods_dirs/clm/Nvmovement/include_user_mods | 1 + .../testdefs/testmods_dirs/clm/Nvmovement/user_nl_clm | 1 + 3 files changed, 11 insertions(+) create mode 100644 cime_config/testdefs/testmods_dirs/clm/Nvmovement/include_user_mods create mode 100644 cime_config/testdefs/testmods_dirs/clm/Nvmovement/user_nl_clm diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index 478dc59bdd..d139ced9c3 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -592,6 +592,15 @@ + + + + + + + + + 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. From 8ffd2714ab1274c8f65ffd9482d9e0747ff5d9ae Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Wed, 6 Aug 2025 15:58:12 -0600 Subject: [PATCH 10/29] Remove namelist abort when both use_nvmovement and matrixcn are .true. --- bld/CLMBuildNamelist.pm | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 898b5e16ca..8ca4ef7d93 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -1006,12 +1006,9 @@ sub setup_cmdl_bgc { $log->fatal_error("$var can only be on with CENTURYKoven2013 soil decomposition"); } - # Set use_nvmovement and check that it is not true at the same time as use_soil_matrixcn + # 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->get_value('use_soil_matrixcn')) ) { - $log->fatal_error("$var can only be on with use_soil_matrixcn = .false."); - } } # end bgc From 275dc23b770c2eb7d20ae143659802d01d7aa8f0 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Wed, 6 Aug 2025 16:00:46 -0600 Subject: [PATCH 11/29] Update testlist_clm.xml for Nvmovement tests with and without matrixcn --- cime_config/testdefs/testlist_clm.xml | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index d139ced9c3..ce240131c9 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -592,13 +592,22 @@ - + + + + + + + + + + - + From f3b438135643eaac50c3319c2d86b3412c62bef8 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Wed, 6 Aug 2025 16:57:16 -0600 Subject: [PATCH 12/29] Draft ChangeLog/ChangeSum --- doc/ChangeLog | 91 +++++++++++++++++++++++++++++++++++++++++++++++++++ doc/ChangeSum | 3 ++ 2 files changed, 94 insertions(+) diff --git a/doc/ChangeLog b/doc/ChangeLog index 6b19a2747b..a3b30eb801 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,95 @@ =============================================================== +Tag name: ctsm5.3.0xx +Originator(s): jinmuluo (Jinmu Luo, Cornell University) +Date: Wed 06 Aug 2025 04:19:47 PM MDT +One-line Summary: New vertical movement scheme for soil nitrate + +Purpose and description of changes +---------------------------------- + + This tag introduces a new vertical movement scheme for soil nitrate but sets 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: + None + +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 - + + regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing): + + derecho ----- OK + izumi ------- + + fates tests: (give name of baseline if different from CTSM tagname, normally fates baselines are fates--) + derecho ----- + izumi ------- + + any other testing (give details below): + + ctsm_sci + derecho ---- + +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.067 +=============================================================== +=============================================================== +Tag name: ctsm5.3.066 +=============================================================== +=============================================================== Tag name: ctsm5.3.065 Originator(s): erik (Erik Kluzek,UCAR/TSS,303-497-1326) Date: Mon 28 Jul 2025 09:30:35 AM MDT diff --git a/doc/ChangeSum b/doc/ChangeSum index 20b0a8f534..f657354b3a 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,5 +1,8 @@ Tag Who Date Summary ============================================================================================================================ + ctsm5.3.0xx jinmuluo xx/xx/2025 New vertical movement scheme for soil nitrate + ctsm5.3.067 / /2025 + ctsm5.3.066 / /2025 ctsm5.3.065 erik 07/28/2025 Merge b4b-dev to master ctsm5.3.064 slevis 07/24/2025 Add time dimension to 1d_wt fields in transient runs ctsm5.3.063 samrabin 07/10/2025 Merge b4b-dev to master From 51d6f02188bfeea525999cbfed6742a8ff83576b Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Thu, 7 Aug 2025 15:07:32 -0600 Subject: [PATCH 13/29] Prevent use_nvmovement = .true. when use_nitrif_denitrif = .false. --- bld/CLMBuildNamelist.pm | 3 +++ 1 file changed, 3 insertions(+) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 8ca4ef7d93..eea1572075 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -1009,6 +1009,9 @@ sub setup_cmdl_bgc { # 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 From 9a5285e524042b62d8f8ccab8790b3cd9d25246a Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Tue, 19 Aug 2025 18:11:44 -0600 Subject: [PATCH 14/29] Rename testmod directory Nvmovement --> nvmovement --- .../clm/{Nvmovement => nvmovement}/include_user_mods | 0 .../testmods_dirs/clm/{Nvmovement => nvmovement}/user_nl_clm | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename cime_config/testdefs/testmods_dirs/clm/{Nvmovement => nvmovement}/include_user_mods (100%) rename cime_config/testdefs/testmods_dirs/clm/{Nvmovement => nvmovement}/user_nl_clm (100%) diff --git a/cime_config/testdefs/testmods_dirs/clm/Nvmovement/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/nvmovement/include_user_mods similarity index 100% rename from cime_config/testdefs/testmods_dirs/clm/Nvmovement/include_user_mods rename to cime_config/testdefs/testmods_dirs/clm/nvmovement/include_user_mods diff --git a/cime_config/testdefs/testmods_dirs/clm/Nvmovement/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/nvmovement/user_nl_clm similarity index 100% rename from cime_config/testdefs/testmods_dirs/clm/Nvmovement/user_nl_clm rename to cime_config/testdefs/testmods_dirs/clm/nvmovement/user_nl_clm From 9603a08eb3a8c6ddab7c206d1906bd6f128a80ef Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Tue, 19 Aug 2025 18:15:37 -0600 Subject: [PATCH 15/29] Change Nvmovement to nvmovement in test-suite tests --- cime_config/testdefs/testlist_clm.xml | 4 ++-- doc/ChangeLog | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index ce240131c9..9ee275804a 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -592,7 +592,7 @@ - + @@ -601,7 +601,7 @@ - + diff --git a/doc/ChangeLog b/doc/ChangeLog index a3b30eb801..6056a5ab19 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -45,9 +45,9 @@ Changes to documentation: 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 + 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: ---------------- From 8646005b3e2ce18e0f5157aa063b4bbe693b0c2b Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Tue, 19 Aug 2025 18:20:55 -0600 Subject: [PATCH 16/29] Update draft ChangeLog/Sum --- doc/ChangeLog | 10 ++-------- doc/ChangeSum | 4 +--- 2 files changed, 3 insertions(+), 11 deletions(-) diff --git a/doc/ChangeLog b/doc/ChangeLog index 6056a5ab19..2eeabdc8af 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,7 +1,7 @@ =============================================================== -Tag name: ctsm5.3.0xx +Tag name: ctsm5.3.071 Originator(s): jinmuluo (Jinmu Luo, Cornell University) -Date: Wed 06 Aug 2025 04:19:47 PM MDT +Date: Tue 19 Aug 2025 06:19:28 PM MDT One-line Summary: New vertical movement scheme for soil nitrate Purpose and description of changes @@ -82,12 +82,6 @@ Other details Pull Requests that document the changes (include PR ids): https://github.com/ESCOMP/ctsm/pull/2992 -=============================================================== -=============================================================== -Tag name: ctsm5.3.067 -=============================================================== -=============================================================== -Tag name: ctsm5.3.066 =============================================================== =============================================================== Tag name: ctsm5.3.065 diff --git a/doc/ChangeSum b/doc/ChangeSum index f657354b3a..fdc5a0890a 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,8 +1,6 @@ Tag Who Date Summary ============================================================================================================================ - ctsm5.3.0xx jinmuluo xx/xx/2025 New vertical movement scheme for soil nitrate - ctsm5.3.067 / /2025 - ctsm5.3.066 / /2025 + ctsm5.3.071 jinmuluo xx/xx/2025 New vertical movement scheme for soil nitrate ctsm5.3.065 erik 07/28/2025 Merge b4b-dev to master ctsm5.3.064 slevis 07/24/2025 Add time dimension to 1d_wt fields in transient runs ctsm5.3.063 samrabin 07/10/2025 Merge b4b-dev to master From dff95b97ad15d0fa111419adb78de59522f9126d Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Tue, 19 Aug 2025 18:25:04 -0600 Subject: [PATCH 17/29] Update comment as recommended in code review --- src/biogeochem/CNDriverMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/biogeochem/CNDriverMod.F90 b/src/biogeochem/CNDriverMod.F90 index e442708c7c..d14e619b03 100644 --- a/src/biogeochem/CNDriverMod.F90 +++ b/src/biogeochem/CNDriverMod.F90 @@ -1075,7 +1075,7 @@ subroutine CNDriverLeaching(bounds, & call t_stopf('SoilNitrogenMovementMod') end if - ! Mineral nitrogen dynamics (deposition, fixation, leaching[only if use_nvmoment=.false.]) + ! Mineral nitrogen dynamics: deposition, fixation. If use_nvmoment false, also leaching. call t_startf('SoilBiogeochemNLeaching') call SoilBiogeochemNLeaching(bounds, num_bgc_soilc, filter_bgc_soilc, & From 04d32d329d560693b0d7124bf206685034ceec0d Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Tue, 19 Aug 2025 18:30:00 -0600 Subject: [PATCH 18/29] Update/correct the nvmovement tests as suggested in the code review --- cime_config/testdefs/testlist_clm.xml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index 9ee275804a..be06600162 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -592,16 +592,16 @@ - + - + - + From 0cebccd151b3981af34739084c2553b0d8be4b6d Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Wed, 20 Aug 2025 15:18:37 -0600 Subject: [PATCH 19/29] Change h2o to H2O in comments for consistency with existing examples --- src/biogeophys/SoilHydrologyType.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/biogeophys/SoilHydrologyType.F90 b/src/biogeophys/SoilHydrologyType.F90 index 138044154f..9be9073bc0 100644 --- a/src/biogeophys/SoilHydrologyType.F90 +++ b/src/biogeophys/SoilHydrologyType.F90 @@ -51,8 +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] + 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 @@ -197,12 +197,12 @@ subroutine InitHistory(this, bounds, use_aquifer_layer) 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', & + 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', & + 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') From b84826c901b32bc420cccc6f79cc368252bf264e Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Wed, 20 Aug 2025 15:22:26 -0600 Subject: [PATCH 20/29] Correct dissove to dissolve in comment --- src/soilbiogeochem/SoilNitrogenMovementMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 index 68649dedea..d1051fba91 100644 --- a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 +++ b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 @@ -74,7 +74,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst real(r8) :: Ldis = 0.1_r8 ! dispersion length (m), Jury et al., 1991 real(r8) :: pcl1, pcl2, qflx1, qflx2, theta, thetasat, nerror real(r8) :: Dw = 1.7e-9_r8 ! Molecular diffusivity of NO3- in water, m2/s - real(r8) :: disp = 1.0_r8 ! dissove percentage + real(r8) :: disp = 1.0_r8 ! dissolve percentage real(r8) :: afunc ! A function in Patankar 1980, figure 5.6 real(r8) :: pcl ! Peclet number in Patankar 1980, foumula 5.18 real(r8) :: dz_node(1:nlevdecomp+1) ! difference between nodes From 1a31282a7169468209ada5d3cfb766e774dbc482 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Wed, 20 Aug 2025 15:40:06 -0600 Subject: [PATCH 21/29] Split variable assignments to one line per assignment --- src/soilbiogeochem/SoilNitrogenMovementMod.F90 | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 index d1051fba91..d38d8a33c7 100644 --- a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 +++ b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 @@ -120,7 +120,8 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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 + 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) @@ -146,12 +147,16 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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 + 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) = disp * smin_no3_vr(c,j) - a_tri(c,j) = 0._r8; b_tri(c,j) = 1._r8; c_tri(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) swliq(c,j) = 1.0_r8 ! change swliq into 1 to be used in the update session below elseif ( j == 1) then @@ -169,7 +174,9 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst ! 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; c_tri(c,j) = 0._r8; b_tri(c,j) = 1._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 From 57aeebf00b8dcefdec8a4811b6797499f2e20f9d Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Wed, 20 Aug 2025 16:07:18 -0600 Subject: [PATCH 22/29] Replace 0.001 with mmh2o_to_m3h2o_per_m2 --- src/main/clm_varcon.F90 | 1 + src/soilbiogeochem/SoilNitrogenMovementMod.F90 | 12 ++++++------ 2 files changed, 7 insertions(+), 6 deletions(-) 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/soilbiogeochem/SoilNitrogenMovementMod.F90 b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 index d38d8a33c7..170f4d47a6 100644 --- a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 +++ b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 @@ -51,7 +51,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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 + 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 @@ -128,7 +128,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst end if ! here we refer the j as the interface of j + zj/2 total_D(c,j) = theta*(Dw * theta**(7/3) * thetasat**(-2)) - total_D(c,j) = total_D(c,j) + Ldis*0.001*abs(qout(c,j)) + total_D(c,j) = total_D(c,j) + Ldis * mmh2o_to_m3h2o_per_m2 * abs(qout(c,j)) else !no gradient for the last layer total_D(c,j) = total_D(c,j-1) @@ -162,8 +162,8 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst elseif ( j == 1) then ! topmost soil layer, flux only interacts with the layer below it conc_trcr(c,j) = disp * smin_no3_vr(c,j)/swliq(c,j) - qflx1 = qin(c,j) * 0.001_r8 ! mmH2O/sec to m3H2O/m2/sec - qflx2 = qout(c,j) * 0.001_r8 + qflx1 = qin(c,j) * mmh2o_to_m3h2o_per_m2 ! mm H2O/s to m3 H2O/m2/s + qflx2 = qout(c,j) * mmh2o_to_m3h2o_per_m2 pcl1 = qflx1*dz_node(j)/total_D(c,j) pcl2 = qflx2*dz_node(j)/total_D(c,j) a_tri(c,j) = 0._r8 @@ -181,8 +181,8 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst else ! Active layers from second one to bedrock-1, concentration should be in gN/m3Water conc_trcr(c,j) = disp * smin_no3_vr(c,j)/swliq(c,j) - qflx1 = qin(c,j) * 0.001_r8 ! mmH2O/sec to m3H2O/m2/sec - qflx2 = qout(c,j) * 0.001_r8 + qflx1 = qin(c,j) * mmh2o_to_m3h2o_per_m2 ! mm H2O/s to m3 H2O/m2/s + qflx2 = qout(c,j) * mmh2o_to_m3h2o_per_m2 pcl1 = qflx1*dz_node(j-1)/total_D(c,j-1) pcl2 = qflx2*dz_node(j)/total_D(c,j) a_tri(c,j) = -total_D(c,j-1)/dz_node(j-1) * afunc(pcl1) - max(qflx1, 0._r8) From 169523b4b281cd549ca29d63c04974e00fc6d374 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Wed, 20 Aug 2025 16:16:27 -0600 Subject: [PATCH 23/29] Make comment clearer for wafc, wafc2 --- src/soilbiogeochem/SoilNitrogenMovementMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 index 170f4d47a6..ff3f7556e1 100644 --- a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 +++ b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 @@ -70,7 +70,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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 ! H2Oliq/(H2Oliq + H2Oice) + real(r8) :: wafc, wafc2 ! Fraction of water that is liquid by mass real(r8) :: Ldis = 0.1_r8 ! dispersion length (m), Jury et al., 1991 real(r8) :: pcl1, pcl2, qflx1, qflx2, theta, thetasat, nerror real(r8) :: Dw = 1.7e-9_r8 ! Molecular diffusivity of NO3- in water, m2/s From 87ed1ea144149b4b8df9081a3edef395235533b9 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Fri, 22 Aug 2025 14:34:06 -0600 Subject: [PATCH 24/29] Fix comment about use_nvmovement in namelist_definition_ctsm.xml --- bld/namelist_files/namelist_definition_ctsm.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 9316947a56..bbece0a7d4 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -995,7 +995,7 @@ How LUNA and Photosynthesis (if needed) will get Leaf nitrogen content group="clm_inparm" value=".false."> 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_soil_matrixcn is true + use_nvmovement cannot be true while use_nitrif_denitrif is false Date: Fri, 22 Aug 2025 15:51:11 -0600 Subject: [PATCH 25/29] Add/correct a couple of comments in the code --- src/soilbiogeochem/SoilNitrogenMovementMod.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 index ff3f7556e1..b3094e5549 100644 --- a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 +++ b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 @@ -80,7 +80,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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], 1 for nonewater + 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_D(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 @@ -89,6 +89,10 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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 afunc(pcl) = max(0._r8, ( 1._r8 - 0.1_r8 * abs(pcl) )**5 ) associate(& From 7d946d94f802022b38be900f99d0b5af28efce55 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Fri, 22 Aug 2025 15:54:46 -0600 Subject: [PATCH 26/29] Replace disp (dissolve percentage) with dissolve_frac --- src/soilbiogeochem/SoilNitrogenMovementMod.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 index b3094e5549..ec263f793c 100644 --- a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 +++ b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 @@ -74,7 +74,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst real(r8) :: Ldis = 0.1_r8 ! dispersion length (m), Jury et al., 1991 real(r8) :: pcl1, pcl2, qflx1, qflx2, theta, thetasat, nerror real(r8) :: Dw = 1.7e-9_r8 ! Molecular diffusivity of NO3- in water, m2/s - real(r8) :: disp = 1.0_r8 ! dissolve percentage + real(r8) :: dissolve_frac = 1.0_r8 ! dissolve fraction real(r8) :: afunc ! A function in Patankar 1980, figure 5.6 real(r8) :: pcl ! Peclet number in Patankar 1980, foumula 5.18 real(r8) :: dz_node(1:nlevdecomp+1) ! difference between nodes @@ -157,7 +157,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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) = disp * smin_no3_vr(c,j) + 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 @@ -165,7 +165,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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) = disp * smin_no3_vr(c,j)/swliq(c,j) + conc_trcr(c,j) = dissolve_frac * smin_no3_vr(c,j)/swliq(c,j) qflx1 = qin(c,j) * mmh2o_to_m3h2o_per_m2 ! mm H2O/s to m3 H2O/m2/s qflx2 = qout(c,j) * mmh2o_to_m3h2o_per_m2 pcl1 = qflx1*dz_node(j)/total_D(c,j) @@ -184,7 +184,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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) = disp * smin_no3_vr(c,j)/swliq(c,j) + conc_trcr(c,j) = dissolve_frac * smin_no3_vr(c,j)/swliq(c,j) qflx1 = qin(c,j) * mmh2o_to_m3h2o_per_m2 ! mm H2O/s to m3 H2O/m2/s qflx2 = qout(c,j) * mmh2o_to_m3h2o_per_m2 pcl1 = qflx1*dz_node(j-1)/total_D(c,j-1) @@ -216,7 +216,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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-disp) + conc_trcr(c,j)*swliq(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 @@ -231,7 +231,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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)*disp + conc_trcr(c,j)*swliq(c,j) + 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 From c473c49139bf385b2f224c4b0e567f680717a9cf Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Fri, 22 Aug 2025 16:26:37 -0600 Subject: [PATCH 27/29] Rename variables pcl1, pcl2, qflx1, qflx2, and Dw as recommended --- .../SoilNitrogenMovementMod.F90 | 35 +++++++++---------- 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 index ec263f793c..abc47d5156 100644 --- a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 +++ b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 @@ -72,11 +72,13 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst real(r8) :: dtime ! land model time step (sec) real(r8) :: wafc, wafc2 ! Fraction of water that is liquid by mass real(r8) :: Ldis = 0.1_r8 ! dispersion length (m), Jury et al., 1991 - real(r8) :: pcl1, pcl2, qflx1, qflx2, theta, thetasat, nerror - real(r8) :: Dw = 1.7e-9_r8 ! Molecular diffusivity of NO3- in water, m2/s + 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) :: afunc ! A function in Patankar 1980, figure 5.6 real(r8) :: pcl ! Peclet number in Patankar 1980, foumula 5.18 + real(r8) :: pcl_in, pcl_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 @@ -131,7 +133,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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_D(c,j) = theta*(Dw * theta**(7/3) * thetasat**(-2)) + total_D(c,j) = theta * (no3_diffusivity_in_water * theta**(7/3) * thetasat**(-2)) total_D(c,j) = total_D(c,j) + Ldis * mmh2o_to_m3h2o_per_m2 * abs(qout(c,j)) else !no gradient for the last layer @@ -166,13 +168,11 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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) - qflx1 = qin(c,j) * mmh2o_to_m3h2o_per_m2 ! mm H2O/s to m3 H2O/m2/s - qflx2 = qout(c,j) * mmh2o_to_m3h2o_per_m2 - pcl1 = qflx1*dz_node(j)/total_D(c,j) - pcl2 = qflx2*dz_node(j)/total_D(c,j) + qflx_out = qout(c,j) * mmh2o_to_m3h2o_per_m2 + pcl_out = qflx_out * dz_node(j) / total_D(c,j) a_tri(c,j) = 0._r8 - c_tri(c,j) = -total_D(c,j)/dz_node(j) * afunc(pcl2) - max(-qflx2, 0._r8) - b_tri(c,j) = total_D(c,j)/dz_node(j) * afunc(pcl2) + max(qflx2, 0._r8) + swliq(c,j)/dtime*dzsoi_decomp(j) + c_tri(c,j) = -total_D(c,j) / dz_node(j) * afunc(pcl_out) - max(-qflx_out, 0._r8) + b_tri(c,j) = total_D(c,j) / dz_node(j) * afunc(pcl_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 @@ -185,14 +185,14 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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) - qflx1 = qin(c,j) * mmh2o_to_m3h2o_per_m2 ! mm H2O/s to m3 H2O/m2/s - qflx2 = qout(c,j) * mmh2o_to_m3h2o_per_m2 - pcl1 = qflx1*dz_node(j-1)/total_D(c,j-1) - pcl2 = qflx2*dz_node(j)/total_D(c,j) - a_tri(c,j) = -total_D(c,j-1)/dz_node(j-1) * afunc(pcl1) - max(qflx1, 0._r8) - c_tri(c,j) = -total_D(c,j)/dz_node(j) * afunc(pcl2) - max(-qflx2, 0._r8) - b_tri(c,j) = total_D(c,j-1)/dz_node(j-1) * afunc(pcl1) + max(-qflx1, 0._r8) + & - total_D(c,j)/dz_node(j) * afunc(pcl2) + max(qflx2, 0._r8) + swliq(c,j)/dtime*dzsoi_decomp(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 + pcl_in = qflx_in * dz_node(j-1) / total_D(c,j-1) + pcl_out = qflx_out * dz_node(j) / total_D(c,j) + a_tri(c,j) = -total_D(c,j-1) / dz_node(j-1) * afunc(pcl_in) - max(qflx_in, 0._r8) + c_tri(c,j) = -total_D(c,j) / dz_node(j) * afunc(pcl_out) - max(-qflx_out, 0._r8) + b_tri(c,j) = total_D(c,j-1) / dz_node(j-1) * afunc(pcl_in) + max(-qflx_in, 0._r8) + & + total_D(c,j) / dz_node(j) * afunc(pcl_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 @@ -222,7 +222,6 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst do fc = 1, num_bgc_soilc c = filter_bgc_soilc(fc) - nerror = mass_old(c) - mass_new(c) ! 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 From 09dbc062fa65a449c23903285c634ab830c87f21 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Mon, 25 Aug 2025 15:30:24 -0600 Subject: [PATCH 28/29] Rename variables as recommended in code review --- .../SoilNitrogenMovementMod.F90 | 36 +++++++++---------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 index abc47d5156..cc352fd827 100644 --- a/src/soilbiogeochem/SoilNitrogenMovementMod.F90 +++ b/src/soilbiogeochem/SoilNitrogenMovementMod.F90 @@ -71,19 +71,19 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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) :: Ldis = 0.1_r8 ! dispersion length (m), Jury et al., 1991 + 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) :: afunc ! A function in Patankar 1980, figure 5.6 - real(r8) :: pcl ! Peclet number in Patankar 1980, foumula 5.18 - real(r8) :: pcl_in, pcl_out ! temporary Peclet numbers + 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_D(bounds%begc:bounds%endc,1:nlevdecomp+1) ! Total diffusivity + 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 @@ -95,7 +95,7 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst ! 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 - afunc(pcl) = max(0._r8, ( 1._r8 - 0.1_r8 * abs(pcl) )**5 ) + 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] @@ -133,11 +133,11 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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_D(c,j) = theta * (no3_diffusivity_in_water * theta**(7/3) * thetasat**(-2)) - total_D(c,j) = total_D(c,j) + Ldis * mmh2o_to_m3h2o_per_m2 * abs(qout(c,j)) + 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_D(c,j) = total_D(c,j-1) + total_diffusivity(c,j) = total_diffusivity(c,j-1) end if end do ! Loop for columns end do ! Loop for depths @@ -169,10 +169,10 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst ! 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 - pcl_out = qflx_out * dz_node(j) / total_D(c,j) + peclet_num_out = qflx_out * dz_node(j) / total_diffusivity(c,j) a_tri(c,j) = 0._r8 - c_tri(c,j) = -total_D(c,j) / dz_node(j) * afunc(pcl_out) - max(-qflx_out, 0._r8) - b_tri(c,j) = total_D(c,j) / dz_node(j) * afunc(pcl_out) + max(qflx_out, 0._r8) + swliq(c,j) / dtime * dzsoi_decomp(j) + 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 @@ -187,12 +187,12 @@ subroutine SoilNitrogenMovement(bounds, num_bgc_soilc, filter_bgc_soilc, waterst 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 - pcl_in = qflx_in * dz_node(j-1) / total_D(c,j-1) - pcl_out = qflx_out * dz_node(j) / total_D(c,j) - a_tri(c,j) = -total_D(c,j-1) / dz_node(j-1) * afunc(pcl_in) - max(qflx_in, 0._r8) - c_tri(c,j) = -total_D(c,j) / dz_node(j) * afunc(pcl_out) - max(-qflx_out, 0._r8) - b_tri(c,j) = total_D(c,j-1) / dz_node(j-1) * afunc(pcl_in) + max(-qflx_in, 0._r8) + & - total_D(c,j) / dz_node(j) * afunc(pcl_out) + max(qflx_out, 0._r8) + swliq(c,j) / dtime * dzsoi_decomp(j) + 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 From f9195c5759f720c7c7573981f276d6ede6e128ed Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Tue, 26 Aug 2025 14:56:32 -0600 Subject: [PATCH 29/29] Final ChangeLog/Sum --- doc/ChangeLog | 19 +++++-------------- doc/ChangeSum | 2 +- 2 files changed, 6 insertions(+), 15 deletions(-) diff --git a/doc/ChangeLog b/doc/ChangeLog index dee5bceccf..69dc644440 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,13 +1,14 @@ =============================================================== Tag name: ctsm5.3.072 Originator(s): jinmuluo (Jinmu Luo, Cornell University) -Date: Mon 25 Aug 2025 11:04:53 AM MDT +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 ---------------------------------- - This tag introduces a new vertical movement scheme for soil nitrate but sets its use to .false. by default. Details about the new scheme in the PR #2992. + 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 -------------------------------------------------------------- @@ -56,25 +57,15 @@ Testing summary: build-namelist tests (if CLMBuildNamelist.pm has changed): - derecho - + 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 ------- - - fates tests: (give name of baseline if different from CTSM tagname, normally fates baselines are fates--) - derecho ----- - izumi ------- - - any other testing (give details below): - - ctsm_sci - derecho ---- + izumi ------- OK Answer changes -------------- - Changes answers relative to baseline: No Other details diff --git a/doc/ChangeSum b/doc/ChangeSum index f380869205..b8340937d7 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,6 +1,6 @@ Tag Who Date Summary ============================================================================================================================ - ctsm5.3.072 jinmuluo 08/xx/2025 New vertical movement scheme for soil nitrate + 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