chem/module_aer_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
  10 MODULE module_aer_drydep
  11
  12 IMPLICIT NONE
  13
  14
  15 CONTAINS
  16
  17
  18 !-----------------------------------------------------------------------
  19         subroutine aer_drydep_driver(                                      &
  20                 id, ktau, dtstep, config_flags, aer_mech_id,               &
  21                 gmt, julday,                                               &
  22                 t_phy, rho_phy, p_phy,                                     &
  23                 alt, p8w, t8w, dz8w, z, z_at_w,                            &
  24                 ust, aer_res, ivgtyp, vegfra, pbl, rmol, znt,              &
  25                 moist, chem, ddvel,                                        &
  26                 h2oai, h2oaj, numgas,                                      &
  27                 ids,ide, jds,jde, kds,kde,                                 &
  28                 ims,ime, jms,jme, kms,kme,                                 &
  29                 its,ite, jts,jte, kts,kte                                  )
  30 !
  31 ! This routines allow different aerosol dry-deposition velocity parameterizations
  32 ! to be used with various MOSAIC and SORGAM chem_opt packages.
  33 !
  34 ! This routine is called when config_flags%aer_drydep_opt >= 100.
  35 ! However, values other than 101,111,201,301,302,311,312 cause an error.
  36 !
  37 ! For MOSAIC packages, calls are made to mosaic_aer_drydep_prep (to compute
  38 ! wet/ambient size and density), then to one of the dry-deposition velocity
  39 ! routines, then to mosaic_aer_drydep_load_ddvel.
  40 ! For SORGAM packages, calls are made to sorgam_aer_drydep_prep (to compute
  41 ! wet/ambient size and density), then to one of the dry-deposition velocity
  42 ! routines, then to sorgam_aer_drydep_load_ddvel.
  43 !
  44 ! config_flags%aer_drydep_opt == 101 -- call sorgam_aer_drydepvel_2, which
  45 !    uses the parameterization of subr vdvg_2 of module_aerosols_sorgam.F
  46 !
  47 ! config_flags%aer_drydep_opt == 111 -- call sorgam_aer_drydepvel_1, which
  48 !    uses the parameterization of subr vdvg of module_aerosols_sorgam.F
  49 !
  50 ! config_flags%aer_drydep_opt == 201 -- call mosaic_aer_drydepvel_1, which
  51 !    uses the parameterization of subr aerosol_depvel_2 of module_mosaic_drydep.F
  52 !
  53 ! config_flags%aer_drydep_opt == 301,302,311,312 -- call zhang2001_aer_drydepvel_1,
  54 !    which uses the parameterization of zhang et al. (2001, atmos. environ., p. 549-560).
  55 !    301 & 302 - uses the parameterization as described in the 2001 paper
  56 !    311 & 312 - has some modifications that increase the interception term (thus
  57 !        increasing depvel for particles in the accumulation size range)
  58 !        and reduce the brownian diffusion and impaction terms (thus reducing
  59 !        depvel for particles in the Aitken and coarse size ranges)
  60 !    301 & 311 - iseason is based on julian date, as in module_dep_simple.F
  61 !    302 & 312 - iseason is always 1 (summer)
  62 !
  63
  64         use module_configure, only:  grid_config_rec_type, num_moist, num_chem
  65         use module_state_description, only:  param_first_scalar
  66
  67         use module_data_mosaic_asect
  68         use module_data_mosaic_other
  69
  70         use modal_aero_data, only:  ntot_amode_cam_mam => ntot_amode
  71         implicit none
  72
  73 !   subr arguments
  74         integer, intent(in) ::   &
  75                 id, ktau, julday, aer_mech_id,   &
  76                 numgas
  77
  78         integer, intent(in) ::   &
  79                 ids, ide, jds, jde, kds, kde,   &
  80                 ims, ime, jms, jme, kms, kme,   &
  81                 its, ite, jts, jte, kts, kte
  82
  83         real, intent(in) :: dtstep, gmt
  84
  85         real, intent(in),   &
  86                 dimension( ims:ime, kms:kme, jms:jme ) :: &
  87                 t_phy, rho_phy, p_phy, &
  88                 alt, p8w, t8w, dz8w, z, z_at_w
  89
  90         integer, intent(in),   &
  91                 dimension( ims:ime, jms:jme ) :: &
  92                 ivgtyp
  93
  94         real, intent(in),   &
  95                 dimension( ims:ime, jms:jme ) :: &
  96                 ust, vegfra, pbl, rmol, znt
  97
  98         real, intent(in),   &
  99                 dimension( its:ite, jts:jte ) :: &
 100                 aer_res
 101
 102         real, intent(in),   &
 103                 dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
 104                 moist
 105
 106         real, intent(inout),   &
 107                 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
 108                 chem
 109
 110         real, intent(inout),   &
 111                 dimension( its:ite, jts:jte, 1:num_chem ) :: &
 112                 ddvel
 113
 114         real, intent(inout),   &
 115                 dimension( ims:ime, kms:kme, jms:jme ) :: &
 116                 h2oai, h2oaj
 117
 118         type(grid_config_rec_type), intent(in) :: config_flags
 119
 120
 121 !   local variables
 122         integer :: i, ibin, iflag_coarse, iflag_nac, &
 123                    iok_aer_drydep_opt, ioptaa, iphase, iseason, isurftype, itype, iwetsurf
 124         integer :: itmpa
 125         integer :: j
 126         integer :: n, nbin_add, numaer
 127
 128         real :: airdensbox, airkinvisc
 129         real :: depresist_a, depresist_d0, depresist_d3
 130         real :: depresist_unstabpblfact
 131         real :: freepath
 132         real :: presbox
 133         real :: tempbox
 134         real :: tmpa, tmp_lnsg, tmp_num
 135         real :: ustar
 136         real :: vsettl_0, vsettl_3
 137         real :: wetdgnum, wetdens, wstar
 138         real :: wetdgnum_si, wetdens_si
 139
 140         real, allocatable :: alnsg_add(:)
 141         real, allocatable :: depvel0_add(:), depvel3_add(:)
 142         real, allocatable :: wetdgnum_add(:), wetdens_add(:)
 143         real, allocatable :: yextra_add(:,:)
 144
 145         character(len=160) :: txtaa
 146
 147
 148         write(*,'(a,2(1x,i10))') &
 149                 '*** aer_drydep_driver - aer_mech_id, aer_drydep_opt = ', &
 150                 aer_mech_id, config_flags%aer_drydep_opt
 151
 152 !   allocate working arrays
 153         if ((aer_mech_id == 1) .or.   &
 154             (aer_mech_id == 2)) then
 155             nbin_add = 3
 156
 157         else if (aer_mech_id == 3) then
 158             nbin_add = sum( nsize_aer(1:ntype_aer) )
 159         else if (aer_mech_id == 4) then
 160             nbin_add = ntot_amode_cam_mam
 161         else
 162             write(txtaa,'(a,2(1x,i10))') &
 163                 '*** aer_drydep_driver - bad aer_mech_id, aer_drydep_opt = ', &
 164                 aer_mech_id, config_flags%aer_drydep_opt
 165             write(*,'(a)') trim(txtaa)
 166             call wrf_error_fatal( txtaa )
 167         end if
 168
 169         allocate( alnsg_add(    nbin_add ) )
 170         allocate( depvel0_add(  nbin_add ) )
 171         allocate( depvel3_add(  nbin_add ) )
 172         allocate( wetdgnum_add( nbin_add ) )
 173         allocate( wetdens_add(  nbin_add ) )
 174         allocate( yextra_add(   nbin_add, 10 ) )
 175
 176
 177 main_j_loop:   &
 178         do j = jts, jte
 179 main_i_loop:   &
 180         do i = its, ite
 181
 182
 183 !   get aerosol properties at lowest model layer
 184         if ((aer_mech_id == 1) .or.   &
 185             (aer_mech_id == 2)) then
 186             call sorgam_aer_drydep_prep(                                   &
 187                 id, ktau, dtstep, config_flags, aer_mech_id,               &
 188                 gmt, julday,                                               &
 189                 t_phy, rho_phy, p_phy, alt, p8w, t8w,                      &
 190                 moist, chem,                                               &
 191                 h2oai, h2oaj,                                              &
 192                 nbin_add, alnsg_add, wetdgnum_add, wetdens_add,            &
 193                 yextra_add,                                                &
 194                 i, j,                                                      &
 195                 ids,ide, jds,jde, kds,kde,                                 &
 196                 ims,ime, jms,jme, kms,kme,                                 &
 197                 its,ite, jts,jte, kts,kte                                  )
 198
 199         else if (aer_mech_id == 3) then
 200             call mosaic_aer_drydep_prep(                                   &
 201                 id, ktau, dtstep, config_flags, aer_mech_id,               &
 202                 gmt, julday,                                               &
 203                 t_phy, rho_phy, p_phy, alt, p8w, t8w,                      &
 204                 moist, chem,                                               &
 205                 nbin_add, alnsg_add, wetdgnum_add, wetdens_add,            &
 206                 yextra_add,                                                &
 207                 i, j,                                                      &
 208                 ids,ide, jds,jde, kds,kde,                                 &
 209                 ims,ime, jms,jme, kms,kme,                                 &
 210                 its,ite, jts,jte, kts,kte                                  )
 211        else if (aer_mech_id == 4) then
 212             call cam_mam_aer_drydep_prep(                                  &
 213                 id, ktau, dtstep, config_flags, aer_mech_id,               &
 214                 gmt, julday,                                               &
 215                 t_phy, rho_phy, p_phy, alt, p8w, t8w,                      &
 216                 moist, chem,                                               &
 217                 nbin_add, alnsg_add, wetdgnum_add, wetdens_add,            &
 218                 yextra_add,                                                &
 219                 i, j,                                                      &
 220                 ids,ide, jds,jde, kds,kde,                                 &
 221                 ims,ime, jms,jme, kms,kme,                                 &
 222                 its,ite, jts,jte, kts,kte                                  )
 223
 224         else
 225             txtaa = '*** aer_drydep_driver - fatal error 200'
 226             call wrf_error_fatal( txtaa )
 227         end if
 228
 229
 230 !   calculate deposition velocities
 231         do ibin = 1, nbin_add
 232
 233         ustar = ust(i,j)
 234         depresist_a = aer_res(i,j)
 235         tempbox = t_phy(i,kts,j)
 236         presbox = p_phy(i,kts,j)
 237
 238         wetdgnum = wetdgnum_add(ibin)
 239         wetdens  = wetdens_add(ibin)
 240         tmp_lnsg = alnsg_add(ibin)
 241
 242         iok_aer_drydep_opt = 0
 243
 244         if ((config_flags%aer_drydep_opt >= 100) .and.   &
 245             (config_flags%aer_drydep_opt <= 199)) then
 246             airkinvisc = 0.0
 247             freepath = 0.0
 248             ustar = max( 1.e-1, ust(i,j) )   ! this done in sorgam_depdriver
 249             tmpa = .001*p_phy(i,kts,j)
 250             presbox = 1.e3*tmpa
 251
 252             airdensbox = 1.0/alt(i,kts,j)
 253
 254             iflag_coarse = 0
 255             if (aer_mech_id <= 2) then
 256                 if (ibin == 3) iflag_coarse = 1
 257             else
 258                 ! yextra_add(ibin,9) = "reference" diameter (m)
 259                 if (yextra_add(ibin,9) >= 1.0e-6) iflag_coarse = 1
 260             end if
 261
 262             if (config_flags%aer_drydep_opt == 101) then
 263                 iok_aer_drydep_opt = 1
 264
 265                 call sorgam_aer_drydepvel_2(   &
 266                     iflag_coarse,   &
 267                     wetdgnum, tmp_lnsg, wetdens,   &
 268                     tempbox, presbox, airdensbox,   &
 269                     ustar, depresist_a,   &
 270                     pbl(i,j), znt(i,j), rmol(i,j),   &
 271                     airkinvisc, freepath,   &
 272                     depvel0_add(ibin), depvel3_add(ibin),   &
 273                     depresist_d0, depresist_d3,   &
 274                     vsettl_0, vsettl_3 )
 275
 276             else if (config_flags%aer_drydep_opt == 111) then
 277                 iok_aer_drydep_opt = 1
 278                 wstar = 0.0
 279
 280                 call sorgam_aer_drydepvel_1(   &
 281                     iflag_coarse,   &
 282                     wetdgnum, tmp_lnsg, wetdens,   &
 283                     tempbox, presbox, airdensbox,   &
 284                     ustar, wstar, depresist_a,   &
 285                     airkinvisc, freepath,   &
 286                     depvel0_add(ibin), depvel3_add(ibin),   &
 287                     depresist_d0, depresist_d3,   &
 288                     vsettl_0, vsettl_3 )
 289
 290             end if
 291
 292             depresist_unstabpblfact = 0.0
 293 !           depresist_d0 = 0.0
 294 !           depresist_d3 = 0.0
 295
 296         else if ((config_flags%aer_drydep_opt >= 200) .and.   &
 297                  (config_flags%aer_drydep_opt <= 299)) then
 298             airkinvisc = 0.0
 299             freepath = 0.0
 300             depresist_unstabpblfact = 0.0
 301
 302             if ((config_flags%aer_drydep_opt == 201) .or.   &
 303                 (config_flags%aer_drydep_opt == 211)) then
 304                 iok_aer_drydep_opt = 1
 305                 airdensbox = rho_phy(i,kts,j)
 306                 ioptaa = 1
 307                 if (config_flags%aer_drydep_opt == 211) ioptaa = 2
 308
 309                 call mosaic_aer_drydepvel_1(   &
 310                     ioptaa, wetdgnum, tmp_lnsg, wetdens,   &
 311                     tempbox, airdensbox, ustar, depresist_a,   &
 312                     airkinvisc, freepath, depresist_unstabpblfact,   &
 313                     depvel0_add(ibin), depvel3_add(ibin),   &
 314                     depresist_d0, depresist_d3,   &
 315                     vsettl_0, vsettl_3 )
 316
 317             end if
 318
 319         else if ((config_flags%aer_drydep_opt >= 300) .and.   &
 320                  (config_flags%aer_drydep_opt <= 399)) then
 321             airkinvisc = 0.0
 322             freepath = 0.0
 323             depresist_unstabpblfact = 0.0
 324
 325             if ((config_flags%aer_drydep_opt == 301) .or. &
 326                 (config_flags%aer_drydep_opt == 302) .or. &
 327                 (config_flags%aer_drydep_opt == 311) .or. &
 328                 (config_flags%aer_drydep_opt == 312)) then
 329                 iok_aer_drydep_opt = 1
 330                 airdensbox = rho_phy(i,kts,j)
 331                 isurftype = ivgtyp(i,j)
 332                 ! should set iwetsurf=1 when rainfall>0 or relhum>~95%
 333                 iwetsurf = 0
 334                 ! this is how iseason is set in module_dep_simple
 335                 ! could do a better job using leaf area index
 336                 iseason = 1
 337                 if (julday < 90 .or. julday > 270) iseason = 2
 338                 if ((config_flags%aer_drydep_opt == 302) .or. &
 339                     (config_flags%aer_drydep_opt == 312)) iseason = 1
 340                 ioptaa = 1
 341                 if ((config_flags%aer_drydep_opt == 311) .or. &
 342                     (config_flags%aer_drydep_opt == 312)) ioptaa = 2
 343
 344                 call zhang2001_aer_drydepvel_1(   &
 345                     wetdgnum, tmp_lnsg, wetdens,   &
 346                     tempbox, airdensbox, ustar, depresist_a,   &
 347                     ioptaa, iseason, isurftype, 2, iwetsurf,   &
 348                     airkinvisc, freepath,   &
 349                     1, 1,   &
 350                     depvel0_add(ibin), depvel3_add(ibin),   &
 351                     depresist_d0, depresist_d3,   &
 352                     vsettl_0, vsettl_3 )
 353
 354             end if
 355
 356         end if
 357
 358         if (iok_aer_drydep_opt <= 0) then
 359             write(txtaa,'(a,2(1x,i10))') &
 360                 '*** aer_drydep_driver - bad aer_mech_id, aer_drydep_opt = ', &
 361                 aer_mech_id, config_flags%aer_drydep_opt
 362             write(*,'(a)') trim(txtaa)
 363             call wrf_error_fatal( txtaa )
 364         end if
 365
 366
 367         end do   ! ibin
 368
 369
 370 !   load deposition velocities into ddvel array
 371         if ((aer_mech_id == 1) .or.   &
 372             (aer_mech_id == 2)) then
 373             call sorgam_aer_drydep_load_ddvel(                             &
 374                 id, ktau, dtstep, config_flags, aer_mech_id,               &
 375                 numgas, ddvel,                                             &
 376                 nbin_add, depvel0_add, depvel3_add,                        &
 377                 i, j,                                                      &
 378                 ids,ide, jds,jde, kds,kde,                                 &
 379                 ims,ime, jms,jme, kms,kme,                                 &
 380                 its,ite, jts,jte, kts,kte                                  )
 381
 382         else if (aer_mech_id == 3) then
 383             call mosaic_aer_drydep_load_ddvel(                             &
 384                 id, ktau, dtstep, config_flags, aer_mech_id,               &
 385                 ddvel,                                                     &
 386                 nbin_add, depvel0_add, depvel3_add,                        &
 387                 i, j,                                                      &
 388                 ids,ide, jds,jde, kds,kde,                                 &
 389                 ims,ime, jms,jme, kms,kme,                                 &
 390                 its,ite, jts,jte, kts,kte                                  )
 391         else if (aer_mech_id == 4) then
 392             call cam_mam_aer_drydep_load_ddvel(                            &
 393                 id, ktau, dtstep, config_flags, aer_mech_id,               &
 394                 ddvel,                                                     &
 395                 nbin_add, depvel0_add, depvel3_add,                        &
 396                 i, j,                                                      &
 397                 ids,ide, jds,jde, kds,kde,                                 &
 398                 ims,ime, jms,jme, kms,kme,                                 &
 399                 its,ite, jts,jte, kts,kte                                  )
 400
 401         else
 402             txtaa = '*** aer_drydep_driver - fatal error 400'
 403             call wrf_error_fatal( txtaa )
 404         end if
 405
 406
 407         end do main_i_loop
 408         end do main_j_loop
 409
 410
 411 !   deallocate working arrays
 412         deallocate( alnsg_add    )
 413         deallocate( depvel0_add  )
 414         deallocate( depvel3_add  )
 415         deallocate( wetdgnum_add )
 416         deallocate( wetdens_add  )
 417         deallocate( yextra_add   )
 418
 419
 420 80000   continue
 421         return
 422         end subroutine aer_drydep_driver
 423
 424
 425
 426 !-----------------------------------------------------------------------
 427         subroutine mosaic_aer_drydep_prep(                                 &
 428                 id, ktau, dtstep, config_flags, aer_mech_id,               &
 429                 gmt, julday,                                               &
 430                 t_phy, rho_phy, p_phy, alt, p8w, t8w,                      &
 431                 moist, chem,                                               &
 432                 nbin_add, alnsg_add, wetdgnum_add, wetdens_add,            &
 433                 yextra_add,                                                &
 434                 i, j,                                                      &
 435                 ids,ide, jds,jde, kds,kde,                                 &
 436                 ims,ime, jms,jme, kms,kme,                                 &
 437                 its,ite, jts,jte, kts,kte                                  )
 438
 439         use module_configure, only:  num_moist, num_chem, &
 440                 grid_config_rec_type
 441         use module_state_description, only:  param_first_scalar
 442
 443         use module_data_mosaic_asect
 444         use module_data_mosaic_other, only:  lunerr, lunout, pi
 445
 446         implicit none
 447
 448 !   subr arguments
 449         integer, intent(in) ::   &
 450                 id, ktau, julday, aer_mech_id, nbin_add,   &
 451                 i, j
 452
 453         integer, intent(in) ::   &
 454                 ids, ide, jds, jde, kds, kde,   &
 455                 ims, ime, jms, jme, kms, kme,   &
 456                 its, ite, jts, jte, kts, kte
 457
 458         real, intent(in) :: dtstep, gmt
 459
 460         real, intent(in),   &
 461                 dimension( ims:ime, kms:kme, jms:jme ) :: &
 462                 t_phy, rho_phy, p_phy, alt, p8w, t8w
 463
 464         real, intent(in),   &
 465                 dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
 466                 moist
 467
 468         real, intent(in),   &
 469                 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
 470                 chem
 471
 472         real, intent(inout),   &
 473                 dimension( 1:nbin_add ) :: &
 474                 alnsg_add, wetdgnum_add, wetdens_add
 475
 476         real, intent(inout),   &
 477                 dimension( 1:nbin_add, 10 ) :: &
 478                 yextra_add
 479
 480         type(grid_config_rec_type), intent(in) :: config_flags
 481
 482
 483 !   local variables
 484         integer :: ibin, ijcount, iok, iphase, itype
 485         integer :: itmpa
 486         integer :: k
 487         integer :: l, ll, l1, n
 488         integer :: p1st
 489
 490         real, parameter :: densdefault = 2.0
 491         real, parameter :: smallmassaa = 1.0e-20
 492         real, parameter :: smallmassbb = 1.0e-30
 493         real, parameter :: piover6 = pi/6.0
 494         real, parameter :: onethird = 1.0/3.0
 495
 496         real :: drydens, drydp, drymass, dryvol
 497         real :: rnum
 498         real :: tmpa, tmp_lnsg
 499         real :: wetdgnum, wetdens, wetdp, wetmass, wetvol
 500         real :: wetdgnum_si, wetdens_si
 501
 502         character(len=160) :: txtaa
 503
 504
 505         n = sum( nsize_aer(1:ntype_aer) )
 506         if (n /= nbin_add) then
 507             write(txtaa,'(a,2(1x,i10))') &
 508                 '*** mosaic_aer_drydep_prep - bad nbin_add, sum(nsize)', &
 509                 nbin_add, n
 510             write(*,'(a)') trim(txtaa)
 511             call wrf_error_fatal( txtaa )
 512             stop
 513         end if
 514
 515
 516 !   set some variables to their wrf-chem "standard" values
 517         p1st = param_first_scalar
 518         lunerr = -1
 519         lunout = -1
 520
 521         ijcount = 0
 522
 523         k = kts
 524
 525         ijcount = ijcount + 1
 526
 527 !   *** for now, just calc vdep_aer for iphase = ai_phase
 528         iphase = ai_phase
 529
 530 !   calculate vdep_aer for each type and size
 531         ibin = 0
 532         do 2900 itype = 1, ntype_aer
 533         do 2800 n = 1, nsize_aer(itype)
 534         ibin = ibin + 1
 535
 536 !   calculate dry and wet mass, volume, density, and mass-mean diameters
 537         dryvol = 0.0
 538         drymass = 0.0
 539         do ll = 1, ncomp_aer(itype)
 540             l1 = massptr_aer(ll,n,itype,iphase)
 541             tmpa = max( chem(i,k,j,l1), 0.0 )
 542             drymass = drymass + tmpa
 543             dryvol = dryvol + tmpa/dens_aer(ll,itype)
 544         end do
 545
 546         l1 = waterptr_aer(n,itype)
 547         tmpa = max( chem(i,k,j,l1), 0.0 )
 548         wetmass = drymass + tmpa
 549         wetvol = dryvol + tmpa/dens_water_aer
 550
 551         l1 = numptr_aer(n,itype,iphase)
 552         rnum = max( chem(i,k,j,l1), 0.0 )
 553
 554         drymass = drymass*28.966e-9 ! ug/kg-air --> g/mol-air
 555         wetmass = wetmass*28.966e-9
 556         dryvol  = dryvol *28.966e-9 ! "u-cm3"/kg-air --> cm3/mol-air
 557         wetvol  = wetvol *28.966e-9
 558         rnum    = rnum*28.966e-3    ! #/kg-air  --> #/mol-air
 559
 560         if (drymass <= smallmassbb) then
 561             drydp = dcen_sect(n,itype)
 562             drydens = densdefault
 563             wetdp = drydp
 564             wetdens = drydens
 565
 566         else
 567             if (drymass <= smallmassaa) then
 568                 wetmass = drymass
 569                 wetvol = dryvol
 570             end if
 571             drydens = drymass/dryvol
 572             wetdens = wetmass/wetvol
 573
 574             if (rnum >= dryvol/volumlo_sect(n,itype)) then
 575                 drydp = dlo_sect(n,itype)
 576             else if (rnum <= dryvol/volumhi_sect(n,itype)) then
 577                 drydp = dhi_sect(n,itype)
 578             else
 579                 drydp = (dryvol/(piover6*rnum))**onethird
 580             end if
 581
 582             if (abs(wetvol) > (1000.*abs(dryvol))) then
 583               tmpa=10.0
 584             else
 585               tmpa=abs(wetvol/dryvol)**onethird
 586               tmpa=max(1.0,min(tmpa,10.0))
 587             endif
 588             wetdp = drydp*tmpa
 589
 590         end if
 591
 592
 593 !   "within section" size distribution currently assumed to be
 594 !   monodisperse, but could assume otherwise
 595         tmp_lnsg = log( 1.0 )
 596         wetdgnum = wetdp * exp( -1.5*tmp_lnsg*tmp_lnsg )
 597
 598         wetdgnum_add(ibin) = wetdgnum * 1.0e-2 ! cm --> m
 599         wetdens_add(ibin) = wetdens * 1.0e3    ! g/cm3 --> kg/m3
 600         alnsg_add(ibin) = tmp_lnsg
 601
 602
 603         yextra_add(ibin,1) = drymass/28.966e-9 ! g/mol-air --> ug/kg-air = "n-kg"/kg-air
 604         yextra_add(ibin,5) = wetmass/28.966e-9
 605         yextra_add(ibin,2) = dryvol/28.966e-6  ! cm3/mol-air --> "n-m3"/kg-air
 606         yextra_add(ibin,6) = wetvol/28.966e-6
 607         yextra_add(ibin,3) = drydens*1.0e3     ! g/cm3 --> kg/m3
 608         yextra_add(ibin,7) = wetdens*1.0e3
 609         yextra_add(ibin,4) = drydp*1.0e-2      ! cm --> m
 610         yextra_add(ibin,8) = wetdp*1.0e-2
 611
 612         yextra_add(ibin,9) = dcen_sect(n,itype)*1.0e-2   ! cm --> m
 613         yextra_add(ibin,10) = rnum/28.966e-3    ! #/mol-air  --> #/kg-air
 614
 615
 616 2800    continue
 617 2900    continue
 618
 619
 620         return
 621         end subroutine mosaic_aer_drydep_prep
 622
 623
 624
 625 !-----------------------------------------------------------------------
 626         subroutine mosaic_aer_drydep_load_ddvel(                           &
 627                 id, ktau, dtstep, config_flags, aer_mech_id,               &
 628                 ddvel,                                                     &
 629                 nbin_add, depvel0_add, depvel3_add,                        &
 630                 i, j,                                                      &
 631                 ids,ide, jds,jde, kds,kde,                                 &
 632                 ims,ime, jms,jme, kms,kme,                                 &
 633                 its,ite, jts,jte, kts,kte                                  )
 634
 635         use module_configure, only:  num_moist, num_chem, &
 636                 grid_config_rec_type
 637         use module_state_description, only:  param_first_scalar
 638
 639         use module_data_mosaic_asect
 640
 641         implicit none
 642
 643 !   subr arguments
 644         integer, intent(in) ::   &
 645                 id, ktau, aer_mech_id, nbin_add,   &
 646                 i, j
 647
 648         integer, intent(in) ::   &
 649                 ids, ide, jds, jde, kds, kde,   &
 650                 ims, ime, jms, jme, kms, kme,   &
 651                 its, ite, jts, jte, kts, kte
 652
 653         real, intent(in) :: dtstep
 654
 655         real, intent(in),   &
 656                 dimension( 1:nbin_add ) :: &
 657                 depvel0_add, depvel3_add
 658
 659         real, intent(out),   &
 660                 dimension( its:ite, jts:jte, 1:num_chem ) :: &
 661                 ddvel
 662
 663         type(grid_config_rec_type), intent(in) :: config_flags
 664
 665
 666 !   local variables
 667         integer :: ibin, iphase, itype
 668         integer :: l, ll, n
 669
 670         character(len=160) :: txtaa
 671
 672
 673         n = sum( nsize_aer(1:ntype_aer) )
 674         if (n /= nbin_add) then
 675             write(txtaa,'(a,2(1x,i10))') &
 676                 '*** mosaic_aer_drydep_load_ddvel - bad nbin_add, sum(nsize)', &
 677                 nbin_add, n
 678             write(*,'(a)') trim(txtaa)
 679             call wrf_error_fatal( txtaa )
 680             stop
 681         end if
 682
 683
 684 ! copy from vdep_aer array to ddvel array
 685         do iphase = 1, nphase_aer
 686         ibin = 0
 687         do itype = 1, ntype_aer
 688         do n = 1, nsize_aer(itype)
 689             ibin = ibin + 1
 690             do ll = -2, ncomp_plustracer_aer(itype)
 691                 if (ll == -2) then
 692                     l = numptr_aer(n,itype,iphase)
 693                 else if (ll == -1) then
 694                     l = -1
 695                     if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
 696                 else if (ll == 0) then
 697                     l = -1
 698                     if (iphase == ai_phase) l = hyswptr_aer(n,itype)
 699                 else
 700                     l = massptr_aer(ll,n,itype,iphase)
 701                 end if
 702                 if (iphase > 1) l = -1
 703
 704                 if ((l >= param_first_scalar) .and. (l <= num_chem)) then
 705                     ddvel(i,j,l) = depvel3_add(ibin)
 706                 end if
 707             end do
 708         end do
 709         end do
 710         end do
 711
 712
 713         return
 714         end subroutine mosaic_aer_drydep_load_ddvel
 715
 716
 717
 718 !------------------------------------------------------------------------
 719         subroutine mosaic_aer_drydepvel_1(   &
 720             ioptaa, dgnum, alnsg, aerodens,   &
 721             temp, airdens, ustar, depresist_a,   &
 722             airkinvisc, freepath, depresist_unstabpblfact,   &
 723             depvel_0, depvel_3,   &
 724             depresist_d0, depresist_d3,   &
 725             vsettl_0, vsettl_3 )
 726 !
 727 !   computes the surface layer resistance term and the
 728 !   gravitational settling velocity term for the 3rd moment
 729 !   of a log-normal aerosol mode
 730 !
 731 !   input parameters
 732 !       ioptaa = 2 - apply a sticking fraction for larger particles
 733 !              = other - sticking fraction = 1 always
 734 !       dgnum - geometric mean diameter for aerosol number (m)
 735 !       alnsg - natural logarithm of the geometric standard deviation
 736 !               for aerosol number
 737 !       aerodens - aerosol density (dgnum and aerodens are for the
 738 !               actual wet distribution)
 739 !       temp - temperature (K)
 740 !       airdens - air density (kg/m^3)
 741 !       airkinvisc - air kinematic viscosity (m^2/s)
 742 !       freepath - air molecular freepath (m)
 743 !       ustar - friction velocity (m/s)
 744 !       depresist_a = aerodynamic resistance (s/m)
 745 !       depresist_unstabpblfact = weseley et al. 1985 factor for increasing
 746 !           depvel in unstable pbl -- either
 747 !               1. + (-0.3*zi/L)**0.667 OR
 748 !               1. + 0.24*((wstar/ustar)**2)
 749 !   output parameters
 750 !       depvel_0 - deposition velocity for 0th (number) moment (m/s)
 751 !       depvel_3 - deposition velocity for 3rd (mass  ) moment (m/s)
 752 !       depresist_d0/_d3 - surface layer resistance for 0th/3rd moment (s/m)
 753 !       vsettl_0/_3 - gravitational settling velocity for 0th/3rd moment (m/s)
 754 !
 755
 756         implicit none
 757
 758         integer, intent(in) ::   &
 759             ioptaa
 760         real, intent(in) ::   &
 761             dgnum, alnsg, aerodens,   &
 762             temp, airdens, ustar, depresist_a
 763         real, intent(inout) ::   &
 764             airkinvisc, freepath, depresist_unstabpblfact
 765         real, intent(out) ::   &
 766             depvel_0, depvel_3,   &
 767             depresist_d0, depresist_d3,   &
 768             vsettl_0, vsettl_3
 769
 770         real :: aerodiffus_0, schmidt_0, stokes_0, facdepresist_d0
 771         real :: aerodiffus_3, schmidt_3, stokes_3, facdepresist_d3
 772         common / aerosol_depvel_cmn01 /   &
 773                 aerodiffus_0, schmidt_0, stokes_0, facdepresist_d0,   &
 774                 aerodiffus_3, schmidt_3, stokes_3, facdepresist_d3
 775
 776         real :: aerodiffus_dgnum, alnsg2,   &
 777             stickfrac,   &
 778             tmpa, tmpb, tmp_dvm,   &
 779             vsettl_dgnum, xknudsen, xknudsenfact
 780
 781         real, parameter :: pi = 3.1415926536
 782 !   gravity = gravitational acceleration in m/s^2
 783         real, parameter :: gravity = 9.80616
 784 !   boltzmann constant in erg/deg-K
 785         real, parameter :: boltzmann = 1.3807e-23
 786
 787
 788         if (airkinvisc <= 0.0) then
 789 !   air kinematic viscosity (m^2/s)
 790             airkinvisc = ( 1.8325e-5 * (416.16/(temp+120.0)) *   &
 791                                 ((temp/296.16)**1.5) ) / airdens
 792         end if
 793         if (freepath <= 0.0) then
 794 !   air molecular freepath (m)
 795             freepath = 7.39758e-2 * airkinvisc / sqrt(temp)
 796         end if
 797         if (depresist_unstabpblfact <= 0) then
 798 !   enhancement factor for unstable pbl
 799             depresist_unstabpblfact = 1.0
 800         end if
 801
 802         xknudsen = 2.*freepath/dgnum
 803         xknudsenfact = xknudsen*1.246
 804         alnsg2 = alnsg*alnsg
 805
 806         vsettl_dgnum = (gravity*aerodens*dgnum*dgnum)/   &
 807                 (18.*airkinvisc*airdens)
 808         vsettl_0 = vsettl_dgnum *   &
 809                 ( exp(2.*alnsg2) + xknudsenfact*exp(0.5*alnsg2) )
 810         vsettl_3 = vsettl_dgnum *   &
 811                 ( exp(8.*alnsg2) + xknudsenfact*exp(3.5*alnsg2) )
 812
 813         aerodiffus_dgnum = (boltzmann*temp)/   &
 814                 (3.*pi*airkinvisc*airdens*dgnum)
 815         aerodiffus_0 = aerodiffus_dgnum *   &
 816                 ( exp(+0.5*alnsg2) + xknudsenfact*exp(+2.*alnsg2) )
 817         aerodiffus_3 = aerodiffus_dgnum *   &
 818                 ( exp(-2.5*alnsg2) + xknudsenfact*exp(-4.*alnsg2) )
 819
 820         schmidt_0 = airkinvisc/aerodiffus_0
 821         schmidt_3 = airkinvisc/aerodiffus_3
 822
 823         stokes_0 = ustar*ustar*vsettl_0/(gravity*airkinvisc)
 824         stokes_3 = ustar*ustar*vsettl_3/(gravity*airkinvisc)
 825
 826         tmp_dvm = dgnum * exp(1.5*alnsg2) ! volume mean diameter
 827
 828         tmpa = (schmidt_0**(-0.66666666)) +   &
 829                 (10.**(-3./max(0.03,stokes_0)))
 830 !       tmpb = tmpa*ustar*(1. + 0.24*wstaroverustar*wstaroverustar)
 831         tmpb = tmpa*ustar*depresist_unstabpblfact
 832         if (ioptaa == 2) then
 833             if (tmp_dvm >= 5.0e-6) then
 834                 stickfrac = exp( -2.0*sqrt(max(0.0,stokes_0)) )
 835                 tmpb = tmpb*stickfrac
 836             end if
 837         end if
 838         depresist_d0 = 1./max( tmpb, 1.e-22 )
 839         facdepresist_d0 = tmpa
 840
 841         tmpa = (schmidt_3**(-0.66666666)) +   &
 842                 (10.**(-3./max(0.03,stokes_3)))
 843 !       tmpb = tmpa*ustar*(1. + 0.24*wstaroverustar*wstaroverustar)
 844         tmpb = tmpa*ustar*depresist_unstabpblfact
 845         if (ioptaa == 2) then
 846             if (tmp_dvm >= 5.0e-6) then
 847                 stickfrac = exp( -2.0*sqrt(max(0.0,stokes_3)) )
 848                 tmpb = tmpb*stickfrac
 849             end if
 850         end if
 851         depresist_d3 = 1./max( tmpb, 1.e-22 )
 852         facdepresist_d3 = tmpa
 853
 854 !   compute overall deposition velocity (binkowski/shankar eqn a33)
 855         tmpa = depresist_a + depresist_d3 +   &
 856                 depresist_a*depresist_d3*vsettl_3
 857         depvel_3 = vsettl_3 + (1./tmpa)
 858
 859         tmpa = depresist_a + depresist_d0 +   &
 860                 depresist_a*depresist_d0*vsettl_0
 861         depvel_0 = vsettl_0 + (1./tmpa)
 862
 863         return
 864         end subroutine mosaic_aer_drydepvel_1
 865
 866
 867
 868 !-----------------------------------------------------------------------
 869         subroutine sorgam_aer_drydep_prep(                                 &
 870                 id, ktau, dtstep, config_flags, aer_mech_id,               &
 871                 gmt, julday,                                               &
 872                 t_phy, rho_phy, p_phy, alt, p8w, t8w,                      &
 873                 moist, chem,                                               &
 874                 h2oai, h2oaj,                                              &
 875                 nbin_add, alnsg_add, wetdgnum_add, wetdens_add,            &
 876                 yextra_add,                                                &
 877                 i, j,                                                      &
 878                 ids,ide, jds,jde, kds,kde,                                 &
 879                 ims,ime, jms,jme, kms,kme,                                 &
 880                 its,ite, jts,jte, kts,kte                                  )
 881
 882         use module_configure, only:  num_moist, num_chem, &
 883                 grid_config_rec_type
 884         use module_state_description
 885
 886         use module_data_sorgam
 887         use module_aerosols_sorgam, only:  modpar
 888
 889         implicit none
 890
 891 !   subr arguments
 892         integer, intent(in) ::   &
 893                 id, ktau, julday, aer_mech_id,   &
 894                 nbin_add, i, j
 895
 896         integer, intent(in) ::   &
 897                 ids, ide, jds, jde, kds, kde,   &
 898                 ims, ime, jms, jme, kms, kme,   &
 899                 its, ite, jts, jte, kts, kte
 900
 901         real, intent(in) :: dtstep, gmt
 902
 903         real, intent(in),   &
 904                 dimension( ims:ime, kms:kme, jms:jme ) :: &
 905                 t_phy, rho_phy, p_phy, alt, p8w, t8w
 906
 907         real, intent(in),   &
 908                 dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
 909                 moist
 910
 911         real, intent(in),   &
 912                 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
 913                 chem
 914
 915         real, intent(in),   &
 916                 dimension( ims:ime, kms:kme, jms:jme ) :: &
 917                 h2oai, h2oaj
 918
 919         real, intent(inout),   &
 920                 dimension( 1:nbin_add ) :: &
 921                 alnsg_add, wetdgnum_add, wetdens_add
 922
 923         real, intent(inout),   &
 924                 dimension( 1:nbin_add, 10 ) :: &
 925                 yextra_add
 926
 927         type(grid_config_rec_type), intent(in) :: config_flags
 928
 929
 930 !
 931 ! local variables and parameters
 932 !
 933
 934 !...blksize set to one in column model ciarev02
 935       integer, parameter :: blksize=1
 936
 937 !...number of aerosol species
 938 !  number of species (gas + aerosol)
 939       integer, parameter :: nspcsda=l1ae  !bs
 940
 941       integer, parameter :: numcells=1
 942       ! actual number of cells in arrays ( default is 1 in box model)
 943       !                   kept to 1 in current version of column model
 944
 945 !bs total # of cond. vapors in ctm
 946       real cblk(blksize,nspcsda) ! main array of variables
 947                                    ! particles [ug/m^3/s]
 948
 949       real blkprs(blksize)         ! pressure in cb
 950       real blkta(blksize)          ! temperature in k
 951       real blkdens(blksize)        ! air density in kg/m3
 952
 953 ! *** output:
 954 ! *** atmospheric properties
 955       real xlm( blksize )           ! atmospheric mean free path [ m ]
 956       real amu( blksize )           ! atmospheric dynamic viscosity [ kg/m s ]
 957
 958 ! *** followng is for future version
 959       real vsed( blksize, naspcssed) ! sedimentation velocity [ m s**-1 ]
 960       real vdep( blksize, naspcsdep) ! deposition velocity [ m s**-1 ]
 961
 962 ! *** modal diameters: [ m ]
 963       real dgnuc( blksize )         ! nuclei mode geometric mean diameter  [ m ]
 964       real dgacc( blksize )         ! accumulation geometric mean diameter [ m ]
 965       real dgcor( blksize )         ! coarse mode geometric mean diameter  [ m ]
 966
 967 ! *** aerosol properties:
 968 ! *** modal mass concentrations [ ug m**3 ]
 969       real pmassn( blksize )        ! mass concentration in aitken mode
 970       real pmassa( blksize )        ! mass concentration in accumulation mode
 971       real pmassc( blksize )        ! mass concentration in coarse mode
 972
 973 ! *** average modal particle densities  [ kg/m**3 ]
 974       real pdensn( blksize )        ! average particle density in nuclei mode
 975       real pdensa( blksize )        ! average particle density in accumulation mode
 976       real pdensc( blksize )        ! average particle density in coarse mode
 977
 978 ! *** average modal knudsen numbers
 979       real knnuc ( blksize )        ! nuclei mode  knudsen number
 980       real knacc ( blksize )        ! accumulation knudsen number
 981       real kncor ( blksize )        ! coarse mode  knudsen number
 982
 983       integer :: k,l
 984
 985       real :: convfac2
 986       real, dimension( ims:ime, jms:jme ) :: aer_res, ust   ! TEMPORARY
 987       real, dimension( kts:kte ) :: p   ! TEMPORARY
 988
 989       character(len=160) :: txtaa
 990
 991
 992       if (3 /= nbin_add) then
 993           write(txtaa,'(a,2(1x,i10))') &
 994             '*** sorgam_aer_drydep_prep - bad nbin_add', &
 995             nbin_add
 996           write(*,'(a)') trim(txtaa)
 997           call wrf_error_fatal( txtaa )
 998           stop
 999       end if
1000
1001
1002 !     print *,'in sorgdepdriver ',its,ite,jts,jte
1003
1004       cblk=epsilc
1005
1006 !     do jj = jts, jte
1007 !     do ii = its, ite
1008 !        do k = kts, kte-1
1009 !           p(k)  = .001*p_phy(ii,k,jj)
1010 !        end do
1011 !     end do
1012 !     end do
1013
1014       k=kts
1015
1016       p(k)  = .001*p_phy(i,k,j)
1017       convfac2 = 1./alt(i,k,j)
1018       blkdens(blksize) = convfac2
1019       blkta(blksize)   = t_phy(i,k,j)
1020       blkprs(blksize)  = 1.e3*p(k)
1021
1022
1023       cblk(1,vso4aj   ) =   max( epsilc, chem(i,k,j,p_so4aj)*convfac2 )
1024       cblk(1,vso4ai   ) =   max( epsilc, chem(i,k,j,p_so4ai)*convfac2 )
1025       cblk(1,vnh4aj   ) =   max( epsilc, chem(i,k,j,p_nh4aj)*convfac2 )
1026       cblk(1,vnh4ai   ) =   max( epsilc, chem(i,k,j,p_nh4ai)*convfac2 )
1027       cblk(1,vno3aj   ) =   max( epsilc, chem(i,k,j,p_no3aj)*convfac2 )
1028       cblk(1,vno3ai   ) =   max( epsilc, chem(i,k,j,p_no3ai)*convfac2 )
1029       if (p_naaj >= param_first_scalar) &
1030          cblk(1,vnaaj ) =   max( epsilc, chem(i,k,j,p_naaj)*convfac2 )
1031       if (p_naai >= param_first_scalar) &
1032          cblk(1,vnaai ) =   max( epsilc, chem(i,k,j,p_naai)*convfac2 )
1033       if (p_claj >= param_first_scalar) &
1034          cblk(1,vclaj ) =   max( epsilc, chem(i,k,j,p_claj)*convfac2 )
1035       if (p_clai >= param_first_scalar) &
1036          cblk(1,vclai ) =   max( epsilc, chem(i,k,j,p_clai)*convfac2 )
1037       cblk(1,vorgaro1j) =   max( epsilc, chem(i,k,j,p_orgaro1j)*convfac2 )
1038       cblk(1,vorgaro1i) =   max( epsilc, chem(i,k,j,p_orgaro1i)*convfac2 )
1039       cblk(1,vorgaro2j) =   max( epsilc, chem(i,k,j,p_orgaro2j)*convfac2 )
1040       cblk(1,vorgaro2i) =   max( epsilc, chem(i,k,j,p_orgaro2i)*convfac2 )
1041       cblk(1,vorgalk1j) =   max( epsilc, chem(i,k,j,p_orgalk1j)*convfac2 )
1042       cblk(1,vorgalk1i) =   max( epsilc, chem(i,k,j,p_orgalk1i)*convfac2 )
1043       cblk(1,vorgole1j) =   max( epsilc, chem(i,k,j,p_orgole1j)*convfac2 )
1044       cblk(1,vorgole1i) =   max( epsilc, chem(i,k,j,p_orgole1i)*convfac2 )
1045       cblk(1,vorgba1j ) =   max( epsilc, chem(i,k,j,p_orgba1j)*convfac2 )
1046       cblk(1,vorgba1i ) =   max( epsilc, chem(i,k,j,p_orgba1i)*convfac2 )
1047       cblk(1,vorgba2j ) =   max( epsilc, chem(i,k,j,p_orgba2j)*convfac2 )
1048       cblk(1,vorgba2i ) =   max( epsilc, chem(i,k,j,p_orgba2i)*convfac2 )
1049       cblk(1,vorgba3j ) =   max( epsilc, chem(i,k,j,p_orgba3j)*convfac2 )
1050       cblk(1,vorgba3i ) =   max( epsilc, chem(i,k,j,p_orgba3i)*convfac2 )
1051       cblk(1,vorgba4j ) =   max( epsilc, chem(i,k,j,p_orgba4j)*convfac2 )
1052       cblk(1,vorgba4i ) =   max( epsilc, chem(i,k,j,p_orgba4i)*convfac2 )
1053       cblk(1,vorgpaj  ) =   max( epsilc, chem(i,k,j,p_orgpaj)*convfac2 )
1054       cblk(1,vorgpai  ) =   max( epsilc, chem(i,k,j,p_orgpai)*convfac2 )
1055       cblk(1,vecj     ) =   max( epsilc, chem(i,k,j,p_ecj)*convfac2 )
1056       cblk(1,veci     ) =   max( epsilc, chem(i,k,j,p_eci)*convfac2 )
1057       cblk(1,vp25aj   ) =   max( epsilc, chem(i,k,j,p_p25j)*convfac2 )
1058       cblk(1,vp25ai   ) =   max( epsilc, chem(i,k,j,p_p25i)*convfac2 )
1059       cblk(1,vantha   ) =   max( epsilc, chem(i,k,j,p_antha)*convfac2 )
1060       cblk(1,vseas    ) =   max( epsilc, chem(i,k,j,p_seas)*convfac2 )
1061       cblk(1,vsoila   ) =   max( epsilc, chem(i,k,j,p_soila)*convfac2 )
1062       cblk(1,vnu0     ) =   max( epsilc, chem(i,k,j,p_nu0)*convfac2 )
1063       cblk(1,vac0     ) =   max( epsilc, chem(i,k,j,p_ac0)*convfac2 )
1064       cblk(1,vcorn    ) =   max( epsilc, chem(i,k,j,p_corn)*convfac2 )
1065
1066       cblk(1,vh2oaj   ) =   h2oaj(i,k,j)
1067       cblk(1,vh2oai   ) =   h2oai(i,k,j)
1068
1069
1070       call modpar(  blksize, nspcsda, numcells, &
1071            cblk,                                &
1072            blkta, blkprs,                       &
1073            pmassn, pmassa, pmassc,              &
1074            pdensn, pdensa, pdensc,              &
1075            xlm, amu,                            &
1076            dgnuc, dgacc, dgcor,                 &
1077            knnuc, knacc, kncor                  )
1078
1079       wetdgnum_add(1) = dgnuc(1)
1080       wetdgnum_add(2) = dgacc(1)
1081       wetdgnum_add(3) = dgcor(1)
1082
1083       wetdens_add(1) = pdensn(1)
1084       wetdens_add(2) = pdensa(1)
1085       wetdens_add(3) = pdensc(1)
1086
1087       alnsg_add(1) = xxlsgn
1088       alnsg_add(2) = xxlsga
1089       alnsg_add(3) = xxlsgc
1090
1091       yextra_add(1:3,1:8) = 0.0
1092
1093       yextra_add(1,5) = pmassn(1)/blkdens(1)   ! ug-AP/kg-air = "n-kg-AP"/kg-air
1094       yextra_add(2,5) = pmassa(1)/blkdens(1)
1095       yextra_add(3,5) = pmassc(1)/blkdens(1)
1096
1097       yextra_add(1,6) = cblk(1,vnu3 )*(pirs/6.0)*1.0e9/blkdens(1)   ! "n-m3-AP"/kg-air
1098       yextra_add(2,6) = cblk(1,vac3 )*(pirs/6.0)*1.0e9/blkdens(1)
1099       yextra_add(3,6) = cblk(1,vcor3)*(pirs/6.0)*1.0e9/blkdens(1)
1100
1101       yextra_add(1,7) = pdensn(1)   ! kg-AP/m3-AP
1102       yextra_add(2,7) = pdensa(1)
1103       yextra_add(3,7) = pdensc(1)
1104
1105       yextra_add(1,8) = dgnuc(1)*exp(1.5*(alnsg_add(1)**2))   ! m
1106       yextra_add(2,8) = dgacc(1)*exp(1.5*(alnsg_add(2)**2))
1107       yextra_add(3,8) = dgcor(1)*exp(1.5*(alnsg_add(3)**2))
1108
1109       yextra_add(1,9) = dginin
1110       yextra_add(2,9) = dginia
1111       yextra_add(3,9) = dginic
1112
1113       yextra_add(1,10) = cblk(1,vnu0 )/blkdens(1)   ! #/kg-air
1114       yextra_add(2,10) = cblk(1,vac0 )/blkdens(1)   ! #/kg-air
1115       yextra_add(3,10) = cblk(1,vcorn)/blkdens(1)   ! #/kg-air
1116
1117         return
1118         end subroutine sorgam_aer_drydep_prep
1119
1120
1121
1122 !-----------------------------------------------------------------------
1123         subroutine sorgam_aer_drydep_load_ddvel(                           &
1124                 id, ktau, dtstep, config_flags, aer_mech_id,               &
1125                 numgas, ddvel,                                             &
1126                 nbin_add, depvel0_add, depvel3_add,                        &
1127                 i, j,                                                      &
1128                 ids,ide, jds,jde, kds,kde,                                 &
1129                 ims,ime, jms,jme, kms,kme,                                 &
1130                 its,ite, jts,jte, kts,kte                                  )
1131
1132         use module_configure, only:  num_moist, num_chem, &
1133                 grid_config_rec_type
1134         use module_state_description, only:  param_first_scalar, &
1135                 p_naai, p_naaj, p_clai, p_claj
1136
1137         use module_data_sorgam
1138
1139         implicit none
1140
1141 !   subr arguments
1142         integer, intent(in) ::   &
1143                 id, ktau, aer_mech_id, &
1144                 numgas, nbin_add, i, j
1145
1146         integer, intent(in) ::   &
1147                 ids, ide, jds, jde, kds, kde,   &
1148                 ims, ime, jms, jme, kms, kme,   &
1149                 its, ite, jts, jte, kts, kte
1150
1151         real, intent(in) :: dtstep
1152
1153         real, intent(in),   &
1154                 dimension( 1:nbin_add ) :: &
1155                 depvel0_add, depvel3_add
1156
1157         real, intent(out),   &
1158                 dimension( its:ite, jts:jte, 1:num_chem ) :: &
1159                 ddvel
1160
1161         type(grid_config_rec_type), intent(in) :: config_flags
1162
1163
1164 !   local variables
1165         integer :: l
1166
1167         real :: vgsa(num_chem)
1168
1169         character(len=160) :: txtaa
1170
1171
1172         if (3 /= nbin_add) then
1173             write(txtaa,'(a,2(1x,i10))') &
1174                 '*** sorgam_aer_drydep_load_ddvel - bad nbin_add', &
1175                 nbin_add
1176             write(*,'(a)') trim(txtaa)
1177             call wrf_error_fatal( txtaa )
1178             stop
1179         end if
1180
1181
1182         vgsa( : )  =  0.0
1183
1184 !       vgsa( vso4aj   )  =  vdep(1, vdmacc )
1185 !       vgsa( vso4ai   )  =  vdep(1, vdmnuc )
1186         vgsa( vso4aj   )  =  depvel3_add( 2 )
1187         vgsa( vso4ai   )  =  depvel3_add( 1 )
1188
1189         vgsa( vnh4aj   )  =  vgsa( vso4aj )
1190         vgsa( vnh4ai   )  =  vgsa( vso4ai )
1191         vgsa( vno3aj   )  =  vgsa( vso4aj )
1192         vgsa( vno3ai   )  =  vgsa( vso4ai )
1193         if (p_naaj >= param_first_scalar) vgsa( vnaaj   )  =  vgsa( vso4aj )
1194         if (p_naai >= param_first_scalar) vgsa( vnaai   )  =  vgsa( vso4ai )
1195         if (p_claj >= param_first_scalar) vgsa( vclaj   )  =  vgsa( vso4aj )
1196         if (p_clai >= param_first_scalar) vgsa( vclai   )  =  vgsa( vso4ai )
1197         vgsa( vorgaro1j)  =  vgsa( vso4aj )
1198         vgsa( vorgaro1i)  =  vgsa( vso4ai )
1199         vgsa( vorgaro2j)  =  vgsa( vso4aj )
1200         vgsa( vorgaro2i)  =  vgsa( vso4ai )
1201         vgsa( vorgalk1j)  =  vgsa( vso4aj )
1202         vgsa( vorgalk1i)  =  vgsa( vso4ai )
1203         vgsa( vorgole1j)  =  vgsa( vso4aj )
1204         vgsa( vorgole1i)  =  vgsa( vso4ai )
1205         vgsa( vorgba1j )  =  vgsa( vso4aj )
1206         vgsa( vorgba1i )  =  vgsa( vso4ai )
1207         vgsa( vorgba2j )  =  vgsa( vso4aj )
1208         vgsa( vorgba2i )  =  vgsa( vso4ai )
1209         vgsa( vorgba3j )  =  vgsa( vso4aj )
1210         vgsa( vorgba3i )  =  vgsa( vso4ai )
1211         vgsa( vorgba4j )  =  vgsa( vso4aj )
1212         vgsa( vorgba4i )  =  vgsa( vso4ai )
1213         vgsa( vorgpaj  )  =  vgsa( vso4aj )
1214         vgsa( vorgpai  )  =  vgsa( vso4ai )
1215         vgsa( vecj     )  =  vgsa( vso4aj )
1216         vgsa( veci     )  =  vgsa( vso4ai )
1217         vgsa( vp25aj   )  =  vgsa( vso4aj )
1218         vgsa( vp25ai   )  =  vgsa( vso4ai )
1219
1220 !       vgsa( vantha   )  =  vdep(1, vdmcor )
1221         vgsa( vantha   )  =  depvel3_add( 3 )
1222         vgsa( vseas    )  =  vgsa( vantha )
1223         vgsa( vsoila   )  =  vgsa( vantha )
1224
1225 !       vgsa( vnu0     )  =  vdep(1, vdnnuc )
1226 !       vgsa( vac0     )  =  vdep(1, vdnacc )
1227 !       vgsa( vcorn    )  =  vdep(1, vdncor )
1228         vgsa( vnu0     )  =  depvel0_add( 1 )
1229         vgsa( vac0     )  =  depvel0_add( 2 )
1230         vgsa( vcorn    )  =  depvel0_add( 3 )
1231
1232         do l = 1, (num_chem - numgas )
1233            ddvel(i, j, l+numgas ) =  vgsa( l )
1234         end do
1235
1236
1237         return
1238         end subroutine sorgam_aer_drydep_load_ddvel
1239
1240
1241
1242 !------------------------------------------------------------------------
1243         subroutine sorgam_aer_drydepvel_1(   &
1244             iflag_coarse,   &
1245             dgnum, alnsg, aerodens,   &
1246             temp, pres, airdens, ustar, wstar, depresist_a,   &
1247             airkinvisc, freepath,   &
1248             depvel_0, depvel_3,   &
1249             depresist_d0, depresist_d3,   &
1250             vsettl_0, vsettl_3 )
1251 !
1252 !   computes the surface layer resistance term and the
1253 !   gravitational settling velocity term for the 0th and 3rd moment
1254 !   of a log-normal aerosol mode
1255 !
1256 !   input parameters
1257 !       dgnum - geometric mean diameter for aerosol number (m)
1258 !       alnsg - natural logarithm of the geometric standard deviation
1259 !               for aerosol number
1260 !       aerodens - aerosol density (dgnum and aerodens are for the
1261 !               actual wet distribution)
1262 !       temp - temperature (K)
1263 !       pres - pressure (Pa)
1264 !       airdens - air density (kg/m^3)
1265 !       ustar - friction velocity (m/s)
1266 !       wstar - convective scale velocity (m/s)
1267 !       depresist_a = aerodynamic resistance (s/m)
1268 !       airkinvisc - air kinematic viscosity (m^2/s)
1269 !       freepath - air molecular freepath (m)
1270 !   output parameters
1271 !       depvel_0 - deposition velocity for 0th (number) moment (m/s)
1272 !       depvel_3 - deposition velocity for 3rd (mass  ) moment (m/s)
1273 !       depresist_d0/_d3 - surface layer resistance for 0th/3rd moment (s/m)
1274 !       vsettl_0/_3 - gravitational settling velocity for 0th/3rd moment (m/s)
1275 !
1276
1277 !       use module_data_sorgam, only: boltz, grav, pss0, tss0, threepi, two3
1278
1279         implicit none
1280
1281         integer, intent(in) ::   &
1282             iflag_coarse
1283         real, intent(in) ::   &
1284             dgnum, alnsg, aerodens,   &
1285             temp, pres, airdens, ustar, wstar, depresist_a
1286         real, intent(inout) ::   &
1287             airkinvisc, freepath
1288         real, intent(out) ::   &
1289             depvel_0, depvel_3,   &
1290             depresist_d0, depresist_d3,   &
1291             vsettl_0, vsettl_3
1292
1293
1294 ! following are from module_data_sorgam
1295 ! avogadro's constant [ 1/mol ]
1296         real, parameter :: avo=6.0221367e23
1297 ! universal gas constant [ j/mol-k ]
1298         real, parameter :: rgasuniv=8.314510
1299 ! boltzmann's constant [ j / k]
1300         real, parameter :: boltz=rgasuniv/avo
1301 ! mean gravitational acceleration [ m/sec**2 ]
1302         real, parameter :: grav=9.80622
1303 !  starting standard surface pressure [ pa ]
1304         real, parameter :: pss0=101325.0
1305 !  starting standard surface temperature [ k ]
1306         real, parameter :: tss0=288.15
1307 !  2/3
1308         real, parameter :: two3=2.0/3.0
1309 !  pi
1310         real*8, parameter :: pirs=3.14159265358979324
1311 !  3*pi
1312         real, parameter :: threepi=3.0*pirs
1313
1314
1315         real, parameter :: bhat = 1.246   ! constant from cunningham slip correction.
1316         real :: alnsgy2, amu
1317         real :: blkta, blkprs, blkdens
1318         real :: dconst1, dconst1y, dconst2, dconst3y
1319         real :: dchat0y, dchat3y
1320         real :: dgyyy
1321         real :: ey1, esy04, esy08, esy16, esy20, esy28, esy32, esy64
1322         real :: esym20, esym32
1323         real :: knyyy
1324         real :: nu
1325         real :: pdensy
1326         real :: ra, rd0y, rd3y
1327         real :: sc0y, sc3y, st0y, st3y
1328         real :: ustfac, utscale
1329         real :: vghat0y, vghat3y
1330         real :: xlm
1331         real :: yextra_vdvg(100,3)
1332
1333
1334         blkta = temp
1335         blkprs = pres
1336         blkdens = airdens
1337
1338 !   amu = air dynamic viscosity (kg/m/s)
1339         if (airkinvisc <= 0.0) then
1340             amu = 1.458e-6*blkta*sqrt(blkta)/(blkta+110.4)
1341             airkinvisc = amu/airdens
1342         else
1343             amu = airkinvisc*airdens
1344         end if
1345
1346 !   air molecular freepath (m)
1347         if (freepath <= 0.0) then
1348             freepath = 6.6328e-8*pss0*blkta/(tss0*blkprs)
1349         end if
1350         xlm = freepath
1351
1352         ra = depresist_a
1353
1354         dgyyy = dgnum
1355         pdensy = aerodens
1356
1357         alnsgy2 = alnsg**2
1358         ey1 = exp(0.125*alnsgy2)
1359         esy04 = ey1**4
1360         esy08 = esy04*esy04
1361 !       esy12 = esy04*esy04*esy04
1362         esy16 = esy08*esy08
1363         esy20 = esy16*esy04
1364 !       esy24 = esy12*esy12
1365         esy28 = esy20*esy08
1366         esy32 = esy16*esy16
1367 !       esy36 = esy16*esy20
1368         esy64 = esy32*esy32
1369         esym20 = 1.0/esy20
1370         esym32 = 1.0/esy32
1371
1372         knyyy = 2.0*xlm/dgyyy
1373
1374         dconst1 = boltz * blkta / ( threepi * amu )
1375         dconst1y = dconst1 / dgyyy
1376
1377         dconst2 = grav / ( 18.0 * amu )
1378         dconst3y = dconst2 * pdensy * dgyyy**2
1379
1380         dchat0y = dconst1y * ( esy04  + bhat * knyyy * esy16 )
1381         dchat3y = dconst1y * ( esym20 + bhat * knyyy * esym32 )
1382         vghat0y = dconst3y * ( esy16  + bhat * knyyy * esy04 )
1383         vghat3y = dconst3y * ( esy64  + bhat * knyyy * esy28 )
1384
1385         nu = amu / blkdens
1386         ustfac = ustar * ustar / ( grav * nu)
1387         utscale = ustar + 0.24 * wstar * wstar / ustar
1388
1389         sc0y = nu / dchat0y
1390         if (iflag_coarse > 0) then
1391             rd0y = 1.0 / ( utscale * ( sc0y**(-two3) ) )
1392         else
1393             st0y = max( vghat0y * ustfac , 0.01)
1394             rd0y = 1.0 / ( utscale * ( sc0y**(-two3) + 10.0**(-3.0 / st0y) ) )
1395         end if
1396
1397         depvel_0 = vghat0y + 1.0 / ( ra + rd0y + rd0y * ra * vghat0y )
1398         vsettl_0 = vghat0y
1399         depresist_d0 = rd0y
1400
1401
1402         sc3y = nu / dchat3y
1403         if (iflag_coarse > 0) then
1404             rd3y = 1.0 / ( utscale * ( sc3y**(-two3) ) )
1405         else
1406             st3y = max( vghat3y * ustfac , 0.01)
1407             rd3y = 1.0 / ( utscale * ( sc3y**(-two3) + 10.0**(-3.0 / st3y) ) )
1408         end if
1409
1410         depvel_3 = vghat3y + 1.0 / ( ra + rd3y + rd3y * ra * vghat3y )
1411         vsettl_3 = vghat3y
1412         depresist_d3 = rd3y
1413
1414
1415         return
1416         end subroutine sorgam_aer_drydepvel_1
1417
1418
1419
1420 !------------------------------------------------------------------------
1421         subroutine sorgam_aer_drydepvel_2(   &
1422             iflag_coarse,   &
1423             dgnum, alnsg, aerodens,   &
1424             temp, pres, airdens, ustar, depresist_a,   &
1425             pblh, zntt, rmolm,   &
1426             airkinvisc, freepath,   &
1427             depvel_0, depvel_3,   &
1428             depresist_d0, depresist_d3,   &
1429             vsettl_0, vsettl_3 )
1430 !
1431 !   computes the surface layer resistance term and the
1432 !   gravitational settling velocity term for the 0th and 3rd moment
1433 !   of a log-normal aerosol mode
1434 !
1435 !   input parameters
1436 !       dgnum - geometric mean diameter for aerosol number (m)
1437 !       alnsg - natural logarithm of the geometric standard deviation
1438 !               for aerosol number
1439 !       aerodens - aerosol density (dgnum and aerodens are for the
1440 !               actual wet distribution)
1441 !       temp - temperature (K)
1442 !       pres - pressure (Pa)
1443 !       airdens - air density (kg/m^3)
1444 !       ustar - friction velocity (m/s)
1445 !       wstar - convective scale velocity (m/s)
1446 !       depresist_a = aerodynamic resistance (s/m)
1447 !       airkinvisc - air kinematic viscosity (m^2/s)
1448 !       freepath - air molecular freepath (m)
1449 !   output parameters
1450 !       depvel_0 - deposition velocity for 0th (number) moment (m/s)
1451 !       depvel_3 - deposition velocity for 3rd (mass  ) moment (m/s)
1452 !       depresist_d0/_d3 - surface layer resistance for 0th/3rd moment (s/m)
1453 !       vsettl_0/_3 - gravitational settling velocity for 0th/3rd moment (m/s)
1454 !
1455
1456 !       use module_data_sorgam, only: boltz, grav, pss0, sqrt2, sqrtpi, tss0, threepi, two3
1457 !       use module_data_sorgam, only: ngausdv, wgaus, y_gq
1458
1459         implicit none
1460
1461         integer, intent(in) ::   &
1462             iflag_coarse
1463         real, intent(in) ::   &
1464             dgnum, alnsg, aerodens,   &
1465             temp, pres, airdens, ustar, depresist_a, pblh, zntt, rmolm
1466         real, intent(inout) ::   &
1467             airkinvisc, freepath
1468         real, intent(out) ::   &
1469             depvel_0, depvel_3,   &
1470             depresist_d0, depresist_d3,   &
1471             vsettl_0, vsettl_3
1472
1473
1474         ! idowescor = 1 means do wesley (85) convective correction to pm dry dep velocities; 0 means dont do correction
1475         integer, parameter :: idowescor = 0
1476
1477 ! following are from module_data_sorgam
1478 ! avogadro's constant [ 1/mol ]
1479         real, parameter :: avo=6.0221367e23
1480 ! universal gas constant [ j/mol-k ]
1481         real, parameter :: rgasuniv=8.314510
1482 ! boltzmann's constant [ j / k]
1483         real, parameter :: boltz=rgasuniv/avo
1484 ! mean gravitational acceleration [ m/sec**2 ]
1485         real, parameter :: grav=9.80622
1486 !  starting standard surface pressure [ pa ]
1487         real, parameter :: pss0=101325.0
1488 !  starting standard surface temperature [ k ]
1489         real, parameter :: tss0=288.15
1490         real*8, parameter :: pirs=3.14159265358979324
1491         real, parameter :: sqrt2=1.4142135623731
1492         real, parameter :: sqrtpi=1.7724539
1493         real, parameter :: two3=2.0/3.0
1494         real, parameter :: threepi=3.0*pirs
1495
1496
1497         real, parameter :: bhat = 1.246   ! constant from cunningham slip correction.
1498         real, parameter :: colctr_bigd=2.e-3,colctr_smald=20.e-6
1499                          ! collector diameters in stokes number and interception efficiency (needleleaf forest)
1500
1501         integer :: n
1502
1503         real :: alnsgy2, amu
1504         real :: blkta, blkprs, blkdens
1505         real :: cunq, czh
1506         real :: dconst1, dconst2, dconst3, dconst3y
1507         real :: dgyyy, dq
1508         real :: eff_dif, eff_imp, eff_int
1509         real :: ey1, esy04, esy08, esy16, esy20, esy28, esy32, esy64
1510         real :: knyyy, knq
1511         real :: nu
1512         real :: pdensy
1513         real :: ra, rbcor, rsurfq
1514         real :: scq, stq, sum0, sum0b, sum3, sum3b
1515         real :: tmpa
1516         real :: utscale
1517         real :: vdplim, vsedq
1518         real :: xlm
1519         real :: yextra_vdvg(100,3)
1520
1521         integer, parameter :: ngausdv = 7
1522         real, save :: y_gq(ngausdv), wgaus(ngausdv)
1523
1524
1525         ! define these locally because those in module_data_sorgam
1526         ! are not always initialized
1527         y_gq = (/ -2.651961356835233, &
1528         -1.673551628767471, -0.816287882858965, -0.0, &
1529         0.816287882858965, 1.673551628767471, 2.651961356835233 /)
1530         wgaus = (/ 0.0009717812450995, &
1531         0.05451558281913, 0.4256072526101, 0.8102646175568, &
1532         0.4256072526101, 0.05451558281913, 0.0009717812450995 /)
1533
1534         blkta = temp
1535         blkprs = pres
1536         blkdens = airdens
1537
1538 !   amu = air dynamic viscosity (kg/m/s)
1539         if (airkinvisc <= 0.0) then
1540             amu = 1.458e-6*blkta*sqrt(blkta)/(blkta+110.4)
1541             airkinvisc = amu/airdens
1542         else
1543             amu = airkinvisc*airdens
1544         end if
1545
1546 !   air molecular freepath (m)
1547         if (freepath <= 0.0) then
1548             freepath = 6.6328e-8*pss0*blkta/(tss0*blkprs)
1549         end if
1550         xlm = freepath
1551
1552         ra = depresist_a
1553
1554         dgyyy = dgnum
1555         pdensy = aerodens
1556
1557         alnsgy2 = alnsg**2
1558         ey1 = exp(0.125*alnsgy2)
1559         esy04 = ey1**4
1560         esy08 = esy04*esy04
1561         esy16 = esy08*esy08
1562         esy20 = esy16*esy04
1563         esy28 = esy20*esy08
1564         esy32 = esy16*esy16
1565         esy64 = esy32*esy32
1566
1567         dconst1 = boltz * blkta / ( threepi * amu )
1568         dconst2 = grav / ( 18.0 * amu )
1569         dconst3 =  ustar/(9.*amu*colctr_bigd)
1570
1571         nu = amu / blkdens
1572
1573         utscale =  1.
1574         if (idowescor.eq.1) then
1575         ! wesley (1985) monin-obukov dependence for convective conditions (sam 10/08)
1576             if (rmolm.lt.0.) then
1577                 czh = -1.*pblh*rmolm
1578                 if (czh.gt.30.0) then
1579                     utscale=0.45*(czh**0.6667)
1580                 else
1581                     utscale=1.+((-300.*rmolm)**0.6667)
1582                 endif
1583             endif
1584         endif   ! end of (idowescor.eq.1) test
1585         utscale = ustar*utscale
1586
1587         sum0 = 0. ; sum0b = 0.
1588         sum3 = 0. ; sum3b = 0.
1589
1590         do n = 1,ngausdv
1591         dq = dgyyy*exp(y_gq(n)*sqrt2*alnsg)  ! diameter (m) at quadrature point
1592         knq = 2.*xlm/dq  ! knudsen number at quadrature point
1593         cunq = 1.+knq*(1.257+.4*exp(-1.1/knq))  ! cunningham correction factor; pruppacher and klett (1980) eq (12-16)
1594         vsedq = pdensy*dconst2*cunq*dq*dq  ! gravitational sedimentation velocity m/s
1595         scq = nu*dq/dconst1/cunq  ! schmidt number, brownian diffusion parameter - same as binkowski and shankar
1596         eff_dif = scq**(-two3)    ! efficiency term for diffusion - same as binkowski and shankar
1597         stq = dconst3*pdensy*dq**2  ! stokes number, peters and eiden (1992)
1598         eff_imp = (stq/(0.8+stq))**2   ! efficiency term for impaction - peters and eiden (1992)
1599 !       eff_int = 0.3*dq/(colctr_smald+dq) ! slinn (1982) interception term, 0.3 prefac insures .2 cm/s at .2 micron diam.
1600         eff_int = (0.00116+0.0061*zntt)*dq/1.414e-7 ! mckeen(2008) intercptn trm, val of .00421 @ ustr=0.475, diam=.1414 micrn, stable, needleleaf evergreen
1601         if (iflag_coarse > 0) then
1602            eff_int = min(1.,eff_int)
1603            rbcor=exp(-2.0*(stq**0.5)) ! rebound correction factor used in slinn (1982)
1604         else
1605            rbcor = 1. ! rebound correction factor
1606         endif
1607         vdplim = utscale*(eff_dif+eff_imp+eff_int)*rbcor
1608 !       vdplim = .002*utscale
1609         vdplim = min(vdplim,.02)
1610 !PMA
1611         vdplim = max(vdplim,1.e-20)
1612         rsurfq = ra+1./vdplim
1613         sum0 = sum0 + wgaus(n)*(vsedq + 1./rsurfq)  ! quadrature sum for 0 moment
1614         sum3 = sum3 + wgaus(n)*(vsedq + 1./rsurfq)*dq**3  ! quadrature sum for 3rd moment
1615         sum0b = sum0b + wgaus(n)*(1./rsurfq)  ! quadrature sum for 0 moment
1616         sum3b = sum3b + wgaus(n)*(1./rsurfq)*dq**3  ! quadrature sum for 3rd moment
1617         end do ! n
1618
1619         ! normalize 0 moment vdep quadrature sum to sqrt(pi) (and number =1 per unit volume)
1620         depvel_0 = sum0/sqrtpi
1621         !normalize 3 moment quad. sum to sqrt(pi) and 3rd moment analytic sum
1622         depvel_3 = sum3/(sqrtpi*exp((1.5*sqrt2*alnsg)**2)*dgyyy**3)
1623
1624         sum0b = sum0b/sqrtpi
1625         sum3b = sum3b/(sqrtpi*exp((1.5*sqrt2*alnsg)**2)*dgyyy**3)
1626
1627         tmpa = 1.0/max( sum0b, 1.0e-10 )
1628         depresist_d0 = max( (tmpa - ra), 1.0e-10 )
1629         tmpa = 1.0/max( sum3b, 1.0e-10 )
1630         depresist_d3 = max( (tmpa - ra), 1.0e-10 )
1631
1632
1633 ! calculate settling velocity
1634 !       knyyy = 2.0*xlm/dgyyy
1635 !       dconst3y = dconst2 * pdensy * dgyyy**2
1636 !       vsettl_0 = dconst3y * ( esy16  + bhat * knyyy * esy04 )
1637 !       vsettl_3 = dconst3y * ( esy64  + bhat * knyyy * esy28 )
1638
1639 ! following are consistent with the quadrature calcs
1640         vsettl_0 = max( (depvel_0 - sum0b), 0.0 )
1641         vsettl_3 = max( (depvel_3 - sum3b), 0.0 )
1642
1643
1644         return
1645         end subroutine sorgam_aer_drydepvel_2
1646
1647
1648
1649 !-----------------------------------------------------------------------
1650      subroutine zhang2001_aer_drydepvel_1( &
1651          dg_num_in, lnsig_in, dens_part_in, &
1652          airtemp_in, airdens_in, ustar_in, resist_a_in, &
1653          zhang_optaa, iseason_in, &
1654          isurftype_in, modesurftype, iwetsurf, &
1655          airkinvisc, freepath,   &
1656          ido_mom0, ido_mom3,   &
1657          depvel_0, depvel_3,   &
1658          resist_s0, resist_s3,   &
1659          vsettl_0, vsettl_3 )
1660 !    calculates surface deposition velocity of particles
1661 !    L. Zhang, S. Gong, J. Padro, and L. Barrie
1662 !    A size-seggregated particle dry deposition scheme for an atmospheric aerosol module
1663 !    Atmospheric Environment, 35, 549-560, 2001.
1664 ! use precision
1665 ! use physconst, only: gravit, tmelt, rhodair
1666
1667       implicit none
1668
1669 !     in
1670       real     :: dg_num_in    ! geometric mean diameter of number distribution (m)
1671       real     :: lnsig_in     ! nat. log. of geometric standard deviation of particles
1672       real     :: dens_part_in ! density of particle material (kg/m3)
1673       real     :: airdens_in   ! air density (kg/m3)
1674       real     :: airtemp_in   ! air temperature (K)
1675       integer  :: zhang_optaa  ! 1=as published; 2=modified
1676       integer  :: iseason_in   ! seasonal category (same as wesely)
1677       integer  :: isurftype_in ! surface type id
1678       integer  :: modesurftype ! if <= 1, isurftype is zhang et al surface type
1679                                ! if >= 2, isurftype is usgs surface type
1680       integer  :: iwetsurf     ! 1 if surface is wet, 0 if dry
1681       real     :: resist_a_in  ! aerodynamic resistance (s/m)
1682       real     :: ustar_in     ! friction velocity (m/s)
1683       real     :: airkinvisc   ! kinematic viscosity of air (m2/s)
1684       real     :: freepath     ! mean free path of air molecules (m)
1685       integer  :: ido_mom0, ido_mom3 ! flag for do/skip moment 0 & 3 calculations
1686
1687 !     out
1688       real     :: depvel_0, depvel_3    ! surface deposition velocity for
1689                                         ! 0th (number) and 3rd (mass) moments (m/s)
1690       real     :: resist_s0, resist_s3  ! surface resistances (s/m)
1691       real     :: vsettl_0, vsettl_3    ! gravitational settling velocities (m/s)
1692
1693 !     local
1694       integer  :: iseason, isurftype
1695       integer  :: moment     ! moment of size distribution (0 for number, 2 for surface area, 3 for volume)
1696
1697       real, parameter :: gravit = 9.806
1698       real, parameter :: boltzmann=1.381e-23 ! J/deg/molecule
1699       real, parameter :: pi = 3.1415926536
1700
1701       real     :: airdens
1702       real     :: airtemp
1703       real     :: beta
1704       real     :: cuncorrect     ! Cunningham correction factor
1705       real     :: dens_part
1706       real     :: depvel         ! deposition velocity (m/s)
1707       real     :: dg_num
1708       real     :: diffus_part    ! particle Brownian diffusivity (m2/s)
1709       real     :: dispersion     ! accounts for influence of size dist dispersion on bulk settling velocity
1710                                 ! assuming dp_part is number mode radius * exp(1.5 ln(sigma))
1711       real     :: dp_mom         ! median diameter (m) for number (when moment=0) or volume (when moment=3)
1712       real     :: dp_vm          ! volume mean diameter = (6*volume_mixratio/(pi*number_mixratio))*0.333
1713       real     :: e_brown        ! collection efficiency for Browning diffusion
1714       real     :: e_impact       ! collection efficiency for impaction
1715       real     :: e_intcept      ! collection efficiency for interception
1716       real     :: e0_brown, e0_impact, e0_intcept   ! "tuning factors" for the above
1717       real     :: lnsig
1718       real     :: resist_a
1719       real     :: resist_sfc     ! surface resistance (s/m)
1720       real     :: schmidt        ! schmidt number
1721       real     :: stickfrac      ! fraction of particles sticking to surface
1722       real     :: stokes         ! stokes number
1723       real     :: tmpa
1724       real     :: ustar
1725       real     :: vsettl         ! gravitational settling velocity of particle (m/s)
1726       real     :: viscosity      ! dynamic viscosity of air (kg/m/s)
1727
1728       real, save :: gammax(15)     ! exponent of schmidt number
1729       data gammax/ 0.56, 0.58, 0.56, 0.56, 0.56, &
1730                    0.54, 0.54, 0.54, 0.54, 0.54, &
1731                    0.54, 0.54, 0.50, 0.50, 0.56/
1732
1733       real, save :: alpha(15)      ! parameter for impaction
1734       data alpha/   1.0,   0.6,   1.1,   0.8,   0.8, &
1735                     1.2,   1.2,  50.0,  50.0,   1.3, &
1736                     2.0,  50.0, 100.0, 100.0,   1.5/
1737
1738       real, save :: rad_collector(5,15) ! radius (m) of surface collectors
1739       data rad_collector/ &
1740            0.002, 0.002, 0.002, 0.002, 0.002, &  ! luc 1
1741            0.005, 0.005, 0.005, 0.005, 0.005, &  ! luc 2
1742            0.002, 0.002, 0.005, 0.005, 0.002, &  ! luc 3
1743            0.005, 0.005, 0.010, 0.010, 0.005, &  ! luc 4
1744            0.005, 0.005, 0.005, 0.005, 0.005, &  ! luc 5
1745            0.002, 0.002, 0.005, 0.005, 0.002, &  ! luc 6
1746            0.002, 0.002, 0.005, 0.005, 0.002, &  ! luc 7
1747            9.999, 9.999, 9.999, 9.999, 9.999, &  ! luc 8
1748            9.999, 9.999, 9.999, 9.999, 9.999, &  ! luc 9
1749            0.010, 0.010, 0.010, 0.010, 0.010, &  ! luc 10
1750            0.010, 0.010, 0.010, 0.010, 0.010, &  ! luc 11
1751            9.999, 9.999, 9.999, 9.999, 9.999, &  ! luc 12
1752            9.999, 9.999, 9.999, 9.999, 9.999, &  ! luc 13
1753            9.999, 9.999, 9.999, 9.999, 9.999, &  ! luc 14
1754            0.010, 0.010, 0.010, 0.010, 0.010  /  ! luc 15
1755
1756       integer, save  :: isurftype_zhang_from_usgs(25) ! maps usgs surface type to zhang type
1757       data isurftype_zhang_from_usgs / &
1758                ! zhang type         <---  usgs type                    --->  wesely type
1759         15, &  ! urban                    1: urban and built-up land         1  urban
1760          7, &  ! crops, mixed farming     2: dryland cropland & pasture      2  agricultural land
1761          7, &  !                          3: irrigated crop & pasture        2
1762          7, &  !                          4: mix dry/irrig crop & pasture    2
1763          7, &  !                          5: cropland/grassland mosaic       2
1764         10, &  ! shrubs & woodland        6: cropland/woodland mosaic        4
1765          6, &  ! grass                    7: grassland                       3  range land
1766          6, &  !                          8: shrubland                       3
1767          6, &  !                          9: mixed shrubland/grassland       3
1768          6, &  !                         10: savanna                         3,    (always summer)
1769          4, &  ! deciduous broadleaf     11: deciduous broadleaf forest      4  deciduous forest
1770          3, &  ! deciduous needle-lf     12: deciduous needleleaf forest     5,    (autumn and winter modi)
1771          2, &  ! evergreen broadleaf     13: evergreen broadleaf forest      4,    (always summer)
1772          1, &  ! evergreen needle-lf     14: evergreen needleleaf forest     5  coniferous forest
1773          5, &  ! mixed broad/needle-lf   15: mixed forest                    6  mixed forest including wetland
1774         14, &  ! ocean                   16: water bodies                    7  water, salt and fresh
1775         11, &  ! wetland w plants        17: herbaceous wetland              9
1776         10, &  !                         18: wooded wetland                  6
1777          8, &  ! desert                  19: barren or sparsely vegetated    8  barren land, mostly desert
1778          9, &  ! tundra                  20: herbaceous tundra               9  nonforested wetland
1779         10, &  !                         21: wooded tundra                   6
1780          9, &  !                         22: mixed tundra                    6
1781          8, &  !                         23: bare ground tundra              8
1782         12, &  ! ice cap & glacier       24: snow or ice                     -,    (always winter)
1783          8  /  !                         25: no data                         8
1784                !
1785                ! zhang 13 (inland water) not used                                wesley 10 (mixed agric/range)
1786                !                                                                      & 11 (rocky w low shrubs) not used
1787                !                         (the usgs ---> wesely mapping is used in module_dep_simple
1788                !                                        and is shown here for reference only)
1789
1790
1791 !  reality checks
1792       airtemp   = max( 50.0,   min( 350.0,  airtemp_in ) )
1793       airdens   = max( 0.01,   min( 5.0,    airdens_in ) )
1794       dg_num    = max( 1.0e-9, min( 1.0e-2, dg_num_in ) )
1795       dens_part = max( 0.1e3,  min( 50.0e3, dens_part_in ) )
1796       lnsig     = max( 0.0,    min( log(5.0), lnsig_in ) )
1797       resist_a  = max( 1.0e-2, min( 1.0e6,  resist_a_in ) )
1798       ustar     = max( 1.0e-6, min( 1.0e2,  ustar_in ) )
1799
1800       isurftype = 8
1801       if (modesurftype <= 1) then
1802          if ((isurftype_in >= 1) .and. (isurftype_in <= 15)) &
1803             isurftype = isurftype_in
1804       else
1805          ! map from usgs to zhang surface type
1806          if ((isurftype_in >= 1) .and. (isurftype_in <= 25)) &
1807             isurftype = isurftype_zhang_from_usgs(isurftype_in)
1808       end if
1809
1810       iseason = 1
1811       if ((iseason_in >= 1) .and. (iseason_in <= 5)) &
1812          iseason = iseason_in
1813
1814 !  for zhang_optaa = 2, do a modified version of zhang et al (2001)
1815 !     the interception efficiency has a first-order dependence on particle
1816 !        diameter (as in slinn 1982) rather than second order
1817 !     the "overall" e0 factor of 3.0 is changed to 1.0 for brownian diffusion
1818 !        and impaction, and 4.0 for interception
1819 !     note 1 -- e0_brown=e0_impact=1.0 reduces the rather large deposition
1820 !        velocities for very small and very large particle sizes
1821 !     note 2 -- e0_intcept=4.0 is equivalent to using e0_intcept=1.0
1822 !        and dividing the radius_collector by 4.0.  this change plus the
1823 !        first-order form for interceptions gives much more interception in
1824 !        the accumulation-mode size range.  the interception in the original
1825 !        zhang et al formulation is (surprisingly) almost negligible!
1826       if (zhang_optaa == 2) then
1827          e0_brown   = 1.0
1828          e0_impact  = 1.0
1829          e0_intcept = 4.0
1830       else
1831          e0_brown   = 3.0
1832          e0_impact  = 3.0
1833          e0_intcept = 3.0
1834       end if
1835
1836
1837 !  following is from subr modpar in module_aerosols_sorgam.F
1838 !     pss0 = 101325.0 ; tss0 = 288.15
1839 !     freepath = 6.6328e-8*pss0*blkta/(tss0*blkprs)
1840       if (freepath <= 0.0) then
1841          freepath = 6.6328e-8*(1.2255/airdens)
1842       end if
1843
1844       if (airkinvisc <= 0.0) then
1845 !  following is from subr modpar in module_aerosols_sorgam.F
1846 !     amu = 1.458e-6*blkta*sqrt(blkta)/(blkta+110.4)
1847          viscosity = 1.458e-6*airtemp*sqrt(airtemp)/(airtemp+110.4)
1848          airkinvisc = viscosity/airdens
1849       else
1850          viscosity = airkinvisc*airdens
1851       end if
1852
1853       dp_vm = dg_num*exp(1.5*lnsig*lnsig)
1854
1855       do moment = 0, 3, 3
1856
1857       if (moment == 0) then
1858          if (ido_mom0 <= 0) then
1859             depvel_0 = 0.0
1860             vsettl_0 = 0.0
1861             resist_s0 = 1.0e30
1862             cycle
1863          end if
1864       else if (moment == 3) then
1865          if (ido_mom3 <= 0) then
1866             depvel_3 = 0.0
1867             vsettl_3 = 0.0
1868             resist_s3 = 1.0e30
1869             cycle
1870          end if
1871       else
1872          cycle
1873       end if
1874
1875 ! use a maximum radius of 100 microns when calculating deposition velocity
1876       dp_mom = dg_num*exp(float(moment)*lnsig*lnsig)
1877       dp_mom = min( 100.0e-6, dp_mom )
1878
1879       tmpa = 2.0*freepath/dp_mom
1880       cuncorrect = 1. + tmpa*(1.257 + 0.4*exp(-1.10/tmpa))
1881       dispersion = exp(2.*lnsig*lnsig)
1882       vsettl = ( dens_part*(dp_mom**2)*gravit*   &
1883                         cuncorrect / (18.0*viscosity) ) * dispersion
1884
1885       diffus_part = boltzmann*airtemp*cuncorrect/   &
1886                          (3.0*pi*viscosity*dp_mom)
1887       schmidt=airkinvisc/diffus_part
1888       e_brown = schmidt**(-gammax(isurftype))
1889
1890       if (isurftype==8  .or. isurftype==9  .or. &
1891           isurftype==12 .or. isurftype==13 .or. isurftype==14) then
1892 !        non-vegetated surface
1893          stokes = vsettl*ustar*ustar/(gravit*airkinvisc)
1894          e_intcept = 0.
1895       else
1896 !        vegetated surface
1897          stokes = vsettl*ustar/(gravit*rad_collector(iseason,isurftype))
1898          if (zhang_optaa == 2) then
1899             tmpa = 0.5*dp_mom/rad_collector(iseason,isurftype)
1900             e_intcept = tmpa
1901          else
1902             tmpa = dp_mom/rad_collector(iseason,isurftype)
1903             e_intcept = 0.5 * tmpa**2
1904          endif
1905       endif
1906
1907       beta = 2.0
1908       e_impact = ( stokes / (alpha(isurftype)+stokes) )**beta
1909
1910 !     stickfrac = 1 if diameter < 5 um or surface is wet
1911       if ((iwetsurf > 0) .or. &
1912           (isurftype == 13) .or. (isurftype == 14) .or. &
1913           (dp_vm < 5.0e-6)) then
1914          stickfrac = 1.0
1915       else
1916          stickfrac = exp(-sqrt(stokes))
1917       end if
1918
1919 !     e0 = 3.0
1920 !     tmpa = e0*ustar*(e_brown+e_intcept+e_impact)*stickfrac
1921       tmpa = ustar*(e0_brown*e_brown + e0_intcept*e_intcept + e0_impact*e_impact)*stickfrac
1922       resist_sfc = 1./max(tmpa,1.0e-30)
1923
1924 !     depvel = vsettl+1./(resist_a+resist_sfc)
1925 ! use this formula rather than zhang 2001 simple formula
1926       depvel = vsettl + 1./(resist_a + resist_sfc + resist_a*resist_sfc*vsettl)
1927
1928       if (moment == 0) then
1929          depvel_0 = depvel
1930          vsettl_0 = vsettl
1931          resist_s0 = resist_sfc
1932       else
1933          depvel_3 = depvel
1934          vsettl_3 = vsettl
1935          resist_s3 = resist_sfc
1936       end if
1937
1938       end do ! moment
1939
1940       return
1941       end subroutine zhang2001_aer_drydepvel_1
1942
1943 !-----------------------------------------------------------------------
1944         subroutine cam_mam_aer_drydep_prep(                                 &
1945                 id, ktau, dtstep, config_flags, aer_mech_id,               &
1946                 gmt, julday,                                               &
1947                 t_phy, rho_phy, p_phy, alt, p8w, t8w,                      &
1948                 moist, chem,                                               &
1949                 nbin_add, alnsg_add, wetdgnum_add, wetdens_add,            &
1950                 yextra_add,                                                &
1951                 i, j,                                                      &
1952                 ids,ide, jds,jde, kds,kde,                                 &
1953                 ims,ime, jms,jme, kms,kme,                                 &
1954                 its,ite, jts,jte, kts,kte                                  )
1955
1956         use module_configure, only:  num_moist, num_chem, &
1957                 grid_config_rec_type
1958         use module_state_description, only:  param_first_scalar
1959
1960         use module_data_cam_mam_asect
1961
1962         implicit none
1963
1964 !   subr arguments
1965         integer, intent(in) ::   &
1966                 id, ktau, julday, aer_mech_id, nbin_add,   &
1967                 i, j
1968
1969         integer, intent(in) ::   &
1970                 ids, ide, jds, jde, kds, kde,   &
1971                 ims, ime, jms, jme, kms, kme,   &
1972                 its, ite, jts, jte, kts, kte
1973
1974         real, intent(in) :: dtstep, gmt
1975
1976         real, intent(in),   &
1977                 dimension( ims:ime, kms:kme, jms:jme ) :: &
1978                 t_phy, rho_phy, p_phy, alt, p8w, t8w
1979
1980         real, intent(in),   &
1981                 dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
1982                 moist
1983
1984         real, intent(in),   &
1985                 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
1986                 chem
1987
1988         real, intent(inout),   &
1989                 dimension( 1:nbin_add ) :: &
1990                 alnsg_add, wetdgnum_add, wetdens_add
1991
1992         real, intent(inout),   &
1993                 dimension( 1:nbin_add, 10 ) :: &
1994                 yextra_add
1995
1996         type(grid_config_rec_type), intent(in) :: config_flags
1997
1998
1999 !   local variables
2000         integer :: ibin, iok, iphase, itype
2001         integer :: itmpa
2002         integer :: k
2003         integer :: l, ll, l1, n
2004         integer :: p1st
2005
2006         real, parameter :: densdefault = 2.0
2007         real, parameter :: smallmassaa = 1.0e-20
2008         real, parameter :: smallmassbb = 1.0e-30
2009         real, parameter :: pi = 3.1415926536
2010         real, parameter :: piover6 = pi/6.0
2011         real, parameter :: onethird = 1.0/3.0
2012
2013         real :: drydens, drydp, drymass, dryvol
2014         real :: rnum
2015         real :: tmpa, tmp_lnsg
2016         real :: wetdgnum, wetdens, wetdp, wetmass, wetvol
2017         real :: wetdgnum_si, wetdens_si
2018
2019         character(len=160) :: txtaa
2020
2021
2022         n = sum( nsize_aer(1:ntype_aer) )
2023         if (n /= nbin_add) then
2024             write(txtaa,'(a,2(1x,i10))') &
2025                 '*** cam_mam_aer_drydep_prep - bad nbin_add, sum(nsize)', &
2026                 nbin_add, n
2027             write(*,'(a)') trim(txtaa)
2028             call wrf_error_fatal( txtaa )
2029             stop
2030         end if
2031
2032
2033 !   set some variables to their wrf-chem "standard" values
2034         p1st = param_first_scalar
2035
2036         k = kts
2037
2038 !   *** for now, just calc vdep_aer for iphase = ai_phase
2039         iphase = ai_phase
2040
2041 !   calculate vdep_aer for each type and size
2042         ibin = 0
2043         do 2900 itype = 1, ntype_aer
2044         do 2800 n = 1, nsize_aer(itype)
2045         ibin = ibin + 1
2046
2047 !   calculate dry and wet mass, volume, density, and mass-mean diameters
2048         dryvol = 0.0
2049         drymass = 0.0
2050         do ll = 1, ncomp_aer(itype)
2051             l1 = massptr_aer(ll,n,itype,iphase)
2052             tmpa = max( chem(i,k,j,l1), 0.0 )
2053             drymass = drymass + tmpa
2054             dryvol = dryvol + tmpa/dens_aer(ll,itype)
2055         end do
2056
2057         l1 = waterptr_aer(n,itype)
2058         tmpa = max( chem(i,k,j,l1), 0.0 )
2059         wetmass = drymass + tmpa
2060         wetvol = dryvol + tmpa/dens_water_aer
2061
2062         l1 = numptr_aer(n,itype,iphase)
2063         rnum = max( chem(i,k,j,l1), 0.0 )
2064
2065         drymass = drymass*28.966e-9 ! ug/kg-air --> g/mol-air
2066         wetmass = wetmass*28.966e-9
2067         dryvol  = dryvol *28.966e-9 ! "u-cm3"/kg-air --> cm3/mol-air
2068         wetvol  = wetvol *28.966e-9
2069         rnum    = rnum*28.966e-3    ! #/kg-air  --> #/mol-air
2070
2071         if (drymass <= smallmassbb) then
2072             drydp = dcen_sect(n,itype)
2073             drydens = densdefault
2074             wetdp = drydp
2075             wetdens = drydens
2076
2077         else
2078             if (drymass <= smallmassaa) then
2079                 wetmass = drymass
2080                 wetvol = dryvol
2081             end if
2082             drydens = drymass/dryvol
2083             wetdens = wetmass/wetvol
2084
2085             if (rnum >= dryvol/volumlo_sect(n,itype)) then
2086                 drydp = dlo_sect(n,itype)
2087             else if (rnum <= dryvol/volumhi_sect(n,itype)) then
2088                 drydp = dhi_sect(n,itype)
2089             else
2090                 drydp = (dryvol/(piover6*rnum))**onethird
2091             end if
2092
2093             if (abs(wetvol) > (1000.*abs(dryvol))) then
2094               tmpa=10.0
2095             else
2096               tmpa=abs(wetvol/dryvol)**onethird
2097               tmpa=max(1.0,min(tmpa,10.0))
2098             endif
2099             wetdp = drydp*tmpa
2100
2101         end if
2102
2103
2104 !   "within section" size distribution currently assumed to be
2105 !   monodisperse, but could assume otherwise
2106         tmp_lnsg = log( sigmag_aer(n,itype) )
2107         wetdgnum = wetdp * exp( -1.5*tmp_lnsg*tmp_lnsg )
2108
2109         wetdgnum_add(ibin) = wetdgnum * 1.0e-2 ! cm --> m
2110         wetdens_add(ibin) = wetdens * 1.0e3    ! g/cm3 --> kg/m3
2111         alnsg_add(ibin) = tmp_lnsg
2112
2113
2114         yextra_add(ibin,1) = drymass/28.966e-9 ! g/mol-air --> ug/kg-air = "n-kg"/kg-air
2115         yextra_add(ibin,5) = wetmass/28.966e-9
2116         yextra_add(ibin,2) = dryvol/28.966e-6  ! cm3/mol-air --> "n-m3"/kg-air
2117         yextra_add(ibin,6) = wetvol/28.966e-6
2118         yextra_add(ibin,3) = drydens*1.0e3     ! g/cm3 --> kg/m3
2119         yextra_add(ibin,7) = wetdens*1.0e3
2120         yextra_add(ibin,4) = drydp*1.0e-2      ! cm --> m
2121         yextra_add(ibin,8) = wetdp*1.0e-2
2122
2123         yextra_add(ibin,9) = dcen_sect(n,itype)*1.0e-2   ! cm --> m
2124         yextra_add(ibin,10) = rnum/28.966e-3    ! #/mol-air  --> #/kg-air
2125
2126
2127 2800    continue
2128 2900    continue
2129
2130
2131         return
2132         end subroutine cam_mam_aer_drydep_prep
2133
2134
2135
2136 !-----------------------------------------------------------------------
2137         subroutine cam_mam_aer_drydep_load_ddvel(                           &
2138                 id, ktau, dtstep, config_flags, aer_mech_id,               &
2139                 ddvel,                                                     &
2140                 nbin_add, depvel0_add, depvel3_add,                        &
2141                 i, j,                                                      &
2142                 ids,ide, jds,jde, kds,kde,                                 &
2143                 ims,ime, jms,jme, kms,kme,                                 &
2144                 its,ite, jts,jte, kts,kte                                  )
2145
2146         use module_configure, only:  num_moist, num_chem, &
2147                 grid_config_rec_type
2148         use module_state_description, only:  param_first_scalar
2149
2150         use module_data_cam_mam_asect
2151
2152         implicit none
2153
2154 !   subr arguments
2155         integer, intent(in) ::   &
2156                 id, ktau, aer_mech_id, nbin_add,   &
2157                 i, j
2158
2159         integer, intent(in) ::   &
2160                 ids, ide, jds, jde, kds, kde,   &
2161                 ims, ime, jms, jme, kms, kme,   &
2162                 its, ite, jts, jte, kts, kte
2163
2164         real, intent(in) :: dtstep
2165
2166         real, intent(in),   &
2167                 dimension( 1:nbin_add ) :: &
2168                 depvel0_add, depvel3_add
2169
2170         real, intent(out),   &
2171                 dimension( its:ite, jts:jte, 1:num_chem ) :: &
2172                 ddvel
2173
2174         type(grid_config_rec_type), intent(in) :: config_flags
2175
2176
2177 !   local variables
2178         integer :: ibin, iphase, itype
2179         integer :: l, ll, n
2180
2181         character(len=160) :: txtaa
2182
2183
2184         n = sum( nsize_aer(1:ntype_aer) )
2185         if (n /= nbin_add) then
2186             write(txtaa,'(a,2(1x,i10))') &
2187                 '*** cam_mam_aer_drydep_load_ddvel - bad nbin_add, sum(nsize)', &
2188                 nbin_add, n
2189             write(*,'(a)') trim(txtaa)
2190             call wrf_error_fatal( txtaa )
2191             stop
2192         end if
2193
2194
2195 ! copy from vdep_aer array to ddvel array
2196         do iphase = 1, nphase_aer
2197         ibin = 0
2198         do itype = 1, ntype_aer
2199         do n = 1, nsize_aer(itype)
2200             ibin = ibin + 1
2201             do ll = -2, ncomp_plustracer_aer(itype)
2202                 if (iphase /= ai_phase) cycle
2203
2204                 if (ll == -2) then
2205                     l = numptr_aer(n,itype,iphase)
2206                 else if (ll == -1) then
2207                     l = -1
2208                     if (iphase == ai_phase) l = waterptr_aer(n,itype)
2209                 else if (ll == 0) then
2210                     l = -1
2211                     if (iphase == ai_phase) l = hyswptr_aer(n,itype)
2212                 else
2213                     l = massptr_aer(ll,n,itype,iphase)
2214                 end if
2215
2216                 if ((l >= param_first_scalar) .and. (l <= num_chem)) then
2217                     if (l == numptr_aer(n,itype,iphase)) then
2218                         ddvel(i,j,l) = depvel0_add(ibin)
2219                     else
2220                         ddvel(i,j,l) = depvel3_add(ibin)
2221                     end if
2222                 end if
2223             end do
2224         end do
2225         end do
2226         end do
2227
2228
2229         return
2230         end subroutine cam_mam_aer_drydep_load_ddvel
2231
2232
2233
2234 !-----------------------------------------------------------------------
2235 END MODULE module_aer_drydep