chem/KPP/mechanisms/mozart/mo_usrrxt.F90

   1 ! lke: corrections to aerosol rxns and clean-up
   2
   3       module mo_usrrxt
   4
   5       private
   6       public :: usrrxt_inti, usrrxt
   7
   8       save
   9
  10       integer :: usr1_ndx, usr2_ndx, usr3_ndx, usr5_ndx, usr6_ndx, usr7_ndx
  11       integer :: usr8_ndx, usr9_ndx, usr11_ndx, usr12_ndx, usr14_ndx, usr15_ndx
  12       integer :: usr16_ndx, usr17_ndx, usr21_ndx, usr22_ndx
  13       integer :: usr23_ndx, usr24_ndx, usr25_ndx, usr26_ndx, usr27_ndx, usr28_ndx
  14       integer :: usr17a_ndx
  15       integer :: so4_ndx, cb2_ndx, oc2_ndx, soa_ndx, nit_ndx
  16       logical :: has_aerosols
  17
  18       contains
  19
  20       subroutine usrrxt_inti
  21 !-----------------------------------------------------------------
  22 !        ... intialize the user reaction constants module
  23 !-----------------------------------------------------------------
  24
  25       use mo_chem_utls, only : get_rxt_ndx, get_spc_ndx
  26
  27       implicit none
  28
  29       usr1_ndx   = get_rxt_ndx( 'usr1' )
  30       usr2_ndx   = get_rxt_ndx( 'usr2' )
  31       usr3_ndx   = get_rxt_ndx( 'usr3' )
  32       usr5_ndx   = get_rxt_ndx( 'usr5' )
  33       usr6_ndx   = get_rxt_ndx( 'usr6' )
  34       usr7_ndx   = get_rxt_ndx( 'usr7' )
  35       usr8_ndx   = get_rxt_ndx( 'usr8' )
  36       usr9_ndx   = get_rxt_ndx( 'usr9' )
  37       usr11_ndx  = get_rxt_ndx( 'usr11' )
  38       usr12_ndx  = get_rxt_ndx( 'usr12' )
  39       usr14_ndx  = get_rxt_ndx( 'usr14' )
  40       usr15_ndx  = get_rxt_ndx( 'usr15' )
  41       usr16_ndx  = get_rxt_ndx( 'usr16' )
  42       usr17_ndx  = get_rxt_ndx( 'usr17' )
  43       usr21_ndx  = get_rxt_ndx( 'usr21' )
  44       usr22_ndx  = get_rxt_ndx( 'usr22' )
  45       usr17a_ndx = get_rxt_ndx( 'usr17a' )
  46       usr23_ndx  = get_rxt_ndx( 'usr23' )
  47       usr24_ndx  = get_rxt_ndx( 'usr24' )
  48       usr25_ndx  = get_rxt_ndx( 'usr25' )
  49       usr26_ndx  = get_rxt_ndx( 'usr26' )
  50       usr27_ndx  = get_rxt_ndx( 'usr27' )
  51       usr28_ndx  = get_rxt_ndx( 'usr28' )
  52       so4_ndx    = get_spc_ndx( 'SO4' )
  53       cb2_ndx    = get_spc_ndx( 'CB2' )
  54       oc2_ndx    = get_spc_ndx( 'OC2' )
  55       soa_ndx    = get_spc_ndx( 'SOA' )
  56       nit_ndx    = get_spc_ndx( 'NH4NO3' )
  57
  58       write(*,*) ' '
  59       write(*,*) 'usrrxt_inti: diagnostics '
  60       write(*,'(10i5)') usr1_ndx, usr2_ndx, usr3_ndx, usr5_ndx, usr6_ndx, usr7_ndx, &
  61                         usr8_ndx, usr9_ndx, usr11_ndx, usr12_ndx, usr14_ndx, usr15_ndx, &
  62                         usr16_ndx, usr17_ndx, usr21_ndx, usr22_ndx
  63
  64       end subroutine usrrxt_inti
  65
  66       subroutine usrrxt( rxt, temp, invariants, h2ovmr,  ps, &
  67                          pmid, m, sulfate, vmr, mmr, &
  68                          relhum, strato_sad, ltrop, lat, ip, plonl, sad_total )
  69
  70 !-----------------------------------------------------------------
  71 !        ... set the user specified reaction rates
  72 !-----------------------------------------------------------------
  73
  74       use chem_mods,    only : nfs, rxntot
  75       use mo_setinv,    only : has_h2o, h2o_ndx
  76       use mo_grid,      only : plev, pcnstm1, pcnst
  77       use mo_constants, only : pi, avo => avogadro, boltz
  78
  79       implicit none
  80
  81 !-----------------------------------------------------------------
  82 !        ... dummy arguments
  83 !-----------------------------------------------------------------
  84       integer, intent(in) :: lat, ip
  85       integer, intent(in) :: plonl
  86       integer, intent(in) :: ltrop(plonl)                ! tropopause vertical index
  87       real, intent(in)    :: temp(plonl,plev)            ! temperature (K)
  88       real, intent(in)    :: m(plonl,plev)               ! total atm density (/cm^3)
  89       real, intent(in)    :: sulfate(plonl,plev)         ! sulfate aerosol (mol/mol)
  90       real, intent(in)    :: strato_sad(plonl,plev)      ! stratospheric aerosol sad (1/cm)
  91       real, intent(in)    :: h2ovmr(plonl,plev)          ! water vapor (mol/mol)
  92       real, intent(in)    :: relhum(plonl,plev)          ! relative humidity
  93       real, intent(in)    :: pmid(plonl,plev)            ! midpoint pressure (Pa)
  94       real, intent(in)    :: ps(plonl)                   ! surface pressure (Pa)
  95       real, intent(in)    :: invariants(plonl,plev,nfs)  ! invariants density (/cm^3)
  96       real, intent(in)    :: vmr(plonl,plev,pcnstm1)     ! species concentrations (mol/mol)
  97       real, intent(in)    :: mmr(plonl,plev,pcnst  )     ! species concentrations (kg/kg)
  98       real, intent(inout) :: rxt(plonl,plev,rxntot)      ! gas phase rates
  99       real, intent(inout) :: sad_total(plonl,plev)       ! total surface area density (cm2/cm3)
 100
 101 !-----------------------------------------------------------------
 102 !        ... local variables
 103 !-----------------------------------------------------------------
 104
 105       real, parameter :: dg = 0.1                        ! mole diffusion =0.1 cm2/s (Dentener, 1993)
 106       real, parameter :: mw_so4 = 98.e-3                 ! so4 molecular wt (kg/mole)
 107
 108 !-----------------------------------------------------------------
 109 !       ... parameters for log-normal distribution by number
 110 ! references:
 111 !   Chin et al., JAS, 59, 461, 2003
 112 !   Liao et al., JGR, 108(D1), 4001, 2003
 113 !   Martin et al., JGR, 108(D3), 4097, 2003
 114 !-----------------------------------------------------------------
 115       real, parameter :: rm_sulf  = 6.95e-6        ! mean radius of sulfate particles (cm) (Chin)
 116       real, parameter :: sd_sulf  = 2.03           ! standard deviation of radius for sulfate (Chin)
 117       real, parameter :: rho_sulf = 1.7e3          ! density of sulfate aerosols (kg/m3) (Chin)
 118
 119       real, parameter :: rm_orgc  = 2.12e-6        ! mean radius of organic carbon particles (cm) (Chin)
 120       real, parameter :: sd_orgc  = 2.20           ! standard deviation of radius for OC (Chin)
 121       real, parameter :: rho_orgc = 1.8e3          ! density of OC aerosols (kg/m3) (Chin)
 122
 123       real, parameter :: rm_bc    = 1.18e-6        ! mean radius of soot/BC particles (cm) (Chin)
 124       real, parameter :: sd_bc    = 2.00           ! standard deviation of radius for BC (Chin)
 125       real, parameter :: rho_bc   = 1.0e3          ! density of BC aerosols (kg/m3) (Chin)
 126
 127 !-----------------------------------------------------------------
 128 !       ... reaction probabilities for heterogeneous reactions
 129 !-----------------------------------------------------------------
 130       real, parameter :: gamma_n2o5 = 0.10         ! from Jacob, Atm Env, 34, 2131, 2000
 131       real, parameter :: gamma_ho2  = 0.20         !
 132       real, parameter :: gamma_no2  = 0.0001       !
 133       real, parameter :: gamma_no3  = 0.001        !
 134
 135 !-----------------------------------------------------------------
 136 !       ... table for hygroscopic growth effect on radius (Chin et al)
 137 !           (no growth effect for mineral dust)
 138 !-----------------------------------------------------------------
 139       real, dimension(7) :: table_rh, table_rfac_sulf, table_rfac_bc, table_rfac_oc, table_rfac_ss
 140
 141       data table_rh(1:7)        / 0.0, 0.5, 0.7, 0.8, 0.9, 0.95, 0.99/
 142       data table_rfac_sulf(1:7) / 1.0, 1.4, 1.5, 1.6, 1.8, 1.9,  2.2/
 143       data table_rfac_oc(1:7)   / 1.0, 1.2, 1.4, 1.5, 1.6, 1.8,  2.2/
 144       data table_rfac_bc(1:7)   / 1.0, 1.0, 1.0, 1.2, 1.4, 1.5,  1.9/
 145       data table_rfac_ss(1:7)   / 1.0, 1.6, 1.8, 2.0, 2.4, 2.9,  4.8/
 146
 147       integer  ::  i, k
 148       real     ::  tp(plonl)                       ! 300/t
 149       real     ::  tinv(plonl)                     ! 1/t
 150       real     ::  ko(plonl)
 151       real     ::  kinf(plonl)
 152       real     ::  fc(plonl)
 153       real     ::  sqrt_t(plonl)                   ! sqrt( temp )
 154       real     ::  exp_fac(plonl)                  ! vector exponential
 155
 156       real     ::  v, rho_air, n, n_exp, r_rd, r_sd
 157       real     ::  dm_sulf, dm_sulf_wet, log_sd_sulf, sfc_sulf, sfc_nit
 158       real     ::  dm_orgc, dm_orgc_wet, log_sd_orgc, sfc_oc, sfc_soa
 159       real     ::  dm_bc, dm_bc_wet, log_sd_bc, sfc_bc
 160       real     ::  rxt_sulf, rxt_nit, rxt_oc, rxt_soa
 161       real     ::  c_n2o5, c_ho2, c_no2, c_no3
 162       real     ::  s_exp
 163
 164       integer  ::  irh, rh_l, rh_u
 165       real     ::  factor, rfac_sulf, rfac_oc, rfac_bc, rfac_ss
 166
 167 !-----------------------------------------------------------------
 168 !       ... o + o2 + m --> o3 + m
 169 !-----------------------------------------------------------------
 170 level_loop : &
 171       do k = 1,plev
 172          tinv(:)           = 1. / temp(:,k)
 173          tp(:)             = 300. * tinv(:)
 174          sqrt_t(:)         = sqrt( temp(:,k) )
 175          if( usr1_ndx > 0 ) then
 176             rxt(:,k,usr1_ndx) = 6.e-34 * tp(:)**2.4
 177          end if
 178 #ifdef IBM
 179 !-----------------------------------------------------------------
 180 !       ... n2o5 + m --> no2 + no3 + m
 181 !-----------------------------------------------------------------
 182          if( usr3_ndx > 0 ) then
 183             if( usr2_ndx > 0 ) then
 184                call vexp( exp_fac, -10990.*tinv, plonl )
 185                rxt(:,k,usr3_ndx) = rxt(:,k,usr2_ndx) * 3.333e26 * exp_fac(:)
 186             else
 187                rxt(:,k,usr3_ndx) = 0.
 188             end if
 189          end if
 190
 191 !-----------------------------------------------------------------
 192 !       set rates for:
 193 !       ... hno3 + oh --> no3 + h2o
 194 !           ho2no2 + m --> ho2 + no2 + m
 195 !           co + oh --> co2 + ho2
 196 !-----------------------------------------------------------------
 197          if( usr5_ndx > 0 ) then
 198             call vexp( exp_fac, 1335.*tinv, plonl )
 199             ko(:) = m(:,k) * 6.5e-34 * exp_fac(:)
 200             call vexp( exp_fac, 2199.*tinv, plonl )
 201             ko(:) = ko(:) / (1. + ko(:)/(2.7e-17*exp_fac(:)))
 202             call vexp( exp_fac, 460.*tinv, plonl )
 203             rxt(:,k,usr5_ndx) = ko(:) + 2.4e-14*exp_fac(:)
 204          end if
 205          if( usr7_ndx > 0 ) then
 206             if( usr6_ndx > 0 ) then
 207                call vexp( exp_fac, -10900.*tinv, plonl )
 208                rxt(:,k,usr7_ndx) = rxt(:,k,usr6_ndx) * exp_fac(:) / 2.1e-27
 209             else
 210                rxt(:,k,usr7_ndx) = 0.
 211             end if
 212          end if
 213          if( usr8_ndx > 0 ) then
 214             rxt(:,k,usr8_ndx) = 1.5e-13 * (1. + 6.e-7*boltz*m(:,k)*temp(:,k))
 215          end if
 216
 217 !-----------------------------------------------------------------
 218 !       ... ho2 + ho2 --> h2o2
 219 !       note: this rate involves the water vapor number density
 220 !-----------------------------------------------------------------
 221          if( usr9_ndx > 0 ) then
 222             if( has_h2o ) then
 223                call vexp( exp_fac, 600.*tinv, plonl )
 224                ko(:)   = 2.3e-13 * exp_fac(:)
 225                call vexp( exp_fac, 1000.*tinv, plonl )
 226                kinf(:) = 1.7e-33 * m(:,k) * exp_fac(:)
 227                call vexp( exp_fac, 2200.*tinv, plonl )
 228                fc(:)   = 1. + 1.4e-21 * invariants(:,k,h2o_ndx) * exp_fac(:)
 229                rxt(:,k,usr9_ndx) = (ko(:) + kinf(:)) * fc(:)
 230             else
 231                rxt(:,k,usr9_ndx) = 0.
 232             end if
 233          end if
 234
 235 !-----------------------------------------------------------------
 236 !       ... mco3 + no2 -> mpan
 237 !-----------------------------------------------------------------
 238          if( usr14_ndx > 0 ) then
 239             rxt(:,k,usr14_ndx) = 1.1e-11 * tp(:) / m(:,k)
 240          end if
 241
 242 !-----------------------------------------------------------------
 243 !       ... pan + m --> ch3co3 + no2 + m
 244 !-----------------------------------------------------------------
 245          call vexp( exp_fac, -14000.*tinv, plonl )
 246          if( usr12_ndx > 0 ) then
 247             if( usr11_ndx > 0 ) then
 248                rxt(:,k,usr12_ndx) = rxt(:,k,usr11_ndx) * 1.111e28 * exp_fac(:)
 249             else
 250                rxt(:,k,usr12_ndx) = 0.
 251             end if
 252          end if
 253
 254 !-----------------------------------------------------------------
 255 !       ... mpan + m --> mco3 + no2 + m
 256 !-----------------------------------------------------------------
 257          if( usr15_ndx > 0 ) then
 258             if( usr14_ndx > 0 ) then
 259                rxt(:,k,usr15_ndx) = rxt(:,k,usr14_ndx) * 1.111e28 * exp_fac(:)
 260             else
 261                rxt(:,k,usr15_ndx) = 0.
 262             end if
 263          end if
 264
 265 !-----------------------------------------------------------------
 266 !       ... xooh + oh -> h2o + oh
 267 !-----------------------------------------------------------------
 268          if( usr21_ndx > 0 ) then
 269             call vexp( exp_fac, 253.*tinv, plonl )
 270             rxt(:,k,usr21_ndx) = temp(:,k)**2 * 7.69e-17 * exp_fac(:)
 271          end if
 272
 273 !-----------------------------------------------------------------
 274 !       ... ch3coch3 + oh -> ro2 + h2o
 275 !-----------------------------------------------------------------
 276          if( usr22_ndx > 0 ) then
 277             call vexp( exp_fac, -2000.*tinv, plonl )
 278             rxt(:,k,usr22_ndx) = 3.82e-11 * exp_fac(:) + 1.33e-13
 279          end if
 280
 281 !-----------------------------------------------------------------
 282 !       ... DMS + OH  --> .5 * SO2
 283 !-----------------------------------------------------------------
 284          if( usr24_ndx > 0 ) then
 285             call vexp( exp_fac, 7460.*tinv, plonl )
 286             ko(:) = 1. + 5.5e-31 * exp_fac * m(:,k) * 0.21
 287             call vexp( exp_fac, 7810.*tinv, plonl )
 288             rxt(:,k,usr24_ndx) = 1.7e-42 * exp_fac * m(:,k) * 0.21 / ko(:)
 289          end if
 290
 291 #else
 292 !-----------------------------------------------------------------
 293 !       ... n2o5 + m --> no2 + no3 + m
 294 !-----------------------------------------------------------------
 295          if( usr3_ndx > 0 ) then
 296             if( usr2_ndx > 0 ) then
 297                rxt(:,k,usr3_ndx) = rxt(:,k,usr2_ndx) * 3.333e26 * exp( -10990.*tinv(:) )
 298             else
 299                rxt(:,k,usr3_ndx) = 0.
 300             end if
 301          end if
 302
 303 !-----------------------------------------------------------------
 304 !       set rates for:
 305 !       ... hno3 + oh --> no3 + h2o
 306 !           ho2no2 + m --> ho2 + no2 + m
 307 !           co + oh --> co2 + ho2
 308 !-----------------------------------------------------------------
 309          if( usr5_ndx > 0 ) then
 310             ko(:) = m(:,k) * 6.5e-34 * exp( 1335.*tinv(:) )
 311             ko(:) = ko(:) / (1. + ko(:)/(2.7e-17*exp( 2199.*tinv(:) )))
 312             rxt(:,k,usr5_ndx) = ko(:) + 2.4e-14*exp( 460.*tinv(:) )
 313          end if
 314          if( usr7_ndx > 0 ) then
 315             if( usr6_ndx > 0 ) then
 316                rxt(:,k,usr7_ndx) = rxt(:,k,usr6_ndx) * exp( -10900.*tinv(:) )/ 2.1e-27
 317             else
 318                rxt(:,k,usr7_ndx) = 0.
 319             end if
 320          end if
 321          if( usr8_ndx > 0 ) then
 322             rxt(:,k,usr8_ndx) = 1.5e-13 * (1. + 6.e-7*boltz*m(:,k)*temp(:,k))
 323          end if
 324
 325 !-----------------------------------------------------------------
 326 !       ... ho2 + ho2 --> h2o2
 327 !       note: this rate involves the water vapor number density
 328 !-----------------------------------------------------------------
 329          if( usr9_ndx > 0 ) then
 330             if( has_h2o ) then
 331                ko(:)   = 2.3e-13 * exp( 600.*tinv(:) )
 332                kinf(:) = 1.7e-33 * m(:,k) * exp( 1000.*tinv(:) )
 333                fc(:)   = 1. + 1.4e-21 * invariants(:,k,h2o_ndx) * exp( 2200.*tinv(:) )
 334                rxt(:,k,usr9_ndx) = (ko(:) + kinf(:)) * fc(:)
 335             else
 336                rxt(:,k,usr9_ndx) = 0.
 337             end if
 338          end if
 339
 340 !-----------------------------------------------------------------
 341 !       ... mco3 + no2 -> mpan
 342 !-----------------------------------------------------------------
 343          if( usr14_ndx > 0 ) then
 344             rxt(:,k,usr14_ndx) = 1.1e-11 * tp(:) / m(:,k)
 345          end if
 346
 347 !-----------------------------------------------------------------
 348 !       ... pan + m --> ch3co3 + no2 + m
 349 !-----------------------------------------------------------------
 350          exp_fac(:) = exp( -14000.*tinv(:) )
 351          if( usr12_ndx > 0 ) then
 352             if( usr11_ndx > 0 ) then
 353                rxt(:,k,usr12_ndx) = rxt(:,k,usr11_ndx) * 1.111e28 * exp_fac(:)
 354             else
 355                rxt(:,k,usr12_ndx) = 0.
 356             end if
 357          end if
 358
 359 !-----------------------------------------------------------------
 360 !       ... mpan + m --> mco3 + no2 + m
 361 !-----------------------------------------------------------------
 362          if( usr15_ndx > 0 ) then
 363             if( usr14_ndx > 0 ) then
 364                rxt(:,k,usr15_ndx) = rxt(:,k,usr14_ndx) * 1.111e28 * exp_fac(:)
 365             else
 366                rxt(:,k,usr15_ndx) = 0.
 367             end if
 368          end if
 369
 370 !-----------------------------------------------------------------
 371 !       ... xooh + oh -> h2o + oh
 372 !-----------------------------------------------------------------
 373          if( usr21_ndx > 0 ) then
 374             rxt(:,k,usr21_ndx) = temp(:,k)**2 * 7.69e-17 * exp( 253.*tinv(:) )
 375          end if
 376
 377 !-----------------------------------------------------------------
 378 !       ... ch3coch3 + oh -> ro2 + h2o
 379 !-----------------------------------------------------------------
 380          if( usr22_ndx > 0 ) then
 381             rxt(:,k,usr22_ndx) = 3.82e-11 * exp( -2000.*tinv(:) ) + 1.33e-13
 382          end if
 383
 384 !-----------------------------------------------------------------
 385 !       ... DMS + OH  --> .5 * SO2
 386 !-----------------------------------------------------------------
 387          if( usr24_ndx > 0 ) then
 388             ko(:) = 1. + 5.5e-31 * exp( 7460.*tinv(:) ) * m(:,k) * 0.21
 389             rxt(:,k,usr24_ndx) = 1.7e-42 * exp( 7810.*tinv(:) ) * m(:,k) * 0.21 / ko(:)
 390          end if
 391 #endif
 392
 393 !-----------------------------------------------------------------
 394 !       ... SO2 + OH  --> SO4  (REFERENCE?? - not Liao)
 395 !-----------------------------------------------------------------
 396          if( usr23_ndx > 0 ) then
 397             fc(:) = 3.0e-31 *(300.*tinv(:))**3.3
 398             ko(:) = fc(:)*m(:,k)/(1. + fc(:)*m(:,k)/1.5e-12)
 399             rxt(:,k,usr23_ndx) = ko(:)*.6**(1. + (log10(fc(:)*m(:,k)/1.5e-12))**2.)**(-1.)
 400          end if
 401
 402 !-----------------------------------------------------------------
 403 !       ... NH3 --> NH4
 404 !-----------------------------------------------------------------
 405          if( usr25_ndx > 0 )  then
 406             rxt(:,k,usr25_ndx) = 0.
 407          end if
 408
 409
 410 !-----------------------------------------------------------------
 411 !       ... exponent for calculating number density
 412 !-----------------------------------------------------------------
 413          n_exp = exp( -4.5*log(sd_sulf)*log(sd_sulf) )
 414
 415          dm_sulf = 2. * rm_sulf
 416          dm_orgc = 2. * rm_orgc
 417          dm_bc   = 2. * rm_bc
 418
 419          log_sd_sulf = log(sd_sulf)
 420          log_sd_orgc = log(sd_orgc)
 421          log_sd_bc = log(sd_bc)
 422
 423 long_loop : &
 424          do i = 1,plonl
 425 !-------------------------------------------------------------------------
 426 !       ... air density (kg/m3)
 427 !-------------------------------------------------------------------------
 428               rho_air = pmid(i,k)/(temp(i,k)*287.04)
 429
 430 !-------------------------------------------------------------------------
 431 !       ... aerosol growth interpolated from M.Chin's table
 432 !-------------------------------------------------------------------------
 433                if (relhum(i,k) >= table_rh(7)) then
 434                  rfac_sulf = table_rfac_sulf(7)
 435                  rfac_oc = table_rfac_oc(7)
 436                  rfac_bc = table_rfac_bc(7)
 437                else
 438                  do irh = 2,7
 439                    if (relhum(i,k) <= table_rh(irh)) then
 440                      exit
 441                    end if
 442                  end do
 443                  rh_l = irh-1
 444                  rh_u = irh
 445
 446                  factor = (relhum(i,k) - table_rh(rh_l))/(table_rh(rh_u) - table_rh(rh_l))
 447
 448                  rfac_sulf = table_rfac_sulf(rh_l) + factor*(table_rfac_sulf(rh_u) - table_rfac_sulf(rh_l))
 449                  rfac_oc = table_rfac_oc(rh_u) + factor*(table_rfac_oc(rh_u) - table_rfac_oc(rh_l))
 450                  rfac_bc = table_rfac_bc(rh_u) + factor*(table_rfac_bc(rh_u) - table_rfac_bc(rh_l))
 451                end if
 452
 453                dm_sulf_wet = dm_sulf * rfac_sulf
 454                dm_orgc_wet = dm_orgc * rfac_oc
 455                dm_bc_wet = dm_bc * rfac_bc
 456
 457 !-------------------------------------------------------------------------
 458 !       ... sulfate aerosols
 459 !-------------------------------------------------------------------------
 460 !-------------------------------------------------------------------------
 461 !       ... use ubvals climatology for stratospheric sulfate surface area density
 462 !-------------------------------------------------------------------------
 463              if( k < ltrop(i) ) then
 464                  sfc_sulf = strato_sad(i,k)
 465              else
 466
 467               if( so4_ndx > 0 ) then
 468 !-------------------------------------------------------------------------
 469 ! convert mass mixing ratio of aerosol to cm3/cm3 (cm^3_aerosol/cm^3_air)
 470 ! v=volume density (m^3/m^3)
 471 ! rho_aer=density of aerosol (kg/m^3)
 472 ! v=m*rho_air/rho_aer   [kg/kg * (kg/m3)_air/(kg/m3)_aer]
 473 !-------------------------------------------------------------------------
 474                 v = mmr(i,k,so4_ndx) * rho_air/rho_sulf
 475 !-------------------------------------------------------------------------
 476 ! calculate the number density of aerosol (aerosols/cm3)
 477 ! assuming a lognormal distribution
 478 ! n  = (aerosols/cm3)
 479 ! dm = geometric mean diameter
 480 !
 481 ! because only the dry mass of the aerosols is known, we
 482 ! use the mean dry radius
 483 !-------------------------------------------------------------------------
 484                 n  = v * (6./pi)*(1./(dm_sulf**3))*n_exp
 485 !-------------------------------------------------------------------------
 486 ! find surface area of aerosols using dm_wet, log_sd
 487 !  (increase of sd due to RH is negligible)
 488 ! and number density calculated above as distribution
 489 ! parameters
 490 ! sfc = surface area of wet aerosols (cm^2/cm^3)
 491 !-------------------------------------------------------------------------
 492                 s_exp    = exp(2.*log_sd_sulf*log_sd_sulf)
 493                 sfc_sulf = n * pi * (dm_sulf_wet**2) * s_exp
 494
 495               else
 496 !-------------------------------------------------------------------------
 497 !  if so4 not simulated, use off-line sulfate and calculate as above
 498 !  convert sulfate vmr to volume density of aerosol (cm^3_aerosol/cm^3_air)
 499 !-------------------------------------------------------------------------
 500                 v = sulfate(i,k) * m(i,k) * mw_so4 / (avo * rho_sulf) *1.e6
 501                 n  = v * (6./pi)*(1./(dm_sulf**3))*n_exp
 502                 s_exp    = exp(2.*log_sd_sulf*log_sd_sulf)
 503                 sfc_sulf = n * pi * (dm_sulf_wet**2) * s_exp
 504
 505               end if
 506              end if
 507
 508 !-------------------------------------------------------------------------
 509 ! ammonium nitrate (follow same procedure as sulfate, using size and density of sulfate)
 510 !-------------------------------------------------------------------------
 511               if( nit_ndx > 0 ) then
 512                 v = mmr(i,k,nit_ndx) * rho_air/rho_sulf
 513                 n  = v * (6./pi)*(1./(dm_sulf**3))*n_exp
 514                 s_exp   = exp(2.*log_sd_sulf*log_sd_sulf)
 515                 sfc_nit = n * pi * (dm_sulf_wet**2) * s_exp
 516               else
 517                 sfc_nit = 0.
 518               end if
 519
 520 !-------------------------------------------------------------------------
 521 ! hydrophylic organic carbon (follow same procedure as sulfate)
 522 !-------------------------------------------------------------------------
 523               if( oc2_ndx > 0 ) then
 524                 v = mmr(i,k,oc2_ndx) * rho_air/rho_orgc
 525                 n  = v * (6./pi)*(1./(dm_orgc**3))*n_exp
 526                 s_exp    = exp(2.*log_sd_orgc*log_sd_orgc)
 527                 sfc_oc   = n * pi * (dm_orgc_wet**2) * s_exp
 528               else
 529                 sfc_oc = 0.
 530               end if
 531
 532 !-------------------------------------------------------------------------
 533 ! secondary organic carbon (follow same procedure as sulfate)
 534 !-------------------------------------------------------------------------
 535               if( soa_ndx > 0 ) then
 536                 v = mmr(i,k,soa_ndx) * rho_air/rho_orgc
 537                 n  = v * (6./pi)*(1./(dm_orgc**3))*n_exp
 538                 s_exp     = exp(2.*log_sd_orgc*log_sd_orgc)
 539                 sfc_soa   = n * pi * (dm_orgc_wet**2) * s_exp
 540               else
 541                 sfc_soa = 0.
 542               end if
 543
 544 !-------------------------------------------------------------------------
 545 ! black carbon (follow same procedure as sulfate)
 546 !-------------------------------------------------------------------------
 547               if( cb2_ndx > 0 ) then
 548                 v = mmr(i,k,cb2_ndx) * rho_air/rho_bc
 549                 n  = v * (6./pi)*(1./(dm_bc**3))*n_exp
 550                 s_exp     = exp(2.*log_sd_bc*log_sd_bc)
 551                 sfc_bc   = n * pi * (dm_bc_wet**2) * s_exp
 552               else
 553                 sfc_bc = 0.
 554               end if
 555
 556 !-------------------------------------------------------------------------
 557 !       ... add up total surface area density for output
 558 !-------------------------------------------------------------------------
 559               sad_total(i,k) = sfc_sulf + sfc_nit + sfc_oc + sfc_soa + sfc_bc
 560
 561
 562 !-------------------------------------------------------------------------
 563 !  Heterogeneous reaction rates for uptake of a gas on an aerosol:
 564 !    rxt = sfc / ( (rad_aer/Dg_gas) + (4/(c_gas*gamma_gas)))
 565 !-------------------------------------------------------------------------
 566 !-------------------------------------------------------------------------
 567 !       ... n2o5 -> 2 hno3  (on sulfate, nh4no3, oc2, soa)
 568 !-------------------------------------------------------------------------
 569               if( usr16_ndx > 0 ) then
 570                 c_n2o5 = 1.40e3 * sqrt_t(i)         ! mean molecular speed of n2o5
 571
 572                 rxt_sulf = sfc_sulf / (0.5*dm_sulf_wet/dg + (4./(c_n2o5*gamma_n2o5)))
 573                 rxt_nit = sfc_nit / (0.5*dm_sulf_wet/dg + (4./(c_n2o5*gamma_n2o5)))
 574                 rxt_oc = sfc_oc / (0.5*dm_orgc_wet/dg + (4./(c_n2o5*gamma_n2o5)))
 575                 rxt_soa = sfc_soa / (0.5*dm_orgc_wet/dg + (4./(c_n2o5*gamma_n2o5)))
 576
 577                 rxt(i,k,usr16_ndx) = rxt_sulf + rxt_nit + rxt_oc + rxt_soa
 578
 579               end if
 580 !-------------------------------------------------------------------------
 581 !       ... no3 -> hno3  (on sulfate, nh4no3, oc, soa)
 582 !-------------------------------------------------------------------------
 583               if( usr17_ndx > 0 ) then
 584                 c_no3 = 1.85e3 * sqrt_t(i)         ! mean molecular speed of no3
 585
 586                 rxt_sulf = sfc_sulf / (0.5*dm_sulf_wet/dg + (4./(c_no3*gamma_no3)))
 587                 rxt_nit = sfc_nit / (0.5*dm_sulf_wet/dg + (4./(c_no3*gamma_no3)))
 588                 rxt_oc = sfc_oc / (0.5*dm_orgc_wet/dg + (4./(c_no3*gamma_no3)))
 589                 rxt_soa = sfc_soa / (0.5*dm_orgc_wet/dg + (4./(c_no3*gamma_no3)))
 590
 591                 rxt(i,k,usr17_ndx) = rxt_sulf + rxt_nit + rxt_oc + rxt_soa
 592               end if
 593 !-------------------------------------------------------------------------
 594 !       ... no2 -> 0.5 * (ho+no+hno3)  (on sulfate, nh4no3, oc2, soa)
 595 !-------------------------------------------------------------------------
 596               if( usr17a_ndx > 0 ) then
 597                 c_no2 = 2.15e3 * sqrt_t(i)         ! mean molecular speed of no2
 598
 599                 rxt_sulf = sfc_sulf / (0.5*dm_sulf_wet/dg + (4./(c_no2*gamma_no2)))
 600                 rxt_nit = sfc_nit / (0.5*dm_sulf_wet/dg + (4./(c_no2*gamma_no2)))
 601                 rxt_oc = sfc_oc / (0.5*dm_orgc_wet/dg + (4./(c_no2*gamma_no2)))
 602                 rxt_soa = sfc_soa / (0.5*dm_orgc_wet/dg + (4./(c_no2*gamma_no2)))
 603
 604                 rxt(i,k,usr17a_ndx) = rxt_sulf + rxt_nit + rxt_oc + rxt_soa
 605               end if
 606 !-------------------------------------------------------------------------
 607 !       ... ho2 -> 0.5 * h2o2  (on sulfate, nh4no3, oc2, soa)
 608 !-------------------------------------------------------------------------
 609               if( usr26_ndx > 0 ) then
 610                 c_ho2 = 2.53e3 * sqrt_t(i)         ! mean molecular speed of ho2
 611
 612                 rxt_sulf = sfc_sulf / (0.5*dm_sulf_wet/dg + (4./(c_ho2*gamma_ho2)))
 613                 rxt_nit = sfc_nit / (0.5*dm_sulf_wet/dg + (4./(c_ho2*gamma_ho2)))
 614                 rxt_oc = sfc_oc / (0.5*dm_orgc_wet/dg + (4./(c_ho2*gamma_ho2)))
 615                 rxt_soa = sfc_soa / (0.5*dm_orgc_wet/dg + (4./(c_ho2*gamma_ho2)))
 616
 617                 rxt(i,k,usr26_ndx) = rxt_sulf + rxt_nit + rxt_oc + rxt_soa
 618               end if
 619             end do long_loop
 620       end do level_loop
 621
 622       end subroutine usrrxt
 623
 624       end module mo_usrrxt