Merging soil NO into the new vertical movement scheme for soil nitrate

bld/namelist_files/namelist_definition_ctsm.xml

@@ -996,6 +996,11 @@ Method for distributing soil thickness across hillslope columns
 Toggle to turn on the plant hydraulic stress model
 </entry>
 
+<entry id="use_soil_nox" type="logical" category="physics"
+        group="clm_inparm" valid_values="" value=".false.">
+Toggle to turn on soil NOx in the N cycle
+</entry>
+
 <entry id="lnc_opt" type="logical" category="clm_nitrogen"
        group="clm_nitrogen" value=".false.">
 How LUNA and Photosynthesis (if needed) will get Leaf nitrogen content

src/biogeochem/CNBalanceCheckMod.F90

@@ -10,7 +10,7 @@ module CNBalanceCheckMod
   use shr_log_mod                     , only : errMsg => shr_log_errMsg
   use decompMod                       , only : bounds_type, subgrid_level_gridcell, subgrid_level_column
   use abortutils                      , only : endrun
-  use clm_varctl                      , only : iulog, use_nitrif_denitrif, use_fates_bgc
+  use clm_varctl                      , only : iulog, use_nitrif_denitrif, use_fates_bgc, use_soil_nox
   use clm_time_manager                , only : get_step_size_real
   use CNVegNitrogenFluxType           , only : cnveg_nitrogenflux_type
   use CNVegNitrogenStateType          , only : cnveg_nitrogenstate_type
@@ -547,6 +547,14 @@ subroutine NBalanceCheck(this, bounds, num_soilc, filter_soilc, &
          smin_no3_leached    => soilbiogeochem_nitrogenflux_inst%smin_no3_leached_col    , & ! Input:  [real(r8) (:) ]  (gN/m2/s) soil mineral NO3 pool loss to leaching 
          smin_no3_runoff     => soilbiogeochem_nitrogenflux_inst%smin_no3_runoff_col     , & ! Input:  [real(r8) (:) ]  (gN/m2/s) soil mineral NO3 pool loss to runoff   
          f_n2o_nit           => soilbiogeochem_nitrogenflux_inst%f_n2o_nit_col           , & ! Input:  [real(r8) (:) ]  (gN/m2/s) flux of N2o from nitrification 
+
+         f_n2o_denit         => soilbiogeochem_nitrogenflux_inst%f_n2o_denit_col         , & ! Input:  [real(r8) (:) ]  (gN/m2/s) flux of N2o from denitrification
+         f_n2_denit         => soilbiogeochem_nitrogenflux_inst%f_n2_denit_col           , & ! Input:  [real(r8) (:) ]  (gN/m2/s) flux of N2 from denitrification
+         f_nox_denit         => soilbiogeochem_nitrogenflux_inst%f_nox_denit_col         , & ! Input:[real(r8) (:) ]  (gN/m2/s) flux of NOx from denitrification
+         f_nox_nit         => soilbiogeochem_nitrogenflux_inst%f_nox_nit_col           , & ! Input:[real(r8) (:) ]  (gN/m2/s) flux of NOx from nitrification
+         f_nox_denit_atmos    => soilbiogeochem_nitrogenflux_inst%f_nox_denit_atmos_col  , & ! Input:[real(r8) (:) ]  (gN/m2/s) flux of NOx from denitrification
+         f_nox_nit_atmos      => soilbiogeochem_nitrogenflux_inst%f_nox_nit_atmos_col    , & ! Input:  [real(r8) (:) ]  (gN/m2/s) flux of NOx from nitrification
+
          som_n_leached       => soilbiogeochem_nitrogenflux_inst%som_n_leached_col       , & ! Input:  [real(r8) (:) ]  (gN/m2/s) total SOM N loss from vertical transport
 
          col_fire_nloss      => cnveg_nitrogenflux_inst%fire_nloss_col                   , & ! Input:  [real(r8) (:) ]  (gN/m2/s) total column-level fire N loss 
@@ -599,9 +607,15 @@ subroutine NBalanceCheck(this, bounds, num_soilc, filter_soilc, &
          col_ninputs_partial(c) = col_ninputs(c)
 
          ! calculate total column-level outputs
-
-         col_noutputs(c) = denit(c)
-
+         if (use_soil_nox) then
+            !denit=f_n2o_denit + f_n2_denit + f_nox_denit
+            col_noutputs(c) = f_n2_denit(c)+f_n2o_denit(c)
+            ! only f_nox_denit_atmos leaves the system, NOx denit trapped in the canopy stays in the system.
+            col_noutputs(c) = col_noutputs(c)+f_nox_denit_atmos(c) 
+         else
+            col_noutputs(c) = denit(c)
+         endif
+
          if( .not.col%is_fates(c) ) then
 
             col_noutputs(c) = col_noutputs(c) + col_fire_nloss(c) + gru_conv_nflux(c)
@@ -626,8 +640,10 @@ subroutine NBalanceCheck(this, bounds, num_soilc, filter_soilc, &
             col_noutputs(c) = col_noutputs(c) + sminn_leached(c)
          else
             col_noutputs(c) = col_noutputs(c) + f_n2o_nit(c)
-
             col_noutputs(c) = col_noutputs(c) + smin_no3_leached(c) + smin_no3_runoff(c)
+            if (use_soil_nox) then
+               col_noutputs(c) = col_noutputs(c) + f_nox_nit_atmos(c)
+            endif
          end if
 
          col_noutputs(c) = col_noutputs(c) - som_n_leached(c)
@@ -652,7 +668,11 @@ subroutine NBalanceCheck(this, bounds, num_soilc, filter_soilc, &
          if (abs(col_errnb(c)) > this%nwarning) then
             write(iulog,*) 'nbalance warning at c =', c, col_errnb(c), col_endnb(c)
             write(iulog,*)'inputs,ffix,nfix,ndep = ',ffix_to_sminn(c)*dt,nfix_to_sminn(c)*dt,ndep_to_sminn(c)*dt
-            write(iulog,*)'outputs,lch,roff,dnit = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,f_n2o_nit(c)*dt
+            write(iulog,*)'outputs,lch,roff,dnit = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,denit(c)*dt
+            if (use_soil_nox) then
+               write(iulog,*)'outputs,fnoxnit,fnoxdenit,fn2 = ',f_nox_nit(c)*dt, f_nox_denit(c)*dt,f_n2_denit(c)*dt
+               write(iulog,*)'outputs,fn2onit,fn2odenit= ',f_n2o_nit(c)*dt, f_n2o_denit(c)*dt
+            end if
          end if
 
       end do ! end of columns loop
@@ -673,9 +693,9 @@ subroutine NBalanceCheck(this, bounds, num_soilc, filter_soilc, &
             write(iulog,*)'inputs,ffix,nfix,ndep = ',ffix_to_sminn(c)*dt,nfix_to_sminn(c)*dt,ndep_to_sminn(c)*dt
          end if
          if(col%is_fates(c))then
-            write(iulog,*)'outputs,lch,roff,dnit,plnt = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,f_n2o_nit(c)*dt,sminn_to_plant(c)*dt
+            write(iulog,*)'outputs,lch,roff,dnit,plnt = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,denit(c)*dt,sminn_to_plant(c)*dt
          else
-            write(iulog,*)'outputs,lch,roff,dnit    = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,f_n2o_nit(c)*dt
+            write(iulog,*)'outputs,lch,roff,dnit    = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,denit(c)*dt
          end if
          call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, msg=errMsg(sourcefile, __LINE__))
       end if

src/biogeochem/CNDriverMod.F90

@@ -45,6 +45,7 @@ module CNDriverMod
   use CLMFatesInterfaceMod            , only : hlm_fates_interface_type
   use CropReprPoolsMod                , only : nrepr
   use SoilHydrologyType               , only: soilhydrology_type
+  use DryDepVelocity                  , only : drydepvel_type
   !
   ! !PUBLIC TYPES:
   implicit none
@@ -103,7 +104,8 @@ subroutine CNDriverNoLeaching(bounds,
        wateratm2lndbulk_inst, canopystate_inst, soilstate_inst, temperature_inst,          &
        soil_water_retention_curve, crop_inst, ch4_inst,                                    &
        dgvs_inst, photosyns_inst, saturated_excess_runoff_inst, energyflux_inst,           &
-       nutrient_competition_method, cnfire_method, dribble_crophrv_xsmrpool_2atm)
+       nutrient_competition_method, cnfire_method, dribble_crophrv_xsmrpool_2atm,          &
+       drydepvel_inst)
     !
     ! !DESCRIPTION:
     ! The core CN code is executed here. Calculates fluxes for maintenance
@@ -207,6 +209,7 @@ subroutine CNDriverNoLeaching(bounds,
     class(fire_method_type)                 , intent(inout) :: cnfire_method
     logical                                 , intent(in)    :: dribble_crophrv_xsmrpool_2atm
     type(hlm_fates_interface_type)          , intent(inout) :: clm_fates
+    type(drydepvel_type)                    , intent(inout)    :: drydepvel_inst
     !
     ! !LOCAL VARIABLES:
     real(r8):: cn_decomp_pools(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_pools)
@@ -479,7 +482,7 @@ subroutine CNDriverNoLeaching(bounds,
                                      cnveg_carbonflux_inst,cnveg_nitrogenstate_inst,cnveg_nitrogenflux_inst,   &
                                      soilbiogeochem_carbonflux_inst,&
                                      soilbiogeochem_state_inst,soilbiogeochem_nitrogenstate_inst,              &
-                                     soilbiogeochem_nitrogenflux_inst,canopystate_inst)
+                                     soilbiogeochem_nitrogenflux_inst,canopystate_inst, drydepvel_inst)
      call t_stopf('soilbiogeochemcompetition')
 
     ! distribute the available N between the competing patches  on the basis of 

src/biogeochem/CNVegetationFacade.F90

@@ -97,6 +97,7 @@ module CNVegetationFacade
   use SoilWaterRetentionCurveMod      , only : soil_water_retention_curve_type
   use CLMFatesInterfaceMod            , only : hlm_fates_interface_type
   use SoilHydrologyType               , only : soilhydrology_type
+  use DryDepVelocity                  , only : drydepvel_type
   !
   implicit none
   private
@@ -128,7 +129,7 @@ module CNVegetationFacade
      type(cn_balance_type)          :: cn_balance_inst
      class(fire_method_type), allocatable :: cnfire_method
      type(dgvs_type)                :: dgvs_inst
-
+     type(drydepvel_type)           :: drydepvel_inst   
      ! Control variables
      logical, private :: reseed_dead_plants              ! Flag to indicate if should reseed dead plants when starting up the model
      logical, private :: dribble_crophrv_xsmrpool_2atm = .False. ! Flag to indicate if should harvest xsmrpool to the atmosphere
@@ -963,7 +964,7 @@ subroutine EcosystemDynamicsPreDrainage(this, bounds, &
        wateratm2lndbulk_inst, canopystate_inst, soilstate_inst, temperature_inst, &
        soil_water_retention_curve, crop_inst, ch4_inst, &
        photosyns_inst, saturated_excess_runoff_inst, energyflux_inst,          &
-       nutrient_competition_method, fireemis_inst)
+       nutrient_competition_method, fireemis_inst, drydepvel_inst)
     !
     ! !DESCRIPTION:
     ! Do the main science for biogeochemistry that needs to be done before hydrology-drainage
@@ -1020,6 +1021,7 @@ subroutine EcosystemDynamicsPreDrainage(this, bounds, &
     class(nutrient_competition_method_type) , intent(inout) :: nutrient_competition_method
     type(fireemis_type)                     , intent(inout) :: fireemis_inst
     type(hlm_fates_interface_type)          , intent(inout) :: clm_fates
+    type(drydepvel_type)                    , intent(inout) :: drydepvel_inst  
     !
     ! !LOCAL VARIABLES:
 
@@ -1054,7 +1056,8 @@ subroutine EcosystemDynamicsPreDrainage(this, bounds, &
          wateratm2lndbulk_inst, canopystate_inst, soilstate_inst, temperature_inst, &
          soil_water_retention_curve, crop_inst, ch4_inst, &
          this%dgvs_inst, photosyns_inst, saturated_excess_runoff_inst, energyflux_inst,          &
-         nutrient_competition_method, this%cnfire_method, this%dribble_crophrv_xsmrpool_2atm)
+         nutrient_competition_method, this%cnfire_method, this%dribble_crophrv_xsmrpool_2atm,    &
+         drydepvel_inst)
 
     ! fire carbon emissions 
     call CNFireEmisUpdate(bounds, num_bgc_vegp, filter_bgc_vegp, &

src/biogeochem/DryDepVelocity.F90

@@ -84,6 +84,7 @@ Module DryDepVelocity
 
      real(r8), pointer, public  :: velocity_patch (:,:) ! Dry Deposition Velocity
      real(r8), pointer, private :: rs_drydep_patch (:)  ! Stomatal resistance associated with dry deposition velocity for Ozone
+     real(r8), pointer, public :: crf_drydep_patch (:)  ! Canopy reduction factor for NO associated with dry deposition velocity
 
    contains
 
@@ -136,6 +137,7 @@ subroutine InitAllocate(this, bounds)
     if ( n_drydep > 0 )then
        allocate(this%velocity_patch(begp:endp, n_drydep));  this%velocity_patch(:,:) = nan
        allocate(this%rs_drydep_patch(begp:endp))         ;  this%rs_drydep_patch(:)  = nan
+       allocate(this%crf_drydep_patch(begp:endp))        ;  this%crf_drydep_patch(:)  = nan 
     end if
 
   end subroutine InitAllocate
@@ -180,6 +182,11 @@ subroutine InitHistory(this, bounds)
          avgflag='A', long_name='Stomatal Resistance Associated with Ozone Dry Deposition Velocity', &
          ptr_patch=this%rs_drydep_patch, default='inactive' )
 
+    this%crf_drydep_patch(begp:endp)= spval
+    call hist_addfld1d ( fname='CRF_DRYDEP_NO', units='unitless',  &
+         avgflag='A', long_name='Canopy Reduction Factor Associated with NO Dry Deposition Velocity', &
+         ptr_patch=this%crf_drydep_patch, default='inactive')
+
   end subroutine InitHistory
 
   !-----------------------------------------------------------------------
@@ -253,7 +260,10 @@ subroutine depvel_compute( bounds, &
 
     !mvm 11/30/2013
     real(r8) :: rlu_lai      ! constant to calculate rlu over bulk canopy
-
+    real(r8) :: ratio_stai_lai     !ratio of stomata area index to lai [unitless]
+    real(r8) :: ks     !factor to adjust SAI for soil canopy reduction [unitless] Yan et al (2005) https://doi.org/10.1029/2004GB002276
+    real(r8) :: kc     !factor to adjust LAI for soil canopy reduction [unitless] Yan et al (2005) https://doi.org/10.1029/2004GB002276
+
     logical  :: has_dew
     logical  :: has_rain
     real(r8), parameter :: rain_threshold      = 1.e-7_r8  ! of the order of 1cm/day expressed in m/s
@@ -268,7 +278,7 @@ subroutine depvel_compute( bounds, &
     real(r8) :: rc    !combined surface resistance
     real(r8) :: cts   !correction to flu rcl and rgs for frost
     real(r8) :: rlux_o3 !to calculate O3 leaf resistance in dew/rain conditions
-
+     
     ! constants
     real(r8), parameter :: slope = 0._r8      ! Used to calculate  rdc in (lower canopy resistance)
     integer, parameter :: wveg_unset = -1     ! Unset Wesley vegetation type
@@ -307,7 +317,9 @@ subroutine depvel_compute( bounds, &
          annlai     =>    canopystate_inst%annlai_patch         , & ! Input:  [real(r8) (:,:) ] 12 months of monthly lai from input data set
 
          velocity   =>    drydepvel_inst%velocity_patch         , & ! Output: [real(r8) (:,:) ] cm/sec
-         rs_drydep  =>    drydepvel_inst%rs_drydep_patch          & ! Output: [real(r8) (:)   ]  stomatal resistance associated with Ozone dry deposition velocity (s/m)
+         rs_drydep  =>    drydepvel_inst%rs_drydep_patch        ,  & ! Output: [real(r8) (:)   ]  stomatal resistance associated with Ozone dry deposition velocity (s/m)
+         crf_drydep  =>   drydepvel_inst%crf_drydep_patch        &    ! Output: [real(r8) (:)   ]  canopy reduction factor for NO associated with dry deposition velocity [unitless]
+
          )
 
       !_________________________________________________________________
@@ -654,6 +666,52 @@ subroutine depvel_compute( bounds, &
                     (1._r8/(rdc+rclx(ispec)))+(1._r8/(rac(index_season,wesveg)+rgsx(ispec))))
                rc = max( 10._r8, rc)
                !
+
+               ! CANOPY REDUCTION FACTOR FOR soil NO fluxes
+               ! (Probably there may be a cleaner way to code this)
+               if ( drydep_list(ispec) == 'NO') then
+                  !factors for canopy reduction
+                  ks=11.6_r8
+                  kc=0.32_r8
+
+             !no canopy reduction factor (ie zero update) over bare land or elai < 0. 
+                  if( clmveg == noveg) then
+                     crf_drydep(pi)=1._r8                 
+                  else
+                     ratio_stai_lai = 0._r8  ! Initialize
+
+                     !needleaf forest (coniferous forest)              
+                     if (clmveg == ndllf_evr_tmp_tree .or. clmveg == ndllf_evr_brl_tree .or.  clmveg == ndllf_dcd_brl_tree ) ratio_stai_lai = 0.03_r8 
+
+          !broadleaf evergreen (rainforest)
+                     if (clmveg == nbrdlf_evr_trp_tree .or. clmveg == nbrdlf_evr_tmp_tree ) ratio_stai_lai = 0.015_r8 
+
+          !broadleaf decidouos 
+                     if (clmveg == nbrdlf_dcd_trp_tree .or. clmveg == nbrdlf_dcd_tmp_tree .or. clmveg == nbrdlf_dcd_brl_tree) then
+                        if (index_season == 1 .or. index_season == 5) then 
+                           ! spring/summer
+                           ratio_stai_lai = 0.015_r8 
+                        else 
+                           !winter/fall
+                           ratio_stai_lai = 0.005_r8
+                        endif
+                     endif
+
+                     !shrubs
+                     if (clmveg == nbrdlf_evr_shrub  .or. clmveg == nbrdlf_dcd_tmp_shrub .or.  clmveg == nbrdlf_dcd_brl_shrub ) ratio_stai_lai = 0.005_r8 
+
+                     !grass
+                     if (clmveg == nc3_arctic_grass .or. clmveg == nc3_nonarctic_grass .or. clmveg == nc4_grass ) ratio_stai_lai = 0.005_r8 
+
+                     !crops
+                     if (clmveg == nc3crop .or. clmveg == nc3irrig ) ratio_stai_lai = 0.008_r8 
+                     if (is_prognostic_crop(clmveg)) ratio_stai_lai = 0.008_r8
+
+                     crf_drydep(pi)=(exp(-ks*elai(pi)*ratio_stai_lai)+exp(-kc*elai(pi)))/2
+                  endif !no veg 
+               endif ! drydep_list(ispec)
+
+
                ! assume no surface resistance for SO2 over water
                !
                if ( drydep_list(ispec) == 'SO2' .and. wesveg == 7 ) then

src/biogeophys/WaterStateType.F90

@@ -40,6 +40,9 @@ module WaterStateType
      real(r8), pointer :: h2osoi_vol_col         (:,:) ! col volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3]  (nlevgrnd)
      real(r8), pointer :: h2osoi_vol_prs_grc     (:,:) ! grc volumetric soil water prescribed (0<=h2osoi_vol<=watsat) [m3/m3]  (nlevgrnd)
      real(r8), pointer :: h2osfc_col             (:)   ! col surface water (mm H2O)
+     real(r8), pointer :: h2osoi_vol_old_col     (:,:) ! col volumetric soil water for previous time step (0<=h2osoi_vol<=watsat) [m3/m3]  (nlevgrnd)
+     real(r8), pointer :: ldry_col               (:,:) ! col dry period rain pulses for NOx from nitrification as state variable [hours]
+
      real(r8), pointer :: snocan_patch           (:)   ! patch canopy snow water (mm H2O)
      real(r8), pointer :: liqcan_patch           (:)   ! patch canopy liquid water (mm H2O)
 
@@ -142,6 +145,15 @@ subroutine InitAllocate(this, bounds, tracer_vars)
          container = tracer_vars, &
          bounds = bounds, subgrid_level = subgrid_level_column, &
          dim2beg = -nlevsno+1, dim2end = nlevmaxurbgrnd)
+    call AllocateVar2d(var = this%h2osoi_vol_old_col, name = 'h2osoi_vol_old_col', &
+         container = tracer_vars, &
+         bounds = bounds, subgrid_level = subgrid_level_column, &
+         dim2beg = 1, dim2end = nlevmaxurbgrnd)
+    call AllocateVar2d(var = this%ldry_col, name = 'ldry_col', &
+         container = tracer_vars, &
+         bounds = bounds, subgrid_level = subgrid_level_column, &
+         dim2beg = 1, dim2end = nlevmaxurbgrnd)       
+
     call AllocateVar1d(var = this%snocan_patch, name = 'snocan_patch', &
          container = tracer_vars, &
          bounds = bounds, subgrid_level = subgrid_level_patch)