#include <misc.h>
#include <preproc.h>
module DUSTMod 2,8
#if (defined DUST)
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: DUSTMod
!
! !DESCRIPTION:
! Routines in this module calculate Dust mobilization and dry deposition for dust.
! Simulates dust mobilization due to wind from the surface into the
! lowest atmospheric layer. On output flx_mss_vrt_dst(ndst) is the surface dust
! emission (kg/m**2/s) [ + = to atm].
! Calculates the turbulent component of dust dry deposition, (the turbulent deposition
! velocity through the lowest atmospheric layer). CAM will calculate the settling
! velocity through the whole atmospheric column. The two calculations will determine
! the dust dry deposition flux to the surface.
!
! !USES:
use shr_kind_mod, only: r8 => shr_kind_r8
use clmtype
use clm_varpar , only : dst_src_nbr, ndst, sz_nbr
use clm_varcon , only : grav, istsoil
use clm_varctl , only : iulog
use abortutils , only : endrun
use subgridAveMod, only: p2l_1d
use clm_varcon, only: spval
!
! !PUBLIC TYPES
implicit none
save
!
! !PUBLIC MEMBER FUNCTIONS:
!
public Dustini ! Initialize variables used in subroutine Dust
public DustEmission ! Dust mobilization
public DustDryDep ! Turbulent dry deposition for dust
!
! !REVISION HISTORY
! Created by Sam Levis, updated to clm2.1 by Mariana Vertenstein
! Source: C. Zender's dust model
!
!EOP
!
! Data private to this module
!
private
real(r8) ovr_src_snk_mss(dst_src_nbr,ndst)
real(r8) tmp1 !Factor in saltation computation (named as in Charlie's code)
real(r8) dmt_vwr(ndst) ![m] Mass-weighted mean diameter resolved
real(r8) stk_crc(ndst) ![frc] Correction to Stokes settling velocity
real(r8) dns_aer ![kg m-3] Aerosol density
!------------------------------------------------------------------------
contains
!------------------------------------------------------------------------
!BOP
!
! !IROUTINE: DustEmission
!
! !INTERFACE:
subroutine DustEmission (lbp, ubp, lbc,ubc,lbl,ubl,num_nolakep, filter_nolakep) 1,4
!
! !DESCRIPTION:
! Dust mobilization. This code simulates dust mobilization due to wind
! from the surface into the lowest atmospheric layer
! On output flx_mss_vrt_dst(ndst) is the surface dust emission
! (kg/m**2/s) [ + = to atm]
! Source: C. Zender's dust model
!
! !USES
use clm_atmlnd , only : clm_a2l
use shr_const_mod, only : SHR_CONST_RHOFW
!
! !ARGUMENTS:
implicit none
integer, intent(in) :: lbp, ubp,lbc,ubc,ubl,lbl ! pft bounds
integer, intent(in) :: num_nolakep ! number of column non-lake points in pft filter
integer, intent(in) :: filter_nolakep(num_nolakep) ! pft filter for non-lake points
!
! !LOCAL VARIABLES
!
! local pointers to implicit in arguments
!
integer , pointer :: pcolumn(:) ! pft's column index
integer , pointer :: plandunit(:) ! pft's landunit index
integer , pointer :: pgridcell(:) ! pft's gridcell index
integer , pointer :: ityplun(:) ! landunit type
real(r8), pointer :: tlai(:) ! one-sided leaf area index, no burying by snow
real(r8), pointer :: tsai(:) ! one-sided stem area index, no burying by snow
real(r8), pointer :: frac_sno(:) ! fraction of ground covered by snow (0 to 1)
real(r8), pointer :: gwc_thr(:) ! threshold gravimetric soil moisture based on clay content
real(r8), pointer :: forc_rho(:) ! density (kg/m**3)
real(r8), pointer :: fv(:) ! friction velocity (m/s) (for dust model)
real(r8), pointer :: u10(:) ! 10-m wind (m/s) (created for dust model)
real(r8), pointer :: mbl_bsn_fct(:) ! basin factor
real(r8), pointer :: mss_frc_cly_vld(:) ! [frc] Mass fraction clay limited to 0.20
real(r8), pointer :: h2osoi_vol(:,:) ! volumetric soil water (0<=h2osoi_vol<=watsat)
real(r8), pointer :: h2osoi_liq(:,:) ! liquid soil water (kg/m2)
real(r8), pointer :: h2osoi_ice(:,:) ! frozen soil water (kg/m2)
real(r8), pointer :: watsat(:,:) ! saturated volumetric soil water
! local pointers to implicit out arguments
!
real(r8), pointer :: flx_mss_vrt_dst(:,:) ! surface dust emission (kg/m**2/s)
real(r8), pointer :: flx_mss_vrt_dst_tot(:) ! total dust flux into atmosphere
! !REVISION HISTORY
! Created by Sam Levis
! Migrated to new data structures by Peter Thornton and Mariana Vertenstein
! !Created by Peter Thornton and Mariana Vertenstein
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
integer :: fp,p,c,l,g,m,n ! indices
real(r8) :: liqfrac ! fraction of total water that is liquid
real(r8) :: wnd_frc_rat ! [frc] Wind friction threshold over wind friction
real(r8) :: wnd_frc_slt_dlt ! [m s-1] Friction velocity increase from saltatn
real(r8) :: wnd_rfr_dlt ! [m s-1] Reference windspeed excess over threshld
real(r8) :: dst_slt_flx_rat_ttl
real(r8) :: flx_mss_hrz_slt_ttl
real(r8) :: flx_mss_vrt_dst_ttl(lbp:ubp)
real(r8) :: frc_thr_wet_fct
real(r8) :: frc_thr_rgh_fct
real(r8) :: wnd_frc_thr_slt
real(r8) :: wnd_rfr_thr_slt
real(r8) :: wnd_frc_slt
real(r8) :: lnd_frc_mbl(lbp:ubp)
real(r8) :: bd
real(r8) :: gwc_sfc
real(r8) :: ttlai(lbp:ubp)
real(r8) :: tlai_lu(lbl:ubl)
!
! 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), pointer :: wtlunit(:) ! weight of pft relative to landunit
real(r8) :: sumwt(lbl:ubl) ! sum of weights
logical :: found ! temporary for error check
integer :: index
!------------------------------------------------------------------------
! Assign local pointers to derived type scalar members (gridcell-level)
forc_rho => clm_a2l%forc_rho
! Assign local pointers to derived type scalar members (landunit-level)
ityplun => clm3%g%l%itype
! Assign local pointers to derived type scalar members (column-level)
frac_sno => clm3%g%l%c%cps%frac_sno
gwc_thr => clm3%g%l%c%cps%gwc_thr
mbl_bsn_fct => clm3%g%l%c%cps%mbl_bsn_fct
mss_frc_cly_vld => clm3%g%l%c%cps%mss_frc_cly_vld
h2osoi_vol => clm3%g%l%c%cws%h2osoi_vol
h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq
h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice
watsat => clm3%g%l%c%cps%watsat
! Assign local pointers to derived type scalar members (pft-level)
pgridcell => clm3%g%l%c%p%gridcell
plandunit => clm3%g%l%c%p%landunit
pcolumn => clm3%g%l%c%p%column
tlai => clm3%g%l%c%p%pps%tlai
tsai => clm3%g%l%c%p%pps%tsai
fv => clm3%g%l%c%p%pps%fv
u10 => clm3%g%l%c%p%pps%u10
flx_mss_vrt_dst => clm3%g%l%c%p%pdf%flx_mss_vrt_dst
flx_mss_vrt_dst_tot => clm3%g%l%c%p%pdf%flx_mss_vrt_dst_tot
!local pointers from subgridAveMod/p2l_1d
wtlunit => clm3%g%l%c%p%wtlunit
ttlai(:) = 0._r8
! make lai average at landunit level
do fp = 1,num_nolakep
p = filter_nolakep(fp)
ttlai(p) = tlai(p)+tsai(p)
enddo
tlai_lu(:) = spval
sumwt(:) = 0._r8
do p = lbp,ubp
if (ttlai(p) /= spval .and. wtlunit(p) /= 0._r8) then
c = pcolumn(p)
l = plandunit(p)
if (sumwt(l) == 0._r8) tlai_lu(l) = 0._r8
tlai_lu(l) = tlai_lu(l) + ttlai(p) * wtlunit(p)
sumwt(l) = sumwt(l) + wtlunit(p)
end if
end do
found = .false.
do l = lbl,ubl
if (sumwt(l) > 1.0_r8 + 1.e-6_r8) then
found = .true.
index = l
exit
else if (sumwt(l) /= 0._r8) then
tlai_lu(l) = tlai_lu(l)/sumwt(l)
end if
end do
if (found) then
write(iulog,*) 'p2l_1d error: sumwt is greater than 1.0 at l= ',index
call endrun()
end if
! Loop through pfts
! initialize variables which get passed to the atmosphere
flx_mss_vrt_dst(lbp:ubp,:)=0._r8
!dir$ concurrent
!cdir nodep
do fp = 1,num_nolakep
p = filter_nolakep(fp)
c = pcolumn(p)
l = plandunit(p)
! the following code from subr. lnd_frc_mbl_get was adapted for lsm use
! purpose: return fraction of each gridcell suitable for dust mobilization
! the "bare ground" fraction of the current sub-gridscale cell decreases
! linearly from 1 to 0 as VAI(=tlai+tsai) increases from 0 to vai_mbl_thr
! if ice sheet, wetland, or lake, no dust allowed
if (ityplun(l) == istsoil) then
if (tlai_lu(l) < vai_mbl_thr) then
lnd_frc_mbl(p) = 1.0_r8 - (tlai_lu(l))/vai_mbl_thr
else
lnd_frc_mbl(p) = 0.0_r8
endif
lnd_frc_mbl(p) = lnd_frc_mbl(p) * (1.0_r8 - frac_sno(c))
else
lnd_frc_mbl(p) = 0.0_r8
end if
end do
!dir$ concurrent
!cdir nodep
do fp = 1,num_nolakep
p = filter_nolakep(fp)
if (lnd_frc_mbl(p)>1.0_r8 .or. lnd_frc_mbl(p)<0.0_r8) then
write(iulog,*)'Error dstmbl: pft= ',p,' lnd_frc_mbl(p)= ',lnd_frc_mbl(p)
call endrun
end if
end do
! reset history output variables before next if-statement to avoid output = inf
!dir$ concurrent
!cdir nodep
do fp = 1,num_nolakep
p = filter_nolakep(fp)
flx_mss_vrt_dst_tot(p) = 0.0_r8
end do
do n = 1, ndst
!dir$ concurrent
!cdir nodep
do fp = 1,num_nolakep
p = filter_nolakep(fp)
flx_mss_vrt_dst(p,n) = 0.0_r8
end do
end do
!dir$ concurrent
!cdir nodep
do fp = 1,num_nolakep
p = filter_nolakep(fp)
c = pcolumn(p)
l = plandunit(p)
g = pgridcell(p)
! only perform the following calculations if lnd_frc_mbl is non-zero
if (lnd_frc_mbl(p) > 0.0_r8) then
! the following comes from subr. frc_thr_rgh_fct_get
! purpose: compute factor by which surface roughness increases threshold
! friction velocity (currently a constant)
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
! adjust threshold velocity for inhibition by moisture
! 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
! 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) ) )
! the following lines come from subr. dst_mbl
! purpose: adjust threshold friction velocity to acct for moisture and
! roughness. The ratio tmp1 / sqrt(forc_rho) comes from
! subr. wnd_frc_thr_slt_get which computes dry threshold
! friction velocity for saltation
wnd_frc_thr_slt = tmp1 / sqrt(forc_rho(g)) * 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)
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)
! 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
! 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(g) * (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
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
end if ! lnd_frc_mbl > 0.0
end do
! the following comes from subr. flx_mss_vrt_dst_prt in C. Zender's code
! purpose: partition total vertical mass flux of dust into transport bins
do n = 1, ndst
do m = 1, dst_src_nbr
!dir$ concurrent
!cdir nodep
do fp = 1,num_nolakep
p = filter_nolakep(fp)
if (lnd_frc_mbl(p) > 0.0_r8) then
flx_mss_vrt_dst(p,n) = flx_mss_vrt_dst(p,n) + ovr_src_snk_mss(m,n) * flx_mss_vrt_dst_ttl(p)
end if
end do
end do
end do
do n = 1, ndst
!dir$ concurrent
!cdir nodep
do fp = 1,num_nolakep
p = filter_nolakep(fp)
if (lnd_frc_mbl(p) > 0.0_r8) then
flx_mss_vrt_dst_tot(p) = flx_mss_vrt_dst_tot(p) + flx_mss_vrt_dst(p,n)
end if
end do
end do
end subroutine DustEmission
!------------------------------------------------------------------------
!BOP
!
! !IROUTINE: subroutine DustDryDep
!
! !INTERFACE:
!
subroutine DustDryDep (lbp, ubp) 2,3
!
! !DESCRIPTION:
!
! Determine Turbulent dry deposition for dust. Calculate the turbulent
! component of dust dry deposition, (the turbulent deposition velocity
! through the lowest atmospheric layer. CAM will calculate the settling
! velocity through the whole atmospheric column. The two calculations
! will determine the dust dry deposition flux to the surface.
! Note: Same process should occur over oceans. For the coupled CCSM,
! we may find it more efficient to let CAM calculate the turbulent dep
! velocity over all surfaces. This would require passing the
! aerodynamic resistance, ram(1), and the friction velocity, fv, from
! the land to the atmosphere component. In that case, dustini need not
! calculate particle diamter (dmt_vwr) and particle density (dns_aer).
! Source: C. Zender's dry deposition code
!
! !USES
!
use shr_const_mod, only : SHR_CONST_PI, SHR_CONST_RDAIR, SHR_CONST_BOLTZ
use clm_atmlnd , only : clm_a2l
!
! !ARGUMENTS:
!
implicit none
integer, intent(in) :: lbp, ubp ! pft bounds
!
! !LOCAL VARIABLES
!
! local pointers to implicit in arguments
!
integer , pointer :: pgridcell(:) ! pft's gridcell index
real(r8), pointer :: pwtgcell(:) ! weight of pft relative to corresponding gridcell
real(r8), pointer :: forc_t(:) ! atm temperature (K)
real(r8), pointer :: forc_pbot(:) ! atm pressure (Pa)
real(r8), pointer :: forc_rho(:) ! atm density (kg/m**3)
real(r8), pointer :: fv(:) ! friction velocity (m/s)
real(r8), pointer :: ram1(:) ! aerodynamical resistance (s/m)
real(r8), pointer :: vlc_trb(:,:) ! Turbulent deposn velocity (m/s)
real(r8), pointer :: vlc_trb_1(:) ! Turbulent deposition velocity 1
real(r8), pointer :: vlc_trb_2(:) ! Turbulent deposition velocity 2
real(r8), pointer :: vlc_trb_3(:) ! Turbulent deposition velocity 3
real(r8), pointer :: vlc_trb_4(:) ! Turbulent deposition velocity 4
!
! !REVISION HISTORY
! Created by Sam Levis
!
!
! !LOCAL VARIABLES
!EOP
!
integer :: p,g,m,n ! indices
real(r8) :: vsc_dyn_atm(lbp:ubp) ! [kg m-1 s-1] Dynamic viscosity of air
real(r8) :: vsc_knm_atm(lbp:ubp) ! [m2 s-1] Kinematic viscosity of atmosphere
real(r8) :: shm_nbr_xpn ! [frc] Sfc-dep exponent for aerosol-diffusion dependence on Schmidt number
real(r8) :: shm_nbr ! [frc] Schmidt number
real(r8) :: stk_nbr ! [frc] Stokes number
real(r8) :: mfp_atm ! [m] Mean free path of air
real(r8) :: dff_aer ! [m2 s-1] Brownian diffusivity of particle
real(r8) :: rss_trb ! [s m-1] Resistance to turbulent deposition
real(r8) :: slp_crc(lbp:ubp,ndst) ! [frc] Slip correction factor
real(r8) :: vlc_grv(lbp:ubp,ndst) ! [m s-1] Settling velocity
real(r8) :: rss_lmn(lbp:ubp,ndst) ! [s m-1] Quasi-laminar layer resistance
real(r8) :: tmp ! temporary
! constants
real(r8),parameter::shm_nbr_xpn_lnd=-2._r8/3._r8 ![frc] shm_nbr_xpn over land
!------------------------------------------------------------------------
! Assign local pointers to derived type members (gridcell-level)
forc_pbot => clm_a2l%forc_pbot
forc_rho => clm_a2l%forc_rho
forc_t => clm_a2l%forc_t
! Assign local pointers to derived type members (pft-level)
pgridcell => clm3%g%l%c%p%gridcell
pwtgcell => clm3%g%l%c%p%wtgcell
fv => clm3%g%l%c%p%pps%fv
ram1 => clm3%g%l%c%p%pps%ram1
vlc_trb => clm3%g%l%c%p%pdf%vlc_trb
vlc_trb_1 => clm3%g%l%c%p%pdf%vlc_trb_1
vlc_trb_2 => clm3%g%l%c%p%pdf%vlc_trb_2
vlc_trb_3 => clm3%g%l%c%p%pdf%vlc_trb_3
vlc_trb_4 => clm3%g%l%c%p%pdf%vlc_trb_4
!dir$ concurrent
!cdir nodep
do p = lbp,ubp
if (pwtgcell(p)>0._r8) then
g = pgridcell(p)
! from subroutine dst_dps_dry (consider adding sanity checks from line 212)
! when code asks to use midlayer density, pressure, temperature,
! I use the data coming in from the atmosphere, ie forc_t, forc_pbot, forc_rho
! Quasi-laminar layer resistance: call rss_lmn_get
! Size-independent thermokinetic properties
vsc_dyn_atm(p) = 1.72e-5_r8 * ((forc_t(g)/273.0_r8)**1.5_r8) * 393.0_r8 / &
(forc_t(g)+120.0_r8) ![kg m-1 s-1] RoY94 p. 102
mfp_atm = 2.0_r8 * vsc_dyn_atm(p) / & ![m] SeP97 p. 455
(forc_pbot(g)*sqrt(8.0_r8/(SHR_CONST_PI*SHR_CONST_RDAIR*forc_t(g))))
vsc_knm_atm(p) = vsc_dyn_atm(p) / forc_rho(g) ![m2 s-1] Kinematic viscosity of air
do m = 1, ndst
slp_crc(p,m) = 1.0_r8 + 2.0_r8 * mfp_atm * &
(1.257_r8+0.4_r8*exp(-1.1_r8*dmt_vwr(m)/(2.0_r8*mfp_atm))) / &
dmt_vwr(m) ![frc] Slip correction factor SeP97 p. 464
vlc_grv(p,m) = (1.0_r8/18.0_r8) * dmt_vwr(m) * dmt_vwr(m) * dns_aer * &
grav * slp_crc(p,m) / vsc_dyn_atm(p) ![m s-1] Stokes' settling velocity SeP97 p. 466
vlc_grv(p,m) = vlc_grv(p,m) * stk_crc(m) ![m s-1] Correction to Stokes settling velocity
end do
end if
end do
do m = 1, ndst
!dir$ concurrent
!cdir nodep
do p = lbp,ubp
if (pwtgcell(p)>0._r8) then
g = pgridcell(p)
stk_nbr = vlc_grv(p,m) * fv(p) * fv(p) / (grav * vsc_knm_atm(p)) ![frc] SeP97 p.965
dff_aer = SHR_CONST_BOLTZ * forc_t(g) * slp_crc(p,m) / & ![m2 s-1]
(3.0_r8*SHR_CONST_PI * vsc_dyn_atm(p) * dmt_vwr(m)) !SeP97 p.474
shm_nbr = vsc_knm_atm(p) / dff_aer ![frc] SeP97 p.972
shm_nbr_xpn = shm_nbr_xpn_lnd ![frc]
! fxm: Turning this on dramatically reduces
! deposition velocity in low wind regimes
! Schmidt number exponent is -2/3 over solid surfaces and
! -1/2 over liquid surfaces SlS80 p. 1014
! if (oro(i)==0.0) shm_nbr_xpn=shm_nbr_xpn_ocn else shm_nbr_xpn=shm_nbr_xpn_lnd
! [frc] Surface-dependent exponent for aerosol-diffusion dependence on Schmidt #
tmp = shm_nbr**shm_nbr_xpn + 10.0_r8**(-3.0_r8/stk_nbr)
rss_lmn(p,m) = 1.0_r8 / (tmp * fv(p)) ![s m-1] SeP97 p.972,965
end if
end do
end do
! Lowest layer: Turbulent deposition (CAM will calc. gravitational dep)
do m = 1, ndst
!dir$ concurrent
!cdir nodep
do p = lbp,ubp
if (pwtgcell(p)>0._r8) then
rss_trb = ram1(p) + rss_lmn(p,m) + ram1(p) * rss_lmn(p,m) * vlc_grv(p,m) ![s m-1]
vlc_trb(p,m) = 1.0_r8 / rss_trb ![m s-1]
end if
end do
end do
!dir$ concurrent
!cdir nodep
do p = lbp,ubp
if (pwtgcell(p)>0._r8) then
vlc_trb_1(p) = vlc_trb(p,1)
vlc_trb_2(p) = vlc_trb(p,2)
vlc_trb_3(p) = vlc_trb(p,3)
vlc_trb_4(p) = vlc_trb(p,4)
end if
end do
end subroutine DustDryDep
!------------------------------------------------------------------------
!BOP
!
! !IROUTINE: subroutine Dustini
!
! !INTERFACE:
!
subroutine Dustini() 2,14
!
! !DESCRIPTION:
!
! Compute source efficiency factor from topography
! Initialize other variables used in subroutine Dust:
! ovr_src_snk_mss(m,n) and tmp1.
! Define particle diameter and density needed by atm model
! as well as by dry dep model
! Source: Paul Ginoux (for source efficiency factor)
! Modifications by C. Zender and later by S. Levis
! Rest of subroutine from C. Zender's dust model
!
! !USES
use shr_const_mod, only: SHR_CONST_PI, SHR_CONST_RDAIR
use decompMod, only : get_proc_bounds
!
! !ARGUMENTS:
implicit none
!
! !REVISION HISTORY
! Created by Samual Levis
!
! !LOCAL VARIABLES
!
! local pointers to implicit in arguments
!
real(r8), pointer :: mbl_bsn_fct(:) !basin factor
!
!
! !LOCAL VARIABLES
!EOP
!
integer :: fc,c,l,m,n ! indices
real(r8) :: ovr_src_snk_frc
real(r8) :: sqrt2lngsdi ! [frc] Factor in erf argument
real(r8) :: lndmaxjovrdmdni ! [frc] Factor in erf argument
real(r8) :: lndminjovrdmdni ! [frc] Factor in erf argument
real(r8) :: ryn_nbr_frc_thr_prx_opt ! [frc] Threshold friction Reynolds number approximation for optimal size
real(r8) :: ryn_nbr_frc_thr_opt_fnc ! [frc] Threshold friction Reynolds factor for saltation calculation
real(r8) :: icf_fct ! Interpartical cohesive forces factor for saltation calc
real(r8) :: dns_fct ! Density ratio factor for saltation calculation
real(r8) :: dmt_min(ndst) ! [m] Size grid minimum
real(r8) :: dmt_max(ndst) ! [m] Size grid maximum
real(r8) :: dmt_ctr(ndst) ! [m] Diameter at bin center
real(r8) :: dmt_dlt(ndst) ! [m] Width of size bin
real(r8) :: slp_crc(ndst) ! [frc] Slip correction factor
real(r8) :: vlm_rsl(ndst) ! [m3 m-3] Volume concentration resolved
real(r8) :: vlc_stk(ndst) ! [m s-1] Stokes settling velocity
real(r8) :: vlc_grv(ndst) ! [m s-1] Settling velocity
real(r8) :: ryn_nbr_grv(ndst) ! [frc] Reynolds number at terminal velocity
real(r8) :: cff_drg_grv(ndst) ! [frc] Drag coefficient at terminal velocity
real(r8) :: tmp ! temporary
real(r8) :: ln_gsd ! [frc] ln(gsd)
real(r8) :: gsd_anl ! [frc] Geometric standard deviation
real(r8) :: dmt_vma ! [m] Mass median diameter analytic She84 p.75 Tabl.1
real(r8) :: dmt_nma ! [m] Number median particle diameter
real(r8) :: lgn_dst ! Lognormal distribution at sz_ctr
real(r8) :: eps_max ! [frc] Relative accuracy for convergence
real(r8) :: eps_crr ! [frc] Current relative accuracy
real(r8) :: itr_idx ! [idx] Counting index
real(r8) :: dns_mdp ! [kg m-3] Midlayer density
real(r8) :: mfp_atm ! [m] Mean free path of air
real(r8) :: vsc_dyn_atm ! [kg m-1 s-1] Dynamic viscosity of air
real(r8) :: vsc_knm_atm ! [kg m-1 s-1] Kinematic viscosity of air
real(r8) :: vlc_grv_old ! [m s-1] Previous gravitational settling velocity
real(r8) :: series_ratio ! Factor for logarithmic grid
real(r8) :: lngsdsqrttwopi_rcp ! Factor in lognormal distribution
real(r8) :: sz_min(sz_nbr) ! [m] Size Bin minima
real(r8) :: sz_max(sz_nbr) ! [m] Size Bin maxima
real(r8) :: sz_ctr(sz_nbr) ! [m] Size Bin centers
real(r8) :: sz_dlt(sz_nbr) ! [m] Size Bin widths
! constants
real(r8) :: dmt_vma_src(dst_src_nbr) = & ! [m] Mass median diameter
(/ 0.832e-6_r8 , 4.82e-6_r8 , 19.38e-6_r8 /) ! BSM96 p. 73 Table 2
real(r8) :: gsd_anl_src(dst_src_nbr) = & ! [frc] Geometric std deviation
(/ 2.10_r8 , 1.90_r8 , 1.60_r8 /) ! BSM96 p. 73 Table 2
real(r8) :: mss_frc_src(dst_src_nbr) = & ! [frc] Mass fraction
(/ 0.036_r8, 0.957_r8, 0.007_r8 /) ! BSM96 p. 73 Table 2
real(r8) :: dmt_grd(5) = & ! [m] Particle diameter grid
(/ 0.1e-6_r8, 1.0e-6_r8, 2.5e-6_r8, 5.0e-6_r8, 10.0e-6_r8 /)
real(r8), parameter :: dmt_slt_opt = 75.0e-6_r8 ! [m] Optim diam for saltation
real(r8), parameter :: dns_slt = 2650.0_r8 ! [kg m-3] Density of optimal saltation particles
! declare erf intrinsic function
real(r8) :: dum ! dummy variable for erf test
#if (defined AIX)
#define ERF erf
#else
#define ERF derf
real(r8) derf
#endif
integer :: begp, endp ! per-proc beginning and ending pft indices
integer :: begc, endc ! per-proc beginning and ending column indices
integer :: begl, endl ! per-proc beginning and ending landunit indices
integer :: begg, endg ! per-proc gridcell ending gridcell indices
!------------------------------------------------------------------------
! Assign local pointers to derived type scalar members (column-level)
mbl_bsn_fct => clm3%g%l%c%cps%mbl_bsn_fct
! Sanity check on erf: erf() in SGI /usr/lib64/mips4/libftn.so is bogus
dum = 1.0_r8
if (abs(0.8427_r8-ERF(dum))/0.8427_r8>0.001_r8) then
write(iulog,*) 'erf(1.0) = ',ERF(dum)
write(iulog,*) 'Dustini: Error function error'
call endrun
end if
dum = 0.0_r8
if (ERF(dum) /= 0.0_r8) then
write(iulog,*) 'erf(0.0) = ',ERF(dum)
write(iulog,*) 'Dustini: Error function error'
call endrun
end if
! the following comes from (1) szdstlgn.F subroutine ovr_src_snk_frc_get
! and (2) dstszdst.F subroutine dst_szdst_ini
! purpose(1): given one set (the "source") of lognormal distributions,
! and one set of bin boundaries (the "sink"), compute and return
! the overlap factors between the source and sink distributions
! purpose(2): set important statistics of size distributions
do m = 1, dst_src_nbr
sqrt2lngsdi = sqrt(2.0_r8) * log(gsd_anl_src(m))
do n = 1, ndst
lndmaxjovrdmdni = log(dmt_grd(n+1)/dmt_vma_src(m))
lndminjovrdmdni = log(dmt_grd(n )/dmt_vma_src(m))
ovr_src_snk_frc = 0.5_r8 * (ERF(lndmaxjovrdmdni/sqrt2lngsdi) - &
ERF(lndminjovrdmdni/sqrt2lngsdi))
ovr_src_snk_mss(m,n) = ovr_src_snk_frc * mss_frc_src(m)
end do
end do
! The following code from subroutine wnd_frc_thr_slt_get was placed
! here because tmp1 needs to be defined just once
ryn_nbr_frc_thr_prx_opt = 0.38_r8 + 1331.0_r8 * (100.0_r8*dmt_slt_opt)**1.56_r8
if (ryn_nbr_frc_thr_prx_opt < 0.03_r8) then
write(iulog,*) 'dstmbl: ryn_nbr_frc_thr_prx_opt < 0.03'
call endrun
else if (ryn_nbr_frc_thr_prx_opt < 10.0_r8) then
ryn_nbr_frc_thr_opt_fnc = -1.0_r8 + 1.928_r8 * (ryn_nbr_frc_thr_prx_opt**0.0922_r8)
ryn_nbr_frc_thr_opt_fnc = 0.1291_r8 * 0.1291_r8 / ryn_nbr_frc_thr_opt_fnc
else
ryn_nbr_frc_thr_opt_fnc = 1.0_r8 - 0.0858_r8 * exp(-0.0617_r8*(ryn_nbr_frc_thr_prx_opt-10.0_r8))
ryn_nbr_frc_thr_opt_fnc = 0.120_r8 * 0.120_r8 * ryn_nbr_frc_thr_opt_fnc * ryn_nbr_frc_thr_opt_fnc
end if
icf_fct = 1.0_r8 + 6.0e-07_r8 / (dns_slt * grav * (dmt_slt_opt**2.5_r8))
dns_fct = dns_slt * grav * dmt_slt_opt
tmp1 = sqrt(icf_fct * dns_fct * ryn_nbr_frc_thr_opt_fnc)
! Set basin factor to 1 for now
call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp)
do c = begc, endc
l = clm3%g%l%c%landunit(c)
if (.not. clm3%g%l%lakpoi(l)) then
mbl_bsn_fct(c) = 1.0_r8
end if
end do
! Introducing particle diameter. Needed by atm model and by dry dep model.
! Taken from Charlie Zender's subroutines dst_psd_ini, dst_sz_rsl,
! grd_mk (dstpsd.F90) and subroutine lgn_evl (psdlgn.F90)
! Charlie allows logarithmic or linear option for size distribution
! however, he hardwires the distribution to logarithmic in his code
! therefore, I take his logarithmic code only
! furthermore, if dst_nbr == 4, he overrides the automatic grid calculation
! he currently works with dst_nbr = 4, so I only take the relevant code
! if ndst ever becomes different from 4, must add call grd_mk (dstpsd.F90)
! as done in subroutine dst_psd_ini
! note that here ndst = dst_nbr
! Override automatic grid with preset grid if available
if (ndst == 4) then
do n = 1, ndst
dmt_min(n) = dmt_grd(n) ![m] Max diameter in bin
dmt_max(n) = dmt_grd(n+1) ![m] Min diameter in bin
dmt_ctr(n) = 0.5_r8 * (dmt_min(n)+dmt_max(n)) ![m] Diameter at bin ctr
dmt_dlt(n) = dmt_max(n)-dmt_min(n) ![m] Width of size bin
end do
else
write(iulog,*) 'Dustini error: ndst must equal to 4 with current code'
call endrun !see more comments above
end if !end if ndst == 4
! Bin physical properties
gsd_anl = 2.0_r8 ! [frc] Geometric std dev PaG77 p. 2080 Table1
ln_gsd = log(gsd_anl)
dns_aer = 2.5e+3_r8 ! [kg m-3] Aerosol density
! Set a fundamental statistic for each bin
dmt_vma = 3.5000e-6_r8 ! [m] Mass median diameter analytic She84 p.75 Table1
! Compute analytic size statistics
! Convert mass median diameter to number median diameter (call vma2nma)
dmt_nma = dmt_vma * exp(-3.0_r8*ln_gsd*ln_gsd) ! [m]
! Compute resolved size statistics for each size distribution
! In C. Zender's code call dst_sz_rsl
do n = 1, ndst
series_ratio = (dmt_max(n)/dmt_min(n))**(1.0_r8/sz_nbr)
sz_min(1) = dmt_min(n)
do m = 2, sz_nbr ! Loop starts at 2
sz_min(m) = sz_min(m-1) * series_ratio
end do
! Derived grid values
do m = 1, sz_nbr-1 ! Loop ends at sz_nbr-1
sz_max(m) = sz_min(m+1) ! [m]
end do
sz_max(sz_nbr) = dmt_max(n) ! [m]
! Final derived grid values
do m = 1, sz_nbr
sz_ctr(m) = 0.5_r8 * (sz_min(m)+sz_max(m))
sz_dlt(m) = sz_max(m)-sz_min(m)
end do
lngsdsqrttwopi_rcp = 1.0_r8 / (ln_gsd*sqrt(2.0_r8*SHR_CONST_PI))
dmt_vwr(n) = 0.0_r8 ! [m] Mass wgted diameter resolved
vlm_rsl(n) = 0.0_r8 ! [m3 m-3] Volume concentration resolved
do m = 1, sz_nbr
! Evaluate lognormal distribution for these sizes (call lgn_evl)
tmp = log(sz_ctr(m)/dmt_nma) / ln_gsd
lgn_dst = lngsdsqrttwopi_rcp * exp(-0.5_r8*tmp*tmp) / sz_ctr(m)
! Integrate moments of size distribution
dmt_vwr(n) = dmt_vwr(n) + sz_ctr(m) * &
SHR_CONST_PI / 6.0_r8 * (sz_ctr(m)**3.0_r8) * & ![m3] Volume
lgn_dst * sz_dlt(m) ![# m-3] Number concentrn
vlm_rsl(n) = vlm_rsl(n) + &
SHR_CONST_PI / 6.0_r8 * (sz_ctr(m)**3.0_r8) * & ![m3] Volume
lgn_dst * sz_dlt(m) ![# m-3] Number concentrn
end do
dmt_vwr(n) = dmt_vwr(n) / vlm_rsl(n) ![m] Mass weighted diameter resolved
end do
! calculate correction to Stokes' settling velocity (subroutine stk_crc_get)
eps_max = 1.0e-4_r8
dns_mdp = 100000._r8 / (295.0_r8*SHR_CONST_RDAIR) ![kg m-3] const prs_mdp & tpt_vrt
! Size-independent thermokinetic properties
vsc_dyn_atm = 1.72e-5_r8 * ((295.0_r8/273.0_r8)**1.5_r8) * 393.0_r8 / &
(295.0_r8+120.0_r8) ![kg m-1 s-1] RoY94 p.102 tpt_mdp=295.0
mfp_atm = 2.0_r8 * vsc_dyn_atm / & !SeP97 p. 455 constant prs_mdp, tpt_mdp
(100000._r8*sqrt(8.0_r8/(SHR_CONST_PI*SHR_CONST_RDAIR*295.0_r8)))
vsc_knm_atm = vsc_dyn_atm / dns_mdp ![m2 s-1] Kinematic viscosity of air
do m = 1, ndst
slp_crc(m) = 1.0_r8 + 2.0_r8 * mfp_atm * &
(1.257_r8+0.4_r8*exp(-1.1_r8*dmt_vwr(m)/(2.0_r8*mfp_atm))) / &
dmt_vwr(m) ! [frc] Slip correction factor SeP97 p.464
vlc_stk(m) = (1.0_r8/18.0_r8) * dmt_vwr(m) * dmt_vwr(m) * dns_aer * &
grav * slp_crc(m) / vsc_dyn_atm ! [m s-1] SeP97 p.466
end do
! For Reynolds number flows Re < 0.1 Stokes' velocity is valid for
! vlc_grv SeP97 p. 466 (8.42). For larger Re, inertial effects become
! important and empirical drag coefficients must be employed
! Implicit equation for Re, Cd, and Vt is SeP97 p. 467 (8.44)
! Using Stokes' velocity rather than iterative solution with empirical
! drag coefficient causes 60% errors for D = 200 um SeP97 p. 468
! Iterative solution for drag coefficient, Reynolds number, and terminal veloc
do m = 1, ndst
! Initialize accuracy and counter
eps_crr = eps_max + 1.0_r8 ![frc] Current relative accuracy
itr_idx = 0 ![idx] Counting index
! Initial guess for vlc_grv is exact for Re < 0.1
vlc_grv(m) = vlc_stk(m) ![m s-1]
do while(eps_crr > eps_max)
! Save terminal velocity for convergence test
vlc_grv_old = vlc_grv(m) ![m s-1]
ryn_nbr_grv(m) = vlc_grv(m) * dmt_vwr(m) / vsc_knm_atm !SeP97 p.460
! Update drag coefficient based on new Reynolds number
if (ryn_nbr_grv(m) < 0.1_r8) then
cff_drg_grv(m) = 24.0_r8 / ryn_nbr_grv(m) !Stokes' law Sep97 p.463 (8.32)
else if (ryn_nbr_grv(m) < 2.0_r8) then
cff_drg_grv(m) = (24.0_r8/ryn_nbr_grv(m)) * &
(1.0_r8 + 3.0_r8*ryn_nbr_grv(m)/16.0_r8 + &
9.0_r8*ryn_nbr_grv(m)*ryn_nbr_grv(m)* &
log(2.0_r8*ryn_nbr_grv(m))/160.0_r8) !Sep97 p.463 (8.32)
else if (ryn_nbr_grv(m) < 500.0_r8) then
cff_drg_grv(m) = (24.0_r8/ryn_nbr_grv(m)) * &
(1.0_r8 + 0.15_r8*ryn_nbr_grv(m)**0.687_r8) !Sep97 p.463 (8.32)
else if (ryn_nbr_grv(m) < 2.0e5_r8) then
cff_drg_grv(m) = 0.44_r8 !Sep97 p.463 (8.32)
else
write(iulog,'(a,es9.2)') "ryn_nbr_grv(m) = ",ryn_nbr_grv(m)
write(iulog,*)'Dustini error: Reynolds number too large in stk_crc_get()'
call endrun
end if
! Update terminal velocity based on new Reynolds number and drag coeff
! [m s-1] Terminal veloc SeP97 p.467 (8.44)
vlc_grv(m) = sqrt(4.0_r8 * grav * dmt_vwr(m) * slp_crc(m) * dns_aer / &
(3.0_r8*cff_drg_grv(m)*dns_mdp))
eps_crr = abs((vlc_grv(m)-vlc_grv_old)/vlc_grv(m)) !Relative convergence
if (itr_idx == 12) then
! Numerical pingpong may occur when Re = 0.1, 2.0, or 500.0
! due to discontinuities in derivative of drag coefficient
vlc_grv(m) = 0.5_r8 * (vlc_grv(m)+vlc_grv_old) ! [m s-1]
end if
if (itr_idx > 20) then
write(iulog,*) 'Dustini error: Terminal velocity not converging ',&
' in stk_crc_get(), breaking loop...'
goto 100 !to next iteration
end if
itr_idx = itr_idx + 1
end do !end while
100 continue !Label to jump to when iteration does not converge
end do !end loop over size
! Compute factors to convert Stokes' settling velocities to
! actual settling velocities
do m = 1, ndst
stk_crc(m) = vlc_grv(m) / vlc_stk(m)
end do
end subroutine Dustini
#endif
end module DUSTMod