chem/module_mosaic_drydep.F

   1 !**********************************************************************************
   2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
   3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
   4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
   5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
   6 !
   7 ! MOSAIC module: see module_mosaic_driver.F for references and terms of use
   8 !**********************************************************************************
   9         module module_mosaic_drydep
  10
  11
  12 ! rce 2005-feb-18 - several fixes for indices of dlo/hi_sect, volumlo/hi_sect,
  13 !     which are now (isize,itype)
  14
  15 ! rce 2004-dec-03 - many changes associated with the new aerosol "pointer"
  16 !     variables in module_data_mosaic_asect
  17
  18
  19         contains
  20
  21
  22 !-----------------------------------------------------------------------
  23         subroutine mosaic_drydep_driver(                                   &
  24                 id, curr_secs, ktau, dtstep, config_flags,                     &
  25                 gmt, julday,                                                   &
  26                 t_phy, rho_phy, p_phy,                                         &
  27                 ust, aer_res,                                                  &
  28                 moist, chem, ddvel,                                            &
  29                 ids,ide, jds,jde, kds,kde,                                     &
  30                 ims,ime, jms,jme, kms,kme,                                     &
  31                 its,ite, jts,jte, kts,kte                                      )
  32
  33         use module_configure, only:  grid_config_rec_type, num_moist, num_chem
  34         use module_state_description, only:  param_first_scalar
  35
  36         use module_data_mosaic_asect
  37         use module_data_mosaic_other
  38         use module_mosaic_driver, only:  mapaer_tofrom_host
  39         use module_peg_util, only:  peg_error_fatal
  40
  41         implicit none
  42
  43 !   subr arguments
  44         integer, intent(in) ::   &
  45                 id, ktau, julday
  46
  47         integer, intent(in) ::   &
  48                 ids, ide, jds, jde, kds, kde,   &
  49                 ims, ime, jms, jme, kms, kme,   &
  50                 its, ite, jts, jte, kts, kte
  51
  52     real(kind=8), intent(in) :: curr_secs
  53
  54         real, intent(in) :: dtstep, gmt
  55
  56         real, intent(in),   &
  57                 dimension( ims:ime, kms:kme, jms:jme ) :: &
  58                 t_phy, rho_phy, p_phy
  59
  60         real, intent(in),   &
  61                 dimension( ims:ime, jms:jme ) :: &
  62                 ust
  63
  64         real, intent(in),   &
  65                 dimension( its:ite, jts:jte ) :: &
  66                 aer_res
  67
  68         real, intent(in),   &
  69                 dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
  70                 moist
  71
  72         real, intent(inout),   &
  73                 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
  74                 chem
  75
  76         real, intent(inout),   &
  77                 dimension( its:ite, jts:jte, 1:num_chem ) :: &
  78                 ddvel
  79
  80
  81         type(grid_config_rec_type), intent(in) :: config_flags
  82
  83
  84 !   local variables
  85         integer idum, jdum
  86         integer it, jt, kt
  87         integer iphase, itype
  88         integer ktmaps, ktmape
  89         integer ll, l1, n
  90         integer levdbg_err, levdbg_info
  91
  92         integer idiagaa_dum, ijcount_dum
  93
  94         real dum
  95         real vdep_aer(maxd_asize,maxd_atype,maxd_aphase)
  96
  97         character*100 msg
  98
  99
 100 !!rcetestdd diagnostics --------------------------------------------------
 101 !       print 93010, ' '
 102 !       print 93010, 'rcetestdd diagnostics from mosaic_drydep_driver'
 103 !       print 93010, 'id, chem_opt, ktau, julday           ',   &
 104 !                     id, config_flags%chem_opt, ktau, julday
 105 !       print 93010, 'ims/e, j, k', ims, ime, jms, jme, kms, kme
 106 !       print 93010, 'its/e, j, k', its, ite, jts, jte, kts, kte
 107 !
 108 !       do jdum = 0, 2
 109 !       do idum = 0, 2
 110 !       jt = jts + ((jte-jts)/2)*jdum
 111 !       it = its + ((ite-its)/2)*idum
 112 !if (idum .eq. 2) it = ite
 113 !if (jdum .eq. 2) jt = jte
 114 !       kt = kts
 115 !       print 93050, 'it, jt, t, p, ust, aer_res', it, jt,   &
 116 !               t_phy(it,kt,jt), p_phy(it,kt,jt), ust(it,jt), aer_res(it,jt)
 117 !       end do
 118 !       end do
 119 !
 120 !93010  format( a, 8(1x,i6) )
 121 !93020  format( a, 8(1p,e14.6) )
 122 !93050  format( a, 2(1x,i4), 8(1p,e14.6) )
 123 !!rcetestdd diagnostics --------------------------------------------------
 124
 125
 126 !   ktmaps,ktmape = first/last wrf kt for which depvel is calculated
 127 !   ktot = number of levels at which depvel is calculated
 128         ktmaps = kts
 129         ktmape = kts
 130         ktot = 1
 131
 132
 133 !   set some variables to their wrf-chem "standard" values
 134         lunerr = -1
 135         lunout = -1
 136         levdbg_err = 0
 137         levdbg_info = 15
 138
 139         iymdcur = 20030822
 140         ihmscur = 0
 141         ihmscur = nint( mod( real(gmt*3600.,8)+curr_secs, 86400.0_8 ) )
 142
 143         t = curr_secs
 144         ncorecnt = ktau - 1
 145
 146 !   reset some variables to "box test" values
 147 !   (*** aboxtest_get_extra_args is for "box testing" only
 148 !   and should be not be called from wrf-chem ***)
 149 !       call aboxtest_get_extra_args( 20,   &
 150 !               iymdcur, ihmscur,   &
 151 !               aboxtest_units_convert, aboxtest_rh_method,   &
 152 !               aboxtest_map_method, aboxtest_gases_fixed,  &
 153 !               lunerr, lunout, t, dum )
 154
 155
 156 !   set "pegasus" grid size variables
 157         itot = ite
 158         jtot = jte
 159         nsubareas = 1
 160
 161         ijcount_dum = 0
 162
 163
 164         do 2920 jt = jts, jte
 165         do 2910 it = its, ite
 166
 167         ijcount_dum = ijcount_dum + 1
 168
 169
 170         call mapaer_tofrom_host( 0,                       &
 171                 ims,ime, jms,jme, kms,kme,                    &
 172                 its,ite, jts,jte, kts,kte,                    &
 173                 it,      jt,      ktmaps,ktmape,              &
 174                 num_moist, num_chem, moist, chem,             &
 175                 t_phy, p_phy, rho_phy                         )
 176
 177
 178 !   compute deposition velocities
 179         idiagaa_dum = 1
 180         idiagaa_dum = 0
 181         if ((jt.ne.jts) .and. (jt.ne.jte) .and.   &
 182                         (jt.ne.(jts+jte)/2)) idiagaa_dum = 0
 183         if ((it.ne.its) .and. (it.ne.ite) .and.   &
 184                         (it.ne.(its+ite)/2)) idiagaa_dum = 0
 185
 186         call mosaic_drydep_1clm( idiagaa_dum, it, jt,   &
 187                 ust(it,jt), aer_res(it,jt), vdep_aer )
 188
 189
 190 ! rce 23-nov-2004 - change to using the "..._aer" species pointers
 191         do iphase = 1, nphase_aer
 192         do itype = 1, ntype_aer
 193         do n = 1, nsize_aer(itype)
 194             do ll = -2, ncomp_plustracer_aer(itype)
 195                 if (ll .eq. -2) then
 196                     l1 = numptr_aer(n,itype,iphase)
 197                 else if (ll .eq. -1) then
 198                     l1 = -1
 199                     if (iphase .eq. ai_phase) l1 = waterptr_aer(n,itype)
 200                 else if (ll .eq. 0) then
 201                     l1 = -1
 202                     if (iphase .eq. ai_phase) l1 = hyswptr_aer(n,itype)
 203                 else
 204                     l1 = massptr_aer(ll,n,itype,iphase)
 205                 end if
 206                 if (l1 .ge. param_first_scalar) then
 207                     ddvel(it,jt,l1) = vdep_aer(n,itype,iphase)
 208                 end if
 209             end do
 210         end do
 211         end do
 212         end do
 213
 214
 215 2910    continue
 216 2920    continue
 217 !       print 93010, 'leaving mosaic_drydep_driver'
 218
 219
 220         return
 221         end subroutine mosaic_drydep_driver
 222
 223
 224 !-----------------------------------------------------------------------
 225         subroutine mosaic_drydep_1clm( idiagaa, it, jt,   &
 226                 ustar_in, depresist_a_in, vdep_aer )
 227
 228         use module_configure, only:  grid_config_rec_type
 229
 230         use module_data_mosaic_asect
 231         use module_data_mosaic_other
 232         use module_mosaic_driver, only:  mapaer_tofrom_host
 233         use module_peg_util, only:  peg_error_fatal
 234
 235         implicit none
 236
 237 !   subr arguments
 238         integer, intent(in) :: idiagaa, it, jt
 239
 240 !   friction velocity (m/s)
 241         real, intent(in) :: ustar_in
 242 !   aerodynamic resistance (s/m)
 243         real, intent(in) :: depresist_a_in
 244
 245 !   deposition velocities (m/s)
 246         real, intent(inout) :: vdep_aer(maxd_asize,maxd_atype,maxd_aphase)
 247
 248 !   local variables
 249         real, parameter :: densdefault = 2.0
 250         real, parameter :: smallmassaa = 1.0e-20
 251         real, parameter :: smallmassbb = 1.0e-30
 252         real, parameter :: piover6 = pi/6.0
 253         real, parameter :: onethird = 1.0/3.0
 254
 255         integer iphase, itype, k, ll, l1, m, n
 256
 257         real airdens, airkinvisc
 258         real depresist_a, depresist_unstabpblfact
 259         real depresist_d0, depresist_d3
 260         real depvel_a0, depvel_a3
 261         real drydens, drydp, drymass, dryvol
 262         real dum, dumalnsg, dumfact, dummass
 263         real freepath
 264         real rnum
 265         real temp
 266         real ustar
 267         real vsettl_0, vsettl_3
 268         real wetdgnum, wetdens, wetdp, wetmass, wetvol
 269
 270
 271 !       if (idiagaa>0) print *, ' '
 272         k = 1
 273         m = 1
 274
 275 !   temperature (K)
 276         temp = rsub(ktemp,k,m)
 277 !   air density (g/cm^3)
 278 !       airdens = cairclm(1)*xmwair
 279         airdens = cairclm(1)*28.966
 280 !   air kinematic viscosity (cm^2/s)
 281         airkinvisc = ( 1.8325e-4 * (416.16/(temp+120.0)) *   &
 282                                 ((temp/296.16)**1.5) ) / airdens
 283 !   air molecular freepath (cm)
 284         freepath = 7.39758e-4 * airkinvisc / sqrt(temp)
 285 !   friction velocity (cm/s)
 286         ustar = ustar_in * 100.0
 287 !   aerodynamic resistance (s/cm)
 288         depresist_a = depresist_a_in * 0.01
 289
 290 !   enhancement factor for unstable pbl
 291         depresist_unstabpblfact = 1.0
 292
 293
 294 !   initialize vdep_aer
 295         vdep_aer(:,:,:) = 0.0
 296
 297 !   *** for now, just calc vdep_aer for iphase = ai_phase
 298         iphase = ai_phase
 299
 300 !   calculate vdep_aer
 301         do 2900 itype = 1, ntype_aer
 302         do 2800 n = 1, nsize_aer(itype)
 303
 304         dryvol = 0.0
 305         drymass = 0.0
 306         do ll = 1, ncomp_aer(itype)
 307             l1 = massptr_aer(ll,n,itype,iphase)
 308             dummass = rsub(l1,k,m)*mw_aer(ll,itype)
 309             drymass = drymass + dummass
 310             dryvol = dryvol + dummass/dens_aer(ll,itype)
 311         end do
 312
 313         l1 = waterptr_aer(n,itype)
 314         dummass = rsub(l1,k,m)*mw_water_aer
 315         wetmass = drymass + dummass
 316         wetvol = dryvol + dummass/dens_water_aer
 317
 318         l1 = numptr_aer(n,itype,iphase)
 319         rnum = rsub(l1,k,m)
 320
 321         if (drymass .le. smallmassbb) then
 322             drydp = dcen_sect(n,itype)
 323             drydens = densdefault
 324             wetdp = drydp
 325             wetdens = drydens
 326             goto 1900
 327         end if
 328
 329 !jdf    if (drymass .le. smallmassbb) then
 330         if (drymass .le. smallmassaa) then
 331             wetmass = drymass
 332             wetvol = dryvol
 333         end if
 334         drydens = drymass/dryvol
 335         wetdens = wetmass/wetvol
 336
 337
 338         if (rnum .ge. dryvol/volumlo_sect(n,itype)) then
 339             drydp = dlo_sect(n,itype)
 340         else if (rnum .le. dryvol/volumhi_sect(n,itype)) then
 341             drydp = dhi_sect(n,itype)
 342         else
 343             drydp = (dryvol/(piover6*rnum))**onethird
 344         end if
 345
 346 !jdf    dumfact = (wetvol/dryvol)**onethird
 347 !jdf    dumfact = min( dumfact, 10.0 )
 348         if(abs(wetvol).gt.(1000.*abs(dryvol))) then
 349           dumfact=10.0
 350         else
 351           dumfact=abs(wetvol/dryvol)**onethird
 352           dumfact=max(1.0,min(dumfact,10.0))
 353         endif
 354 !jdf
 355         wetdp = drydp*dumfact
 356
 357 1900    continue
 358
 359
 360 !
 361 !   get surface resistance and settling velocity for mass (moment 3)
 362 !   and number (moment 0)
 363 !
 364         dumalnsg = log( 1.0 )
 365         wetdgnum = wetdp * exp( -1.5*dumalnsg*dumalnsg )
 366         call aerosol_depvel_2(   &
 367             wetdgnum, dumalnsg, wetdens,   &
 368             temp, airdens, airkinvisc, freepath,   &
 369             ustar, depresist_unstabpblfact,   &
 370             depresist_d0, vsettl_0,   &
 371             depresist_d3, vsettl_3 )
 372
 373 !
 374 !   compute overall deposition velocity (binkowski/shankar eqn a33)
 375 !
 376         dum = depresist_a + depresist_d3 +   &
 377                     depresist_a*depresist_d3*vsettl_3
 378         depvel_a3 = vsettl_3 + (1./dum)
 379
 380         dum = depresist_a + depresist_d0 +   &
 381                     depresist_a*depresist_d0*vsettl_0
 382         depvel_a0 = vsettl_0 + (1./dum)
 383
 384 !   cm/s --> m/s
 385         vdep_aer(n,itype,iphase) = depvel_a3 * 0.01
 386
 387
 388 !
 389 !   diagnostic output
 390 !
 391         if (idiagaa>0) print 9310, it, jt, n, itype, iphase,   &
 392                 dcen_sect(n,itype), drydp, wetdp,   &
 393                 drydens, wetdens, vdep_aer(n,itype,iphase),   &
 394                 vsettl_3, depresist_d3, depresist_a
 395 9310    format( 'aerdep', 2i4, 3i3, 1p, 3e10.2,   &
 396                 2x, 0p, 2f5.2, 2x, 1p, 4e10.2 )
 397
 398
 399 2800    continue
 400 2900    continue
 401
 402
 403         return
 404         end subroutine mosaic_drydep_1clm
 405
 406
 407 !------------------------------------------------------------------------
 408 !------------------------------------------------------------------------
 409         subroutine aerosol_depvel_2(   &
 410             dgnum, alnsg, aerodens,   &
 411             temp, airdens, airkinvisc, freepath,   &
 412             ustar, depresist_unstabpblfact,   &
 413             depresist_d0, vsettl_0,   &
 414             depresist_d3, vsettl_3 )
 415 !
 416 !   computes the surface layer resistance term and the
 417 !   gravitational settling velocity term for the 3rd moment
 418 !   of a log-normal aerosol mode
 419 !
 420 !   input parameters
 421 !       dgnum - geometric mean diameter for aerosol number (cm)
 422 !       alnsg - natural logarithm of the geometric standard deviation
 423 !               for aerosol number
 424 !       aerodens - aerosol density (dgnum and aerodens are for the
 425 !               actual wet distribution)
 426 !       temp - temperature (K)
 427 !       airdens - air density (g/cm^3)
 428 !       airkinvisc - air kinematic viscosity (cm^2/s)
 429 !       freepath - air molecular freepath (cm)
 430 !       ustar - friction velocity (cm/s)
 431 !       depresist_unstabpblfact = weseley et al. 1985 factor for increasing
 432 !           depvel in unstable pbl -- either
 433 !               1. + (-0.3*zi/L)**0.667 OR
 434 !               1. + 0.24*((wstar/ustar)**2)
 435 !   output parameters
 436 !       depresist_d3 - surface layer resistance for 3rd (mass) moment (s/cm)
 437 !       vsettl_3 - gravitational settling velocity for 3rd moment (cm/s)
 438 !       depresist_d0 - surface layer resistance for 0th (number) moment (s/cm)
 439 !       vsettl_0 - gravitational settling velocity for 0th moment (cm/s)
 440 !
 441
 442         implicit none
 443
 444         real dgnum, alnsg, aerodens,   &
 445             temp, airdens, airkinvisc, freepath,   &
 446             ustar, depresist_unstabpblfact,   &
 447             depresist_d0, vsettl_0,   &
 448             depresist_d3, vsettl_3
 449
 450         real aerodiffus_0, schmidt_0, stokes_0, facdepresist_d0
 451         real aerodiffus_3, schmidt_3, stokes_3, facdepresist_d3
 452         common / aerosol_depvel_cmn01 /   &
 453                 aerodiffus_0, schmidt_0, stokes_0, facdepresist_d0,   &
 454                 aerodiffus_3, schmidt_3, stokes_3, facdepresist_d3
 455
 456         real xknudsen, xknudsenfact, alnsg2, duma, dumb,   &
 457         vsettl_dgnum, aerodiffus_dgnum
 458
 459         real pi
 460         parameter (pi = 3.1415926536)
 461 !   gravity = gravitational acceleration in cm/s^2
 462         real gravity
 463         parameter (gravity = 980.616)
 464 !   boltzmann constant in erg/deg-K
 465         real boltzmann
 466         parameter (boltzmann = 1.3807e-16)
 467
 468         xknudsen = 2.*freepath/dgnum
 469         xknudsenfact = xknudsen*1.246
 470         alnsg2 = alnsg*alnsg
 471
 472         vsettl_dgnum = (gravity*aerodens*dgnum*dgnum)/   &
 473                 (18.*airkinvisc*airdens)
 474         vsettl_0 = vsettl_dgnum *   &
 475                 ( exp(2.*alnsg2) + xknudsenfact*exp(0.5*alnsg2) )
 476         vsettl_3 = vsettl_dgnum *   &
 477                 ( exp(8.*alnsg2) + xknudsenfact*exp(3.5*alnsg2) )
 478
 479         aerodiffus_dgnum = (boltzmann*temp)/   &
 480                 (3.*pi*airkinvisc*airdens*dgnum)
 481         aerodiffus_0 = aerodiffus_dgnum *   &
 482                 ( exp(+0.5*alnsg2) + xknudsenfact*exp(+2.*alnsg2) )
 483         aerodiffus_3 = aerodiffus_dgnum *   &
 484                 ( exp(-2.5*alnsg2) + xknudsenfact*exp(-4.*alnsg2) )
 485
 486         schmidt_0 = airkinvisc/aerodiffus_0
 487         schmidt_3 = airkinvisc/aerodiffus_3
 488
 489         stokes_0 = ustar*ustar*vsettl_0/(gravity*airkinvisc)
 490         stokes_3 = ustar*ustar*vsettl_3/(gravity*airkinvisc)
 491
 492         duma = (schmidt_0**(-0.66666666)) +   &
 493                 (10.**(-3./max(0.03,stokes_0)))
 494 !       dumb = duma*ustar*(1. + 0.24*wstaroverustar*wstaroverustar)
 495         dumb = duma*ustar*depresist_unstabpblfact
 496         depresist_d0 = 1./max( dumb, 1.e-20 )
 497         facdepresist_d0 = duma
 498
 499         duma = (schmidt_3**(-0.66666666)) +   &
 500                 (10.**(-3./max(0.03,stokes_3)))
 501 !       dumb = duma*ustar*(1. + 0.24*wstaroverustar*wstaroverustar)
 502         dumb = duma*ustar*depresist_unstabpblfact
 503         depresist_d3 = 1./max( dumb, 1.e-20 )
 504         facdepresist_d3 = duma
 505
 506         return
 507         end subroutine aerosol_depvel_2
 508
 509
 510 !-----------------------------------------------------------------------
 511         end module module_mosaic_drydep