Move the dust emission source function soil erodibility for the Zender scheme from CAM to CTSM

bld/CLMBuildNamelist.pm

@@ -1560,6 +1560,7 @@ sub process_namelist_inline_logic {
   setup_logic_snowpack($opts,  $nl_flags, $definition, $defaults, $nl);
   setup_logic_fates($opts,  $nl_flags, $definition, $defaults, $nl);
   setup_logic_misc($opts, $nl_flags, $definition, $defaults, $nl);
+  setup_logic_zender_soilerod($opts, $nl_flags, $definition, $defaults, $nl);
 
   #########################################
   # namelist group: atm2lnd_inparm
@@ -4245,6 +4246,17 @@ sub setup_logic_misc {
 
 #-------------------------------------------------------------------------------
 
+sub setup_logic_zender_soilerod {
+  #
+  # Handle the Zender soil eroditability file
+  #
+  my ($opts, $nl_flags, $definition, $defaults, $nl) = @_;
+  add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldfilename_zendersoilerod' );
+  add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_meshfile_zendersoilerod' );
+}
+
+#-------------------------------------------------------------------------------
+
 sub write_output_files {
   my ($opts, $nl_flags, $defaults, $nl) = @_;
 
@@ -4293,6 +4305,7 @@ sub write_output_files {
   push @groups, "ch4finundated";
   push @groups, "soilbgc_decomp";
   push @groups, "clm_canopy_inparm";
+  push @groups, "zendersoilerod";
   if (remove_leading_and_trailing_quotes($nl->get_value('snow_cover_fraction_method')) eq 'SwensonLawrence2012') {
      push @groups, "scf_swenson_lawrence_2012_inparm";
   }

bld/namelist_files/namelist_defaults_ctsm.xml

@@ -2537,6 +2537,16 @@ lnd/clm2/surfdata_map/release-clm5.0.30/surfdata_ne0np4.CONUS.ne30x8_hist_78pfts
 <stream_meshfile_ch4finundated finundation_method="ZWT_inversion"
 >lnd/clm2/paramdata/finundated_inversiondata_0.9x1_ESMFmesh_cdf5_130621.nc</stream_meshfile_ch4finundated>
 
+<!-- ===========================================  -->
+<!-- Defaults for dust namelist items             -->
+<!-- Defaults for Zender's soil erodibility data  -->
+<!-- ===========================================  -->
+
+<stream_fldfilename_zendersoilerod
+>lnd/clm2/dustemisdata/dst_source2x2tunedcam6-2x2-forCLM_cdf5_c230312.nc</stream_fldfilename_zendersoilerod>
+<stream_meshfile_zendersoilerod
+>lnd/clm2/dustemisdata/dust_2x2_ESMFmesh_cdf5_c230730.nc</stream_meshfile_zendersoilerod>
+
 <!-- =========================================  -->
 <!-- Defaults for different BGC/decomp modes    -->
 <!-- =========================================  -->

bld/namelist_files/namelist_definition_ctsm.xml

@@ -1693,6 +1693,25 @@ Mapping method from Nitrogen deposition input file to the model resolution
     copy     = copy using the same indices
 </entry>
 
+<!-- ========================================================================================  -->
+<!-- Dust namelist and Zender soil erodibility stream                                          -->
+<!-- ========================================================================================  -->
+
+<!-- <entry id="dust_emis_fact" type="real" category="dataset" -->
+<!--        group="dust_nl" valid_values="" > -->
+<!-- Tuning parameter for dust emissions. -->
+<!-- </entry> -->
+
+<entry id="stream_fldfilename_zendersoilerod" type="char*256" category="datasets"
+       input_pathname="abs" group="zendersoilerod" valid_values="" >
+Filename of input stream data for Zender's soil erodibility source function
+</entry>
+
+<entry id="stream_meshfile_zendersoilerod" type="char*256" category="datasets"
+       input_pathname="abs" group="zendersoilerod" valid_values="" >
+mesh filename of input stream data for Zender's soil erodibility source function
+</entry>
+
 <!-- ========================================================================================  -->
 <!-- finundated stream -->
 <!-- ========================================================================================  -->

src/biogeochem/DUSTMod.F90

@@ -30,6 +30,7 @@ module DUSTMod
   use LandunitType         , only : lun
   use ColumnType           , only : col
   use PatchType            , only : patch
+  use ZenderSoilErodStreamType,  only : zendersoilerodstream_type
   !  
   ! !PUBLIC TYPES
   implicit none
@@ -59,7 +60,8 @@ module DUSTMod
      real(r8), pointer, private :: vlc_trb_2_patch           (:)   ! turbulent deposition velocity 2(m/s)
      real(r8), pointer, private :: vlc_trb_3_patch           (:)   ! turbulent deposition velocity 3(m/s)
      real(r8), pointer, private :: vlc_trb_4_patch           (:)   ! turbulent deposition velocity 4(m/s)
-     real(r8), pointer, private :: mbl_bsn_fct_col           (:)   ! basin factor
+     type(zendersoilerodstream_type), private :: zendersoilerodstream ! Zender soil erodibility stream data
+     real(r8), pointer, private :: mbl_bsn_fct_col           (:)   ! [dimensionless] basin factor, or soil erodibility, time-constant
 
    contains
 
@@ -78,13 +80,16 @@ module DUSTMod
 contains
 
   !------------------------------------------------------------------------
-  subroutine Init(this, bounds)
+  !subroutine Init(this, bounds)
+  subroutine Init(this, bounds, NLFilename)
 
     class(dust_type) :: this
     type(bounds_type), intent(in) :: bounds  
+    character(len=*),  intent(in) :: NLFilename
 
     call this%InitAllocate (bounds)
     call this%InitHistory  (bounds)
+    call this%zendersoilerodstream%Init( bounds, NLFilename )
     call this%InitCold     (bounds)
     call this%InitDustVars (bounds)
 
@@ -112,7 +117,7 @@ subroutine InitAllocate(this, bounds)
     allocate(this%vlc_trb_2_patch           (begp:endp))        ; this%vlc_trb_2_patch           (:)   = nan 
     allocate(this%vlc_trb_3_patch           (begp:endp))        ; this%vlc_trb_3_patch           (:)   = nan
     allocate(this%vlc_trb_4_patch           (begp:endp))        ; this%vlc_trb_4_patch           (:)   = nan
-    allocate(this%mbl_bsn_fct_col           (begc:endc))        ; this%mbl_bsn_fct_col     (:)   = nan
+    allocate(this%mbl_bsn_fct_col           (begc:endc))        ; this%mbl_bsn_fct_col           (:)   = nan
 
   end subroutine InitAllocate
 
@@ -129,9 +134,11 @@ subroutine InitHistory(this, bounds)
     !
     ! !LOCAL VARIABLES:
     integer :: begp,endp
+    integer :: begc,endc
     !------------------------------------------------------------------------
 
     begp = bounds%begp; endp = bounds%endp
+    begc = bounds%begc; endc = bounds%endp
 
     this%flx_mss_vrt_dst_tot_patch(begp:endp) = spval
     call hist_addfld1d (fname='DSTFLXT', units='kg/m2/s',  &
@@ -158,6 +165,11 @@ subroutine InitHistory(this, bounds)
          avgflag='A', long_name='turbulent deposition velocity 4', &
          ptr_patch=this%vlc_trb_4_patch, default='inactive')
 
+    !this%mbl_bsn_fct_col(begc:endc) = spval
+    !call hist_addfld1d (fname='LND_MBL', units='fraction',  &
+    !     avgflag='A', long_name='Soil erodibility factor', &
+    !     ptr_col=this%mbl_bsn_fct_col)
+
   end subroutine InitHistory
 
   !-----------------------------------------------------------------------
@@ -169,17 +181,48 @@ subroutine InitCold(this, bounds)
     !
     ! !LOCAL VARIABLES:
     integer :: c,l
+    !integer ier
     !-----------------------------------------------------------------------
 
+    !#### dmleung added 11 Mar 2023 ####################################
+    ! read in global tuning factor from namelist
+    !namelist /dust_nl/ dust_emis_fact
+    ! Default values for namelist
+    !dust_emis_fact      = 1
+    ! Read soilm_streams namelist
+    !if (masterproc) then
+    !   open( newunit=nu_nml, file=trim(NLFilename), status='old', iostat=nml_error )
+    !   call find_nlgroup_name(nu_nml, 'dust_nl', status=nml_error)
+    !   if (nml_error == 0) then
+    !      read(nu_nml, nml=dust_nl,iostat=nml_error)
+    !      if (nml_error /= 0) then
+    !         call endrun(subname // ':: ERROR reading dust_nl namelist')
+    !      end if
+    !   else
+    !      call endrun(subname // ':: ERROR finding dust_nl namelist')
+    !   end if
+    !   close(nu_nml)
+    !endif
+    !call shr_mpi_bcast(dust_emis_fact, mpicom)
+    !call mpi_bcast (dust_emis_fact, 1, MPI_REAL8,0, mpicom, ier)
+    !if (masterproc) then
+    !   write(iulog,*) '  dust_emis_fact  = ',dust_emis_fact
+    !end if
+
     ! Set basin factor to 1 for now
 
-    do c = bounds%begc, bounds%endc
-       l = col%landunit(c)
+    !if (dust_emission_scheme == 'Leung2023')
+    !   do c = bounds%begc, bounds%endc
+    !      l = col%landunit(c)
 
-       if (.not.lun%lakpoi(l)) then
-          this%mbl_bsn_fct_col(c) = 1.0_r8
-       end if
-    end do
+    !      if (.not.lun%lakpoi(l)) then
+    !         this%mbl_bsn_fct_col(c) = 1.0_r8
+    !      end if
+    !   end do
+    !else if (dust_emission_scheme == 'Zender2003')
+       call this%zendersoilerodstream%CalcDustSource( bounds, &
+                                this%mbl_bsn_fct_col(bounds%begc:bounds%endc) )
+    !end if
 
   end subroutine InitCold
 
@@ -262,7 +305,7 @@ subroutine DustEmission (bounds, &
          fv                  => frictionvel_inst%fv_patch            , & ! Input:  [real(r8) (:)   ]  friction velocity (m/s) (for dust model)          
          u10                 => frictionvel_inst%u10_patch           , & ! Input:  [real(r8) (:)   ]  10-m wind (m/s) (created for dust model)          
 
-         mbl_bsn_fct         => dust_inst%mbl_bsn_fct_col            , & ! Input:  [real(r8) (:)   ]  basin factor                                      
+         mbl_bsn_fct         => dust_inst%mbl_bsn_fct_col            , & ! Input:  [real(r8) (:)   ]  basin factor
          flx_mss_vrt_dst     => dust_inst%flx_mss_vrt_dst_patch      , & ! Output: [real(r8) (:,:) ]  surface dust emission (kg/m**2/s)               
          flx_mss_vrt_dst_tot => dust_inst%flx_mss_vrt_dst_tot_patch    & ! Output: [real(r8) (:)   ]  total dust flux back to atmosphere (pft)
          )
@@ -436,6 +479,7 @@ subroutine DustEmission (bounds, &
 
             dst_slt_flx_rat_ttl = 100.0_r8 * exp( log(10.0_r8) * (13.4_r8 * mss_frc_cly_vld(c) - 6.0_r8) )
             flx_mss_vrt_dst_ttl(p) = flx_mss_hrz_slt_ttl * dst_slt_flx_rat_ttl
+            !flx_mss_vrt_dst_ttl(p) = flx_mss_hrz_slt_ttl * dst_slt_flx_rat_ttl / dust_emis_fact
 
          end if   ! lnd_frc_mbl > 0.0