A new physically based dust emission scheme with more aeolian physics

bld/namelist_files/namelist_definition_clm4_5.xml

@@ -2572,4 +2572,13 @@ use case.)
 
 </entry>
 
+<!-- added by dmleung to CLM5, 17 Dec 2021 -->
+<!-- roughness factor f_eff   -dmleung, 22 Apr 2020 -->
+<entry id="rough_fct" type="char*256" input_pathname="abs" category="datasets"
+       group="clm_inparm" valid_values="" >
+Full pathname of time-invariant roughness factor dataset for calculating drag partition of dust emission.
+</entry>
+
+
+
 </namelist_definition>

src/biogeophys/SoilStateInitTimeConstMod.F90

@@ -100,7 +100,8 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
     use clm_varcon          , only : zsoi, dzsoi, zisoi, spval
     use clm_varcon          , only : secspday, pc, mu, denh2o, denice, grlnd
     use clm_varctl          , only : use_cn, use_lch4, use_fates
-    use clm_varctl          , only : iulog, fsurdat, paramfile, soil_layerstruct
+    use clm_varctl          , only : iulog, fsurdat, paramfile, soil_layerstruct, &
+                                     rough_fct    ! -dmleung added to CESM2/CLM5 17 Dec 2021
     use landunit_varcon     , only : istdlak, istwet, istsoil, istcrop, istice_mec
     use column_varcon       , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv 
     use fileutils           , only : getfil
@@ -151,6 +152,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
     real(r8) ,pointer  :: sand3d (:,:)                  ! read in - soil texture: percent sand (needs to be a pointer for use in ncdio)
     real(r8) ,pointer  :: clay3d (:,:)                  ! read in - soil texture: percent clay (needs to be a pointer for use in ncdio)
     real(r8) ,pointer  :: organic3d (:,:)               ! read in - organic matter: kg/m3 (needs to be a pointer for use in ncdio)
+    real(r8) ,pointer  :: roughfct2d(:)                 ! read in - time-invariant roughness factor,   dmleung added 17 Dec 2021
     character(len=256) :: locfn                         ! local filename
     integer            :: ipedof  
     integer            :: begp, endp
@@ -297,6 +299,25 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
 
     call ncd_pio_closefile(ncid)
 
+    !############ subsection for input data for new dust emission scheme ##############
+
+    ! dmleung added to CESM2/CLM5 17 Dec 2021
+    allocate(roughfct2d(begg:endg))                                   ! dmleung, 16 Jul 2020
+    ! here to read roughness factor file, 16 Jul 2020
+    !write(iulog,*) 'Attempting to read roughness factor data, by dmleung .....'
+    call getfil (rough_fct, locfn, 0)
+    call ncd_pio_openfile (ncid, locfn, 0)
+    call ncd_io(ncid=ncid, varname='F_eff', flag='read', data=roughfct2d, dim1name=grlnd, readvar=readvar)
+    !write(iulog,*) 'initialize pft level roughness factor from roughfct2d(g) to roughfct(p)'  
+
+    do p = begp,endp
+       g = patch%gridcell(p)
+       soilstate_inst%roughfct_patch(p) = roughfct2d(g)
+    end do
+
+    call ncd_pio_closefile(ncid)
+    !##################################################################################
+
     ! --------------------------------------------------------------------
     ! get original soil depths to be used in interpolation of sand and clay
     ! --------------------------------------------------------------------
@@ -614,7 +635,8 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename)
     do c = begc,endc
        g = col%gridcell(c)
 
-       soilstate_inst%gwc_thr_col(c) = 0.17_r8 + 0.14_r8 * clay3d(g,1) * 0.01_r8
+       soilstate_inst%gwc_thr_col(c) = 0.01_r8*(0.17_r8*clay3d(g,1) + 0.0014_r8*clay3d(g,1)*clay3d(g,1))  !Fecan et al. (1999) -jfk, dmleung coded 27 Nov 2021 for CLM clay fraction
+       !soilstate_inst%gwc_thr_col(c) = 0.17_r8 + 0.14_r8 * clay3d(g,1) * 0.01_r8 -dmleung commented 27 Nov 2021
        soilstate_inst%mss_frc_cly_vld_col(c) = min(clay3d(g,1) * 0.01_r8, 0.20_r8)
     end do
 

src/biogeophys/SoilStateType.F90

@@ -24,12 +24,15 @@ module SoilStateType
 
      ! sand/ clay/ organic matter
      real(r8), pointer :: sandfrac_patch       (:)   ! patch sand fraction
-     real(r8), pointer :: clayfrac_patch       (:)   ! patch clay fraction
+     real(r8), pointer :: clayfrac_patch       (:)   ! patch clay fraction 
      real(r8), pointer :: mss_frc_cly_vld_col  (:)   ! col mass fraction clay limited to 0.20
      real(r8), pointer :: cellorg_col          (:,:) ! col organic matter for gridcell containing column (1:nlevsoi)
      real(r8), pointer :: cellsand_col         (:,:) ! sand value for gridcell containing column (1:nlevsoi)
      real(r8), pointer :: cellclay_col         (:,:) ! clay value for gridcell containing column (1:nlevsoi)
      real(r8), pointer :: bd_col               (:,:) ! col bulk density of dry soil material [kg/m^3] (CN)
+     !####################### for new dust emission scheme  -dmleung ############################
+     real(r8), pointer :: roughfct_patch       (:)   ! roughness factor, 17 dec 2021
+     !###########################################################################################
 
      ! hydraulic properties
      real(r8), pointer :: hksat_col            (:,:) ! col hydraulic conductivity at saturation (mm H2O /s) 
@@ -131,6 +134,9 @@ subroutine InitAllocate(this, bounds)
     allocate(this%cellclay_col         (begc:endc,nlevsoi))             ; this%cellclay_col         (:,:) = nan 
     allocate(this%bd_col               (begc:endc,nlevgrnd))            ; this%bd_col               (:,:) = nan
 
+    !################ dmleung added 14 Dec 2021 ########################
+    allocate(this%roughfct_patch       (begp:endp))                     ; this%roughfct_patch       (:)   = nan
+    !###################################################################
     allocate(this%hksat_col            (begc:endc,nlevgrnd))            ; this%hksat_col            (:,:) = spval
     allocate(this%hksat_min_col        (begc:endc,nlevgrnd))            ; this%hksat_min_col        (:,:) = spval
     allocate(this%hk_l_col             (begc:endc,nlevgrnd))            ; this%hk_l_col             (:,:) = nan   

src/cpl/lnd_import_export.F90

@@ -206,7 +206,13 @@ subroutine lnd_import( bounds, x2l, glc_present, atm2lnd_inst, glc2lnd_inst)
 
        ! Check that solar, specific-humidity and LW downward aren't negative
        if ( atm2lnd_inst%forc_lwrad_not_downscaled_grc(g) <= 0.0_r8 )then
-          call endrun( sub//' ERROR: Longwave down sent from the atmosphere model is negative or zero' )
+          !####added by dmleung to avoid negative forc_lwrad_not_downscaled_grc(g) 17 Dec 2021###########
+          write(iulog,*) 'FLDS <= 0; FLDS = ', atm2lnd_inst%forc_lwrad_not_downscaled_grc(g)
+          write(iulog,*) 'gridcell index = ', g
+          atm2lnd_inst%forc_lwrad_not_downscaled_grc(g) = 0.1_r8
+          write(iulog,*) 'FLDS set to 0.1'
+          !call endrun( sub//' ERROR: Longwave down sent from the atmosphere model is negative or zero' )
+          !##### should not be here but I am not sure why sometimes LW radiation is negative yet #######
        end if
        if ( (atm2lnd_inst%forc_solad_grc(g,1) < 0.0_r8) .or.  (atm2lnd_inst%forc_solad_grc(g,2) < 0.0_r8) &
        .or. (atm2lnd_inst%forc_solai_grc(g,1) < 0.0_r8) .or.  (atm2lnd_inst%forc_solai_grc(g,2) < 0.0_r8) ) then

src/main/clm_varctl.F90

@@ -92,6 +92,7 @@ module clm_varctl
   character(len=fname_len), public :: nrevsn     = ' '        ! restart data file name for branch run
   character(len=fname_len), public :: fsnowoptics  = ' '      ! snow optical properties file name
   character(len=fname_len), public :: fsnowaging   = ' '      ! snow aging parameters file name
+  character(len=fname_len), public :: rough_fct   = ' '       ! roughness factor, dmleung, added 17 Dec 2021
 
   !----------------------------------------------------------
   ! Flag to read ndep rather than obtain it from coupler

src/main/controlMod.F90

@@ -262,6 +262,9 @@ subroutine control_init( )
          use_grainproduct, use_snicar_frc, use_vancouver, use_mexicocity, use_noio, &
          use_nguardrail
 
+    ! dust emission, dmleung 14 Dec 2021
+    namelist /clm_inparm/  &
+         rough_fct    !-dmleung 17 Dec 2021
 
     ! ----------------------------------------------------------------------
     ! Default values
@@ -606,6 +609,9 @@ subroutine control_spmd()
     call mpi_bcast (fsnowoptics, len(fsnowoptics),  MPI_CHARACTER, 0, mpicom, ier)
     call mpi_bcast (fsnowaging,  len(fsnowaging),   MPI_CHARACTER, 0, mpicom, ier)
 
+    ! initialize input data for new dust emission module dmleung 14 Dec 2021
+    call mpi_bcast (rough_fct, len(rough_fct), MPI_CHARACTER, 0, mpicom, ier)! added by dmleung, 17 Dec 2021
+
     ! Irrigation
     call mpi_bcast(irrigate, 1, MPI_LOGICAL, 0, mpicom, ier)
 
@@ -823,6 +829,13 @@ subroutine control_print ()
     else
        write(iulog,*) '   land frac data = ',trim(fatmlndfrc)
     end if
+    !##### for dmleung's input data for new dust emission module #####
+    if (rough_fct == ' ') then    ! -dmleung, 17 Dec 2021
+       write(iulog,*) '   rough_fct surface dataset not set'
+    else
+       write(iulog,*) '   surface data   = ',trim(rough_fct)
+    end if
+    !#################################################################
     if (use_cn) then
        if (suplnitro /= suplnNon)then
           write(iulog,*) '   Supplemental Nitrogen mode is set to run over Patches: ', &

src/main/lnd2atmMod.F90

@@ -115,6 +115,14 @@ subroutine lnd2atm_minimal(bounds, &
          p2c_scale_type='unity', c2l_scale_type= 'urbanf', l2g_scale_type='unity')
 
     do g = bounds%begg,bounds%endg
+       !#####added by dmleung to avoid negative eflx_lwrad_out_grc(g) 11 Dec 2021 ###################
+       if (lnd2atm_inst%eflx_lwrad_out_grc(g) < 0) then
+          write(iulog,*) 'FIRE = ', lnd2atm_inst%eflx_lwrad_out_grc(g)
+          write(iulog,*) 'gridcell index = ', g
+          !write(iulog,*) 'sb constant = ', sb
+          lnd2atm_inst%eflx_lwrad_out_grc(g) = 100_r8
+       end if
+       !##### should not be here but I am not sure why sometimes LW radiation is negative yet #######
        lnd2atm_inst%t_rad_grc(g) = sqrt(sqrt(lnd2atm_inst%eflx_lwrad_out_grc(g)/sb))
     end do