From e23347add05545305756e51bc18461cd162a654a Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Fri, 15 Apr 2022 14:42:33 -0600 Subject: [PATCH 01/13] adding three new namelist items: clay fraction, roughness factor, and LULC area fractions --- src/main/controlMod.F90 | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index d46df04a5d..b8d0a3576a 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -259,6 +259,10 @@ 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/ & + clay_frc, & !-dmleung 14 Dec 2021 + rough_fct, lulc_frc !-dmleung 17 Dec 2021 ! ---------------------------------------------------------------------- ! Default values @@ -603,6 +607,11 @@ 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 (clay_frc, len(clay_frc), MPI_CHARACTER, 0, mpicom, ier) ! added by dmleung, 14 Dec 2021 + call mpi_bcast (rough_fct, len(rough_fct), MPI_CHARACTER, 0, mpicom, ier)! added by dmleung, 17 Dec 2021 + call mpi_bcast (lulc_frc, len(lulc_frc), MPI_CHARACTER, 0, mpicom, ier) ! added by dmleung, 17 Dec 2021 + ! Irrigation call mpi_bcast(irrigate, 1, MPI_LOGICAL, 0, mpicom, ier) @@ -817,6 +826,23 @@ subroutine control_print () else write(iulog,*) ' land frac data = ',trim(fatmlndfrc) end if + !##### for dmleung's input data for new dust emission module ##### + if (clay_frc == ' ') then ! -dml, 2 Mar 2021, added 14 Dec 2021 + write(iulog,*) ' clay_frc surface dataset not set' + else + write(iulog,*) ' surface data = ',trim(clay_frc) + end if + if (rough_fct == ' ') then ! -dmleung, 17 Dec 2020 + write(iulog,*) ' rough_fct surface dataset not set' + else + write(iulog,*) ' surface data = ',trim(rough_fct) + end if + if (lulc_frc == ' ') then ! -dmleung, 17 Dec 2021 + write(iulog,*) ' lulc_frc surface dataset not set' + else + write(iulog,*) ' surface data = ',trim(lulc_frc) + end if + !################################################################# if (use_cn) then if (suplnitro /= suplnNon)then write(iulog,*) ' Supplemental Nitrogen mode is set to run over Patches: ', & From 92497fbbb517dc68e06b48f633ef99951f2c4e91 Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Fri, 15 Apr 2022 14:50:04 -0600 Subject: [PATCH 02/13] Define the three namelist variables for reading in the datasets. --- .../namelist_definition_clm4_5.xml | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/bld/namelist_files/namelist_definition_clm4_5.xml b/bld/namelist_files/namelist_definition_clm4_5.xml index 09af1a1756..806e39e8e2 100644 --- a/bld/namelist_files/namelist_definition_clm4_5.xml +++ b/bld/namelist_files/namelist_definition_clm4_5.xml @@ -2482,4 +2482,25 @@ use case.) + + + +Full pathname of time-invariant clay fraction dataset for calculating dust emission. + + + + + +Full pathname of time-invariant roughness factor dataset for calculating drag partition of dust emission. + + + + +Full pathname of time-invariant GLCNMO LULC area fraction for calculating hybrid drag partition. + + + From a7a95d2decd1ef1a500b72c5920c67a5a7e60335 Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Mon, 18 Apr 2022 12:38:52 -0600 Subject: [PATCH 03/13] Declare & initiate pointer variables for the read-in datasets --- src/biogeophys/SoilStateType.F90 | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/biogeophys/SoilStateType.F90 b/src/biogeophys/SoilStateType.F90 index 18baac4976..ea0755f3bf 100644 --- a/src/biogeophys/SoilStateType.F90 +++ b/src/biogeophys/SoilStateType.F90 @@ -30,6 +30,13 @@ module SoilStateType 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 :: clayfrc_patch (:) ! alternative clay fraction, 14 dec 2021 + real(r8), pointer :: roughfct_patch (:) ! roughness factor, 17 dec 2021 + real(r8), pointer :: rockfrc_patch (:) ! GLCNMO rock/bare area fraction, 17 dec 2021 + real(r8), pointer :: vegefrc_patch (:) ! GLCNMO shrub/crop area fraction, 17 dec 2021 + real(r8), pointer :: sparfrc_patch (:) ! GLCNMO sparse vegetation area fraction, 17 dec 2021 + !########################################################################################### ! hydraulic properties real(r8), pointer :: hksat_col (:,:) ! col hydraulic conductivity at saturation (mm H2O /s) @@ -131,6 +138,13 @@ 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%clayfrc_patch (begp:endp)) ; this%clayfrc_patch (:) = nan + allocate(this%roughfct_patch (begp:endp)) ; this%roughfct_patch (:) = nan + allocate(this%rockfrc_patch (begp:endp)) ; this%rockfrc_patch (:) = nan + allocate(this%vegefrc_patch (begp:endp)) ; this%vegefrc_patch (:) = nan + allocate(this%sparfrc_patch (begp:endp)) ; this%sparfrc_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 From fdf0192a6ffd184a6df473209e5404133f704891 Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Mon, 18 Apr 2022 12:44:54 -0600 Subject: [PATCH 04/13] Use the pio library to read the data files. Also, change the gravimetric soil moisture threshold to the equation by Fecan et al. (1999) --- src/biogeophys/SoilStateInitTimeConstMod.F90 | 63 +++++++++++++++++++- 1 file changed, 61 insertions(+), 2 deletions(-) diff --git a/src/biogeophys/SoilStateInitTimeConstMod.F90 b/src/biogeophys/SoilStateInitTimeConstMod.F90 index a157e16bb6..d9aae77d91 100644 --- a/src/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/biogeophys/SoilStateInitTimeConstMod.F90 @@ -100,7 +100,9 @@ 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, & + clay_frc, & ! -dmleung added clay_frc to CESM2/CLM5 14 Dec 2021 + rough_fct, lulc_frc ! -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 +153,11 @@ 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 :: clayfrc2d(:) ! read in - clay percentage (called fraction here) from SoilGrids, dmleung added 14 Dec 2021 + real(r8) ,pointer :: roughfct2d(:) ! read in - time-invariant roughness factor, dmleung added 17 Dec 2021 + real(r8) ,pointer :: rockfrc2d(:) ! read in - rock area fraction from GLCNMO, dmleung added 17 Dec 2021 + real(r8) ,pointer :: vegefrc2d(:) ! read in - shrub/crop area fraction from GLCNMO, dmleung added 17 Dec 2021 + real(r8) ,pointer :: sparfrc2d(:) ! read in - sparse vegetation area fraction from GLCNMO, dmleung added 17 Dec 2021 character(len=256) :: locfn ! local filename integer :: ipedof integer :: begp, endp @@ -297,6 +304,56 @@ 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 14 Dec 2021 + allocate(clayfrc2d(begg:endg)) ! 2 Mar 2021 + ! read SoilGrids clay fraction, 2 Mar 2021 + write(iulog,*) 'Attempting to read alternative (e.g., SoilGrids) clay fraction data, by dmleung .....' + call getfil (clay_frc, locfn, 0) + call ncd_pio_openfile (ncid, locfn, 0) + call ncd_io(ncid=ncid, varname='f_clay', flag='read', data=clayfrc2d, dim1name=grlnd, readvar=readvar) + + do p = begp,endp + g = patch%gridcell(p) + soilstate_inst%clayfrc_patch(p) = clayfrc2d(g)/100_r8 + end do + + call ncd_pio_closefile(ncid) + + ! 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,*) 'dimension of roughness = ', shape(roughfct2d) ! -dml + 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 + + ! dmleung added to CESM2/CLM5 17 Dec 2021 + allocate(rockfrc2d(begg:endg));allocate(vegefrc2d(begg:endg));allocate(sparfrc2d(begg:endg)) ! 7 Jul 2021 + write(iulog,*) 'Attempting to read GLCNMO LULC area fraction data, by dmleung .....' + call getfil (lulc_frc, locfn, 0) + call ncd_pio_openfile (ncid, locfn, 0) + call ncd_io(ncid=ncid, varname='A_r', flag='read', data=rockfrc2d, dim1name=grlnd, readvar=readvar) + call ncd_io(ncid=ncid, varname='A_v', flag='read', data=vegefrc2d, dim1name=grlnd, readvar=readvar) + call ncd_io(ncid=ncid, varname='A_s', flag='read', data=sparfrc2d, dim1name=grlnd, readvar=readvar) + + do p = begp,endp + g = patch%gridcell(p) + soilstate_inst%rockfrc_patch(p) = rockfrc2d(g) + soilstate_inst%vegefrc_patch(p) = vegefrc2d(g) + soilstate_inst%sparfrc_patch(p) = sparfrc2d(g) + end do + + call ncd_pio_closefile(ncid) + !################################################################################## + ! -------------------------------------------------------------------- ! get original soil depths to be used in interpolation of sand and clay ! -------------------------------------------------------------------- @@ -614,7 +671,9 @@ 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.01_r8*(0.17_r8*clayfrc2d(g) + 0.0014_r8*clayfrc2d(g)*clayfrc2d(g)) !Fecan et al. (1999) -jfk, dmleung coded 14 Dec 2021 for SoilGrids 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 From 31ed4fc2c93a1e30cffa56a7b271a11d6415b670 Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Mon, 18 Apr 2022 13:59:26 -0600 Subject: [PATCH 05/13] Add entries to run input data files --- src/main/clm_varctl.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index ef8c25761e..ae9077e4dd 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -92,6 +92,9 @@ 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 :: clay_frc = ' ' ! alternative clay fraction, dmleung, added 14 Dec 2021 + character(len=fname_len), public :: rough_fct = ' ' ! roughness factor, dmleung, added 17 Dec 2021 + character(len=fname_len), public :: lulc_frc = ' ' ! LULC area fraction, dmleung, added 17 Dec 2021 !---------------------------------------------------------- ! Flag to read ndep rather than obtain it from coupler From a80d69ac7c0d7ea5e39d2f788ba998a0ed662b00 Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Mon, 18 Apr 2022 14:01:07 -0600 Subject: [PATCH 06/13] All the changes of dust emission scheme are here. Major changes include: 1. Use Shao et al. (2000) dust emissions threshold scheme instead of Iversen and White (1982) scheme. 2. Use Kok et al. (2014) dust emission scheme instead of Zender et al. (2003) scheme. 3. Change the optimal soil diameter of 75 micrometer to 130 micrometer. 4. Change the threshold LAI value of dust emission (from 0.3) to 1 (or 0.5 is also acceptable). 5. Add the drag partition effect due to rocks (using the roughness factor input data) and plants (using CLM LAI). 6. Remove Owen effect which is less relevant to the dust emission process. 7. Add the dust emission intermittency effect due to the PBL turbulence. Other code changes: 8. Add multiple variables for output. 9. Pull out some major dust emission calculations out of the if loop (if lnd_frc_mbl > 1) so that variables are outputted regardless of the possibility of saltation and dust emission. --- src/biogeochem/DUSTMod.F90 | 438 ++++++++++++++++++++++++++++++++++--- 1 file changed, 408 insertions(+), 30 deletions(-) diff --git a/src/biogeochem/DUSTMod.F90 b/src/biogeochem/DUSTMod.F90 index ad2e5bf6a9..cb57c657e9 100644 --- a/src/biogeochem/DUSTMod.F90 +++ b/src/biogeochem/DUSTMod.F90 @@ -60,7 +60,30 @@ module DUSTMod 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 - + !########### added by dmleung 27 Nov 2021 ######################################################################## + real(r8), pointer, private :: dst_emiss_coeff_patch (:) ! dust emission coefficient (unitless) + real(r8), pointer, private :: wnd_frc_thr_patch (:) ! wet fluid threshold (m/s) + real(r8), pointer, private :: wnd_frc_thr_dry_patch (:) ! dry fluid threshold (m/s) + real(r8), pointer, private :: lnd_frc_mble_patch (:) ! land mobile fraction -dmleung + real(r8), pointer, private :: liq_frac_patch (:) ! liquid fraction of total water + real(r8), pointer, private :: wnd_frc_soil_patch (:) ! soil wind friction velocity (m/s) + real(r8), pointer, private :: gwc_patch (:) ! gravimetric water content (kg/kg) + !########### added by dmleung 2 Dec 2021 ######################################################################### + real(r8), pointer, private :: intrmtncy_fct_patch (:) ! intermittency factor, accounting for turbulence shutting down dust emissions (unitless) + real(r8), pointer, private :: stblty_patch (:) ! stability parameter for checking stability condition (stblty < 0 is unstable atmosphere) + real(r8), pointer, private :: u_mean_slt_patch (:) ! wind speed 0.1 m level of dust saltation (m/s) + real(r8), pointer, private :: u_sd_slt_patch (:) ! sd of wind speed 0.1 m level of dust saltation (m/s) + real(r8), pointer, private :: u_fld_thr_patch (:) ! fluid threshold wind speed 0.1 m level of dust saltation (m/s) + real(r8), pointer, private :: u_impct_thr_patch (:) ! impact threshold wind speed at 0.1 m level of dust saltation (m/s) + real(r8), pointer, private :: thr_crs_rate_patch (:) ! threshold crossing rate (unitless) + real(r8), pointer, private :: prb_crs_fld_thr_patch (:) ! probability of wind speed crossing fluid threshold + real(r8), pointer, private :: prb_crs_impct_thr_patch (:) ! probability of wind speed crossing impact threshold + !########### added by dmleung 17 Dec 2021 ######################################################################## + real(r8), pointer, private :: ustar_patch (:) ! output friction velocity for SP mode (m/s) + !########### added by dmleung 20 Dec 2021 ######################################################################## + real(r8), pointer, private :: ssr_patch (:) ! [dimless] integrated shear stress ratiio, defined by Okin (2008) and then integrated by Caroline Pierre et al. (2014) + real(r8), pointer, private :: lai_patch (:) ! [m2 leaf /m2 land] LAI+SAI for calculating Okin's drag partition, averaged to landunit level + real(r8), pointer, private :: frc_thr_rghn_fct_patch (:) ! [dimless] hybrid drag partition (or called roughness) factor contains procedure , public :: Init @@ -113,7 +136,30 @@ subroutine InitAllocate(this, bounds) 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 - + !#### added by dmleung 27 Nov 2021 ##################################### + allocate(this%dst_emiss_coeff_patch (begp:endp)) ; this%dst_emiss_coeff_patch (:) = nan + allocate(this%wnd_frc_thr_patch (begp:endp)) ; this%wnd_frc_thr_patch (:) = nan + allocate(this%wnd_frc_thr_dry_patch (begp:endp)) ; this%wnd_frc_thr_dry_patch (:) = nan + allocate(this%lnd_frc_mble_patch (begp:endp)) ; this%lnd_frc_mble_patch (:) = nan + allocate(this%wnd_frc_soil_patch (begp:endp)) ; this%wnd_frc_soil_patch (:) = nan + allocate(this%gwc_patch (begp:endp)) ; this%gwc_patch (:) = nan + allocate(this%liq_frac_patch (begp:endp)) ; this%liq_frac_patch (:) = nan + !#### added by dmleung 2 Dec 2021 ###################################### + allocate(this%intrmtncy_fct_patch (begp:endp)) ; this%intrmtncy_fct_patch (:) = nan + allocate(this%stblty_patch (begp:endp)) ; this%stblty_patch (:) = nan + allocate(this%u_mean_slt_patch (begp:endp)) ; this%u_mean_slt_patch (:) = nan + allocate(this%u_sd_slt_patch (begp:endp)) ; this%u_sd_slt_patch (:) = nan + allocate(this%u_fld_thr_patch (begp:endp)) ; this%u_fld_thr_patch (:) = nan + allocate(this%u_impct_thr_patch (begp:endp)) ; this%u_impct_thr_patch (:) = nan + allocate(this%thr_crs_rate_patch (begp:endp)) ; this%thr_crs_rate_patch (:) = nan + allocate(this%prb_crs_fld_thr_patch (begp:endp)) ; this%prb_crs_fld_thr_patch (:) = nan + allocate(this%prb_crs_impct_thr_patch (begp:endp)) ; this%prb_crs_impct_thr_patch (:) = nan + !#### added by dmleung 17 Dec 2021 ###################################### + allocate(this%ustar_patch (begp:endp)) ; this%ustar_patch (:) = nan + !#### added by dmleung 17 Dec 2021 ###################################### + allocate(this%ssr_patch (begp:endp)) ; this%ssr_patch (:) = nan + allocate(this%lai_patch (begp:endp)) ; this%lai_patch (:) = nan + allocate(this%frc_thr_rghn_fct_patch (begp:endp)) ; this%frc_thr_rghn_fct_patch (:) = nan end subroutine InitAllocate !------------------------------------------------------------------------ @@ -158,6 +204,92 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='turbulent deposition velocity 4', & ptr_patch=this%vlc_trb_4_patch, default='inactive') + !#####added by dmleung 27 Nov 2021######################################### + this%dst_emiss_coeff_patch(begp:endp) = spval + call hist_addfld1d (fname='C_d', units='dimensionless', & + avgflag='A', long_name='dust emission coefficient', & + ptr_patch=this%dst_emiss_coeff_patch, set_lake=0._r8, set_urb=0._r8) + this%wnd_frc_thr_patch(begp:endp) = spval + call hist_addfld1d (fname='WND_FRC_FT', units='m/s', & + avgflag='A', long_name='fluid threshold friction velocity', & + ptr_patch=this%wnd_frc_thr_patch, set_lake=0._r8, set_urb=0._r8) + this%wnd_frc_thr_dry_patch(begp:endp) = spval + call hist_addfld1d (fname='WND_FRC_FT_DRY', units='m/s', & + avgflag='A', long_name='dry fluid threshold friction velocity', & + ptr_patch=this%wnd_frc_thr_dry_patch, set_lake=0._r8, set_urb=0._r8) + this%wnd_frc_soil_patch(begp:endp) = spval + call hist_addfld1d (fname='WND_FRC_SOIL', units='m/s', & + avgflag='A', long_name='soil surface wind friction velocity', & + ptr_patch=this%wnd_frc_soil_patch, set_lake=0._r8, set_urb=0._r8) + this%lnd_frc_mble_patch(begp:endp) = spval + call hist_addfld1d (fname='LND_FRC_MBLE', units='dimensionless', & + avgflag='A', long_name='land mobile fraction', & + ptr_patch=this%lnd_frc_mble_patch, set_lake=0._r8, set_urb=0._r8) + this%gwc_patch(begp:endp) = spval + call hist_addfld1d (fname='GWC', units='kg/kg', & + avgflag='A', long_name='gravimetric water content', & + ptr_patch=this%gwc_patch, set_lake=0._r8, set_urb=0._r8) + this%liq_frac_patch(begp:endp) = spval + call hist_addfld1d (fname='LIQ_FRAC', units='dimensionless', & + avgflag='A', long_name='fraction of total water that is liquid', & + ptr_patch=this%liq_frac_patch, set_lake=0._r8, set_urb=0._r8) + !#####added by dmleung 2 Dec 2021 ######################################### + this%u_mean_slt_patch(begp:endp) = spval + call hist_addfld1d (fname='U_S_MEAN', units='m/s', & + avgflag='A', long_name='mean wind velocity at saltation level', & + ptr_patch=this%u_mean_slt_patch, set_lake=0._r8, set_urb=0._r8) + this%u_sd_slt_patch(begp:endp) = spval + call hist_addfld1d (fname='U_S_SIGMA', units='m/s', & + avgflag='A', long_name='sd of wind velocity at saltation level', & + ptr_patch=this%u_sd_slt_patch, set_lake=0._r8, set_urb=0._r8) + this%stblty_patch(begp:endp) = spval + call hist_addfld1d (fname='ZETA', units='', & + avgflag='A', long_name='stability parameter', & + ptr_patch=this%stblty_patch, set_lake=0._r8, set_urb=0._r8) + this%u_fld_thr_patch(begp:endp) = spval + call hist_addfld1d (fname='U_FT', units='m/s', & + avgflag='A', long_name='fluid threshold velocity at saltation level', & + ptr_patch=this%u_fld_thr_patch, set_lake=0._r8, set_urb=0._r8) + this%u_impct_thr_patch(begp:endp) = spval + call hist_addfld1d (fname='U_IT', units='m/s', & + avgflag='A', long_name='impact threshold velocity at saltation level', & + ptr_patch=this%u_impct_thr_patch, set_lake=0._r8, set_urb=0._r8) + this%thr_crs_rate_patch(begp:endp) = spval + call hist_addfld1d (fname='ALPHA', units='', & + avgflag='A', long_name='threshold crossing rate', & + ptr_patch=this%thr_crs_rate_patch, set_lake=0._r8, set_urb=0._r8) + this%prb_crs_fld_thr_patch(begp:endp) = spval + call hist_addfld1d (fname='P_FT', units='', & + avgflag='A', long_name='probability of crossing fluid threshold', & + ptr_patch=this%prb_crs_fld_thr_patch, set_lake=0._r8, set_urb=0._r8) + this%prb_crs_impct_thr_patch(begp:endp) = spval + call hist_addfld1d (fname='P_IT', units='', & + avgflag='A', long_name='probability of crossing impact threshold', & + ptr_patch=this%prb_crs_impct_thr_patch, set_lake=0._r8, set_urb=0._r8) + this%intrmtncy_fct_patch(begp:endp) = spval + call hist_addfld1d (fname='ETA', units='', & + avgflag='A', long_name='intermittency factor', & + ptr_patch=this%intrmtncy_fct_patch, set_lake=0._r8, set_urb=0._r8) + !#####added by dmleung 2 Dec 2021 ######################################### + this%ustar_patch(begp:endp) = spval + call hist_addfld1d (fname='USTAR', units='m/s', & + avgflag='A', long_name='friction velocity', & + ptr_patch=this%ustar_patch, set_lake=0._r8, set_urb=0._r8) + !#####added by dmleung 20 Dec 2021 ######################################## + this%ssr_patch(begp:endp) = spval + call hist_addfld1d (fname='SSR', units='m/s', & + avgflag='A', long_name='Okin-Pierre shear stress ratio', & + ptr_patch=this%ssr_patch, set_lake=0._r8, set_urb=0._r8) + this%lai_patch(begp:endp) = spval + call hist_addfld1d (fname='LAI', units='m/s', & + avgflag='A', long_name='landunit-mean LAI for Okin-Pierre scheme', & + ptr_patch=this%lai_patch, set_lake=0._r8, set_urb=0._r8) + this%frc_thr_rghn_fct_patch(begp:endp) = spval + call hist_addfld1d (fname='FRC_THR_RGHN_FCT', units='dimensionless', & + avgflag='A', long_name='hybrid drag partition (or roughness) factor', & + ptr_patch=this%frc_thr_rghn_fct_patch, set_lake=0._r8, set_urb=0._r8) + !########################################################################## + end subroutine InitHistory !----------------------------------------------------------------------- @@ -223,7 +355,7 @@ subroutine DustEmission (bounds, & real(r8) :: flx_mss_vrt_dst_ttl(bounds%begp:bounds%endp) real(r8) :: frc_thr_wet_fct real(r8) :: frc_thr_rgh_fct - real(r8) :: wnd_frc_thr_slt + !real(r8) :: wnd_frc_thr_slt ! dmleung commented and put below 2 Dec 2021 real(r8) :: wnd_rfr_thr_slt real(r8) :: wnd_frc_slt real(r8) :: lnd_frc_mbl(bounds%begp:bounds%endp) @@ -234,12 +366,38 @@ subroutine DustEmission (bounds, & real(r8) :: sumwt(bounds%begl:bounds%endl) ! sum of weights logical :: found ! temporary for error check integer :: index + !########### added by dmleung 27 Nov 2021 ######################################### + real(r8) :: tmp2 ! calculates the dry fluid threshold using Shao and Lu (2000) scheme; replace the tmp1 (Iversen and White, 1982) that was passed from Dustini to DustEmission; tmp2 will be calculated here 23 May 2020 -dmleung + real(r8) :: wnd_frc_thr_slt_std ! [m/s] The soil threshold friction speed at standard air density (1.2250 kg/m3) -jfk + real(r8) :: frag_expt ! fragmentation exponent, -dmleung 22 Jun 2021 + !########### added by dmleung 2 Dec 2021 for intermittency scheme ################# + real(r8) :: wnd_frc_thr_slt_it ! [m/s] created for impact threshold friction velocity, dmleung 9 Jun 2021 + real(r8) :: wnd_frc_thr_slt ! [m/s] used for wet fluid threshold friction velocity, dmleung 9 Jun 2021 + !########### added by dmleung 20 Dec 2021 for drag partition effect ################# + real(r8) :: K_length ! [dimless] normalized mean interobstacle distance, or called gap length (Okin, 2008) ! ! constants ! real(r8), parameter :: cst_slt = 2.61_r8 ! [frc] Saltation constant real(r8), parameter :: flx_mss_fdg_fct = 5.0e-4_r8 ! [frc] Empir. mass flx tuning eflx_lh_vegt - real(r8), parameter :: vai_mbl_thr = 0.3_r8 ! [m2 m-2] VAI threshold quenching dust mobilization + !real(r8), parameter :: vai_mbl_thr = 0.3_r8 ! [m2 m-2] VAI threshold quenching dust mobilization + !####### added by dmleung 27 Nov 2021 ########################################################################### + real(r8), parameter :: vai_mbl_thr = 1.0_r8 ! [m2 m-2] new VAI threshold; dmleung suggests 1 or 0.5, and the default 0.3 seems a bit too small -dmleung 27 Nov 2021 + real(r8), parameter :: Cd0 = 4.4e-5_r8 ! [dimless] proportionality constant in calculation of dust emission coefficient -jfk + real(r8), parameter :: Ca = 2.7_r8 ! [dimless] proportionality constant in scaling of dust emission exponent -jfk + real(r8), parameter :: Ce = 2.0_r8 ! [dimless] proportionality constant scaling exponential dependence of dust emission coefficient on standardized soil threshold friction speed -jfk + real(r8), parameter :: C_tune = 0.05_r8 ! [dimless] global tuning constant for vertical dust flux; set to produce ~same global dust flux in control sim (I_2000) as old parameterization -jfk + real(r8), parameter :: wnd_frc_thr_slt_std_min = 0.16_r8 ! [m/s] minimum standardized soil threshold friction speed -jfk + real(r8), parameter :: forc_rho_std = 1.2250_r8 ! [kg/m3] density of air at standard pressure (101325) and temperature (293 K) -jfk + real(r8), parameter :: dns_slt = 2650.0_r8 ! [kg m-3] Density of optimal saltation particles, dml 23 May 2020 + !####### added by dmleung 2 Dec 2021 for intermittency ########################################################## + real(r8), parameter :: B_it = 0.82_r8 ! [dimless] ratio = u_star_it / u_star_ft0 (may need to change into a fn of moisture later on) -dml + real(r8), parameter :: k = 0.4_r8 ! [dimless] von Karman constant -dml + !####### added by dmleung 2 Dec 2021 for Okin (2008) drag partition for plants ########################################################## + real(r8), parameter :: f_0 = 0.32_r8 ! [dimless] SSR in the immediate lee of a plant, dimensionless + real(r8), parameter :: c_e = 4.8_r8 ! [dimless] e-folding distance velocity recovery, dimensionless + !################################################################################################################ + !################################################################################################################ !------------------------------------------------------------------------ associate( & @@ -262,7 +420,35 @@ subroutine DustEmission (bounds, & 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) + flx_mss_vrt_dst_tot => dust_inst%flx_mss_vrt_dst_tot_patch , & ! Output: [real(r8) (:) ] total dust flux back to atmosphere (pft) + ! the following are added by dmleung 27 Nov 2021 + dst_emiss_coeff => dust_inst%dst_emiss_coeff_patch , & ! Output dust emission coefficient + wnd_frc_thr => dust_inst%wnd_frc_thr_patch , & ! output impact threshold -dmleung + wnd_frc_thr_dry => dust_inst%wnd_frc_thr_dry_patch , & ! output dry threshold + lnd_frc_mble => dust_inst%lnd_frc_mble_patch , & ! -dmleung, 3 Feb 2020 + wnd_frc_soil => dust_inst%wnd_frc_soil_patch , & ! soil friction velocity u_*s = (u_*)(f_eff) + gwc => dust_inst%gwc_patch , & ! output gravimetric water content + liq_frac => dust_inst%liq_frac_patch , & + ! added by dmleung 8 Jul 2019, recoded 2 Dec 2021 + intrmtncy_fct => dust_inst%intrmtncy_fct_patch , & + stblty => dust_inst%stblty_patch , & + u_mean_slt => dust_inst%u_mean_slt_patch , & + u_sd_slt => dust_inst%u_sd_slt_patch , & + u_fld_thr => dust_inst%u_fld_thr_patch , & + u_impct_thr => dust_inst%u_impct_thr_patch , & + thr_crs_rate => dust_inst%thr_crs_rate_patch , & + prb_crs_fld_thr => dust_inst%prb_crs_fld_thr_patch , & + prb_crs_impct_thr => dust_inst%prb_crs_impct_thr_patch , & + ! added by dmleung 17 Dec 2021 + roughfct => soilstate_inst%roughfct_patch , & + rockfrc => soilstate_inst%rockfrc_patch , & + vegefrc => soilstate_inst%vegefrc_patch , & + sparfrc => soilstate_inst%sparfrc_patch , & + ustar => dust_inst%ustar_patch , & ! Output friction velocity for SP mode + ! added by dmleung 20 Dec 2021 + ssr => dust_inst%ssr_patch , & + lai => dust_inst%lai_patch , & + frc_thr_rghn_fct => dust_inst%frc_thr_rghn_fct_patch & ) ttlai(bounds%begp : bounds%endp) = 0._r8 @@ -340,6 +526,30 @@ subroutine DustEmission (bounds, & do fp = 1,num_nolakep p = filter_nolakep(fp) flx_mss_vrt_dst_tot(p) = 0.0_r8 + ! the following are added by dmleung 27 Nov 2021 + dst_emiss_coeff(p) = 0.0_r8 + wnd_frc_thr(p) = 0.0_r8 + wnd_frc_thr_dry(p) = 0.0_r8 + lnd_frc_mble(p) = 0.0_r8 + wnd_frc_soil(p) = 0.0_r8 + gwc(p) = 0.0_r8 + liq_frac(p) = 0.0_r8 + ! dmleung's edit, 8 Jul 2019; added by dmleung 2 Dec 2021 + u_mean_slt(p) = 0.0_r8 + u_sd_slt(p) = 0.0_r8 + stblty(p) = 0.0_r8 + u_fld_thr(p) = 0.0_r8 + u_impct_thr(p) = 0.0_r8 + thr_crs_rate(p) = 0.0_r8 + prb_crs_fld_thr(p) = 0.0_r8 + prb_crs_impct_thr(p) = 0.0_r8 + intrmtncy_fct(p) = 0.0_r8 + ! dmleung's edit for including friction velcoity for SP mode output, 17 Dec 2021 + ustar(p) = 0.0_r8 + ! dmleung's edit, 20 Dec 2021 + ssr(p) = 0.0_r8 + lai(p) = 0.0_r8 + frc_thr_rghn_fct(p) = 0.0_r8 end do do n = 1, ndst do fp = 1,num_nolakep @@ -354,6 +564,101 @@ subroutine DustEmission (bounds, & l = patch%landunit(p) g = patch%gridcell(p) + !################################################################################################ + ! put dust emission calculation here to output threshold friction velocity for the whole globe, + ! not just when lnd_frc_mbl = 0. Edited by dmleung 27 Nov 2021 + bd = (1._r8-watsat(c,1))*2.7e3_r8 ![kg m-3] Bulk density of dry surface soil + gwc_sfc = h2osoi_vol(c,1)*SHR_CONST_RHOFW/bd ![kg kg-1] Gravimetric H2O cont + if (gwc_sfc > gwc_thr(c)) then + frc_thr_wet_fct = sqrt(1.0_r8 + 1.21_r8 * (100.0_r8*(gwc_sfc - gwc_thr(c)))**0.68_r8) + else + frc_thr_wet_fct = 1.0_r8 + end if + + ! output moisture variables -dmleung, coded Jul 2020, recoded 18 Mar 2021, added to CLM5 27 Nov 2021 + gwc(p) = gwc_sfc ! output surface gravimetric water content + + ! slevis: adding liqfrac here, because related to effects from soil water + liqfrac = max( 0.0_r8, min( 1.0_r8, h2osoi_liq(c,1) / (h2osoi_ice(c,1)+h2osoi_liq(c,1)+1.0e-6_r8) ) ) !-dmleung 27 Nov 2021 + ! output liquid fraction -dmleung 27 Nov 2021 + liq_frac(p) = liqfrac + + !####################################################################################################### + ! calculate Shao & Lu (2000) dust emission threshold scheme here + ! use tmp1 from DUSTini for Iversen and White I&W (1982) (75 um is optimal); use tmp2 for S&L (2000) (107 um is optimal) + ! recoded to CLM5 27 Nov 2021 + !####################################################################################################### + + tmp2 = 1.0_r8*sqrt(0.0123_r8 * (dns_slt*grav*140.0e-6_r8 + 5.0e-4_r8/140.0e-6_r8)) ! calculate S&L (2000) scheme here for threshold + wnd_frc_thr_dry(p) = tmp2 / sqrt(forc_rho(c)) ! output dry fluid threshold + wnd_frc_thr_slt = tmp2 / sqrt(forc_rho(c)) * frc_thr_wet_fct !* frc_thr_rgh_fct ! use as threshold in this module + wnd_frc_thr_slt_it = B_it * tmp2 / sqrt(forc_rho(c)) ! define impact threshold -dml 9 Jun 2021, recoded to CLM5 27 Nov 2021 + + !wnd_frc_thr_dry(p) = tmp1 / sqrt(forc_rho(c)) ! output dry fluid threshold + !wnd_frc_thr_slt = tmp1 / sqrt(forc_rho(c)) * frc_thr_wet_fct !* frc_thr_rgh_fct ! use as threshold in this module + !wnd_frc_thr_slt_it = B_it * tmp1 / sqrt(forc_rho(c)) ! define impact threshold -dmleung 9 Jun 2021, recoded to CLM5 27 Nov 2021 + ! the above formula is true for Iversen and White (1982) and Shao and Lu (2000) scheme -dmleung, 23 Feb 2020, added to CLM5 27 Nov 2021 + wnd_frc_thr(p) = wnd_frc_thr_slt ! output fluid threshold -dmleung + + ! use emission threshold to calculate standardized threshold and dust emission coefficient dmleung 27 Nov 2021 + wnd_frc_thr_slt_std = wnd_frc_thr_slt * sqrt(forc_rho(c) / forc_rho_std) ! standardized soil threshold friction speed -jfk (defined using fluid threshold + dst_emiss_coeff(p) = Cd0 * exp(-Ce * (wnd_frc_thr_slt_std - wnd_frc_thr_slt_std_min) / wnd_frc_thr_slt_std_min) ! save dust emission coefficient here for all grids, -dml, 1 Mar 2021 + + ! framentation exponent dmleung 27 Nov 2021; moved to this block 23 Dec 2021 + frag_expt = (Ca * (wnd_frc_thr_slt_std - wnd_frc_thr_slt_std_min) / wnd_frc_thr_slt_std_min) ! fragmentation exponent, defined in Kok et al. (2014a) -dmleung 27 Nov 2021 + if (frag_expt > 3_r8) then ! set fragmentation exponent to be 3 or 5 at maximum, to avoid local AOD blowup + frag_expt = 3_r8 + end if + + !################ drag partition effect, and soil friction velocity############################ + ! subsection on computing vegetation drag partition and hybrid drag partition factors + ! in our scheme, drag partition effect is applied on the wind instead of the threshold + ! -dmleung, 7 Jul 2021 , coded to CLM5 27 Nov 2021 + !############################################################################################## + ! the following comes from subr. frc_thr_rgh_fct_get + ! purpose: compute factor by which surface roughness increases threshold + ! friction velocity (currently a constant) + + if (lnd_frc_mbl(p) > 0.0_r8 .AND. tlai_lu(l)<=1_r8) then + ! vegetation drag partition equation following Gregory Okin (2008) + Caroline Pierre et al. (2014), dmleung 20 Dec 2021 + lai(p) = tlai_lu(l) ! LAI+SAI averaged to landunit level; saved for output + if (lai(p) < 0.1_r8) then + lai(p) = 0.1_r8 ! setting LAI ~ 0.1 to be a threshold value as computing K involves 1 / LAI + end if + ! calculate Okin's shear stress ratio (which is drag partition factor) using Pierre's equation + K_length = 2_r8 * (1_r8/lai(p) - 1_r8) ! Here LAI has to be non-zero to avoid blowup + ssr(p) = (K_length+f_0*c_e)/(K_length+c_e) + + frc_thr_rgh_fct = (rockfrc(p)*(roughfct(p))**3_r8 + (vegefrc(p)+sparfrc(p))*(ssr(p))**3_r8 )**(0.3333_r8) ! land cover weighted mean using static GLCNMo bare land fraction LC0, dmleung 20 Dec 2021 + + wnd_frc_slt = fv(p) * frc_thr_rgh_fct ! wnd_frc_slt will be used in the dust emission equation -dmleung + + frc_thr_rghn_fct(p) = frc_thr_rgh_fct ! save hybrid drag partition factor, dmleung 20 Dec 2021 + else + wnd_frc_slt = fv(p) ! The value here is not important since once lnd_frc_mbl(p) <= 0.0_r8 there will be no emission. + frc_thr_rghn_fct(p) = 0.0_r8 ! save hybrid drag partition factor, dmleung 20 Dec 2021 + end if + + !##########end of drag partition effect ####################################################### + + !############ Add Owen effect; if not, comment out this block !-dmleung, 27 Nov 2021 ########### + ! the following if-block comes from subr. wnd_frc_slt_get + ! purpose: compute the saltating friction velocity + ! theory: saltation roughens the boundary layer, AKA "Owen's effect" + + !if (u10(p) >= wnd_rfr_thr_slt) then + ! wnd_rfr_dlt = u10(p) - wnd_rfr_thr_slt + ! wnd_frc_slt_dlt = 0.003_r8 * wnd_rfr_dlt * wnd_rfr_dlt + ! wnd_frc_slt = wnd_frc_slt + wnd_frc_slt_dlt ! careful that RHS is now wnd_frc_slt instead of fv(p) + ! ! because wnd_frc_slt takes drag partition effect into account, but fv(p) doesn't. dmleung 27 Nov 2021 + !end if + !########## end of Owen effect ################################################################ + + ! save soil friction velocity and roughness effect before the if-statement, -dml, 1 Mar 2021, coded to CLM5 27 Nov 2021 + wnd_frc_soil(p) = wnd_frc_slt ! save soil friction velocity for CLM output, which has drag partition and Owen effect -dml + ustar(p) = fv(p) ! save friction velocity for SP mode (use_cn=0) since only CN/BGC mode (use_cn=1) has FV output -dmleung 17 Dec 2021 + ! save land mobile fraction + lnd_frc_mble(p) = lnd_frc_mbl(p) ! save land mobile fraction first, before the if-statement, -dml, 1 Mar 2021 ! only perform the following calculations if lnd_frc_mbl is non-zero if (lnd_frc_mbl(p) > 0.0_r8) then @@ -362,7 +667,7 @@ subroutine DustEmission (bounds, & ! purpose: compute factor by which surface roughness increases threshold ! friction velocity (currently a constant) - frc_thr_rgh_fct = 1.0_r8 + !frc_thr_rgh_fct = 1.0_r8 ! the following comes from subr. frc_thr_wet_fct_get ! purpose: compute factor by which soil moisture increases threshold friction velocity @@ -370,17 +675,17 @@ subroutine DustEmission (bounds, & ! modified 4/5/2002 (slevis) to use gravimetric instead of volumetric ! water content - bd = (1._r8-watsat(c,1))*2.7e3_r8 ![kg m-3] Bulk density of dry surface soil - gwc_sfc = h2osoi_vol(c,1)*SHR_CONST_RHOFW/bd ![kg kg-1] Gravimetric H2O cont - if (gwc_sfc > gwc_thr(c)) then - frc_thr_wet_fct = sqrt(1.0_r8 + 1.21_r8 * (100.0_r8*(gwc_sfc - gwc_thr(c)))**0.68_r8) - else - frc_thr_wet_fct = 1.0_r8 - end if + !bd = (1._r8-watsat(c,1))*2.7e3_r8 ![kg m-3] Bulk density of dry surface soil + !gwc_sfc = h2osoi_vol(c,1)*SHR_CONST_RHOFW/bd ![kg kg-1] Gravimetric H2O cont + !if (gwc_sfc > gwc_thr(c)) then + ! frc_thr_wet_fct = sqrt(1.0_r8 + 1.21_r8 * (100.0_r8*(gwc_sfc - gwc_thr(c)))**0.68_r8) + !else + ! frc_thr_wet_fct = 1.0_r8 + !end if ! slevis: adding liqfrac here, because related to effects from soil water - liqfrac = max( 0.0_r8, min( 1.0_r8, h2osoi_liq(c,1) / (h2osoi_ice(c,1)+h2osoi_liq(c,1)+1.0e-6_r8) ) ) + !liqfrac = max( 0.0_r8, min( 1.0_r8, h2osoi_liq(c,1) / (h2osoi_ice(c,1)+h2osoi_liq(c,1)+1.0e-6_r8) ) ) ! the following lines come from subr. dst_mbl ! purpose: adjust threshold friction velocity to acct for moisture and @@ -388,52 +693,125 @@ subroutine DustEmission (bounds, & ! subr. wnd_frc_thr_slt_get which computes dry threshold ! friction velocity for saltation - wnd_frc_thr_slt = tmp1 / sqrt(forc_rho(c)) * frc_thr_wet_fct * frc_thr_rgh_fct + !wnd_frc_thr_slt = tmp1 / sqrt(forc_rho(c)) * frc_thr_wet_fct * frc_thr_rgh_fct ! reset these variables which will be updated in the following if-block - wnd_frc_slt = fv(p) + !wnd_frc_slt = fv(p) flx_mss_hrz_slt_ttl = 0.0_r8 flx_mss_vrt_dst_ttl(p) = 0.0_r8 ! the following line comes from subr. dst_mbl ! purpose: threshold saltation wind speed - wnd_rfr_thr_slt = u10(p) * wnd_frc_thr_slt / fv(p) + wnd_rfr_thr_slt = u10(p) * wnd_frc_thr_slt / fv(p) ! keep and use if I want Z03 scheme -dmleung ! the following if-block comes from subr. wnd_frc_slt_get ! purpose: compute the saltating friction velocity ! theory: saltation roughens the boundary layer, AKA "Owen's effect" - if (u10(p) >= wnd_rfr_thr_slt) then - wnd_rfr_dlt = u10(p) - wnd_rfr_thr_slt - wnd_frc_slt_dlt = 0.003_r8 * wnd_rfr_dlt * wnd_rfr_dlt - wnd_frc_slt = fv(p) + wnd_frc_slt_dlt - end if + !if (u10(p) >= wnd_rfr_thr_slt) then + ! wnd_rfr_dlt = u10(p) - wnd_rfr_thr_slt + ! wnd_frc_slt_dlt = 0.003_r8 * wnd_rfr_dlt * wnd_rfr_dlt + ! wnd_frc_slt = fv(p) + wnd_frc_slt_dlt + !end if ! the following comes from subr. flx_mss_hrz_slt_ttl_Whi79_get ! purpose: compute vertically integrated streamwise mass flux of particles - if (wnd_frc_slt > wnd_frc_thr_slt) then - wnd_frc_rat = wnd_frc_thr_slt / wnd_frc_slt - flx_mss_hrz_slt_ttl = cst_slt * forc_rho(c) * (wnd_frc_slt**3.0_r8) * & - (1.0_r8 - wnd_frc_rat) * (1.0_r8 + wnd_frc_rat) * (1.0_r8 + wnd_frc_rat) / grav + !if (wnd_frc_slt > wnd_frc_thr_slt) then! if want to use fluid threshold for dust emission, uncomment this one, -dmleung 2 Dec 2021 + if (wnd_frc_slt > wnd_frc_thr_slt_it) then! if want to use impact threshold for dust emission, uncomment this one, -dmleung 2 Dec 2021 + + !################### for Zender et al. (2003) scheme -dmleung ########################### + !################ uncomment the below block if want to use Z03 scheme ################### + !wnd_frc_rat = wnd_frc_thr_slt / wnd_frc_slt + !flx_mss_hrz_slt_ttl = cst_slt * forc_rho(c) * (wnd_frc_slt**3.0_r8) * & + ! (1.0_r8 - wnd_frc_rat) * (1.0_r8 + wnd_frc_rat) * (1.0_r8 + wnd_frc_rat) / grav ! the following loop originates from subr. dst_mbl ! purpose: apply land sfc and veg limitations and global tuning factor ! slevis: multiply flx_mss_hrz_slt_ttl by liqfrac to incude the effect ! of frozen soil - flx_mss_hrz_slt_ttl = flx_mss_hrz_slt_ttl * lnd_frc_mbl(p) * mbl_bsn_fct(c) * & - flx_mss_fdg_fct * liqfrac + !flx_mss_hrz_slt_ttl = flx_mss_hrz_slt_ttl * lnd_frc_mbl(p) * mbl_bsn_fct(c) * & + ! flx_mss_fdg_fct * liqfrac + + ! dmleung moved to this block + !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 + !######################################################################################## + + !################### for Kok et al. (2014) scheme -dmleung ############################## + !################ uncomment the below block if want to use K14 scheme ################### + + ! if want to use fluid threshold for dust emission, uncomment this one, -dmleung 27 Nov 2021 + !flx_mss_vrt_dst_ttl(p) = dst_emiss_coeff(p) * mss_frc_cly_vld(c) * forc_rho(c) * ((wnd_frc_slt**2.0_r8 - wnd_frc_thr_slt**2.0_r8) / wnd_frc_thr_slt_std) * (wnd_frc_slt / wnd_frc_thr_slt)**frag_expt ! change forc_rho(g) to forc_rho(c) to avoid passing Nan values to the coupler -Longlei ! if want to use fluid threshold for dust emission, uncomment this one, -dml 27 Nov 2021 + + ! if want to use impact threshold for dust emission, uncomment this one, -dmleung 2 Dec 2021 + flx_mss_vrt_dst_ttl(p) = dst_emiss_coeff(p) * mss_frc_cly_vld(c) * forc_rho(c) * ((wnd_frc_slt**2.0_r8 - wnd_frc_thr_slt_it**2.0_r8) / wnd_frc_thr_slt_std) * (wnd_frc_slt / wnd_frc_thr_slt_it)**frag_expt ! if want to use impact threshold for dust emission, uncomment this one, -dml 2 Dec 2021 + + ! account for bare soil fraction, frozen soil fraction, and apply global tuning parameter (Kok et al. 2014) + flx_mss_vrt_dst_ttl(p) = flx_mss_vrt_dst_ttl(p) * lnd_frc_mbl(p) * C_tune * liqfrac + !######################################################################################## end if ! the following comes from subr. flx_mss_vrt_dst_ttl_MaB95_get ! purpose: diagnose total vertical mass flux of dust from vertically ! integrated streamwise mass flux - 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 + !dst_slt_flx_rat_ttl = 100.0_r8 * exp( log(10.0_r8) * (13.4_r8 * mss_frc_cly_vld(c) - 6.0_r8) ) ! dmleung commented and moved to the previous block + !flx_mss_vrt_dst_ttl(p) = flx_mss_hrz_slt_ttl * dst_slt_flx_rat_ttl + + !############## added by dmleung 2 Dec 2021 ############################################# + ! subsection for intermittency factor calculation + ! need to use with impact threshold and cannot be used with fluid threshold + ! Danny M. Leung, 24 Jun 2019, readded into CLM5 by dmleung 2 Dec 2021 + ! 2 Dec 2021 note: assume no buoyancy contribution to the wind fluctuation (u_sd_slt), so no obul(p) is needed. It is shown to be important for the wind fluctuations contribute little to the intermittency factor. We might add this back in the future revisions. + + ! mean lowpass-filtered wind speed at 0.1 m saltation height (assuming aerodynamic roughness length = 1e-4 m globally for ease; also assuming neutral condition) + u_mean_slt(p) = (wnd_frc_slt/k) * log(0.1_r8 / 1e-4_r8) + + ! sd of lowpass-filtered wind speed + !if (obul(p)==0) then + ! zetaobu = 0 + !else + !zetaobu = zii(p) / obul(p) ! For now zii is a constant of 1000 m in CLM -dml, 24 Aug 2021 + ! zetaobu = 1000_r8 / obul(p) ! For now zii is a constant of 1000 m in CLM -dml, 24 Aug 2021 + !end if + !stblty(p) = zetaobu ! zetaobu get outputted as the Obukhov stability parameter + stblty(p) = 0 ! -dmleung 2 Dec 2021: use 0 for now, assuming no buoyancy contribution. Might uncomment the above lines in future revisions. + if ((12_r8 - 0.5_r8 * stblty(p)) .GE. 0.001_r8) then + u_sd_slt(p) = wnd_frc_slt * (12_r8 - 0.5_r8 * stblty(p))**0.333_r8 + else + u_sd_slt(p) = 0.001_r8 ! should have used 0 theoretically; used 0.001 here to avoid undefined values + end if + + ! threshold velocities + ! Here wnd_frc_thr_slt is the fluid threshold; wnd_frc_thr_dry(p) is the dry fluid threshold; B_it*wnd_frc_thr_dry(p) is the impact threshold, -dml, 1 Mar 2021 + ! fluid threshold wind at 0.1 m saltation height + u_fld_thr(p) = (wnd_frc_thr_slt/k) * log(0.1_r8 / 1e-4_r8) + ! impact threshold wind at 0.1 m saltation height + u_impct_thr(p) = (wnd_frc_thr_slt_it/k) * log(0.1_r8 / 1e-4_r8) ! to avoid model error + + ! threshold crossing rate + thr_crs_rate(p) = (exp((u_fld_thr(p)**2_r8 - u_impct_thr(p)**2_r8 - 2_r8 * u_mean_slt(p) * (u_fld_thr(p) - u_impct_thr(p))) / (2_r8 * u_sd_slt(p)**2_r8)) + 1_r8)**(-1_r8) + + ! probability that lowpass-filtered wind speed does not exceed u_ft + prb_crs_fld_thr(p) = 0.5_r8 * (1_r8 + erf((u_fld_thr(p) - u_mean_slt(p)) / (1.414_r8 * u_sd_slt(p)))) + ! probability that lowpass-filtered wind speed does not exceed u_it + prb_crs_impct_thr(p) = 0.5_r8 * (1_r8 + erf((u_impct_thr(p) - u_mean_slt(p)) / (1.414_r8 * u_sd_slt(p)))) + + ! intermittency factor (from 0 to 1) + intrmtncy_fct(p) = 1_r8 - prb_crs_fld_thr(p) + thr_crs_rate(p) * (prb_crs_fld_thr(p) - prb_crs_impct_thr(p)) + + ! multiply dust emission flux by intermittency factor + if (intrmtncy_fct(p) /= intrmtncy_fct(p)) then ! if intrmtncy_fct(p) is not NaN then multiply by intermittency factor; this statement is needed because dust emission flx_mss_vrt_dst_ttl(p) has to be non NaN (at least zero) to be outputted, dmleung 9 Jun 2021 + flx_mss_vrt_dst_ttl(p) = flx_mss_vrt_dst_ttl(p) ! -dmleung + else + flx_mss_vrt_dst_ttl(p) = flx_mss_vrt_dst_ttl(p) * intrmtncy_fct(p) ! multiply dust flux by intermittency -dmleung + end if + + !############### end my subsection here -dmleung ######################################## end if ! lnd_frc_mbl > 0.0 From 2806c5b11ea6562c9c36513441d5449a5185c3db Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Mon, 18 Apr 2022 16:19:27 -0600 Subject: [PATCH 07/13] Avoid negative downwelling LW radiation from CAM --- src/main/lnd2atmMod.F90 | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/main/lnd2atmMod.F90 b/src/main/lnd2atmMod.F90 index 2a4a00d5cd..46a2ca354d 100644 --- a/src/main/lnd2atmMod.F90 +++ b/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 From 9b26a7bb300784fc8acec01e4061ac3e6e497061 Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Mon, 18 Apr 2022 16:25:30 -0600 Subject: [PATCH 08/13] Avoid negative downwelling LW radiation from CAM --- src/cpl/lnd_import_export.F90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/cpl/lnd_import_export.F90 b/src/cpl/lnd_import_export.F90 index 4923b721af..3e9e9bb0a6 100644 --- a/src/cpl/lnd_import_export.F90 +++ b/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 From 51a538247a3e1ed393ec0702e89fe95129a46a0a Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Mon, 11 Jul 2022 16:30:38 -0600 Subject: [PATCH 09/13] Per Erik Kluzek's request, dmleung commented out all lines relevant to the alternative clay fraction dataset from SoilGrids. CLM will use its own clay dataset. (dmleung edited 11 Jul 2022) --- .../namelist_definition_clm4_5.xml | 9 ++++--- src/biogeophys/SoilStateInitTimeConstMod.F90 | 27 ++++++++++--------- src/biogeophys/SoilStateType.F90 | 6 ++--- src/main/clm_varctl.F90 | 2 +- src/main/controlMod.F90 | 15 ++++++----- 5 files changed, 31 insertions(+), 28 deletions(-) diff --git a/bld/namelist_files/namelist_definition_clm4_5.xml b/bld/namelist_files/namelist_definition_clm4_5.xml index e109f1fd12..8655e17752 100644 --- a/bld/namelist_files/namelist_definition_clm4_5.xml +++ b/bld/namelist_files/namelist_definition_clm4_5.xml @@ -2574,10 +2574,11 @@ use case.) - -Full pathname of time-invariant clay fraction dataset for calculating dust emission. - + + + + diff --git a/src/biogeophys/SoilStateInitTimeConstMod.F90 b/src/biogeophys/SoilStateInitTimeConstMod.F90 index d9aae77d91..275407b068 100644 --- a/src/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/biogeophys/SoilStateInitTimeConstMod.F90 @@ -101,7 +101,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) 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, & - clay_frc, & ! -dmleung added clay_frc to CESM2/CLM5 14 Dec 2021 + !clay_frc, & ! -dmleung added clay_frc to CESM2/CLM5 14 Dec 2021; commented out per Erik Kluzek's request 11 Jul 2022 rough_fct, lulc_frc ! -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 @@ -153,7 +153,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 :: clayfrc2d(:) ! read in - clay percentage (called fraction here) from SoilGrids, dmleung added 14 Dec 2021 + !real(r8) ,pointer :: clayfrc2d(:) ! read in - clay percentage (called fraction here) from SoilGrids, dmleung added 14 Dec 2021; commented out per Erik Kluzek's request; 11 Jul 2022 real(r8) ,pointer :: roughfct2d(:) ! read in - time-invariant roughness factor, dmleung added 17 Dec 2021 real(r8) ,pointer :: rockfrc2d(:) ! read in - rock area fraction from GLCNMO, dmleung added 17 Dec 2021 real(r8) ,pointer :: vegefrc2d(:) ! read in - shrub/crop area fraction from GLCNMO, dmleung added 17 Dec 2021 @@ -306,19 +306,20 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) !############ subsection for input data for new dust emission scheme ############## ! dmleung added to CESM2/CLM5 14 Dec 2021 - allocate(clayfrc2d(begg:endg)) ! 2 Mar 2021 + ! commented out per Erik Kluzek's request; not incorpating into CLM. 11 Jul 2022 + !allocate(clayfrc2d(begg:endg)) ! 2 Mar 2021 ! read SoilGrids clay fraction, 2 Mar 2021 - write(iulog,*) 'Attempting to read alternative (e.g., SoilGrids) clay fraction data, by dmleung .....' - call getfil (clay_frc, locfn, 0) - call ncd_pio_openfile (ncid, locfn, 0) - call ncd_io(ncid=ncid, varname='f_clay', flag='read', data=clayfrc2d, dim1name=grlnd, readvar=readvar) + !write(iulog,*) 'Attempting to read alternative (e.g., SoilGrids) clay fraction data, by dmleung .....' + !call getfil (clay_frc, locfn, 0) + !call ncd_pio_openfile (ncid, locfn, 0) + !call ncd_io(ncid=ncid, varname='f_clay', flag='read', data=clayfrc2d, dim1name=grlnd, readvar=readvar) - do p = begp,endp - g = patch%gridcell(p) - soilstate_inst%clayfrc_patch(p) = clayfrc2d(g)/100_r8 - end do + !do p = begp,endp + ! g = patch%gridcell(p) + ! soilstate_inst%clayfrc_patch(p) = clayfrc2d(g)/100_r8 + !end do - call ncd_pio_closefile(ncid) + !call ncd_pio_closefile(ncid) ! dmleung added to CESM2/CLM5 17 Dec 2021 allocate(roughfct2d(begg:endg)) ! dmleung, 16 Jul 2020 @@ -672,7 +673,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) g = col%gridcell(c) 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.01_r8*(0.17_r8*clayfrc2d(g) + 0.0014_r8*clayfrc2d(g)*clayfrc2d(g)) !Fecan et al. (1999) -jfk, dmleung coded 14 Dec 2021 for SoilGrids clay fraction + !soilstate_inst%gwc_thr_col(c) = 0.01_r8*(0.17_r8*clayfrc2d(g) + 0.0014_r8*clayfrc2d(g)*clayfrc2d(g)) !Fecan et al. (1999) -jfk, dmleung coded 14 Dec 2021 for SoilGrids clay fraction; commented out per Erik Kluzek's request 11 Jul 2022 !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 diff --git a/src/biogeophys/SoilStateType.F90 b/src/biogeophys/SoilStateType.F90 index ea0755f3bf..9b5d25eaea 100644 --- a/src/biogeophys/SoilStateType.F90 +++ b/src/biogeophys/SoilStateType.F90 @@ -24,14 +24,14 @@ 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 :: clayfrc_patch (:) ! alternative clay fraction, 14 dec 2021 + !real(r8), pointer :: clayfrc_patch (:) ! alternative clay fraction, 14 dec 2021; eliminated per Erik Kluzek's request, 11 Jul 2022 real(r8), pointer :: roughfct_patch (:) ! roughness factor, 17 dec 2021 real(r8), pointer :: rockfrc_patch (:) ! GLCNMO rock/bare area fraction, 17 dec 2021 real(r8), pointer :: vegefrc_patch (:) ! GLCNMO shrub/crop area fraction, 17 dec 2021 @@ -139,7 +139,7 @@ subroutine InitAllocate(this, bounds) allocate(this%bd_col (begc:endc,nlevgrnd)) ; this%bd_col (:,:) = nan !################ dmleung added 14 Dec 2021 ######################## - allocate(this%clayfrc_patch (begp:endp)) ; this%clayfrc_patch (:) = nan + !allocate(this%clayfrc_patch (begp:endp)) ; this%clayfrc_patch (:) = nan !eliminated per Erik Kluzek's request, dmelung 11 Jul 2022 allocate(this%roughfct_patch (begp:endp)) ; this%roughfct_patch (:) = nan allocate(this%rockfrc_patch (begp:endp)) ; this%rockfrc_patch (:) = nan allocate(this%vegefrc_patch (begp:endp)) ; this%vegefrc_patch (:) = nan diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index 73e774b8fc..0cbba97f32 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -92,7 +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 :: clay_frc = ' ' ! alternative clay fraction, dmleung, added 14 Dec 2021 + !character(len=fname_len), public :: clay_frc = ' ' ! alternative clay fraction, dmleung, added 14 Dec 2021; eliminate the use of alternative clay fraction per Erik Kluzek's request 11 Jul 2022 character(len=fname_len), public :: rough_fct = ' ' ! roughness factor, dmleung, added 17 Dec 2021 character(len=fname_len), public :: lulc_frc = ' ' ! LULC area fraction, dmleung, added 17 Dec 2021 diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index d516f2d0a1..84d4f45f5d 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -264,7 +264,7 @@ subroutine control_init( ) ! dust emission, dmleung 14 Dec 2021 namelist /clm_inparm/ & - clay_frc, & !-dmleung 14 Dec 2021 + !clay_frc, & !-dmleung 14 Dec 2021; eliminated per Erik Kluzek's request, dmleung 11 Jul 2022 rough_fct, lulc_frc !-dmleung 17 Dec 2021 ! ---------------------------------------------------------------------- @@ -611,7 +611,7 @@ subroutine control_spmd() 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 (clay_frc, len(clay_frc), MPI_CHARACTER, 0, mpicom, ier) ! added by dmleung, 14 Dec 2021 + !call mpi_bcast (clay_frc, len(clay_frc), MPI_CHARACTER, 0, mpicom, ier) ! added by dmleung, 14 Dec 2021; eliminated per Erk Kluzek's request, dmleung 11 Jul 2022 call mpi_bcast (rough_fct, len(rough_fct), MPI_CHARACTER, 0, mpicom, ier)! added by dmleung, 17 Dec 2021 call mpi_bcast (lulc_frc, len(lulc_frc), MPI_CHARACTER, 0, mpicom, ier) ! added by dmleung, 17 Dec 2021 @@ -833,11 +833,12 @@ subroutine control_print () write(iulog,*) ' land frac data = ',trim(fatmlndfrc) end if !##### for dmleung's input data for new dust emission module ##### - if (clay_frc == ' ') then ! -dml, 2 Mar 2021, added 14 Dec 2021 - write(iulog,*) ' clay_frc surface dataset not set' - else - write(iulog,*) ' surface data = ',trim(clay_frc) - end if + ! eliminate the use of alternative clay fraction per Erik Kluzek's request, dmleung 11 Jul 2022 + !if (clay_frc == ' ') then ! -dml, 2 Mar 2021, added 14 Dec 2021 + ! write(iulog,*) ' clay_frc surface dataset not set' + !else + ! write(iulog,*) ' surface data = ',trim(clay_frc) + !end if if (rough_fct == ' ') then ! -dmleung, 17 Dec 2020 write(iulog,*) ' rough_fct surface dataset not set' else From 13c8e51f8e95ff47f596580ad1b861a542f0b823 Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Mon, 11 Jul 2022 17:10:47 -0600 Subject: [PATCH 10/13] dmleung removed all lines relevant to the alternative clay dataset from SoilGrids (edit 11 Jul 2022). --- .../namelist_definition_clm4_5.xml | 8 -------- src/biogeophys/SoilStateInitTimeConstMod.F90 | 18 ------------------ src/biogeophys/SoilStateType.F90 | 2 -- src/main/clm_varctl.F90 | 1 - src/main/controlMod.F90 | 10 +--------- 5 files changed, 1 insertion(+), 38 deletions(-) diff --git a/bld/namelist_files/namelist_definition_clm4_5.xml b/bld/namelist_files/namelist_definition_clm4_5.xml index 8655e17752..710e70bad2 100644 --- a/bld/namelist_files/namelist_definition_clm4_5.xml +++ b/bld/namelist_files/namelist_definition_clm4_5.xml @@ -2572,14 +2572,6 @@ use case.) - - - - - - - shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) - use clm_varpar , only : dst_src_nbr, ndst, sz_nbr + use clm_varpar , only : dst_src_nbr, ndst, sz_nbr, & + natpft_size ! -dmleung added 24 Jul 2022 use clm_varcon , only : grav, spval use landunit_varcon , only : istcrop, istsoil use clm_varctl , only : iulog @@ -30,6 +31,9 @@ module DUSTMod use LandunitType , only : lun use ColumnType , only : col use PatchType , only : patch + use clm_instur , only : wt_lunit, wt_nat_patch ! dmleung added 24 Jul 2022 + use landunit_varcon , only : istsoil, istcrop ! dmleung added 24 Jul 2022 (refering to main/landunit_varcon.F90, for wt_lunit, istsoil=1 is nat veg, istcrop=2 is crop) + use pftconMod , only : noveg ! dmleung added 24 Jul 2022 ! ! !PUBLIC TYPES implicit none @@ -375,6 +379,9 @@ subroutine DustEmission (bounds, & real(r8) :: wnd_frc_thr_slt ! [m/s] used for wet fluid threshold friction velocity, dmleung 9 Jun 2021 !########### added by dmleung 20 Dec 2021 for drag partition effect ################# real(r8) :: K_length ! [dimless] normalized mean interobstacle distance, or called gap length (Okin, 2008) + !########### added by dmleung 22 Jul 2022 for LUH2 land cover #################### + real(r8) :: bare_frc ! LUH2 bare soil land cover fraction + real(r8) :: veg_frc ! LUH2 natural vegetation + crop land cover fraction ! ! constants ! @@ -629,7 +636,14 @@ subroutine DustEmission (bounds, & K_length = 2_r8 * (1_r8/lai(p) - 1_r8) ! Here LAI has to be non-zero to avoid blowup ssr(p) = (K_length+f_0*c_e)/(K_length+c_e) - frc_thr_rgh_fct = (rockfrc(p)*(roughfct(p))**3_r8 + (vegefrc(p)+sparfrc(p))*(ssr(p))**3_r8 )**(0.3333_r8) ! land cover weighted mean using static GLCNMo bare land fraction LC0, dmleung 20 Dec 2021 + !frc_thr_rgh_fct = (rockfrc(p)*(roughfct(p))**3_r8 + (vegefrc(p)+sparfrc(p))*(ssr(p))**3_r8 )**(0.3333_r8) ! land cover weighted mean using static GLCNMo bare land fraction LC0, dmleung 20 Dec 2021i; dmleung commented 24 Jul 2022 + + ! dmleung added calculation of LUH2 bare vs veg fraction within a grid 24 Jul 2022 + bare_frc = wt_lunit(g,istsoil) * wt_nat_patch(g,noveg) + veg_frc = wt_lunit(g,istsoil) * sum(wt_nat_patch(g,(noveg+1):natpft_size)) + wt_lunit(g,istcrop) + + frc_thr_rgh_fct = (bare_frc*(roughfct(p))**3_r8 + veg_frc*(ssr(p))**3_r8 )**(0.3333_r8) ! land cover weighted mean using LUH2 land cover, dmleung 24 Jul 2022 + wnd_frc_slt = fv(p) * frc_thr_rgh_fct ! wnd_frc_slt will be used in the dust emission equation -dmleung From 2ed57e598b2469f409fdb164654c6cf13707013b Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Mon, 1 Aug 2022 15:46:21 -0600 Subject: [PATCH 12/13] 1 Aug 2022 dmleung talked to ekluzek and changed the index of wt_lunit in DUSTMod.F90. It is correct to use natveg_lb:natveg_ub instead of natveg_size to denote the indices of wt_lunit. In calculating the land cover fraction of bare soil, noveg is is used instead of natveg_lb. dmleung also added a history field for diagnostics of the standardized fluid threshold friction velocity (m/s) for the dust emission scheme in DUSTMod.F90. --- src/biogeochem/DUSTMod.F90 | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/src/biogeochem/DUSTMod.F90 b/src/biogeochem/DUSTMod.F90 index 1dee188896..5cea60dc51 100644 --- a/src/biogeochem/DUSTMod.F90 +++ b/src/biogeochem/DUSTMod.F90 @@ -16,7 +16,7 @@ module DUSTMod use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) use clm_varpar , only : dst_src_nbr, ndst, sz_nbr, & - natpft_size ! -dmleung added 24 Jul 2022 + natpft_lb, natpft_ub, natpft_size ! -dmleung added 24 Jul 2022 use clm_varcon , only : grav, spval use landunit_varcon , only : istcrop, istsoil use clm_varctl , only : iulog @@ -88,6 +88,8 @@ module DUSTMod real(r8), pointer, private :: ssr_patch (:) ! [dimless] integrated shear stress ratiio, defined by Okin (2008) and then integrated by Caroline Pierre et al. (2014) real(r8), pointer, private :: lai_patch (:) ! [m2 leaf /m2 land] LAI+SAI for calculating Okin's drag partition, averaged to landunit level real(r8), pointer, private :: frc_thr_rghn_fct_patch (:) ! [dimless] hybrid drag partition (or called roughness) factor + !########### added by dmleung 28 Jul 2022 ######################################################################## + real(r8), pointer, private :: wnd_frc_thr_std_patch (:) ! standardized fluid threshold friction velocity (m/s) contains procedure , public :: Init @@ -164,6 +166,8 @@ subroutine InitAllocate(this, bounds) allocate(this%ssr_patch (begp:endp)) ; this%ssr_patch (:) = nan allocate(this%lai_patch (begp:endp)) ; this%lai_patch (:) = nan allocate(this%frc_thr_rghn_fct_patch (begp:endp)) ; this%frc_thr_rghn_fct_patch (:) = nan + !#### added by dmleung 28 Jul 2022 ###################################### + allocate(this%wnd_frc_thr_std_patch (begp:endp)) ; this%wnd_frc_thr_std_patch (:) = nan end subroutine InitAllocate !------------------------------------------------------------------------ @@ -292,6 +296,11 @@ subroutine InitHistory(this, bounds) call hist_addfld1d (fname='FRC_THR_RGHN_FCT', units='dimensionless', & avgflag='A', long_name='hybrid drag partition (or roughness) factor', & ptr_patch=this%frc_thr_rghn_fct_patch, set_lake=0._r8, set_urb=0._r8) + !#####added by dmleung 28 Jul 2022 ######################################## + this%wnd_frc_thr_std_patch(begp:endp) = spval + call hist_addfld1d (fname='WND_FRC_FT_STD', units='m/s', & + avgflag='A', long_name='standardized fluid threshold friction velocity', & + ptr_patch=this%wnd_frc_thr_std_patch, set_lake=0._r8, set_urb=0._r8) !########################################################################## end subroutine InitHistory @@ -455,7 +464,9 @@ subroutine DustEmission (bounds, & ! added by dmleung 20 Dec 2021 ssr => dust_inst%ssr_patch , & lai => dust_inst%lai_patch , & - frc_thr_rghn_fct => dust_inst%frc_thr_rghn_fct_patch & + frc_thr_rghn_fct => dust_inst%frc_thr_rghn_fct_patch , & + ! added by dmleung 28 Jul 2022 + wnd_frc_thr_std => dust_inst%wnd_frc_thr_std_patch & ) ttlai(bounds%begp : bounds%endp) = 0._r8 @@ -557,6 +568,8 @@ subroutine DustEmission (bounds, & ssr(p) = 0.0_r8 lai(p) = 0.0_r8 frc_thr_rghn_fct(p) = 0.0_r8 + ! dmleung added 28 Jul 2022 + wnd_frc_thr_std(p) = 0.0_r8 end do do n = 1, ndst do fp = 1,num_nolakep @@ -609,6 +622,7 @@ subroutine DustEmission (bounds, & ! use emission threshold to calculate standardized threshold and dust emission coefficient dmleung 27 Nov 2021 wnd_frc_thr_slt_std = wnd_frc_thr_slt * sqrt(forc_rho(c) / forc_rho_std) ! standardized soil threshold friction speed -jfk (defined using fluid threshold + wnd_frc_thr_std(p) = wnd_frc_thr_slt_std ! output standardized fluid threshold -dmleung added 28 Jul 2022 dst_emiss_coeff(p) = Cd0 * exp(-Ce * (wnd_frc_thr_slt_std - wnd_frc_thr_slt_std_min) / wnd_frc_thr_slt_std_min) ! save dust emission coefficient here for all grids, -dml, 1 Mar 2021 ! framentation exponent dmleung 27 Nov 2021; moved to this block 23 Dec 2021 @@ -640,7 +654,7 @@ subroutine DustEmission (bounds, & ! dmleung added calculation of LUH2 bare vs veg fraction within a grid 24 Jul 2022 bare_frc = wt_lunit(g,istsoil) * wt_nat_patch(g,noveg) - veg_frc = wt_lunit(g,istsoil) * sum(wt_nat_patch(g,(noveg+1):natpft_size)) + wt_lunit(g,istcrop) + veg_frc = wt_lunit(g,istsoil) * sum(wt_nat_patch(g,(noveg+1):natpft_ub)) + wt_lunit(g,istcrop) frc_thr_rgh_fct = (bare_frc*(roughfct(p))**3_r8 + veg_frc*(ssr(p))**3_r8 )**(0.3333_r8) ! land cover weighted mean using LUH2 land cover, dmleung 24 Jul 2022 From b998c1e8cb8c52622eea8193dbe1166437721ca0 Mon Sep 17 00:00:00 2001 From: "Danny M. Leung" Date: Tue, 2 Aug 2022 13:20:54 -0600 Subject: [PATCH 13/13] 2 Aug 2022 dmleung updated the parameters within the dry fluid threshold scheme of Shao and Lu (2000) to be consistent with dmleung's dust emission scheme paper. dmleung also removed most of the relevant paths of the existing land cover dataset to stick with the LUH2 land cover within CLM. --- .../namelist_definition_clm4_5.xml | 5 ---- src/biogeochem/DUSTMod.F90 | 5 +--- src/biogeophys/SoilStateInitTimeConstMod.F90 | 26 +++---------------- src/biogeophys/SoilStateType.F90 | 6 ----- src/main/clm_varctl.F90 | 1 - src/main/controlMod.F90 | 8 +----- 6 files changed, 5 insertions(+), 46 deletions(-) diff --git a/bld/namelist_files/namelist_definition_clm4_5.xml b/bld/namelist_files/namelist_definition_clm4_5.xml index 710e70bad2..bd2046ea5f 100644 --- a/bld/namelist_files/namelist_definition_clm4_5.xml +++ b/bld/namelist_files/namelist_definition_clm4_5.xml @@ -2579,11 +2579,6 @@ use case.) Full pathname of time-invariant roughness factor dataset for calculating drag partition of dust emission. - - -Full pathname of time-invariant GLCNMO LULC area fraction for calculating hybrid drag partition. - diff --git a/src/biogeochem/DUSTMod.F90 b/src/biogeochem/DUSTMod.F90 index 5cea60dc51..ec0d702582 100644 --- a/src/biogeochem/DUSTMod.F90 +++ b/src/biogeochem/DUSTMod.F90 @@ -457,9 +457,6 @@ subroutine DustEmission (bounds, & prb_crs_impct_thr => dust_inst%prb_crs_impct_thr_patch , & ! added by dmleung 17 Dec 2021 roughfct => soilstate_inst%roughfct_patch , & - rockfrc => soilstate_inst%rockfrc_patch , & - vegefrc => soilstate_inst%vegefrc_patch , & - sparfrc => soilstate_inst%sparfrc_patch , & ustar => dust_inst%ustar_patch , & ! Output friction velocity for SP mode ! added by dmleung 20 Dec 2021 ssr => dust_inst%ssr_patch , & @@ -609,7 +606,7 @@ subroutine DustEmission (bounds, & ! recoded to CLM5 27 Nov 2021 !####################################################################################################### - tmp2 = 1.0_r8*sqrt(0.0123_r8 * (dns_slt*grav*140.0e-6_r8 + 5.0e-4_r8/140.0e-6_r8)) ! calculate S&L (2000) scheme here for threshold + tmp2 = 1.0_r8*sqrt(0.0123_r8 * (dns_slt*grav*130.0e-6_r8 + 1.65e-4_r8/130.0e-6_r8)) ! calculate S&L (2000) scheme here for threshold; gamma = 1.65e-4 following S&L00, D_p = 127 um ~ 130 um following dmleung's dust paper. As this is a global constant, this line can be put outside the loop to save computational power. wnd_frc_thr_dry(p) = tmp2 / sqrt(forc_rho(c)) ! output dry fluid threshold wnd_frc_thr_slt = tmp2 / sqrt(forc_rho(c)) * frc_thr_wet_fct !* frc_thr_rgh_fct ! use as threshold in this module wnd_frc_thr_slt_it = B_it * tmp2 / sqrt(forc_rho(c)) ! define impact threshold -dml 9 Jun 2021, recoded to CLM5 27 Nov 2021 diff --git a/src/biogeophys/SoilStateInitTimeConstMod.F90 b/src/biogeophys/SoilStateInitTimeConstMod.F90 index 023e41719d..4d13a6d9a7 100644 --- a/src/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/biogeophys/SoilStateInitTimeConstMod.F90 @@ -101,7 +101,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) 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, & - rough_fct, lulc_frc ! -dmleung added to CESM2/CLM5 17 Dec 2021 + 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 @@ -153,9 +153,6 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) 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 - real(r8) ,pointer :: rockfrc2d(:) ! read in - rock area fraction from GLCNMO, dmleung added 17 Dec 2021 - real(r8) ,pointer :: vegefrc2d(:) ! read in - shrub/crop area fraction from GLCNMO, dmleung added 17 Dec 2021 - real(r8) ,pointer :: sparfrc2d(:) ! read in - sparse vegetation area fraction from GLCNMO, dmleung added 17 Dec 2021 character(len=256) :: locfn ! local filename integer :: ipedof integer :: begp, endp @@ -307,34 +304,17 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) ! 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 .....' + !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,*) 'dimension of roughness = ', shape(roughfct2d) ! -dml - write(iulog,*) 'initialize pft level roughness factor from roughfct2d(g) to roughfct(p)' + !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 - ! dmleung added to CESM2/CLM5 17 Dec 2021 - allocate(rockfrc2d(begg:endg));allocate(vegefrc2d(begg:endg));allocate(sparfrc2d(begg:endg)) ! 7 Jul 2021 - write(iulog,*) 'Attempting to read GLCNMO LULC area fraction data, by dmleung .....' - call getfil (lulc_frc, locfn, 0) - call ncd_pio_openfile (ncid, locfn, 0) - call ncd_io(ncid=ncid, varname='A_r', flag='read', data=rockfrc2d, dim1name=grlnd, readvar=readvar) - call ncd_io(ncid=ncid, varname='A_v', flag='read', data=vegefrc2d, dim1name=grlnd, readvar=readvar) - call ncd_io(ncid=ncid, varname='A_s', flag='read', data=sparfrc2d, dim1name=grlnd, readvar=readvar) - - do p = begp,endp - g = patch%gridcell(p) - soilstate_inst%rockfrc_patch(p) = rockfrc2d(g) - soilstate_inst%vegefrc_patch(p) = vegefrc2d(g) - soilstate_inst%sparfrc_patch(p) = sparfrc2d(g) - end do - call ncd_pio_closefile(ncid) !################################################################################## diff --git a/src/biogeophys/SoilStateType.F90 b/src/biogeophys/SoilStateType.F90 index d5fe966f98..3c7a4ba825 100644 --- a/src/biogeophys/SoilStateType.F90 +++ b/src/biogeophys/SoilStateType.F90 @@ -32,9 +32,6 @@ module SoilStateType 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 - real(r8), pointer :: rockfrc_patch (:) ! GLCNMO rock/bare area fraction, 17 dec 2021 - real(r8), pointer :: vegefrc_patch (:) ! GLCNMO shrub/crop area fraction, 17 dec 2021 - real(r8), pointer :: sparfrc_patch (:) ! GLCNMO sparse vegetation area fraction, 17 dec 2021 !########################################################################################### ! hydraulic properties @@ -139,9 +136,6 @@ subroutine InitAllocate(this, bounds) !################ dmleung added 14 Dec 2021 ######################## allocate(this%roughfct_patch (begp:endp)) ; this%roughfct_patch (:) = nan - allocate(this%rockfrc_patch (begp:endp)) ; this%rockfrc_patch (:) = nan - allocate(this%vegefrc_patch (begp:endp)) ; this%vegefrc_patch (:) = nan - allocate(this%sparfrc_patch (begp:endp)) ; this%sparfrc_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 diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index 6235f5cc52..ee7ff62bd5 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -93,7 +93,6 @@ module clm_varctl 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 - character(len=fname_len), public :: lulc_frc = ' ' ! LULC area fraction, dmleung, added 17 Dec 2021 !---------------------------------------------------------- ! Flag to read ndep rather than obtain it from coupler diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index 9fae63250d..2f8d57d26a 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -264,7 +264,7 @@ subroutine control_init( ) ! dust emission, dmleung 14 Dec 2021 namelist /clm_inparm/ & - rough_fct, lulc_frc !-dmleung 17 Dec 2021 + rough_fct !-dmleung 17 Dec 2021 ! ---------------------------------------------------------------------- ! Default values @@ -611,7 +611,6 @@ subroutine control_spmd() ! 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 - call mpi_bcast (lulc_frc, len(lulc_frc), MPI_CHARACTER, 0, mpicom, ier) ! added by dmleung, 17 Dec 2021 ! Irrigation call mpi_bcast(irrigate, 1, MPI_LOGICAL, 0, mpicom, ier) @@ -836,11 +835,6 @@ subroutine control_print () else write(iulog,*) ' surface data = ',trim(rough_fct) end if - if (lulc_frc == ' ') then ! -dmleung, 17 Dec 2021 - write(iulog,*) ' lulc_frc surface dataset not set' - else - write(iulog,*) ' surface data = ',trim(lulc_frc) - end if !################################################################# if (use_cn) then if (suplnitro /= suplnNon)then