MODULE module_uoc_dustwd
!----------------------------------------------------------------------------
! Dust wet deposition module of the University of Cologne, Germany.
! Dust wet deposition scheme developed by E Jung (2004, PhD Thesis).
! Implementation into WRF and modifications by C Frick, 2014 (claudia.frick@uni-koeln.de)
! 2014-2015, updates and modifications, Martina Klose (mklose@uni-koeln.de)
!----------------------------------------------------------------------------
USE module_state_description ! num_chem, p_qr, ...
USE module_model_constants ! rhowater in kg/m3
USE physconst ! rair, ...
USE module_data_gocart_dust
CONTAINS
subroutine uoc_dustwd_driver(precr,chem,p_phy,t_phy, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
dtstepc, &
dustwd_1, dustwd_2, &
dustwd_3, dustwd_4, &
dustwd_5, &
wetdep_1, wetdep_2, &
wetdep_3, wetdep_4, &
wetdep_5, &
dustwdload_1, dustwdload_2, &
dustwdload_3, dustwdload_4, &
dustwdload_5, &
alt, dz8w, epsilc )
IMPLICIT NONE
INTEGER :: debug_level
CHARACTER*(100) :: text
real :: dustold ! help variable
integer :: d, i, j, k, l ! loop
integer, parameter :: nbins=5 ! number of dust bins
integer, parameter :: nbinsa=5 ! number of dust bins
integer :: bins(nbinsa) ! dust bin numbers in chem
! real, dimension(nbins) :: dbin ! max. dimension of dust in the bin
! data dbin/2.5,5.,10.,20./ ! size cut diameter (um) - from qf03
real, dimension(nbins) :: dbinm ! mean dimension of dust (um) in the bin
real, dimension(nbins) :: dbinmm ! mean dimension of dust (m) in the bin
! data dbinm/1.25,3.75,7.5,15./ ! mean size cut diameter (um)
real :: conver9, converi9 ! transformation ug/kg-dryair to kg/kg-dryair and vice versa
real :: conver6, converi6 ! transformation um to m and vice versa
real :: wt ! terminal fall velocity of dust in the bin (m/s)
real :: rmin, rmax ! minimum/maximum raindrop diameter (m)
integer, parameter :: nrbins=30 ! number of raindrop bins
real, dimension(nrbins) :: raind ! raindrop diameter in the middle of the raindrop bin (m)
real, dimension(nrbins+1) :: rend ! raindrop diameter at the end of the raindrop bin (m)
real, dimension(nrbins) :: dsd_rn ! raindrop size distribution
real, dimension(nrbins) :: vt ! raindrop terminal fall velocity (m/s)
real :: z ! help parameter for the determination of the rain bins
real :: visca ! (dynamic) viscosity of air (g/cm s)
real :: tw, viscw ! temperature of water (C) and (dynamic) viscosity of water (g/cm s)
real :: colece ! collection efficiency
real :: scrate, scavn ! scavanging rate (1/s)
real :: delR, delv, rnflx, carea ! help parameter for the determination of the scavanging rate
real :: tair, pair ! current air temperature and pressure
real :: rhoair ! dry air density (kg/m3) - rhoair=pa/(Ra*Ta)
INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ! domain grid
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
REAL, INTENT(IN) :: dtstepc, & ! time step (s)
epsilc
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN) :: precr
! precr - rain precipitation rate at all levels (kg/m2 s)
REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
INTENT(INOUT) :: chem
! chem(i,k,j,p_dust_?) - dust mixing ratio for bin np_dust_? (ug/kg-dryair)
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN) :: p_phy, t_phy
! pressure (Pa) and temperature (K)
real, dimension( ims:ime, kms:kme, jms:jme ) :: rnrate
! precipitation rate (mm/h) - rnrate=precr/3600 using rhowater (1l = 1kg)
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: &
dustwd_1, dustwd_2, dustwd_3, dustwd_4, dustwd_5
! loss in dust mixing ratio due to wet deposition (ug/kg-dryair) (current time step)
REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ) :: wdbins
! loss in dust mixing ratio due to wet deposition - help variable (ug/kg-dryair)
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
dustwdload_1, dustwdload_2, dustwdload_3, dustwdload_4, dustwdload_5
! loss in dustload due to wet deposition (ug/m2) (current time step)
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
wetdep_1, wetdep_2, wetdep_3, wetdep_4, wetdep_5
! loss in dust mixing ratio due to wet deposition in one coulmn for all size bins (ug/m2/s)
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: alt, dz8w
! altitude = 1/rhoair (kg/m3) & vertical grid-spacing delta-z / dz between full levels (m)
den_dust(:) = 2650.
dbinm(:) = 2.0d0*reff_dust(:)
CALL get_wrf_debug_level(debug_level)
rmin = 100.*1.e-6
rmax = 5800.*1.e-6 ! minimum/maximum raindrop diameter (m)
conver6 = 1.e-6 ! transformation um to m
converi6 = 1.e6 ! transformation m to um
conver9 = 1.e-9 ! transformation ug/kg-dryair to kg/kg-dryair
converi9 = 1.e9 ! transformation kg/kg-dryair to ug/kg-dryair
bins(1)=p_dust_1 ! dust bin number in chem
bins(2)=p_dust_2 ! dust bin number in chem
bins(3)=p_dust_3 ! dust bin number in chem
bins(4)=p_dust_4 ! dust bin number in chem
bins(5)=p_dust_5 ! dust bin number in chem
dustwd_1(its:ite,kts:kte,jts:jte) = 0.
dustwd_2(its:ite,kts:kte,jts:jte) = 0.
dustwd_3(its:ite,kts:kte,jts:jte) = 0.
dustwd_4(its:ite,kts:kte,jts:jte) = 0.
dustwd_5(its:ite,kts:kte,jts:jte) = 0.
wetdep_1(its:ite,jts:jte) = 0.
wetdep_2(its:ite,jts:jte) = 0.
wetdep_3(its:ite,jts:jte) = 0.
wetdep_4(its:ite,jts:jte) = 0.
wetdep_5(its:ite,jts:jte) = 0.
dustwdload_1(its:ite,jts:jte) = 0.
dustwdload_2(its:ite,jts:jte) = 0.
dustwdload_3(its:ite,jts:jte) = 0.
dustwdload_4(its:ite,jts:jte) = 0.
dustwdload_5(its:ite,jts:jte) = 0.
wdbins(its:ite,kts:kte,jts:jte,1)=dustwd_1(its:ite,kts:kte,jts:jte) ! dust dust loss in dust bin 1
wdbins(its:ite,kts:kte,jts:jte,2)=dustwd_2(its:ite,kts:kte,jts:jte) ! dust dust loss in dust bin 2
wdbins(its:ite,kts:kte,jts:jte,3)=dustwd_3(its:ite,kts:kte,jts:jte) ! dust dust loss in dust bin 3
wdbins(its:ite,kts:kte,jts:jte,4)=dustwd_4(its:ite,kts:kte,jts:jte) ! dust dust loss in dust bin 4
wdbins(its:ite,kts:kte,jts:jte,5)=dustwd_5(its:ite,kts:kte,jts:jte) ! dust dust loss in dust bin 5
rnrate = precr*3600. ! precipitation rate (mm/h) - rnrate=precr*3600 using rhowater
dbinmm = dbinm*conver6 ! mean dimension of dust (m) in the bin
tw=0 ! water temperature (C)
call dviscw(tw,viscw) ! determine water (dynamic) viscosity viscw
! raindrop classes
do l = 1,nrbins+1
z = alog10(rmin*converi6)+(alog10(rmax*converi6)-alog10(rmin*converi6))*real(l-1)/real(nrbins)
rend(l) = real(int(10.**z))*conver6
enddo
do l = 1,nrbins
z = (alog10(rend(l)*converi6)+alog10(rend(l+1)*converi6))*0.5
raind(l) = real(int(10.**z))*conver6
enddo
! bin 1 bin 2 bin nrbin
! |--------------|------------!-------------!--------------|
! rend(1) rend(2) rend(3) rend(nrbin+1)
! raind(1) raind(2) raind(nrbin)
do j = jts, jte
do k = kts, kte
do i = its, ite
if ( precr(i,k,j).gt.0.) then ! precipitation exists
call dsd_rain(rnrate(i,k,j),raind,nrbins,dsd_rn,3) ! 3 = Gamma Distribution
tair=t_phy(i,k,j)
pair=p_phy(i,k,j)
rhoair = pair/(rair*tair) ! dry air density (kg/m3) - rhoair=pa/(Ra*Ta)
call dvisca(tair,visca)
do l = 1,nrbins
call fallv(raind(l),tair,pair,rhoair,visca,vt(l))
enddo
do d=1,nbins ! d = dust bins // dust classes
if ( chem(i,k,j,bins(d)).gt.0.) then ! dust exists
call w_t(wt,dbinmm(d),den_dust(d),rhoair)
scrate=0.
scavn=0.
do l = 1,nrbins
! calculate collection efficiency for each dust bin (already in this loop) by each rain bin (already in this loop)
call coleff(dbinmm(d),raind(l),den_dust(d),pair,tair,wt,vt(l),rhoair,visca,viscw,colece)
! calculate scavanging rate
delR = (rend(l+1)-rend(l))*1.e2 ! cm
delv = (vt(l)-wt)*1.e2 ! cm/s
delv = amax1(delv,0.)
rnflx = dsd_rn(l)*delv
carea = pi*0.25*(raind(l)*1.e2+dbinmm(d)*1.e2)**2.
scavn = scavn+carea*colece*rnflx*delR
enddo
scrate = scavn
dustold = chem(i,k,j,bins(d))
chem(i,k,j,bins(d))=chem(i,k,j,bins(d))-max(0.,chem(i,k,j,bins(d))*scrate*dtstepc)
chem(i,k,j,bins(d))=max(chem(i,k,j,bins(d)),epsilc)
wdbins(i,k,j,d)=max(epsilc,dustold-chem(i,k,j,bins(d)))
endif
enddo
! mklose: fifths size bin is used now, hence the following is not needed:
! dustold = chem(i,k,j,bins(5))
! chem(i,k,j,bins(5))=chem(i,k,j,bins(1))+chem(i,k,j,bins(2))+chem(i,k,j,bins(3))+chem(i,k,j,bins(4))
! chem(i,k,j,bins(5))=max(chem(i,k,j,bins(5)),epsilc)
! wdbins(i,k,j,5)=max(epsilc,dustold-chem(i,k,j,bins(5)))
endif
enddo
enddo
enddo
do j=jts,jte
do i=its,ite
do k=kts,kte
dustwdload_1(i,j)= max(epsilc,dustwdload_1(i,j) + wdbins(i,k,j,1)/alt(i,k,j) * dz8w(i,k,j))
dustwdload_2(i,j)= max(epsilc,dustwdload_2(i,j) + wdbins(i,k,j,2)/alt(i,k,j) * dz8w(i,k,j))
dustwdload_3(i,j)= max(epsilc,dustwdload_3(i,j) + wdbins(i,k,j,3)/alt(i,k,j) * dz8w(i,k,j))
dustwdload_4(i,j)= max(epsilc,dustwdload_4(i,j) + wdbins(i,k,j,4)/alt(i,k,j) * dz8w(i,k,j))
dustwdload_5(i,j)= max(epsilc,dustwdload_5(i,j) + wdbins(i,k,j,5)/alt(i,k,j) * dz8w(i,k,j))
dustwd_1(i,k,j)=max(epsilc,wdbins(i,k,j,1))
dustwd_2(i,k,j)=max(epsilc,wdbins(i,k,j,2))
dustwd_3(i,k,j)=max(epsilc,wdbins(i,k,j,3))
dustwd_4(i,k,j)=max(epsilc,wdbins(i,k,j,4))
dustwd_5(i,k,j)=max(epsilc,wdbins(i,k,j,5))
enddo
wetdep_1(i,j)= max(epsilc,wdbins(i,kts,j,1)/alt(i,kts,j) * dz8w(i,kts,j) / dtstepc)
wetdep_2(i,j)= max(epsilc,wdbins(i,kts,j,2)/alt(i,kts,j) * dz8w(i,kts,j) / dtstepc)
wetdep_3(i,j)= max(epsilc,wdbins(i,kts,j,3)/alt(i,kts,j) * dz8w(i,kts,j) / dtstepc)
wetdep_4(i,j)= max(epsilc,wdbins(i,kts,j,4)/alt(i,kts,j) * dz8w(i,kts,j) / dtstepc)
wetdep_5(i,j)= max(epsilc,wdbins(i,kts,j,5)/alt(i,kts,j) * dz8w(i,kts,j) / dtstepc)
enddo
enddo
end subroutine uoc_dustwd_driver
!=======================================================================
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
! - further updates M. Klose
subroutine coleff(d,rain,prho,p,t,pv,rv,arho,visca,viscw,ecolec)
! the program to calculate collection efficiency
! it is based on Slinn's empirical equation(1983)
! and theoretical collision efficiencies by
! Mason(1971), Klett and Davis (1973) and Beard and Ochs (1984)
! input
! =====
real, intent (in) :: d ! particle diameter in m
real, intent (in) :: rain ! raindrio diameter in m
real(8), intent (in) :: prho ! particle density in kg/m3
real, intent (in) :: p ! atmosphere pressure in Pa
real, intent (in) :: t ! atmospheric temperature in K
real, intent (in) :: pv ! dust fall velocity in m/s
real, intent (in) :: rv ! rain fall velocity in m/s
real, intent (in) :: arho ! rair density in kg/m3
real, intent (in) :: visca, viscw
! output
! ======
real, intent (out) :: ecolec
! local variables
! ===============
real :: beta
real :: rmin
real(8) :: E_im,E_in,E_br
real :: dd, raind, rhop, pr, ta, vp, vs, rhoa
CHARACTER*(100) :: text
integer :: debug_level
data rhow/1/ ! raindrop density g/cm^3
CALL get_wrf_debug_level(debug_level)
dd = d*1.e6 ! m to um
raind = rain*1.e6 ! m to um
rhop = prho*1.e-3 ! kg/m3 to g/cm3
pr = p*1.e-2 ! Pa to mb or hPa
ta = t-273.15 ! K to C
vp = pv*1.e2 ! m/s to cm/s
vs = rv*1.e2 ! m/s to cm/s
rhoa = arho*1.e-3 ! kg/m3 to g/cm3
!=======================================================================
ecol=0
if(dd.le.2) then
call eslinn(raind,dd,pr,ta,rhop,vp,vs,rhoa,visca,viscw,ecol)
elseif(dd.gt.2 .and. dd.le.40) then
call eimpact(rhop,dd,raind,ecol)
elseif(dd.gt.40) then
ecol=0
endif
ecolec=ecol
! return
end subroutine coleff
!****************************************************************
!****************************************************************
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
! - further updates M. Klose
! xlambda is set fix instead of calculated - cfrick
subroutine eslinn(dr,dp,pair,tair,rhop,vp,vs,rhoa,visca,viscw,eslin)
! input
! =====
real, intent (in) :: dp ! particle diameter in um
real, intent (in) :: dr ! raindrio diameter in um
real, intent (in) :: rhop ! particle density in g/cm3
real, intent (in) :: pair ! atmosphere pressure in hPa
real, intent (in) :: tair ! atmospheric temperature in C
real, intent (in) :: vp ! dust fall velocity in cm/s
real, intent (in) :: vs ! rain fall velocity in cm/s
real, intent (in) :: rhoa ! rair density in g/cm3
real, intent (in) :: visca, viscw
! output
! ======
real, intent (out) :: eslin
real :: pi
CHARACTER*(100) :: text
integer :: debug_level
real :: gnarf, ren
E_br=0
E_in=0
E_im=0
drc=dr*1.e-4 ! rain diameter in cm
dpc=dp*1.e-4 ! particle diameter in cm
pr=pair
ta=tair
pi=acos(-1.0)
xk=1.381e-16 ! Boltzmann constant (g cm2/s2)/K
! xd=3.75e-8
! prs=pr*1.e+3
! xlambda=xk*(ta+273.15)/(pi*sqrt(2.0)*prs*xd*xd)
omg=viscw/visca
ren=0.5*drc*vs*rhoa/visca ! reynolds number
sstar=(1.2+alog(1+ren)/12)/(1+alog(1+ren)) ! S*
xlambda=0.0651*1.e-4 ! mean free path of air in cm
cc=1+(2*xlambda/dpc)*(1.257+0.4*exp(-1.1*dpc/(2*xlambda))) ! Cunningham's correction for small particles
tau=dpc**2*rhop*cc/(18*visca)
stn=2*tau*(vs-vp)/drc ! Stokes number
diff=(xk*(ta+273.15)*cc)/(3*pi*visca*dpc)
scn=visca/(diff*rhoa) ! Schmidt number
phi=dpc/drc
pen=ren*scn ! peclet number
! Brownian diffusion
E_br=4*(1+0.4*sqrt(ren)*scn**(1.0/3.0)+0.16*sqrt(ren)*sqrt(scn))/pen
! Interception
E_in=4*phi*(1/omg+(1+2*sqrt(ren))*phi)
! Impaction
if(stn.ge.sstar) then
E_im=((stn-sstar)/(stn-sstar+2.0/3.0))**1.5
endif
eslin=E_br+E_in
! return
end subroutine eslinn
!=======================================================================
!*******************************************************************************
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
! - further updates M. Klose
subroutine eimpact(rhop,dp,raind,e_im)
USE module_data_uoc_wd ! mklose [17122014]: renamed from we_ematrix
integer, parameter :: lp=17, lr=17
! include 'we_ematrix.dat'
! input
! =====
real, intent (in) :: dp ! particle diameter in um
real, intent (in) :: raind ! raindrio diameter in um
real, intent (in) :: rhop ! particle density in g/cm3
! output
! ======
real, intent (out) :: e_im
real :: ee(lr)
real :: u(lr+1)
real :: enew(lr+1)
real :: e_r(lp)
real :: e_c(lp)
real :: f(lp)
real :: v(lp+1)
real :: dnew(lp+1)
CHARACTER*(100) :: text
integer :: debug_level
CALL get_wrf_debug_level(debug_level)
! rr, rp : radiuse of raindrop, dust particle in um
! edata : dataset of composite collection efficiency
!===============================================================================
if(raind.lt.100 .or. raind.gt.6000) then
text = 'raindrop diameter out of range'
text=trim(trim(text)//" - ERROR - UoC dust wet deposition")
call wrf_debug (debug_level,text)
endif
if(dp.lt.2 .or. dp.gt.40) then
text = 'particle diameter out of range'
text=trim(trim(text)//" - ERROR - UoC dust wet deposition")
call wrf_debug (debug_level,text)
endif
! radius of raindrop, dust particle
rr=raind*0.5
rp=dp*0.5
do i=1,lp
! calculate e_r(1:lp) corresponding to dustr(1:lp) for rr
jj=0
do j=1,lr
ee(j)=edatawd(i,j)
if(rr.ge.rainwd(j)) jj=j
enddo
if(jj.eq.0) then
text = 'error in raindrop radius'
text=trim(trim(text)//" - ERROR - UoC dust wet deposition")
call wrf_debug (debug_level,text)
endif
if(rr.gt.rainwd(jj)) then
lrp1=lr+1
do l=1,jj
u(l)=rainwd(l)
enddo
u(jj+1)=rr
do l=jj+1,lr
u(l+1)=rainwd(l)
enddo
call intrpl(lr,log(rainwd),ee,lrp1,log(u),enew)
e_r(i)=enew(jj+1)
elseif(rr.eq.rainwd(jj)) then
e_r(i)=ee(jj)
endif
enddo
! interpolate e_d for rp from e_r
rmax=sqrt(1/rhop)
do i=1,lp
if(dustwd(i).le.1) then
f(i)=1
elseif(dustwd(i).gt.1 .and. dustwd(i).lt.3.98) then
f(i)=(1-rmax)*(dustwd(i)-3.98)**2/(1-3.98)**2+rmax
elseif(dustwd(i).ge.3.98) then
f(i)=rmax
endif
enddo
do i=1,lp
e_c(i)=e_r(i)*f(i)
enddo
ii=0
e_d=0
do i=1,lp
if(rp.ge.dustwd(i)) ii=i
enddo
if(ii.eq.0) then
text = 'wrong particle radius'
text=trim(trim(text)//" - ERROR - UoC dust wet deposition")
call wrf_debug (debug_level,text)
endif
if(rp.gt.dustwd(ii)) then
do l=1,ii
v(l)=dustwd(l)
enddo
v(ii+1)=rp
do l=ii+1,lp
v(l+1)=dustwd(l)
enddo
lpp1=lp+1
call intrpl(lp,dustwd,e_c,lpp1,v,dnew)
e_d=dnew(ii+1)
elseif(rp.eq.dustwd(ii)) then
e_d=e_c(ii)
endif
e_im=e_d
! return
end subroutine eimpact
!=======================================================================
!***********************************************************************
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
! - further updates M. Klose
subroutine dvisca(tair,dvisc)
! tair : air temperature in K
! visca : dynamic viscosity of air in g/cm/s
integer, parameter :: l=9, n=1
real, dimension(l) :: ta, visca
real, dimension(l) :: x, y
real, dimension(n) :: u, v
data ta/-173,-73,0,20,25,27,127,227,327/
data visca/0.71e-4,1.33e-4,1.72e-4,1.797e-4,1.818e-4, &
& 1.86e-4,2.31e-4,2.71e-4,3.08e-4/
! values from crc handbook of chemistry and physics
do i = 1,l
x(i)=ta(i)+273
y(i)=visca(i)
enddo
u(1)=tair
call intrpl(l,x,y,n,u,v)
dvisc=v(1)
! return
end subroutine dvisca
!
!************************************************************
!************************************************************
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
! - further updates M. Klose
subroutine dviscw(twt,dvisc)
integer, parameter :: l=11, n=1
real, dimension(l) :: tw, viscw
real, dimension(l) :: x, y
real, dimension(l) :: u, v
data tw/0,10,20,30,40,50,60,70,80,90,100/
data viscw/1.7930e-2,1.3070e-2,1.002e-2,0.7977e-2,0.6532e-2, &
& 0.5470e-2,0.4665e-2,0.404e-2,0.3544e-2,0.3145e-2, &
& 0.2818e-2/
! tw : water temperature in degree C
! dynamic viscosity of water in g/cm/s
! values from crc handbook of chemistry and physics
! properties of water in the range 0 - 100c
do i = 1,l
x(i)=tw(i)+273.15
y(i)=viscw(i)
enddo
u(1)=twt+273.15
call intrpl(l,x,y,n,u,v)
dvisc=v(1)
! return
end subroutine dviscw
!*******************************************************************************
!************************************************************
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
! - further updates by M. Klose, 2015
subroutine intrpl(l,x,y,n,u,v)
!************************************************************
! algorithm appeared in comm. acm, vol. 15, no. 10,
! p. 914.
!
! interpolation of a single-valued function
! this subroutine interpolates, from values of the function
! given as ordinates of input data points in an x-y plane
! and for a given set of x values (abscissas), the values of
! a single-valued function y = y(x).
! the input parameters are
! l = number of input data points
! (must be 2 or greater)
! x = array of dimension l storing the x values
! (abscissas) of input data points
! (in ascending order)
! y = array of dimension l storing the y values
! (ordinates) of input data points
! n = number of points at which interpolation of the
! y value (ordinate) is desired
! (must be 1 or greater)
! u = array of dimension n storing the x values
! (abscissas) of desired points
! the output parameter is
! v = array of dimension n where the interpolated y
! values (ordinates) are to be displayed
! declaration statements
integer, intent(in) :: l, n
real, dimension(l) :: x, y
real, dimension(n) :: u, v
equivalence (p0,x3),(q0,y3),(q1,t3)
real :: m1,m2,m3,m4,m5
equivalence (uk,dx),(imn,x2,a1,m1),(imx,x5,a5,m5), &
& (j,sw,sa),(y2,w2,w4,q2),(y5,w3,q3)
real :: a1, a2, a3, a4, a5
integer :: debug_level
CALL get_wrf_debug_level(debug_level)
! preliminary processing
10 l0=l
lm1=l0-1
lm2=lm1-1
lp1=l0+1
n0=n
if(lm2.lt.0) go to 90
if(n0.le.0) go to 91
do 11 i=2,l0
if(x(i-1)-x(i)) 11,95,96
11 continue
ipv=0
! main do-loop
do 80 k=1,n0
uk=u(k)
! routine to locate the desired point
20 if(lm2.eq.0) go to 27
if(uk.ge.x(l0)) go to 26
if(uk.lt.x(1)) go to 25
imn=2
imx=l0
21 i=(imn+imx)/2
if(uk.ge.x(i)) go to 23
22 imx=i
go to 24
23 imn=i+1
24 if(imx.gt.imn) go to 21
i=imx
go to 30
25 i=1
go to 30
26 i=lp1
go to 30
27 i=2
! check if i = ipv
30 if(i.eq.ipv) go to 70
ipv=i
! routines to pick up necessary x and y values and
! to estimate them if necessary
40 j=i
if(j.eq.1) j=2
if(j.eq.lp1) j=l0
x3=x(j-1)
y3=y(j-1)
x4=x(j)
y4=y(j)
a3=x4-x3
m3=(y4-y3)/a3
if(lm2.eq.0) go to 43
if(j.eq.2) go to 41
x2=x(j-2)
y2=y(j-2)
a2=x3-x2
m2=(y3-y2)/a2
if(j.eq.l0) go to 42
41 x5=x(j+1)
y5=y(j+1)
a4=x5-x4
m4=(y5-y4)/a4
if(j.eq.2) m2=m3+m3-m4
go to 45
42 m4=m3+m3-m2
go to 45
43 m2=m3
m4=m3
45 if(j.le.3) go to 46
a1=x2-x(j-3)
m1=(y2-y(j-3))/a1
go to 47
46 m1=m2+m2-m3
47 if(j.ge.lm1) go to 48
a5=x(j+2)-x5
m5=(y(j+2)-y5)/a5
go to 50
48 m5=m4+m4-m3
! numerical differentiation
50 if(i.eq.lp1) go to 52
w2=abs(m4-m3)
w3=abs(m2-m1)
sw=w2+w3
if(sw.ne.0.0) go to 51
w2=0.5
w3=0.5
sw=1.0
51 t3=(w2*m2+w3*m3)/sw
if(i.eq.1) go to 54
52 w3=abs(m5-m4)
w4=abs(m3-m2)
sw=w3+w4
if(sw.ne.0.0) go to 53
w3=0.5
w4=0.5
sw=1.0
53 t4=(w3*m3+w4*m4)/sw
if(i.ne.lp1) go to 60
t3=t4
sa=a2+a3
t4=0.5*(m4+m5-a2*(a2-a3)*(m2-m3)/(sa*sa))
x3=x4
y3=y4
a3=a2
m3=m4
go to 60
54 t4=t3
sa=a3+a4
t3=0.5*(m1+m2-a4*(a3-a4)*(m3-m4)/(sa*sa))
x3=x3-a4
y3=y3-m2*a4
a3=a4
m3=m2
! determination of the coefficients
60 q2=(2.0*(m3-t3)+m3-t4)/a3
q3=(-m3-m3+t3+t4)/(a3*a3)
! computation of the polynomial
70 dx=uk-p0
80 v(k)=q0+dx*(q1+dx*(q2+dx*q3))
return
! error exit
90 call wrf_debug (debug_level,'error 90 in subroutine intrpl - ERROR - UoC dust wet deposition') !write (6,2090)
go to 99
91 call wrf_debug (debug_level,'error 91 in subroutine intrpl - ERROR - UoC dust wet deposition') !write (6,2091)
go to 99
95 call wrf_debug (debug_level,'error 95 in subroutine intrpl - ERROR - UoC dust wet deposition') !write (6,2095)
go to 97
96 call wrf_debug (debug_level,'error 96 in subroutine intrpl - ERROR - UoC dust wet deposition') !write (6,2096)
97 call wrf_debug (debug_level,'error 97 in subroutine intrpl - ERROR - UoC dust wet depositionk') !write (6,2097) i,x(i)
99 call wrf_debug (debug_level,'error 98 in subroutine intrpl - ERROR - UoC dust wet deposition') !write (6,2099) l0,n0
return
! format statements
! 2090 format(1x/22h *** l = 1 or less./)
! 2091 format(1x/22h *** n = 0 or less./)
! 2095 format(1x/27h *** identical x values./)
! 2096 format(1x/33h *** x values out of sequence./)
! 2097 format(6h i =,i7,10x,6hx(i) =,e12.3)
! 2099 format(6h l =,i7,10x,3hn =,i7/&
! & 36h error detected in routine intrpl)
end subroutine intrpl
!************************************************************
!=======================================================================
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
! - further updates M. Klose
subroutine fallv(d,t,p,rhoair,visc,vt)
!=======================================================================
! the program to calculate the termincal velocity of raindrop
! this is based on Beard (1976)
!
! input
! ta : temperature in K
! dropd : raindrop diameter in m
!
! output
! vt : terminal velocity in m/s
!=======================================================================
real :: gg, pr0, t0, visc0, rhow, rhoa, dp, dl, cl, csc, vt, ta, dropd, pr
real :: a0, a1, a2, a3, a4, a5, a6
real :: b0, b1, b2, b3, b4, b5
real :: c1, c2, dan, x, vy, rey, c3, bon, phn, rhoair, p, t, d
CHARACTER*(100) :: text
integer :: debug_level
CALL get_wrf_debug_level(debug_level)
! compute surface tension
! sig : the surface tension in mN/m
dropd = d*1.e6 ! m to um
ta = t-273.15 ! K to C
pr = p*1.e-2 ! Pa to hPa (mb)
call sfctens(ta,sig)
gg=980 ! cm/s^2
pr0=1013.25 ! mb
t0=293.15 ! K
visc0=1.818e-4 ! g/cm/sec
rhow=1 ! g/cm^3
rhoa = rhoair*1.e-3 ! kg/m3 to g/cm3
! compute the terminal velocity
vt=0
dp=dropd
dl=dropd*1.e-4 ! diameter in cm
c1=(rhow-rhoa)*gg/(18*visc)
cl=6.62e-6*(visc/visc0)*(pr0/pr)*sqrt((ta+273.15)/t0) !cm
csc=1+2.51*cl/dl
if(dp.lt.19.) then
vt=c1*csc*dl**2
vt=vt*1.e-2 ! cm/s to m/s
elseif(dp.ge.19. .and. dp.lt.1070.) then
a0=-0.318657e+1
a1=+0.992696
a2=-0.153193e-2
a3=-0.987059e-3
a4=-0.578878e-3
a5=+0.855176e-4
a6=-0.327815e-5
c2=4*rhoa*(rhow-rhoa)*gg/(3*visc*visc)
dan=c2*dl**3
x=alog(dan)
vy=a0+a1*x+a2*x**2+a3*x**3+a4*x**4+a5*x**5+a6*x**6
rey=csc*exp(vy)
vt=visc*rey/(rhoa*dl)
vt=vt*1.e-2 ! cm/s to m/s
elseif(dp.ge.1070. .and. dp.le.7000.) then
b0=-0.500015e+1
b1=+0.523778e+1
b2=-0.204914e+1
b3=+0.475294
b4=-0.542819e-1
b5=+0.238449e-2
c3=4*(rhow-rhoa)*gg/(3*sig)
bon=c3*dl*dl ! Bond Number
phn=sig**3*rhoa**2/(visc**4*gg*(rhow-rhoa)) ! Physical Property Number
y=alog(bon*phn**(1.0/6.0))
vy=b0+b1*y+b2*y**2+b3*y**3+b4*y**4+b5*y**5
rey=phn**(1.0/6.0)*exp(vy)
vt=visc*rey/(rhoa*dl)
vt=vt*1.e-2 ! cm/s to m/s
endif
! return
end subroutine fallv
!=======================================================================
!=======================================================================
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
! - further updates, M. Klose
subroutine sfctens(temp,sig)
! temp : temperature in K
! sig : surface tension in mN/m or dyne/cm
integer, parameter :: n=14
real, dimension(n) :: t, sfctn
data t/0,5,10,15,18,20,25,30,40,50,60,70,80,100/
data sfctn/75.6,74.9,74.22,73.49,73.05,72.75,71.97,71.18,69.56, &
& 67.91,66.18,64.4,62.6,58.9/
real :: ta
kk = 1
xsfc = 0
ta=max(temp-273.15,0.0)
do i = 1, 14
if(i.lt.14) then
if(ta.ge.t(i) .and. ta.lt.t(i+1)) then
kk=i
goto 1
endif
else
if(ta.ge.t(i)) then
kk=14
goto 1
endif
endif
enddo
1 continue
if(kk.lt.14) then
sig = sfctn(kk)+ &
& (sfctn(kk+1)-sfctn(kk))*(ta-t(kk))/(t(kk+1)-t(kk))
else
sig = sfctn(kk)+ &
& (sfctn(kk)-sfctn(kk-1))*(ta-t(kk))/(t(kk)-t(kk-1))
endif
! return
end subroutine sfctens
!***********************************************************************
!***********************************************************************
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
! - further updates M. Klose
subroutine dsd_rain(rnrate,raind,nrbin,dsd_rn,ldsd)
! input
real, intent(in) :: raind(nrbin), rnrate ! in m and mm/h
integer, intent(in) :: ldsd, nrbin
! local
real :: dsd_mp(nrbin)
real :: dsd_ss(nrbin)
real :: dsd_wt(nrbin)
real :: dsd_fl(nrbin)
real :: ng
real :: nt
real :: pi, dr, wc, d0, dcm, const, dumm
! output
real, intent(out) :: dsd_rn(nrbin)
!.......................................................................
d0=0
do nk=1,nrbin
dr=raind(nk)*1.e3 ! m to mm
if(ldsd.eq.1) then
!
! Marshall-Palmer(1948)
dsd_mp(nk)=8.e3*exp(-4.1*dr*rnrate**(-0.21)) ! 1/m^3/mm
dsd_mp(nk)=dsd_mp(nk)*1.e-5 ! 1/cm^3/cm
dsd_rn(nk)=dsd_mp(nk)
!
elseif(ldsd.eq.2) then
! Sekhon and Srivastava
dsd_ss(nk)=0.07*rnrate**0.37*exp(-3.8*dr/rnrate**0.14) ! cm^(-3) cm^(-1)
dsd_rn(nk)=dsd_ss(nk)
!
elseif(ldsd.eq.3) then
! Willis-Tattelman(1989)
wc=0.062*rnrate**0.913 ! water content in g/m^3
d0=0.1571*wc**0.1681 ! median volume diameter in cm
ng=512.85*wc*1.e-6/d0**4*d0**(-2.16)
dcm=dr*0.1 ! diameter in cm
dsd_wt(nk)=ng*dcm**2.16*exp(-5.5880*dcm/d0) ! #/cm^3/cm
dsd_rn(nk)=dsd_wt(nk)
elseif(ldsd.eq.4) then
! Feingold-Levin(1986)
pi=acos(-1.0)
const=sqrt(2*pi)*alog(1.43)
nt=172*rnrate**0.22 ! 1/m^3
dg=0.72*rnrate**0.21
dumm=0.5*alog(dr/dg)**2/alog(1.43)**2
dsd_fl(nk)=nt/const/dr*exp(-dumm) ! 1/m^3/mm
dsd_fl(nk)=dsd_fl(nk)*1.e-5 ! 1/cm^3/cm
dsd_rn(nk)=dsd_fl(nk)
endif
!
enddo
return
end subroutine dsd_rain
!***********************************************************************
!***********************************************************************
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
! - further updates M. Klose
subroutine w_p(d,rho,wt)
! Yaping Shao 28-07-92
! Calculate Terminal velocity for spherical grain
! wt : terminal velocity [m/s]
! d : particle diameter [m]
! rho : particle density [kg/m3]
! Require CDSPH
! Parameter specifications
IMPLICIT NONE ! cfrick
REAL, INTENT(OUT) :: wt
REAL :: lambda
REAL :: rden, err, test, x, xl, xh
REAL :: f, fm, fl, fh
INTEGER :: i, isign
real, intent(in) :: rho, d
real :: cc
REAL, PARAMETER :: visc=1.5E-05
INTEGER :: debug_level
CHARACTER*(100) :: text
! VARIABLE X IS WT, FUNCTION F IS: CD(RE)*WT**2 - 4*D*G*RDEN/3
CALL get_wrf_debug_level(debug_level)
rden=rho ! particle density
lambda=0.0651*1.e-6 ! mean free path of air
cc=1+(2*lambda/d)*(1.257+0.4*exp(-1.1*d/(2*lambda))) ! Cunningham's correction for small particles
err=1E-6
xl=1E-12
xh=1E3
wt = xl
fl = CDSPH(wt*d/visc)*wt**2/cc - 4.0*d*g*rden/3.0
wt = xh
fh = CDSPH(wt*d/visc)*wt**2/cc - 4.0*d*g*rden/3.0
IF (fh*fl.GT.0) call wrf_debug (debug_level,'w_t: F(XL), F(XH) HAVE SAME SIGN - ERROR - UoC dust wet deposition')
IF (fh-fl.EQ.0) call wrf_debug (debug_level,'w_t: F(XL) = F(XH) - ERROR - UoC dust wet deposition')
isign=1
IF (fh.LT.fl) isign=-1
fh=isign*fh
fl=isign*fl
DO 1 i=1,100
x=(xl+xh)/2
wt = x
f = CDSPH(wt*d/visc)*wt**2/cc - 4.0*d*g*rden/3.0
fm=isign*f
IF (fm.GT.fh.OR.fm.LT.fl) call wrf_debug (debug_level,'w_t: F(X) NON-MONOTONIC - ERROR - UoC dust wet deposition')
IF (fm.GT.0) xh=x
IF (fm.LT.0) xl=x
IF (fm.EQ.0) then
return
ENDIF
test=ABS(xh-xl)/ABS(x)
IF (test.LT.ABS(err)) then
return
ENDIF
1 CONTINUE
call wrf_debug (debug_level,'w_t: NO SOLUTION IN 100 LOOPS - ERROR - UoC dust wet deposition')
END subroutine w_p
!***************************************************************
!***********************************************************************
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
! - further updates M. Klose
SUBROUTINE w_t (wt, dmm, rhop, rhoa)
! Yaping Shao 28-07-92
! Calculate Terminal velocity for spherical grain
! wt : terminal velocity [m/s]
! dmm: particle diameter in [m]
! rhop: particle density [kg/m3]
! rhoa: air density [kg/m3]
! Require CDSPH
! Parameter specifications
! include Cunningham's correction - Claudia Frick 24-04-2014
REAL(8), INTENT(IN) :: rhop
REAL, INTENT(IN) :: rhoa
REAL :: wt, dmm
REAL :: RDEN, D, ERR, TEST, X, XL, XH
REAL :: F, FM, FL, FH
REAL, PARAMETER :: VISC=1.5E-05
INTEGER :: I, ISIGN
CHARACTER*(100) :: text
REAL :: lambda
real :: cc
! VARIABLE X IS WT, FUNCTION F IS: CD(RE)*WT**2 - 4*D*G*RDEN/3
D = dmm
RDEN=rhop/rhoa ! = sigma
lambda=0.0651*1.e-6 ! mean free path of air
cc=1+(2*lambda/D)*(1.257+0.4*exp(-1.1*D/(2*lambda))) ! Cunningham's correction for small particles
ERR=1E-6
XL=1E-12
XH=1E3
WT = XL
FL = CDSPH(WT*D/VISC)*WT**2/cc - 4.0*D*G*RDEN/3.0
WT = XH
FH = CDSPH(WT*D/VISC)*WT**2/cc - 4.0*D*G*RDEN/3.0
IF (FH*FL.GT.0) call wrf_debug (debug_level,'w_t: F(XL), F(XH) HAVE SAME SIGN - ERROR - UoC dust wet deposition')
IF (FH-FL.EQ.0) call wrf_debug (debug_level,'w_t: F(XL) = F(XH) - ERROR - UoC dust wet deposition')
ISIGN=1
IF (FH.LT.FL) ISIGN=-1
FH=ISIGN*FH
FL=ISIGN*FL
DO 1 I=1,100
X=(XL+XH)/2
WT = X
F = CDSPH(WT*D/VISC)*WT**2/cc - 4.0*D*G*RDEN/3.0
FM=ISIGN*F
IF (FM.GT.FH.OR.FM.LT.FL) call wrf_debug (debug_level,'w_t: F(X) NON-MONOTONIC - ERROR - UoC dust wet deposition')
IF (FM.GT.0) XH=X
IF (FM.LT.0) XL=X
IF (FM.EQ.0) then
return
ENDIF
TEST=ABS(XH-XL)/ABS(X)
IF (TEST.LT.ABS(ERR)) then
return
ENDIF
1 CONTINUE
call wrf_debug (debug_level,'w_t: NO SOLUTION IN 100 LOOPS - ERROR - UoC dust wet deposition')
END subroutine w_t
!***************************************************************
!---------------------------------------------------------------
! CEMSYS5 - smaller modifications by cfrick for WRF implementation
real FUNCTION CDSPH(RE)
! MRR, 30-SEP-87
! ADAPTED 10-OCT-89 (NO ISTOKE, NO CRASH FOR RE>50000)
! CALCULATES DRAG COEFFICIENT CD FOR A SPHERE, AS A FUNCTION OF
! REYNOLDS NUMBER RE, WHERE:
! DRAG FORCE = (RHO * PI * D**2 * U**2) / 8
! RE = U * D / VISCOSITY
! (D = SPHERE DIAMETER, U = VELOCITY, RHO = FLUID DENSITY)
! ALGORITHM FROM MORSI AND ALEXANDER (1972, JFM 55, 193-208), CHECKED
! AGAINST THEIR TABLE FOR CD(SPHERE).
! FOR RE > 50000, SET CDSPH = 0.48802, THE VALUE AT RE=50000=5*10**4.
! THIS IS OK TO RE=3*10**5 (APPROX), BEYOND WHICH DRAG CRISIS OCCURS
! AND SETS CDSPH TO ABOUT 0.1.
IMPLICIT NONE ! cfrick
real :: RE
IF (RE.LE.0) THEN
WRITE(6,*) RE
STOP 'CDSPH: RE.LE.0'
ELSE IF (RE.LE.0.1) THEN
CDSPH = 24/RE
ELSE IF (RE.LE.1) THEN
CDSPH = 22.73/RE + 0.0903/RE**2 + 3.69
ELSE IF (RE.LE.10) THEN
CDSPH = 29.1667/RE - 3.8889/RE**2 + 1.222
ELSE IF (RE.LE.100) THEN
CDSPH = 46.5/RE - 116.67/RE**2 + 0.6167
ELSE IF (RE.LE.1000) THEN
CDSPH = 98.33/RE - 2778/RE**2 + 0.3644
ELSE IF (RE.LE.5000) THEN
CDSPH = 148.62/RE - 47500/RE**2 + 0.357
ELSE IF (RE.LE.10000) THEN
CDSPH = -490.546/RE + 578700/RE**2 + 0.46
ELSE IF (RE.LE.50000) THEN
CDSPH = -1662.5/RE + 5416700/RE**2 + 0.5191
ELSE
CDSPH = 0.48802
END IF
END FUNCTION CDSPH
!--------------------------------------------------------------------
END MODULE module_uoc_dustwd