chem/module_uoc_dustwd.F

   1 MODULE module_uoc_dustwd
   2
   3
   4 !----------------------------------------------------------------------------
   5 ! Dust wet deposition module of the University of Cologne, Germany.
   6 ! Dust wet deposition scheme developed by E Jung (2004, PhD Thesis).
   7 ! Implementation into WRF and modifications by C Frick, 2014 (claudia.frick@uni-koeln.de)
   8 ! 2014-2015, updates and modifications, Martina Klose (mklose@uni-koeln.de)
   9 !----------------------------------------------------------------------------
  10
  11 USE module_state_description   ! num_chem, p_qr, ...
  12 USE module_model_constants     ! rhowater in kg/m3
  13 USE physconst                  ! rair, ...
  14 USE module_data_gocart_dust
  15
  16  CONTAINS
  17
  18 subroutine uoc_dustwd_driver(precr,chem,p_phy,t_phy,                       &
  19                              ids,ide, jds,jde, kds,kde,                    &
  20                              ims,ime, jms,jme, kms,kme,                    &
  21                              its,ite, jts,jte, kts,kte,                    &
  22                              dtstepc,                                      &
  23                              dustwd_1, dustwd_2,                           &
  24                              dustwd_3, dustwd_4,                           &
  25                              dustwd_5,                                     &
  26                              wetdep_1, wetdep_2,                           &
  27                              wetdep_3, wetdep_4,                           &
  28                              wetdep_5,                                     &
  29                              dustwdload_1, dustwdload_2,                   &
  30                              dustwdload_3, dustwdload_4,                   &
  31                              dustwdload_5,                                 &
  32                              alt, dz8w, epsilc                             )
  33
  34 IMPLICIT NONE
  35
  36  INTEGER :: debug_level
  37  CHARACTER*(100) :: text
  38
  39 real :: dustold                           ! help variable
  40
  41 integer :: d, i, j, k, l                  ! loop
  42
  43 integer, parameter :: nbins=5                 ! number of dust bins
  44 integer, parameter :: nbinsa=5                ! number of dust bins
  45 integer :: bins(nbinsa)                       ! dust bin numbers in chem
  46 ! real, dimension(nbins)    :: dbin              ! max. dimension of dust in the bin
  47 ! data dbin/2.5,5.,10.,20./                      ! size cut diameter (um) - from qf03
  48 real, dimension(nbins)    :: dbinm             ! mean dimension of dust (um) in the bin
  49 real, dimension(nbins)    :: dbinmm            ! mean dimension of dust (m) in the bin
  50 ! data dbinm/1.25,3.75,7.5,15./                  ! mean size cut diameter (um)
  51
  52 real :: conver9, converi9                 ! transformation ug/kg-dryair to kg/kg-dryair and vice versa
  53 real :: conver6, converi6                 ! transformation um to m and vice versa
  54
  55 real :: wt                                ! terminal fall velocity of dust in the bin (m/s)
  56
  57 real :: rmin, rmax                        ! minimum/maximum raindrop diameter (m)
  58
  59 integer, parameter        :: nrbins=30    ! number of raindrop bins
  60 real, dimension(nrbins)   :: raind        ! raindrop diameter in the middle of the raindrop bin (m)
  61 real, dimension(nrbins+1) :: rend         ! raindrop diameter at the end of the raindrop bin (m)
  62 real, dimension(nrbins)   :: dsd_rn       ! raindrop size distribution
  63 real, dimension(nrbins)   :: vt           ! raindrop terminal fall velocity (m/s)
  64 real                      :: z            ! help parameter for the determination of the rain bins
  65
  66 real :: visca                             ! (dynamic) viscosity of air (g/cm s)
  67 real :: tw, viscw                         ! temperature of water (C) and (dynamic) viscosity of water (g/cm s)
  68
  69 real :: colece                            ! collection efficiency
  70
  71 real :: scrate, scavn                     ! scavanging rate (1/s)
  72 real :: delR, delv, rnflx, carea          ! help parameter for the determination of the scavanging rate
  73 real :: tair, pair                        ! current air temperature and pressure
  74 real :: rhoair                            ! dry air density (kg/m3) - rhoair=pa/(Ra*Ta)
  75
  76 INTEGER,  INTENT(IN)    ::  ids,ide, jds,jde, kds,kde,    & ! domain grid
  77                             ims,ime, jms,jme, kms,kme,    &
  78                             its,ite, jts,jte, kts,kte
  79
  80 REAL, INTENT(IN)        :: dtstepc, & ! time step (s)
  81                            epsilc
  82
  83 REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),                        &
  84           INTENT(IN)    ::                                   precr
  85 ! precr - rain precipitation rate at all levels (kg/m2 s)
  86
  87 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),                &
  88           INTENT(INOUT) ::                                   chem
  89 ! chem(i,k,j,p_dust_?) - dust mixing ratio for bin np_dust_? (ug/kg-dryair)
  90
  91 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                           &
  92           INTENT(IN)    ::                                   p_phy, t_phy
  93 ! pressure (Pa) and temperature (K)
  94
  95 real, dimension( ims:ime, kms:kme, jms:jme ) ::              rnrate
  96 ! precipitation rate (mm/h) - rnrate=precr/3600 using rhowater (1l = 1kg)
  97
  98 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::                     &
  99 dustwd_1, dustwd_2, dustwd_3, dustwd_4, dustwd_5
 100 ! loss in dust mixing ratio due to wet deposition (ug/kg-dryair) (current time step)
 101
 102 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ) ::           wdbins
 103 ! loss in dust mixing ratio due to wet deposition - help variable (ug/kg-dryair)
 104
 105 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) ::                     &
 106 dustwdload_1, dustwdload_2, dustwdload_3, dustwdload_4, dustwdload_5
 107 ! loss in dustload due to wet deposition (ug/m2) (current time step)
 108
 109 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) ::                     &
 110 wetdep_1, wetdep_2, wetdep_3, wetdep_4, wetdep_5
 111 ! loss in dust mixing ratio due to wet deposition in one coulmn for all size bins (ug/m2/s)
 112
 113 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: alt, dz8w
 114 ! altitude = 1/rhoair (kg/m3) & vertical grid-spacing delta-z / dz between full levels (m)
 115
 116 den_dust(:) = 2650.
 117
 118 dbinm(:) = 2.0d0*reff_dust(:)
 119
 120 CALL get_wrf_debug_level(debug_level)
 121
 122 rmin =  100.*1.e-6
 123 rmax = 5800.*1.e-6              ! minimum/maximum raindrop diameter (m)
 124
 125  conver6  = 1.e-6               ! transformation um to m
 126  converi6 = 1.e6                ! transformation m to um
 127
 128  conver9  = 1.e-9               ! transformation ug/kg-dryair to kg/kg-dryair
 129  converi9 = 1.e9                ! transformation kg/kg-dryair to ug/kg-dryair
 130
 131 bins(1)=p_dust_1                ! dust bin number in chem
 132 bins(2)=p_dust_2                ! dust bin number in chem
 133 bins(3)=p_dust_3                ! dust bin number in chem
 134 bins(4)=p_dust_4                ! dust bin number in chem
 135 bins(5)=p_dust_5                ! dust bin number in chem
 136
 137 dustwd_1(its:ite,kts:kte,jts:jte) = 0.
 138 dustwd_2(its:ite,kts:kte,jts:jte) = 0.
 139 dustwd_3(its:ite,kts:kte,jts:jte) = 0.
 140 dustwd_4(its:ite,kts:kte,jts:jte) = 0.
 141 dustwd_5(its:ite,kts:kte,jts:jte) = 0.
 142 wetdep_1(its:ite,jts:jte) = 0.
 143 wetdep_2(its:ite,jts:jte) = 0.
 144 wetdep_3(its:ite,jts:jte) = 0.
 145 wetdep_4(its:ite,jts:jte) = 0.
 146 wetdep_5(its:ite,jts:jte) = 0.
 147 dustwdload_1(its:ite,jts:jte) = 0.
 148 dustwdload_2(its:ite,jts:jte) = 0.
 149 dustwdload_3(its:ite,jts:jte) = 0.
 150 dustwdload_4(its:ite,jts:jte) = 0.
 151 dustwdload_5(its:ite,jts:jte) = 0.
 152
 153 wdbins(its:ite,kts:kte,jts:jte,1)=dustwd_1(its:ite,kts:kte,jts:jte)             ! dust dust loss in dust bin 1
 154 wdbins(its:ite,kts:kte,jts:jte,2)=dustwd_2(its:ite,kts:kte,jts:jte)             ! dust dust loss in dust bin 2
 155 wdbins(its:ite,kts:kte,jts:jte,3)=dustwd_3(its:ite,kts:kte,jts:jte)             ! dust dust loss in dust bin 3
 156 wdbins(its:ite,kts:kte,jts:jte,4)=dustwd_4(its:ite,kts:kte,jts:jte)             ! dust dust loss in dust bin 4
 157 wdbins(its:ite,kts:kte,jts:jte,5)=dustwd_5(its:ite,kts:kte,jts:jte)             ! dust dust loss in dust bin 5
 158
 159 rnrate = precr*3600.            ! precipitation rate (mm/h) - rnrate=precr*3600 using rhowater
 160
 161 dbinmm = dbinm*conver6          ! mean dimension of dust (m) in the bin
 162
 163 tw=0                            ! water temperature (C)
 164 call dviscw(tw,viscw)           ! determine water (dynamic) viscosity viscw
 165
 166 ! raindrop classes
 167 do l = 1,nrbins+1
 168  z = alog10(rmin*converi6)+(alog10(rmax*converi6)-alog10(rmin*converi6))*real(l-1)/real(nrbins)
 169  rend(l) = real(int(10.**z))*conver6
 170 enddo
 171 do l = 1,nrbins
 172  z = (alog10(rend(l)*converi6)+alog10(rend(l+1)*converi6))*0.5
 173  raind(l) = real(int(10.**z))*conver6
 174 enddo
 175 !          bin 1        bin 2                     bin nrbin
 176 !    |--------------|------------!-------------!--------------|
 177 ! rend(1)        rend(2)      rend(3)                    rend(nrbin+1)
 178 !         raind(1)      raind(2)                 raind(nrbin)
 179
 180
 181 do j = jts, jte
 182  do k = kts, kte
 183   do i = its, ite
 184    if ( precr(i,k,j).gt.0.) then ! precipitation exists
 185
 186     call dsd_rain(rnrate(i,k,j),raind,nrbins,dsd_rn,3) ! 3 = Gamma Distribution
 187
 188     tair=t_phy(i,k,j)
 189     pair=p_phy(i,k,j)
 190     rhoair = pair/(rair*tair)     ! dry air density (kg/m3) - rhoair=pa/(Ra*Ta)
 191
 192     call dvisca(tair,visca)
 193
 194     do l = 1,nrbins
 195      call fallv(raind(l),tair,pair,rhoair,visca,vt(l))
 196     enddo
 197
 198     do d=1,nbins ! d = dust bins // dust classes
 199
 200      if ( chem(i,k,j,bins(d)).gt.0.) then ! dust exists
 201
 202       call w_t(wt,dbinmm(d),den_dust(d),rhoair)
 203
 204       scrate=0.
 205       scavn=0.
 206       do l = 1,nrbins
 207        ! calculate collection efficiency for each dust bin (already in this loop) by each rain bin (already in this loop)
 208        call coleff(dbinmm(d),raind(l),den_dust(d),pair,tair,wt,vt(l),rhoair,visca,viscw,colece)
 209        ! calculate scavanging rate
 210        delR = (rend(l+1)-rend(l))*1.e2 ! cm
 211        delv = (vt(l)-wt)*1.e2 ! cm/s
 212        delv = amax1(delv,0.)
 213        rnflx = dsd_rn(l)*delv
 214        carea = pi*0.25*(raind(l)*1.e2+dbinmm(d)*1.e2)**2.
 215        scavn = scavn+carea*colece*rnflx*delR
 216       enddo
 217       scrate = scavn
 218
 219       dustold = chem(i,k,j,bins(d))
 220       chem(i,k,j,bins(d))=chem(i,k,j,bins(d))-max(0.,chem(i,k,j,bins(d))*scrate*dtstepc)
 221       chem(i,k,j,bins(d))=max(chem(i,k,j,bins(d)),epsilc)
 222       wdbins(i,k,j,d)=max(epsilc,dustold-chem(i,k,j,bins(d)))
 223
 224      endif
 225     enddo
 226
 227 ! mklose: fifths size bin is used now, hence the following is not needed:
 228 !     dustold = chem(i,k,j,bins(5))
 229 !     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))
 230 !     chem(i,k,j,bins(5))=max(chem(i,k,j,bins(5)),epsilc)
 231 !     wdbins(i,k,j,5)=max(epsilc,dustold-chem(i,k,j,bins(5)))
 232
 233    endif
 234
 235   enddo
 236  enddo
 237 enddo
 238
 239 do j=jts,jte
 240  do i=its,ite
 241   do k=kts,kte
 242    dustwdload_1(i,j)= max(epsilc,dustwdload_1(i,j) + wdbins(i,k,j,1)/alt(i,k,j) * dz8w(i,k,j))
 243    dustwdload_2(i,j)= max(epsilc,dustwdload_2(i,j) + wdbins(i,k,j,2)/alt(i,k,j) * dz8w(i,k,j))
 244    dustwdload_3(i,j)= max(epsilc,dustwdload_3(i,j) + wdbins(i,k,j,3)/alt(i,k,j) * dz8w(i,k,j))
 245    dustwdload_4(i,j)= max(epsilc,dustwdload_4(i,j) + wdbins(i,k,j,4)/alt(i,k,j) * dz8w(i,k,j))
 246    dustwdload_5(i,j)= max(epsilc,dustwdload_5(i,j) + wdbins(i,k,j,5)/alt(i,k,j) * dz8w(i,k,j))
 247    dustwd_1(i,k,j)=max(epsilc,wdbins(i,k,j,1))
 248    dustwd_2(i,k,j)=max(epsilc,wdbins(i,k,j,2))
 249    dustwd_3(i,k,j)=max(epsilc,wdbins(i,k,j,3))
 250    dustwd_4(i,k,j)=max(epsilc,wdbins(i,k,j,4))
 251    dustwd_5(i,k,j)=max(epsilc,wdbins(i,k,j,5))
 252   enddo
 253    wetdep_1(i,j)= max(epsilc,wdbins(i,kts,j,1)/alt(i,kts,j) * dz8w(i,kts,j) / dtstepc)
 254    wetdep_2(i,j)= max(epsilc,wdbins(i,kts,j,2)/alt(i,kts,j) * dz8w(i,kts,j) / dtstepc)
 255    wetdep_3(i,j)= max(epsilc,wdbins(i,kts,j,3)/alt(i,kts,j) * dz8w(i,kts,j) / dtstepc)
 256    wetdep_4(i,j)= max(epsilc,wdbins(i,kts,j,4)/alt(i,kts,j) * dz8w(i,kts,j) / dtstepc)
 257    wetdep_5(i,j)= max(epsilc,wdbins(i,kts,j,5)/alt(i,kts,j) * dz8w(i,kts,j) / dtstepc)
 258  enddo
 259 enddo
 260
 261
 262 end subroutine uoc_dustwd_driver
 263
 264
 265 !=======================================================================
 266 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
 267 !         - further updates M. Klose
 268
 269       subroutine coleff(d,rain,prho,p,t,pv,rv,arho,visca,viscw,ecolec)
 270
 271 !     the program to calculate collection efficiency
 272 !     it is based on Slinn's empirical equation(1983)
 273 !     and theoretical collision efficiencies by
 274 !     Mason(1971), Klett and Davis (1973) and Beard and Ochs (1984)
 275
 276 ! input
 277 ! =====
 278       real, intent (in)    :: d                                             ! particle diameter in m
 279       real, intent (in)    :: rain                                          ! raindrio diameter in m
 280       real(8), intent (in) :: prho                                          ! particle density in kg/m3
 281       real, intent (in)    :: p                                             ! atmosphere pressure in Pa
 282       real, intent (in)    :: t                                             ! atmospheric temperature in K
 283       real, intent (in)    :: pv                                            ! dust fall velocity in m/s
 284       real, intent (in)    :: rv                                            ! rain fall velocity in m/s
 285       real, intent (in)    :: arho                                          ! rair density in kg/m3
 286       real, intent (in)    :: visca, viscw
 287
 288 ! output
 289 ! ======
 290       real, intent (out) ::  ecolec
 291
 292 ! local variables
 293 ! ===============
 294       real    :: beta
 295       real    :: rmin
 296       real(8) :: E_im,E_in,E_br
 297       real    :: dd, raind, rhop, pr, ta, vp, vs, rhoa
 298
 299       CHARACTER*(100) :: text
 300       integer :: debug_level
 301
 302       data rhow/1/                                                       ! raindrop density g/cm^3
 303
 304       CALL get_wrf_debug_level(debug_level)
 305
 306       dd = d*1.e6       ! m to um
 307       raind = rain*1.e6 ! m to um
 308       rhop = prho*1.e-3 ! kg/m3 to g/cm3
 309       pr = p*1.e-2      ! Pa to mb or hPa
 310       ta = t-273.15     ! K to C
 311       vp = pv*1.e2      ! m/s to cm/s
 312       vs = rv*1.e2      ! m/s to cm/s
 313       rhoa = arho*1.e-3 ! kg/m3 to g/cm3
 314
 315       !=======================================================================
 316       ecol=0
 317
 318       if(dd.le.2) then
 319         call eslinn(raind,dd,pr,ta,rhop,vp,vs,rhoa,visca,viscw,ecol)
 320       elseif(dd.gt.2 .and. dd.le.40) then
 321         call eimpact(rhop,dd,raind,ecol)
 322       elseif(dd.gt.40) then
 323         ecol=0
 324       endif
 325
 326       ecolec=ecol
 327
 328 !       return
 329       end subroutine coleff
 330 !****************************************************************
 331
 332
 333 !****************************************************************
 334 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
 335 !         - further updates M. Klose
 336
 337 ! xlambda is set fix instead of calculated - cfrick
 338
 339       subroutine eslinn(dr,dp,pair,tair,rhop,vp,vs,rhoa,visca,viscw,eslin)
 340
 341 ! input
 342 ! =====
 343       real, intent (in) :: dp                                             ! particle diameter in um
 344       real, intent (in) :: dr                                             ! raindrio diameter in um
 345       real, intent (in) :: rhop                                           ! particle density in g/cm3
 346       real, intent (in) :: pair                                           ! atmosphere pressure in hPa
 347       real, intent (in) :: tair                                           ! atmospheric temperature in C
 348       real, intent (in) :: vp                                             ! dust fall velocity in cm/s
 349       real, intent (in) :: vs                                             ! rain fall velocity in cm/s
 350       real, intent (in) :: rhoa                                           ! rair density in g/cm3
 351       real, intent (in) :: visca, viscw
 352
 353 ! output
 354 ! ======
 355       real, intent (out) ::  eslin
 356
 357       real :: pi
 358
 359       CHARACTER*(100) :: text
 360       integer :: debug_level
 361
 362       real :: gnarf, ren
 363
 364       E_br=0
 365       E_in=0
 366       E_im=0
 367
 368       drc=dr*1.e-4                                        ! rain diameter in cm
 369       dpc=dp*1.e-4                                        ! particle diameter in cm
 370
 371       pr=pair
 372       ta=tair
 373
 374       pi=acos(-1.0)
 375       xk=1.381e-16                                        ! Boltzmann constant (g cm2/s2)/K
 376 !      xd=3.75e-8
 377 !      prs=pr*1.e+3
 378 !      xlambda=xk*(ta+273.15)/(pi*sqrt(2.0)*prs*xd*xd)
 379
 380       omg=viscw/visca
 381       ren=0.5*drc*vs*rhoa/visca                           ! reynolds number
 382       sstar=(1.2+alog(1+ren)/12)/(1+alog(1+ren))          ! S*
 383
 384       xlambda=0.0651*1.e-4 ! mean free path of air in cm
 385       cc=1+(2*xlambda/dpc)*(1.257+0.4*exp(-1.1*dpc/(2*xlambda))) ! Cunningham's correction for small particles
 386
 387       tau=dpc**2*rhop*cc/(18*visca)
 388       stn=2*tau*(vs-vp)/drc                               ! Stokes number
 389       diff=(xk*(ta+273.15)*cc)/(3*pi*visca*dpc)
 390       scn=visca/(diff*rhoa)                               ! Schmidt number
 391       phi=dpc/drc
 392       pen=ren*scn                                         ! peclet number
 393
 394 ! Brownian diffusion
 395
 396       E_br=4*(1+0.4*sqrt(ren)*scn**(1.0/3.0)+0.16*sqrt(ren)*sqrt(scn))/pen
 397
 398 ! Interception
 399
 400       E_in=4*phi*(1/omg+(1+2*sqrt(ren))*phi)
 401
 402 ! Impaction
 403
 404       if(stn.ge.sstar) then
 405           E_im=((stn-sstar)/(stn-sstar+2.0/3.0))**1.5
 406       endif
 407
 408       eslin=E_br+E_in
 409
 410 !       return
 411       end subroutine eslinn
 412 !=======================================================================
 413
 414 !*******************************************************************************
 415 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
 416 !         - further updates M. Klose
 417
 418 subroutine eimpact(rhop,dp,raind,e_im)
 419
 420     USE module_data_uoc_wd   ! mklose [17122014]: renamed from we_ematrix
 421
 422       integer, parameter :: lp=17, lr=17
 423
 424 !       include 'we_ematrix.dat'
 425
 426 ! input
 427 ! =====
 428       real, intent (in) :: dp                                             ! particle diameter in um
 429       real, intent (in) :: raind                                          ! raindrio diameter in um
 430       real, intent (in) :: rhop                                           ! particle density in g/cm3
 431
 432 ! output
 433 ! ======
 434       real, intent (out) ::  e_im
 435
 436       real  ::  ee(lr)
 437       real  ::  u(lr+1)
 438       real  ::  enew(lr+1)
 439       real  ::  e_r(lp)
 440       real  ::  e_c(lp)
 441       real  ::  f(lp)
 442       real  ::  v(lp+1)
 443       real  ::  dnew(lp+1)
 444
 445       CHARACTER*(100) :: text
 446       integer :: debug_level
 447
 448       CALL get_wrf_debug_level(debug_level)
 449
 450 ! rr, rp       : radiuse of raindrop, dust particle in um
 451 ! edata        : dataset of composite collection efficiency
 452 !===============================================================================
 453
 454       if(raind.lt.100 .or. raind.gt.6000) then
 455        text = 'raindrop diameter out of range'
 456        text=trim(trim(text)//" - ERROR - UoC dust wet deposition")
 457        call wrf_debug (debug_level,text)
 458       endif
 459       if(dp.lt.2 .or. dp.gt.40) then
 460        text = 'particle diameter out of range'
 461        text=trim(trim(text)//" - ERROR - UoC dust wet deposition")
 462        call wrf_debug (debug_level,text)
 463       endif
 464
 465 ! radius of raindrop, dust particle
 466
 467       rr=raind*0.5
 468       rp=dp*0.5
 469
 470       do i=1,lp
 471 ! calculate e_r(1:lp) corresponding to dustr(1:lp) for rr
 472
 473         jj=0
 474         do j=1,lr
 475           ee(j)=edatawd(i,j)
 476           if(rr.ge.rainwd(j)) jj=j
 477         enddo
 478         if(jj.eq.0) then
 479          text = 'error in raindrop radius'
 480          text=trim(trim(text)//" - ERROR - UoC dust wet deposition")
 481          call wrf_debug (debug_level,text)
 482         endif
 483         if(rr.gt.rainwd(jj)) then
 484            lrp1=lr+1
 485            do l=1,jj
 486              u(l)=rainwd(l)
 487            enddo
 488              u(jj+1)=rr
 489            do l=jj+1,lr
 490              u(l+1)=rainwd(l)
 491            enddo
 492            call intrpl(lr,log(rainwd),ee,lrp1,log(u),enew)
 493            e_r(i)=enew(jj+1)
 494         elseif(rr.eq.rainwd(jj)) then
 495            e_r(i)=ee(jj)
 496         endif
 497       enddo
 498
 499 ! interpolate e_d for rp from e_r
 500
 501       rmax=sqrt(1/rhop)
 502
 503       do i=1,lp
 504         if(dustwd(i).le.1) then
 505            f(i)=1
 506         elseif(dustwd(i).gt.1 .and. dustwd(i).lt.3.98) then
 507            f(i)=(1-rmax)*(dustwd(i)-3.98)**2/(1-3.98)**2+rmax
 508         elseif(dustwd(i).ge.3.98) then
 509            f(i)=rmax
 510         endif
 511       enddo
 512
 513       do i=1,lp
 514         e_c(i)=e_r(i)*f(i)
 515       enddo
 516
 517       ii=0
 518       e_d=0
 519       do i=1,lp
 520         if(rp.ge.dustwd(i)) ii=i
 521       enddo
 522       if(ii.eq.0) then
 523        text = 'wrong particle radius'
 524        text=trim(trim(text)//" - ERROR - UoC dust wet deposition")
 525        call wrf_debug (debug_level,text)
 526       endif
 527       if(rp.gt.dustwd(ii)) then
 528          do l=1,ii
 529            v(l)=dustwd(l)
 530          enddo
 531            v(ii+1)=rp
 532          do l=ii+1,lp
 533            v(l+1)=dustwd(l)
 534          enddo
 535            lpp1=lp+1
 536            call intrpl(lp,dustwd,e_c,lpp1,v,dnew)
 537            e_d=dnew(ii+1)
 538       elseif(rp.eq.dustwd(ii)) then
 539            e_d=e_c(ii)
 540       endif
 541
 542       e_im=e_d
 543
 544 !       return
 545       end subroutine eimpact
 546 !=======================================================================
 547
 548 !***********************************************************************
 549 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
 550 !         - further updates M. Klose
 551
 552       subroutine dvisca(tair,dvisc)
 553
 554 ! tair  : air temperature in K
 555 ! visca : dynamic viscosity of air in g/cm/s
 556
 557       integer, parameter :: l=9, n=1
 558
 559       real, dimension(l) :: ta, visca
 560       real, dimension(l) :: x, y
 561       real, dimension(n) :: u, v
 562
 563       data ta/-173,-73,0,20,25,27,127,227,327/
 564       data visca/0.71e-4,1.33e-4,1.72e-4,1.797e-4,1.818e-4, &
 565      &           1.86e-4,2.31e-4,2.71e-4,3.08e-4/
 566
 567 ! values from crc handbook of chemistry and physics
 568
 569       do i = 1,l
 570          x(i)=ta(i)+273
 571          y(i)=visca(i)
 572       enddo
 573
 574       u(1)=tair
 575
 576       call intrpl(l,x,y,n,u,v)
 577
 578       dvisc=v(1)
 579
 580 !       return
 581       end subroutine dvisca
 582 !
 583 !************************************************************
 584
 585 !************************************************************
 586 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
 587 !         - further updates M. Klose
 588
 589       subroutine dviscw(twt,dvisc)
 590
 591       integer, parameter :: l=11, n=1
 592
 593       real, dimension(l) :: tw, viscw
 594       real, dimension(l) :: x, y
 595       real, dimension(l) :: u, v
 596
 597       data tw/0,10,20,30,40,50,60,70,80,90,100/
 598       data viscw/1.7930e-2,1.3070e-2,1.002e-2,0.7977e-2,0.6532e-2, &
 599      &           0.5470e-2,0.4665e-2,0.404e-2,0.3544e-2,0.3145e-2, &
 600      &           0.2818e-2/
 601
 602 ! tw : water temperature in degree C
 603
 604 ! dynamic viscosity of water in g/cm/s
 605 ! values from crc handbook of chemistry and physics
 606 ! properties of water in the range 0 - 100c
 607
 608       do i = 1,l
 609          x(i)=tw(i)+273.15
 610          y(i)=viscw(i)
 611       enddo
 612
 613       u(1)=twt+273.15
 614
 615       call intrpl(l,x,y,n,u,v)
 616
 617       dvisc=v(1)
 618
 619 !       return
 620       end subroutine dviscw
 621
 622 !*******************************************************************************
 623
 624 !************************************************************
 625 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
 626 !         - further updates by M. Klose, 2015
 627
 628 subroutine  intrpl(l,x,y,n,u,v)
 629 !************************************************************
 630 !     algorithm appeared in comm. acm, vol. 15, no. 10,
 631 !     p. 914.
 632 !
 633 ! interpolation of a single-valued function
 634 ! this subroutine interpolates, from values of the function
 635 ! given as ordinates of input data points in an x-y plane
 636 ! and for a given set of x values (abscissas), the values of
 637 ! a single-valued function y = y(x).
 638 ! the input parameters are
 639 !     l  = number of input data points
 640 !          (must be 2 or greater)
 641 !     x  = array of dimension l storing the x values
 642 !          (abscissas) of input data points
 643 !          (in ascending order)
 644 !     y  = array of dimension l storing the y values
 645 !          (ordinates) of input data points
 646 !     n  = number of points at which interpolation of the
 647 !          y value (ordinate) is desired
 648 !          (must be 1 or greater)
 649 !     u  = array of dimension n storing the x values
 650 !          (abscissas) of desired points
 651 ! the output parameter is
 652 !     v  = array of dimension n where the interpolated y
 653 !          values (ordinates) are to be displayed
 654 ! declaration statements
 655       integer, intent(in) :: l, n
 656       real, dimension(l)  :: x, y
 657       real, dimension(n)  :: u, v
 658       equivalence  (p0,x3),(q0,y3),(q1,t3)
 659       real                :: m1,m2,m3,m4,m5
 660       equivalence  (uk,dx),(imn,x2,a1,m1),(imx,x5,a5,m5), &
 661      &             (j,sw,sa),(y2,w2,w4,q2),(y5,w3,q3)
 662       real    :: a1, a2, a3, a4, a5
 663       integer :: debug_level
 664       CALL get_wrf_debug_level(debug_level)
 665 ! preliminary processing
 666
 667    10 l0=l
 668       lm1=l0-1
 669       lm2=lm1-1
 670       lp1=l0+1
 671       n0=n
 672       if(lm2.lt.0)        go to 90
 673       if(n0.le.0)         go to 91
 674       do 11  i=2,l0
 675         if(x(i-1)-x(i))   11,95,96
 676    11   continue
 677       ipv=0
 678 ! main do-loop
 679       do 80  k=1,n0
 680         uk=u(k)
 681 ! routine to locate the desired point
 682    20   if(lm2.eq.0)      go to 27
 683         if(uk.ge.x(l0))   go to 26
 684         if(uk.lt.x(1))    go to 25
 685         imn=2
 686         imx=l0
 687    21   i=(imn+imx)/2
 688         if(uk.ge.x(i))    go to 23
 689    22   imx=i
 690         go to 24
 691    23   imn=i+1
 692    24   if(imx.gt.imn)    go to 21
 693         i=imx
 694         go to 30
 695    25   i=1
 696         go to 30
 697    26   i=lp1
 698         go to 30
 699    27   i=2
 700 ! check if i = ipv
 701    30   if(i.eq.ipv)      go to 70
 702         ipv=i
 703 ! routines to pick up necessary x and y values and
 704 !          to estimate them if necessary
 705    40   j=i
 706         if(j.eq.1)        j=2
 707         if(j.eq.lp1)      j=l0
 708         x3=x(j-1)
 709         y3=y(j-1)
 710         x4=x(j)
 711         y4=y(j)
 712         a3=x4-x3
 713         m3=(y4-y3)/a3
 714         if(lm2.eq.0)      go to 43
 715         if(j.eq.2)        go to 41
 716         x2=x(j-2)
 717         y2=y(j-2)
 718         a2=x3-x2
 719         m2=(y3-y2)/a2
 720         if(j.eq.l0)       go to 42
 721    41   x5=x(j+1)
 722         y5=y(j+1)
 723         a4=x5-x4
 724         m4=(y5-y4)/a4
 725         if(j.eq.2)        m2=m3+m3-m4
 726         go to 45
 727    42   m4=m3+m3-m2
 728         go to 45
 729    43   m2=m3
 730         m4=m3
 731    45   if(j.le.3)        go to 46
 732         a1=x2-x(j-3)
 733         m1=(y2-y(j-3))/a1
 734         go to 47
 735    46   m1=m2+m2-m3
 736    47   if(j.ge.lm1)      go to 48
 737         a5=x(j+2)-x5
 738         m5=(y(j+2)-y5)/a5
 739         go to 50
 740    48   m5=m4+m4-m3
 741 ! numerical differentiation
 742    50   if(i.eq.lp1)      go to 52
 743         w2=abs(m4-m3)
 744         w3=abs(m2-m1)
 745         sw=w2+w3
 746         if(sw.ne.0.0)     go to 51
 747         w2=0.5
 748         w3=0.5
 749         sw=1.0
 750    51   t3=(w2*m2+w3*m3)/sw
 751         if(i.eq.1)        go to 54
 752    52   w3=abs(m5-m4)
 753         w4=abs(m3-m2)
 754         sw=w3+w4
 755         if(sw.ne.0.0)     go to 53
 756         w3=0.5
 757         w4=0.5
 758         sw=1.0
 759    53   t4=(w3*m3+w4*m4)/sw
 760         if(i.ne.lp1)      go to 60
 761         t3=t4
 762         sa=a2+a3
 763         t4=0.5*(m4+m5-a2*(a2-a3)*(m2-m3)/(sa*sa))
 764         x3=x4
 765         y3=y4
 766         a3=a2
 767         m3=m4
 768         go to 60
 769    54   t4=t3
 770         sa=a3+a4
 771         t3=0.5*(m1+m2-a4*(a3-a4)*(m3-m4)/(sa*sa))
 772         x3=x3-a4
 773         y3=y3-m2*a4
 774         a3=a4
 775         m3=m2
 776 ! determination of the coefficients
 777    60   q2=(2.0*(m3-t3)+m3-t4)/a3
 778         q3=(-m3-m3+t3+t4)/(a3*a3)
 779 ! computation of the polynomial
 780    70   dx=uk-p0
 781    80   v(k)=q0+dx*(q1+dx*(q2+dx*q3))
 782       return
 783 ! error exit
 784    90 call wrf_debug (debug_level,'error 90 in subroutine intrpl - ERROR - UoC dust wet deposition') !write (6,2090)
 785       go to 99
 786    91 call wrf_debug (debug_level,'error 91 in subroutine intrpl - ERROR - UoC dust wet deposition') !write (6,2091)
 787       go to 99
 788    95 call wrf_debug (debug_level,'error 95 in subroutine intrpl - ERROR - UoC dust wet deposition') !write (6,2095)
 789       go to 97
 790    96 call wrf_debug (debug_level,'error 96 in subroutine intrpl - ERROR - UoC dust wet deposition') !write (6,2096)
 791    97 call wrf_debug (debug_level,'error 97 in subroutine intrpl - ERROR - UoC dust wet depositionk') !write (6,2097)  i,x(i)
 792    99 call wrf_debug (debug_level,'error 98 in subroutine intrpl - ERROR - UoC dust wet deposition') !write (6,2099)  l0,n0
 793       return
 794 ! format statements
 795 ! 2090 format(1x/22h  ***   l = 1 or less./)
 796 ! 2091 format(1x/22h  ***   n = 0 or less./)
 797 ! 2095 format(1x/27h  ***   identical x values./)
 798 ! 2096 format(1x/33h  ***   x values out of sequence./)
 799 ! 2097 format(6h   i =,i7,10x,6hx(i) =,e12.3)
 800 ! 2099 format(6h   l =,i7,10x,3hn =,i7/&
 801 !     &       36h error detected in routine    intrpl)
 802       end subroutine intrpl
 803 !************************************************************
 804
 805 !=======================================================================
 806 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
 807 !         - further updates M. Klose
 808
 809       subroutine fallv(d,t,p,rhoair,visc,vt)
 810 !=======================================================================
 811 !      the program to calculate the termincal velocity of raindrop
 812 !      this is based on Beard (1976)
 813 !
 814 !      input
 815 !      ta    : temperature in K
 816 !      dropd : raindrop diameter in m
 817 !
 818 !      output
 819 !      vt    : terminal velocity in m/s
 820 !=======================================================================
 821
 822 real :: gg, pr0, t0, visc0, rhow, rhoa, dp, dl, cl, csc, vt, ta, dropd, pr
 823 real :: a0, a1, a2, a3, a4, a5, a6
 824 real :: b0, b1, b2, b3, b4, b5
 825 real :: c1, c2, dan, x, vy, rey, c3, bon, phn, rhoair, p, t, d
 826
 827  CHARACTER*(100) :: text
 828  integer :: debug_level
 829  CALL get_wrf_debug_level(debug_level)
 830
 831 ! compute surface tension
 832 ! sig : the surface tension  in mN/m
 833
 834         dropd = d*1.e6 ! m to um
 835         ta = t-273.15  ! K to C
 836         pr = p*1.e-2   ! Pa to hPa (mb)
 837
 838         call sfctens(ta,sig)
 839
 840         gg=980                                     ! cm/s^2
 841         pr0=1013.25                               ! mb
 842         t0=293.15                                 ! K
 843         visc0=1.818e-4                            ! g/cm/sec
 844         rhow=1                                    ! g/cm^3
 845
 846         rhoa = rhoair*1.e-3 ! kg/m3 to g/cm3
 847
 848 ! compute the terminal velocity
 849
 850         vt=0
 851
 852         dp=dropd
 853         dl=dropd*1.e-4                            ! diameter in cm
 854
 855         c1=(rhow-rhoa)*gg/(18*visc)
 856         cl=6.62e-6*(visc/visc0)*(pr0/pr)*sqrt((ta+273.15)/t0)  !cm
 857         csc=1+2.51*cl/dl
 858
 859         if(dp.lt.19.) then
 860
 861          vt=c1*csc*dl**2
 862          vt=vt*1.e-2 ! cm/s to m/s
 863
 864         elseif(dp.ge.19. .and. dp.lt.1070.) then
 865
 866          a0=-0.318657e+1
 867          a1=+0.992696
 868          a2=-0.153193e-2
 869          a3=-0.987059e-3
 870          a4=-0.578878e-3
 871          a5=+0.855176e-4
 872          a6=-0.327815e-5
 873
 874          c2=4*rhoa*(rhow-rhoa)*gg/(3*visc*visc)
 875          dan=c2*dl**3
 876          x=alog(dan)
 877          vy=a0+a1*x+a2*x**2+a3*x**3+a4*x**4+a5*x**5+a6*x**6
 878          rey=csc*exp(vy)
 879          vt=visc*rey/(rhoa*dl)
 880          vt=vt*1.e-2 ! cm/s to m/s
 881
 882         elseif(dp.ge.1070. .and. dp.le.7000.) then
 883
 884          b0=-0.500015e+1
 885          b1=+0.523778e+1
 886          b2=-0.204914e+1
 887          b3=+0.475294
 888          b4=-0.542819e-1
 889          b5=+0.238449e-2
 890
 891          c3=4*(rhow-rhoa)*gg/(3*sig)
 892          bon=c3*dl*dl                               ! Bond Number
 893          phn=sig**3*rhoa**2/(visc**4*gg*(rhow-rhoa)) ! Physical Property Number
 894
 895          y=alog(bon*phn**(1.0/6.0))
 896          vy=b0+b1*y+b2*y**2+b3*y**3+b4*y**4+b5*y**5
 897          rey=phn**(1.0/6.0)*exp(vy)
 898          vt=visc*rey/(rhoa*dl)
 899          vt=vt*1.e-2 ! cm/s to m/s
 900
 901         endif
 902
 903 !         return
 904         end subroutine fallv
 905 !=======================================================================
 906
 907 !=======================================================================
 908 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
 909 !         - further updates, M. Klose
 910
 911         subroutine sfctens(temp,sig)
 912
 913 ! temp : temperature in K
 914 ! sig  : surface tension in mN/m or dyne/cm
 915         integer, parameter :: n=14
 916         real, dimension(n) :: t, sfctn
 917         data t/0,5,10,15,18,20,25,30,40,50,60,70,80,100/
 918         data sfctn/75.6,74.9,74.22,73.49,73.05,72.75,71.97,71.18,69.56, &
 919      &             67.91,66.18,64.4,62.6,58.9/
 920
 921         real :: ta
 922
 923         kk = 1
 924         xsfc = 0
 925         ta=max(temp-273.15,0.0)
 926
 927         do i = 1, 14
 928            if(i.lt.14) then
 929               if(ta.ge.t(i) .and. ta.lt.t(i+1)) then
 930                  kk=i
 931                  goto 1
 932               endif
 933            else
 934               if(ta.ge.t(i)) then
 935                  kk=14
 936                  goto 1
 937               endif
 938            endif
 939         enddo
 940
 941 1       continue
 942
 943         if(kk.lt.14) then
 944              sig = sfctn(kk)+ &
 945      &       (sfctn(kk+1)-sfctn(kk))*(ta-t(kk))/(t(kk+1)-t(kk))
 946         else
 947              sig = sfctn(kk)+ &
 948      &       (sfctn(kk)-sfctn(kk-1))*(ta-t(kk))/(t(kk)-t(kk-1))
 949         endif
 950
 951 !         return
 952         end subroutine sfctens
 953 !***********************************************************************
 954
 955 !***********************************************************************
 956 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
 957 !         - further updates M. Klose
 958
 959         subroutine dsd_rain(rnrate,raind,nrbin,dsd_rn,ldsd)
 960
 961 ! input
 962         real, intent(in) :: raind(nrbin), rnrate ! in m and mm/h
 963         integer, intent(in) :: ldsd, nrbin
 964 ! local
 965         real :: dsd_mp(nrbin)
 966         real :: dsd_ss(nrbin)
 967         real :: dsd_wt(nrbin)
 968         real :: dsd_fl(nrbin)
 969         real :: ng
 970         real :: nt
 971         real :: pi, dr, wc, d0, dcm, const, dumm
 972 ! output
 973         real, intent(out) :: dsd_rn(nrbin)
 974 !.......................................................................
 975
 976         d0=0
 977         do nk=1,nrbin
 978
 979           dr=raind(nk)*1.e3                            ! m to mm
 980
 981           if(ldsd.eq.1) then
 982 !
 983 ! Marshall-Palmer(1948)
 984           dsd_mp(nk)=8.e3*exp(-4.1*dr*rnrate**(-0.21))  ! 1/m^3/mm
 985           dsd_mp(nk)=dsd_mp(nk)*1.e-5                   ! 1/cm^3/cm
 986           dsd_rn(nk)=dsd_mp(nk)
 987 !
 988           elseif(ldsd.eq.2) then
 989
 990 ! Sekhon and Srivastava
 991           dsd_ss(nk)=0.07*rnrate**0.37*exp(-3.8*dr/rnrate**0.14) ! cm^(-3) cm^(-1)
 992           dsd_rn(nk)=dsd_ss(nk)
 993 !
 994           elseif(ldsd.eq.3) then
 995
 996 ! Willis-Tattelman(1989)
 997           wc=0.062*rnrate**0.913                        ! water content in g/m^3
 998           d0=0.1571*wc**0.1681                          ! median volume diameter in cm
 999           ng=512.85*wc*1.e-6/d0**4*d0**(-2.16)
1000           dcm=dr*0.1                                    ! diameter in cm
1001           dsd_wt(nk)=ng*dcm**2.16*exp(-5.5880*dcm/d0)   ! #/cm^3/cm
1002           dsd_rn(nk)=dsd_wt(nk)
1003
1004           elseif(ldsd.eq.4) then
1005
1006 ! Feingold-Levin(1986)
1007           pi=acos(-1.0)
1008           const=sqrt(2*pi)*alog(1.43)
1009           nt=172*rnrate**0.22                           ! 1/m^3
1010           dg=0.72*rnrate**0.21
1011           dumm=0.5*alog(dr/dg)**2/alog(1.43)**2
1012           dsd_fl(nk)=nt/const/dr*exp(-dumm)             ! 1/m^3/mm
1013           dsd_fl(nk)=dsd_fl(nk)*1.e-5                   ! 1/cm^3/cm
1014           dsd_rn(nk)=dsd_fl(nk)
1015
1016           endif
1017 !
1018         enddo
1019
1020         return
1021         end subroutine dsd_rain
1022
1023 !***********************************************************************
1024
1025
1026 !***********************************************************************
1027 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
1028 !         - further updates M. Klose
1029
1030 subroutine w_p(d,rho,wt)
1031
1032 ! Yaping Shao 28-07-92
1033
1034 ! Calculate Terminal velocity for spherical grain
1035 ! wt  : terminal velocity    [m/s]
1036 ! d   : particle diameter    [m]
1037 ! rho : particle density     [kg/m3]
1038 ! Require CDSPH
1039 ! Parameter specifications
1040
1041 IMPLICIT NONE ! cfrick
1042
1043       REAL, INTENT(OUT) :: wt
1044       REAL              :: lambda
1045       REAL              :: rden, err, test, x, xl, xh
1046       REAL              :: f, fm, fl, fh
1047       INTEGER           :: i, isign
1048       real, intent(in)  :: rho, d
1049       real              :: cc
1050
1051       REAL, PARAMETER :: visc=1.5E-05
1052       INTEGER :: debug_level
1053
1054       CHARACTER*(100) :: text
1055
1056 ! VARIABLE X IS WT, FUNCTION F IS: CD(RE)*WT**2 - 4*D*G*RDEN/3
1057
1058
1059        CALL get_wrf_debug_level(debug_level)
1060
1061        rden=rho ! particle density
1062        lambda=0.0651*1.e-6 ! mean free path of air
1063        cc=1+(2*lambda/d)*(1.257+0.4*exp(-1.1*d/(2*lambda))) ! Cunningham's correction for small particles
1064
1065        err=1E-6
1066        xl=1E-12
1067        xh=1E3
1068
1069        wt = xl
1070        fl = CDSPH(wt*d/visc)*wt**2/cc - 4.0*d*g*rden/3.0
1071
1072        wt = xh
1073        fh = CDSPH(wt*d/visc)*wt**2/cc - 4.0*d*g*rden/3.0
1074
1075        IF (fh*fl.GT.0) call wrf_debug (debug_level,'w_t: F(XL), F(XH) HAVE SAME SIGN - ERROR - UoC dust wet deposition')
1076        IF (fh-fl.EQ.0) call wrf_debug (debug_level,'w_t: F(XL) = F(XH) - ERROR - UoC dust wet deposition')
1077
1078        isign=1
1079        IF (fh.LT.fl) isign=-1
1080
1081        fh=isign*fh
1082        fl=isign*fl
1083
1084        DO 1 i=1,100
1085          x=(xl+xh)/2
1086          wt = x
1087          f = CDSPH(wt*d/visc)*wt**2/cc - 4.0*d*g*rden/3.0
1088          fm=isign*f
1089          IF (fm.GT.fh.OR.fm.LT.fl) call wrf_debug (debug_level,'w_t: F(X) NON-MONOTONIC - ERROR - UoC dust wet deposition')
1090          IF (fm.GT.0) xh=x
1091          IF (fm.LT.0) xl=x
1092          IF (fm.EQ.0) then
1093           return
1094          ENDIF
1095          test=ABS(xh-xl)/ABS(x)
1096          IF (test.LT.ABS(err)) then
1097           return
1098          ENDIF
1099 1      CONTINUE
1100        call wrf_debug (debug_level,'w_t: NO SOLUTION IN 100 LOOPS - ERROR - UoC dust wet deposition')
1101
1102 END subroutine w_p
1103 !***************************************************************
1104
1105
1106 !***********************************************************************
1107 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
1108 !         - further updates M. Klose
1109
1110 SUBROUTINE w_t (wt, dmm, rhop, rhoa)
1111
1112 ! Yaping Shao 28-07-92
1113
1114 ! Calculate Terminal velocity for spherical grain
1115 ! wt : terminal velocity    [m/s]
1116 ! dmm: particle diameter in [m]
1117 ! rhop: particle density    [kg/m3]
1118 ! rhoa: air density         [kg/m3]
1119 ! Require CDSPH
1120 ! Parameter specifications
1121
1122 ! include Cunningham's correction - Claudia Frick 24-04-2014
1123
1124       REAL(8), INTENT(IN) :: rhop
1125       REAL,    INTENT(IN) :: rhoa
1126
1127       REAL :: wt, dmm
1128       REAL :: RDEN, D, ERR, TEST, X, XL, XH
1129       REAL :: F, FM, FL, FH
1130       REAL, PARAMETER :: VISC=1.5E-05
1131       INTEGER :: I, ISIGN
1132
1133       CHARACTER*(100) :: text
1134       REAL            :: lambda
1135       real            :: cc
1136
1137
1138 ! VARIABLE X IS WT, FUNCTION F IS: CD(RE)*WT**2 - 4*D*G*RDEN/3
1139
1140       D = dmm
1141       RDEN=rhop/rhoa ! = sigma
1142       lambda=0.0651*1.e-6 ! mean free path of air
1143       cc=1+(2*lambda/D)*(1.257+0.4*exp(-1.1*D/(2*lambda))) ! Cunningham's correction for small particles
1144
1145       ERR=1E-6
1146       XL=1E-12
1147       XH=1E3
1148       WT = XL
1149       FL = CDSPH(WT*D/VISC)*WT**2/cc - 4.0*D*G*RDEN/3.0
1150       WT = XH
1151       FH = CDSPH(WT*D/VISC)*WT**2/cc - 4.0*D*G*RDEN/3.0
1152       IF (FH*FL.GT.0) call wrf_debug (debug_level,'w_t: F(XL), F(XH) HAVE SAME SIGN - ERROR - UoC dust wet deposition')
1153       IF (FH-FL.EQ.0) call wrf_debug (debug_level,'w_t: F(XL) = F(XH) - ERROR - UoC dust wet deposition')
1154       ISIGN=1
1155       IF (FH.LT.FL) ISIGN=-1
1156       FH=ISIGN*FH
1157       FL=ISIGN*FL
1158       DO 1 I=1,100
1159         X=(XL+XH)/2
1160         WT = X
1161         F = CDSPH(WT*D/VISC)*WT**2/cc - 4.0*D*G*RDEN/3.0
1162         FM=ISIGN*F
1163         IF (FM.GT.FH.OR.FM.LT.FL) call wrf_debug (debug_level,'w_t: F(X) NON-MONOTONIC - ERROR - UoC dust wet deposition')
1164         IF (FM.GT.0) XH=X
1165         IF (FM.LT.0) XL=X
1166         IF (FM.EQ.0) then
1167          return
1168         ENDIF
1169         TEST=ABS(XH-XL)/ABS(X)
1170         IF (TEST.LT.ABS(ERR)) then
1171          return
1172         ENDIF
1173 1     CONTINUE
1174       call wrf_debug (debug_level,'w_t: NO SOLUTION IN 100 LOOPS - ERROR - UoC dust wet deposition')
1175
1176 END subroutine w_t
1177 !***************************************************************
1178
1179
1180 !---------------------------------------------------------------
1181 ! CEMSYS5 - smaller modifications by cfrick for WRF implementation
1182 real FUNCTION CDSPH(RE)
1183
1184 ! MRR, 30-SEP-87
1185 ! ADAPTED 10-OCT-89 (NO ISTOKE, NO CRASH FOR RE>50000)
1186 ! CALCULATES DRAG COEFFICIENT CD FOR A SPHERE, AS A FUNCTION OF
1187 ! REYNOLDS NUMBER RE, WHERE:
1188 ! DRAG FORCE = (RHO * PI * D**2 * U**2) / 8
1189 !         RE = U * D / VISCOSITY
1190 !         (D = SPHERE DIAMETER, U = VELOCITY, RHO = FLUID DENSITY)
1191 ! ALGORITHM FROM MORSI AND ALEXANDER (1972, JFM 55, 193-208), CHECKED
1192 ! AGAINST THEIR TABLE FOR CD(SPHERE).
1193 ! FOR RE > 50000, SET CDSPH = 0.48802, THE VALUE AT RE=50000=5*10**4.
1194 ! THIS IS OK TO RE=3*10**5 (APPROX), BEYOND WHICH DRAG CRISIS OCCURS
1195 ! AND SETS CDSPH TO ABOUT 0.1.
1196
1197 IMPLICIT NONE ! cfrick
1198
1199 real :: RE
1200
1201 IF (RE.LE.0) THEN
1202  WRITE(6,*) RE
1203  STOP 'CDSPH: RE.LE.0'
1204 ELSE IF (RE.LE.0.1) THEN
1205  CDSPH = 24/RE
1206 ELSE IF (RE.LE.1) THEN
1207  CDSPH = 22.73/RE + 0.0903/RE**2 + 3.69
1208 ELSE IF (RE.LE.10) THEN
1209  CDSPH = 29.1667/RE - 3.8889/RE**2 + 1.222
1210 ELSE IF (RE.LE.100) THEN
1211  CDSPH = 46.5/RE - 116.67/RE**2 + 0.6167
1212 ELSE IF (RE.LE.1000) THEN
1213  CDSPH = 98.33/RE - 2778/RE**2 + 0.3644
1214 ELSE IF (RE.LE.5000) THEN
1215  CDSPH = 148.62/RE - 47500/RE**2 + 0.357
1216 ELSE IF (RE.LE.10000) THEN
1217  CDSPH = -490.546/RE + 578700/RE**2 + 0.46
1218 ELSE IF (RE.LE.50000) THEN
1219  CDSPH = -1662.5/RE + 5416700/RE**2 + 0.5191
1220 ELSE
1221  CDSPH = 0.48802
1222 END IF
1223
1224 END FUNCTION CDSPH
1225 !--------------------------------------------------------------------
1226
1227
1228 END MODULE module_uoc_dustwd