!>\file dust_fengsha_mod.F90
!! This file contains the FENGSHA dust scheme.
module dust_fengsha_mod
!
! This module developed by Barry Baker (NOAA ARL)
! For serious questions contact barry.baker@noaa.gov
!
! 07/16/2019 - Adapted for NUOPC/GOCART, R. Montuoro
! 02/01/2020 - Adapted for FV3/CCPP, Haiqin Li
use rrfs_smoke_data
use machine , only : kind_phys
use dust_data_mod
implicit none
private
public :: gocart_dust_fengsha_driver
contains
subroutine gocart_dust_fengsha_driver(data, dt, &
chem,rho_phy,smois,p8w,ssm, &
isltyp,vegfra,snowh,xland,area,g,emis_dust, &
ust,znt,clay,sand,rdrag,uthr, &
num_emis_dust,num_moist,num_chem,num_soil_layers, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
IMPLICIT NONE
type(smoke_data), intent(inout) :: data
INTEGER, INTENT(IN ) :: &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
num_emis_dust,num_moist,num_chem,num_soil_layers
INTEGER,DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: isltyp
REAL(kind_phys), DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), INTENT(INOUT) :: chem
REAL(kind_phys), DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL, INTENT(INOUT) :: emis_dust
REAL(kind_phys), DIMENSION( ims:ime, num_soil_layers, jms:jme ), INTENT(IN) :: smois
REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: ssm
REAL(kind_phys), DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: vegfra, &
snowh, &
xland, &
area, &
ust, &
znt, &
clay, &
sand, &
rdrag, &
uthr
REAL(kind_phys), DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: &
p8w, &
rho_phy
REAL(kind_phys), INTENT(IN) :: dt,g
! Local variables
integer :: nmx,smx,i,j,k,imx,jmx,lmx
integer,dimension (1,1) :: ilwi
real(kind_phys), DIMENSION (1,1) :: erodtot
REAL(kind_phys), DIMENSION (1,1) :: gravsm
REAL(kind_phys), DIMENSION (1,1) :: drylimit
real(kind_phys), DIMENSION (5) :: tc,bems
real(kind_phys), dimension (1,1) :: airden,airmas,ustar
real(kind_phys), dimension (1) :: dxy
real(kind_phys), dimension (3) :: massfrac
real(kind_phys) :: conver,converi
real(kind_phys) :: R
! threshold values
conver=1.e-9
converi=1.e9
! Number of dust bins
imx=1
jmx=1
lmx=1
nmx=ndust
smx=nsalt
k=kts
do j=jts,jte
do i=its,ite
! Don't do dust over water!!!
ilwi(1,1)=0
if(xland(i,j).lt.1.5)then
ilwi(1,1)=1
! Total concentration at lowest model level. This is still hardcoded for 5 bins.
! if(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) then
! tc(:)=1.e-16*conver
! else
tc(1)=chem(i,kts,j,p_dust_1)*conver
tc(2)=chem(i,kts,j,p_dust_2)*conver
tc(3)=chem(i,kts,j,p_dust_3)*conver
tc(4)=chem(i,kts,j,p_dust_4)*conver
tc(5)=chem(i,kts,j,p_dust_5)*conver
! endif
! Air mass and density at lowest model level.
airmas(1,1)=-(p8w(i,kts+1,j)-p8w(i,kts,j))*area(i,j)/g
airden(1,1)=rho_phy(i,kts,j)
ustar(1,1)=ust(i,j)
dxy(1)=area(i,j)
! Mass fractions of clay, silt, and sand.
massfrac(1)=clay(i,j)
massfrac(2)=1-(clay(i,j)+sand(i,j))
massfrac(3)=sand(i,j)
! Total erodibility.
erodtot(1,1) = ssm(i,j) ! SUM(erod(i,j,:))
! Don't allow roughness lengths greater than 20 cm to be lofted.
! This kludge accounts for land use types like urban areas and
! forests which would otherwise show up as high dust emitters.
! This is a placeholder for a more widely accepted kludge
! factor in the literature, which reduces lofting for rough areas.
! Forthcoming...
IF (znt(i,j) .gt. 0.2) then
ilwi(1,1)=0
endif
! limit where there is lots of vegetation
if (vegfra(i,j) .gt. .17) then
ilwi(1,1) = 0
endif
! limit where there is snow on the ground
if (snowh(i,j) .gt. 0) then
ilwi(1,1) = 0
endif
! Do not allow areas with bedrock, lava, or land-ice to loft
IF (isltyp(i,j) .eq. 15 .or. isltyp(i,j) .eq. 16. .or. &
isltyp(i,j) .eq. 18) then
ilwi(1,1)=0
ENDIF
IF (isltyp(i,j) .eq. 0)then
ilwi(1,1)=0
endif
if(ilwi(1,1) == 0 ) cycle
! Calculate gravimetric soil moisture and drylimit.
gravsm(1,1)=100.*smois(i,1,j)/((1.-maxsmc(isltyp(i,j)))*(2.65*(1.-clay(i,j))+2.50*clay(i,j)))
drylimit(1,1)=14.0*clay(i,j)*clay(i,j)+17.0*clay(i,j)
! get drag partition
! FENGSHA uses the drag partition correction of MacKinnon et al 2004
! doi:10.1016/j.geomorph.2004.03.009
if (dust_calcdrag .ne. 1) then
call fengsha_drag(data,znt(i,j),R)
else
! use the precalculated version derived from ASCAT; Prigent et al. (2012,2015)
! doi:10.1109/TGRS.2014.2338913 & doi:10.5194/amt-5-2703-2012
! pick only valid values
if (rdrag(i,j) > 0.) then
R = real(rdrag(i,j), kind=kind_phys)
else
cycle
endif
endif
! Call dust emission routine.
call source_dust(data, imx, jmx, lmx, nmx, smx, dt, tc, ustar, massfrac, &
erodtot, dxy, gravsm, airden, airmas, &
bems, g, drylimit, dust_alpha, dust_gamma, R, uthr(i,j))
! if(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) then
! dustin(i,j,1:5)=tc(1:5)*converi
! else
chem(i,kts,j,p_dust_1)=tc(1)*converi
chem(i,kts,j,p_dust_2)=tc(2)*converi
chem(i,kts,j,p_dust_3)=tc(3)*converi
chem(i,kts,j,p_dust_4)=tc(4)*converi
chem(i,kts,j,p_dust_5)=tc(5)*converi
! endif
! chem(i,kts,j,p_dust_1)=tc(1)*converi
! chem(i,kts,j,p_dust_2)=tc(2)*converi
! chem(i,kts,j,p_dust_3)=tc(3)*converi
! chem(i,kts,j,p_dust_4)=tc(4)*converi
! chem(i,kts,j,p_dust_5)=tc(5)*converi
! For output diagnostics
emis_dust(i,1,j,p_edust1)=bems(1)
emis_dust(i,1,j,p_edust2)=bems(2)
emis_dust(i,1,j,p_edust3)=bems(3)
emis_dust(i,1,j,p_edust4)=bems(4)
emis_dust(i,1,j,p_edust5)=bems(5)
endif
enddo
enddo
!
end subroutine gocart_dust_fengsha_driver
SUBROUTINE source_dust(data, imx, jmx, lmx, nmx, smx, dt1, tc, ustar, massfrac, &
erod, dxy, gravsm, airden, airmas, bems, g0, drylimit, alpha, &
gamma, R, uthres)
! ****************************************************************************
! * Evaluate the source of each dust particles size bin by soil emission
! *
! * Input:
! * EROD Fraction of erodible grid cell (-)
! * GRAVSM Gravimetric soil moisture (g/g)
! * DRYLIMIT Upper GRAVSM limit for air-dry soil (g/g)
! * ALPHA Constant to fudge the total emission of dust (1/m)
! * GAMMA Tuning constant for erodibility (-)
! * DXY Surface of each grid cell (m2)
! * AIRMAS Mass of air for each grid box (kg)
! * AIRDEN Density of air for each grid box (kg/m3)
! * USTAR Friction velocity (m/s)
! * DT1 Time step (s)
! * NMX Number of dust bins (-)
! * SMX Number of saltation bins (-)
! * IMX Number of I points (-)
! * JMX Number of J points (-)
! * LMX Number of L points (-)
! * R Drag Partition (-)
! * UTHRES FENGSHA Dry Threshold Velocities (m/s)
! *
! * Data:
! * MASSFRAC Fraction of mass in each of 3 soil classes (-)
! * SPOINT Pointer to 3 soil classes (-)
! * DEN_DUST Dust density (kg/m3)
! * DEN_SALT Saltation particle density (kg/m3)
! * REFF_SALT Reference saltation particle diameter (m)
! * REFF_DUST Reference dust particle diameter (m)
! * LO_DUST Lower diameter limits for dust bins (m)
! * UP_DUST Upper diameter limits for dust bins (m)
! * FRAC_SALT Soil class mass fraction for saltation bins (-)
! *
! * Parameters:
! * CMB Constant of proportionality (-)
! * MMD_DUST Mass median diameter of dust (m)
! * GSD_DUST Geometric standard deviation of dust (-)
! * LAMBDA Side crack propagation length (m)
! * CV Normalization constant (-)
! * G0 Gravitational acceleration (m/s2)
! * G Gravitational acceleration in cgs (cm/s2)
! *
! * Working:
! * U_TS0 "Dry" threshold friction velocity (m/s)
! * U_TS Moisture-adjusted threshold friction velocity (m/s)
! * RHOA Density of air in cgs (g/cm3)
! * DEN Dust density in cgs (g/cm3)
! * DIAM Dust diameter in cgs (cm)
! * DMASS Saltation mass distribution (-)
! * DSURFACE Saltation surface area per unit mass (m2/kg)
! * DS_REL Saltation surface area distribution (-)
! * SALT Saltation flux (kg/m/s)
! * DLNDP Dust bin width (-)
! * EMIT Total vertical mass flux (kg/m2/s)
! * EMIT_VOL Total vertical volume flux (m/s)
! * DSRC Mass of emitted dust (kg/timestep/cell)
! *
! * Output:
! * TC Total concentration of dust (kg/kg/timestep/cell)
! * BEMS Source of each dust type (kg/timestep/cell)
! *
! ****************************************************************************
implicit none
type(smoke_data), intent(inout) :: data
INTEGER, INTENT(IN) :: imx,jmx,lmx,nmx,smx
REAL(kind_phys), INTENT(IN) :: dt1
REAL(kind_phys), INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
REAL(kind_phys), INTENT(IN) :: ustar(imx,jmx)
REAL(kind_phys), INTENT(IN) :: massfrac(3)
REAL(kind_phys), INTENT(IN) :: erod(imx,jmx)
REAL(kind_phys), INTENT(IN) :: dxy(jmx)
REAL(kind_phys), INTENT(IN) :: gravsm(imx,jmx)
REAL(kind_phys), INTENT(IN) :: airden(imx,jmx,lmx)
REAL(kind_phys), INTENT(IN) :: airmas(imx,jmx,lmx)
REAL(kind_phys), INTENT(OUT) :: bems(imx,jmx,nmx)
REAL(kind_phys), INTENT(IN) :: g0
REAL(kind_phys), INTENT(IN) :: drylimit(imx,jmx)
!! Sandblasting mass efficiency, aka "fudge factor" (based on Tegen et al,
!! 2006 and Hemold et al, 2007)
!
! REAL, PARAMETER :: alpha=1.8E-8 ! (m^-1)
REAL(kind_phys), INTENT(IN) :: alpha
! Experimental optional exponential tuning constant for erodibility.
! 0 < gamma < 1 -> more relative impact by low erodibility regions.
REAL(kind_phys), INTENT(IN) :: gamma
REAL(kind_phys), INTENT(IN) :: R
REAL(kind_phys), INTENT(IN) :: uthres
REAL(kind_phys) :: den(smx), diam(smx)
REAL(kind_phys) :: dvol(nmx), distr_dust(nmx), dlndp(nmx)
REAL(kind_phys) :: dsurface(smx), ds_rel(smx)
REAL(kind_phys) :: u_ts0, u_ts, dsrc, dmass, dvol_tot
REAL(kind_phys) :: salt,emit, emit_vol, stotal
REAL(kind_phys) :: rhoa, g
INTEGER :: i, j, n
! Sandblasting mass efficiency, beta.
! Beta maxes out for clay fractions above 0.2 = betamax.
REAL(kind_phys), PARAMETER :: betamax=5.25E-4
REAL(kind_phys) :: beta
integer :: styp
! Constant of proportionality from Marticorena et al, 1997 (unitless)
! Arguably more ~consistent~ fudge than alpha, which has many walnuts
! sprinkled throughout the literature. - GC
REAL(kind_phys), PARAMETER :: cmb=1.0
! REAL, PARAMETER :: cmb=2.61 ! from White,1979
! Parameters used in Kok distribution function. Advise not to play with
! these without the expressed written consent of someone who knows what
! they're doing. - GC
REAL(kind_phys), PARAMETER :: mmd_dust=3.4D-6 ! median mass diameter (m)
REAL(kind_phys), PARAMETER :: gsd_dust=3.0 ! geom. std deviation
REAL(kind_phys), PARAMETER :: lambda=12.0D-6 ! crack propagation length (m)
REAL(kind_phys), PARAMETER :: cv=12.62D-6 ! normalization constant
! Calculate saltation surface area distribution from sand, silt, and clay
! mass fractions and saltation bin fraction. This will later become a
! modifier to the total saltation flux. The reasoning here is that the
! size and availability of saltators affects saltation efficiency. Based
! on Eqn. (32) in Marticorena & Bergametti, 1995 (hereon, MB95).
DO n=1,smx
dmass=massfrac(spoint(n))*frac_salt(n)
dsurface(n)=0.75*dmass/(den_salt(n)*reff_salt(n))
ENDDO
! The following equation yields relative surface area fraction. It will only
! work if you are representing the "full range" of all three soil classes.
! For this reason alone, we have incorporated particle sizes that encompass
! the clay class, to account for the its relative area over the basal
! surface, even though these smaller bins would be unlikely to play any large
! role in the actual saltation process. - GC
stotal=SUM(dsurface(:))
DO n=1,smx
ds_rel(n)=dsurface(n)/stotal
ENDDO
! Calculate total dust emission due to saltation of sand sized particles.
! Begin by calculating DRY threshold friction velocity (u_ts0). Next adjust
! u_ts0 for moisture to get threshold friction velocity (u_ts). Then
! calculate saltation flux (salt) where ustar has exceeded u_ts. Finally,
! calculate total dust emission (tot_emit), taking into account erodibility.
! Set DRY threshold friction velocity to input value
u_ts0 = uthres
g = g0*1.0E2
emit=0.0
DO n = 1, smx
den(n) = den_salt(n)*1.0D-3 ! (g cm^-3)
diam(n) = 2.0*reff_salt(n)*1.0D2 ! (cm)
DO i = 1,imx
DO j = 1,jmx
rhoa = airden(i,j,1)*1.0D-3 ! (g cm^-3)
! FENGSHA uses the 13 category soil type from the USDA
! call calc_fengsha_styp(massfrac(1),massfrac(3),massfrac(2),styp)
! Fengsha uses threshold velocities based on dale gilletes data
! call fengsha_utst(styp,uthres,u_ts0)
! Friction velocity threshold correction function based on physical
! properties related to moisture tension. Soil moisture greater than
! dry limit serves to increase threshold friction velocity (making
! it more difficult to loft dust). When soil moisture has not reached
! dry limit, treat as dry
IF (gravsm(i,j) > drylimit(i,j)) THEN
u_ts = MAX(0.0D+0,u_ts0*(sqrt(1.0+1.21*(gravsm(i,j)-drylimit(i,j))**0.68)) / R)
ELSE
u_ts = u_ts0 / R
END IF
! Calculate total vertical mass flux (note beta has units of m^-1)
! Beta acts to tone down dust in areas with so few dust-sized particles that the
! lofting efficiency decreases. Otherwise, super sandy zones would be huge dust
! producers, which is generally not the case. Equation derived from wind-tunnel
! experiments (see MB95).
beta=10**(13.6*massfrac(1)-6.0) ! (unitless)
if (massfrac(1) <= 0.2) then
beta=10**(13.4*massfrac(1)-6.0)
else
beta = 2.E-4
endif
!---------------------------------------------------------------------
! formula of Draxler & Gillette (2001) Atmos. Environ.
! F = K A (r/g) U* ( U*^2 - Ut*^2 )
!
! where:
! F = vertical emission flux [g/m**2-s]
! K = constant 2.0E-04 [1/m]
! A = 0~3.5 mean = 2.8 (fudge factor)
! U* = friction velocity [m/s]
! Ut* = threshold friction velocity [m/s]
!
!--------------------------------------------------------------------
IF (ustar(i,j) .gt. u_ts) then
call fengsha_hflux(data,ustar(i,j),u_ts,beta, salt)
salt = alpha * cmb * ds_rel(n) * airden(i,j,1) / g0 * salt * (erod(i,j)**gamma) * beta
else
salt = 0.
endif
! EROD is taken into account above
emit = emit + salt
END DO
END DO
END DO
! Now that we have the total dust emission, distribute into dust bins using
! lognormal distribution (Dr. Jasper Kok, in press), and
! calculate total mass emitted over the grid box over the timestep.
!
! In calculating the Kok distribution, we assume upper and lower limits to each bin.
! For reff_dust=(/0.73D-6,1.4D-6,2.4D-6,4.5D-6,8.0D-6/) (default),
! lower limits were ASSUMED at lo_dust=(/0.1D-6,1.0D-6,1.8D-6,3.0D-6,6.0D-6/)
! upper limits were ASSUMED at up_dust=(/1.0D-6,1.8D-6,3.0D-6,6.0D-6,10.0D-6/)
! These may be changed within module_data_gocart_dust.F, but make sure it is
! consistent with reff_dust values. These values were taken from the original
! GOCART bin configuration. We use them here to calculate dust bin width, dlndp.
! dVol is the volume distribution. You know...if you were wondering. GC
dvol_tot=0.
DO n=1,nmx
dlndp(n)=LOG(up_dust(n)/lo_dust(n))
dvol(n)=(2.0*reff_dust(n)/cv)*(1.+ERF(LOG(2.0*reff_dust(n)/mmd_dust)/(SQRT(2.)*LOG(gsd_dust))))*&
EXP(-(2.0*reff_dust(n)/lambda)**3.0)*dlndp(n)
dvol_tot=dvol_tot+dvol(n)
! Convert mass flux to volume flux
!emit_vol=emit/den_dust(n) ! (m s^-1)
END DO
DO n=1,nmx
distr_dust(n)=dvol(n)/dvol_tot
!print *,"distr_dust(",n,")=",distr_dust(n)
END DO
! Now distribute total vertical emission into dust bins and update concentration.
DO n=1,nmx
DO i=1,imx
DO j=1,jmx
! Calculate total mass emitted
dsrc = emit*distr_dust(n)*dxy(j)*dt1 ! (kg)
IF (dsrc < 0.0) dsrc = 0.0
! Update dust mixing ratio at first model level.
tc(i,j,1,n) = tc(i,j,1,n) + dsrc / airmas(i,j,1) ! (kg/kg)
! bems(i,j,n) = dsrc ! diagnostic
!bems(i,j,n) = 1000.*dsrc/(dxy(j)*dt1) ! diagnostic (g/m2/s)
bems(i,j,n) = 1.e+9*dsrc/(dxy(j)*dt1) ! diagnostic (ug/m2/s) !lzhang
END DO
END DO
END DO
END SUBROUTINE source_dust
subroutine fengsha_utst(data,styp,uth, ut)
implicit none
type(smoke_data), intent(inout) :: data
integer, intent(in) :: styp
real(kind_phys), dimension(fengsha_maxstypes), intent(in) :: uth
real(kind_phys), intent(out) :: ut
ut = uth(styp)
! real (kind_phys) :: uth(13) = &
! (/ 0.08, & ! Sand - 1
! 0.20, & ! Loamy Sand - 2
! 0.30, & ! Sandy Loam - 3
! 0.30, & ! Silt Loam - 4
! 0.35, & ! Silt - 5
! 0.60, & ! Loam - 6
! 0.30, & ! Sandy Clay Loam - 7
! 0.35, & ! Silty Clay Loam - 8
! 0.45, & ! Clay Loam - 9
! 0.45, & ! Sandy Clay - 10
! 0.45, & ! Silty Clay - 11
! 0.60, & ! Clay - 12
! 9.999 /) ! Other - 13
return
end subroutine fengsha_utst
subroutine calc_fengsha_styp(data, clay, sand, silt, type)
implicit none
type(smoke_data), intent(inout) :: data
!---------------------------------------------------------------
! Function: calculate soil type based on USDA definition.
! Source: USDA soil texture calculator
!
! Defintion of soil types:
!
!
! NOAH 1 2 3 4 5 6 7 8 9 10 11 12
! PX 1 2 3 4 - 5 6 7 8 9 10 11
! Soil "Sand" "Loamy Sand" "Sandy Loam" "Silt Loam" "Silt" "Loam" "Sandy Clay Loam" "Silt Clay Loam" "Clay Loam" "Sandy Clay" "Silty Clay" "Clay"
!---------------------------------------------------------------
REAL(kind_phys), intent(in) :: clay, sand, silt
integer, intent(out) :: type
real(kind_phys) :: cly, snd, slt
type = 0
snd = sand * 100.
cly = clay * 100.
slt = silt * 100.
if (slt+1.5*cly .lt. 15) type = 1 ! snd
if (slt+1.5*cly .ge. 15 .and.slt+1.5*cly .lt. 30) type = 2 ! loamy snd
if (cly .ge. 7 .and. cly .lt. 20 .and. snd .gt. 52 .and. slt+2*cly .ge. 30) type = 3 ! sndy loam (cond 1)
if (cly .lt. 7 .and. slt .lt. 50 .and. slt+2*cly .ge. 30) type = 3 ! sndy loam (cond 2)
if (slt .ge. 50 .and. cly .ge. 12 .and.cly .lt. 27 ) type = 4 ! slt loam (cond 1)
if (slt .ge. 50 .and. slt .lt. 80 .and.cly .lt. 12) type = 4 ! slt loam (cond 2)
if (slt .ge. 80 .and. cly .lt. 12) type = 5 ! slt
if (cly .ge. 7 .and. cly .lt. 27 .and.slt .ge. 28 .and. slt .lt. 50 .and.snd .le. 52) type = 6 ! loam
if (cly .ge. 20 .and. cly .lt. 35 .and.slt .lt. 28 .and. snd .gt. 45) type = 7 ! sndy cly loam
if (cly .ge. 27 .and. cly .lt. 40 .and.snd .lt. 20) type = 8 ! slt cly loam
if (cly .ge. 27 .and. cly .lt. 40 .and.snd .ge. 20 .and. snd .le. 45) type = 9 ! cly loam
if (cly .ge. 35 .and. snd .gt. 45) type = 10 ! sndy cly
if (cly .ge. 40 .and. slt .ge. 40) type = 11 ! slty cly
if (cly .ge. 40 .and. snd .le. 45 .and.slt .lt. 40) type = 12 ! clay
return
end subroutine calc_fengsha_styp
subroutine fengsha_drag(data,z0,R)
implicit none
type(smoke_data), intent(inout) :: data
real(kind_phys), intent(in) :: z0
real(kind_phys), intent(out) :: R
real(kind_phys), parameter :: z0s = 1.0e-04 !Surface roughness for ideal bare surface [m]
! ------------------------------------------------------------------------
! Function: Calculates the MacKinnon et al. 2004 Drag Partition Correction
!
! R = 1.0 - log(z0 / z0s) / log( 0.7 * (12255./z0s) ** 0.8)
!
!--------------------------------------------------------------------------
! Drag partition correction. See MacKinnon et al. (2004),
! doi:10.1016/j.geomorph.2004.03.009
R = 1.0 - log(z0 / z0s) / log( 0.7 * (12255./z0s) ** 0.8)
! Drag partition correction. See Marticorena et al. (1997),
! doi:10.1029/96JD02964
!R = 1.0 - log(z0 / z0s) / log( 0.7 * (10./z0s) ** 0.8)
return
end subroutine fengsha_drag
subroutine fengsha_hflux(data,ust,utst, kvh, salt)
!---------------------------------------------------------------------
! Function: Calculates the Horizontal Saltation Flux, Q, and then
! calculates the vertical flux.
!
! formula of Draxler & Gillette (2001) Atmos. Environ.
! F = K A (r/g) U* ( U*^2 - Ut*^2 )
!
! where:
! F = vertical emission flux [g/m**2-s]
! K = constant 2.0E-04 [1/m]
! A = 0~3.5 mean = 2.8 (fudge factor)
! U* = friction velocity [m/s]
! Ut* = threshold friction velocity [m/s]
!
!--------------------------------------------------------------------
implicit none
type(smoke_data), intent(inout) :: data
real(kind_phys), intent(in) :: ust, & ! friction velocity
utst, & ! threshold friction velocity
kvh ! vertical to horizontal mass flux ratio
real(kind_phys), intent(out) :: salt
real(kind_phys) :: Q
Q = ust * (ust * ust - utst * utst)
salt = Q ! sdep * kvh * Q
return
end subroutine fengsha_hflux
end module dust_fengsha_mod