chem/module_mosaic_newnucb.F

   1         module module_mosaic_newnucb
   2
   3 !-----------------------------------------------------------------------
   4 ! 23-may-2008 - created
   5 !       subr mosaic_newnuc_1box & wexler_nuc_mosaic_1box taken
   6 !           from wrfchem (amt version) module_mosaic_newnuc.F
   7 !       subr subr mer07_veh02_nuc_mosaic_1box, binary_nuc_vehk2002,
   8 !           & ternary_nuc_merik2007 taken from cam3 (modal_aero version)
   9 !           modal_aero_newnuc.F90
  10 !       mosaic_newnuc_1box required significant changes
  11 !       other routines required no or minor changes
  12 !-----------------------------------------------------------------------
  13
  14
  15
  16         use module_data_mosaic_kind, only:  r8
  17         use module_data_mosaic_main, only:  pi !BSINGH - 05/28/2013(RCE updates)
  18         use module_peg_util
  19
  20
  21
  22         implicit none
  23         !BSINGH - 05/28/2013(RCE updates)
  24         ! adjustment factors
  25         real(r8), parameter, public :: adjust_factor_dnaitdt = 1.0_r8  ! applied to final dnait/dt
  26         real(r8), parameter, public :: adjust_factor_bin_tern_ratenucl = 1.0_r8  ! applied to binary/ternary nucleation rate
  27         real(r8), parameter, public :: adjust_factor_pbl_ratenucl = 1.0_r8  ! applied to boundary layer nucleation rate
  28         !BSINGH - 05/28/2013(RCE updates ENDS)
  29
  30         contains
  31
  32
  33
  34 !-----------------------------------------------------------------------
  35         subroutine mosaic_newnuc_1box(   &
  36                 istat_newnuc, idiagbb_in,  &
  37                 itstep, itype_newnuc, dtchem,   &
  38                 temp_box, patm_box, cair_box, rh_box,   &
  39                 dens_nh4so4a_in, mw_so4a, mw_nh4a, mw_h2o, mw_air,   &
  40                 fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam, fact_gasmr,   &
  41                 rbox0, rbox,   &
  42                 adrydens_box, awetdens_box,   &
  43                 adrydpav_box, awetdpav_box, adryqmas_box )
  44 !
  45 !   calculates new particle nucleation for a single grid box
  46 !       over timestep dtchem
  47 !   works with mosaic sectional aerosol packages
  48 !
  49 !   uses the following nucleation parameterizations
  50 !       mnewnuc_flag1 = 1 -- h2so4-nh3-h2o ternary nuc. of merikanto et al. (2007)
  51 !       mnewnuc_flag1 = 2 -- h2so4-h2o binary nuc. of vehkamaki et al. (2002)
  52 !       mnewnuc_flag1 = 3 -- h2so4-h2o binary nuc. of wexler et al. (1994)
  53
  54 !BSINGH - 05/28/2013(RCE updates)
  55 !       mnewnuc_flag1 = 11 -- empirical first  order boundary layer nucleation
  56 !       mnewnuc_flag1 = 12 -- empirical second order boundary layer nucleation
  57 !BSINGH - 05/28/2013(RCE updates ENDS)
  58         use module_data_mosaic_main, only:  ntot_used, piover6, third
  59         use module_data_mosaic_boxmod, only:  kh2so4, knh3, kso2
  60         use module_data_mosaic_aero, only:  mnewnuc_flag1
  61         use module_data_mosaic_asecthp, only:   &
  62                 ai_phase, dens_aer, dens_water_aer,   &
  63                 dcen_sect, dlo_sect, dhi_sect,   &
  64                 hyswptr_aer, lptr_nh4_aer, lptr_so4_aer, lunerr,   &
  65                 massptr_aer, maxd_asize, maxd_atype,   &
  66                 ncomp_aer, nsize_aer, numptr_aer, smallmassbb,   &
  67                 volumcen_sect, volumlo_sect, volumhi_sect, waterptr_aer
  68
  69
  70 !   subr arguments
  71         integer,  intent(inout) :: istat_newnuc   ! =0 if no problems
  72         integer,  intent(in)    :: idiagbb_in     ! controls diagnostics
  73         integer,  intent(in)    :: itstep         ! host code time step index
  74         integer,  intent(in)    :: itype_newnuc   ! "type" that gets the new particles
  75         real(r8), intent(in)    :: dtchem         ! time step for nucleation
  76         real(r8), intent(in)    :: temp_box       ! air temp (K)
  77         real(r8), intent(in)    :: patm_box       ! air pressure (atm)
  78         real(r8), intent(in)    :: cair_box       ! air molar density (mol/cm^3)
  79         real(r8), intent(in)    :: rh_box         ! relative humidity (0-1)
  80
  81         real(r8), intent(in) :: dens_nh4so4a_in   ! nh4so4 density (units same as dens_aer)
  82         real(r8), intent(in) :: mw_so4a, mw_nh4a, mw_h2o, mw_air   ! molec wghts (g/mol)
  83
  84         real(r8), intent(in) :: fact_apmassmr, fact_apnumbmr, &
  85                                 fact_apdens, fact_apdiam, fact_gasmr
  86 !   the module_data_mosaic_asecthp densities are such that
  87 !       dens_aer*fact_apdens give (g/cm^3)
  88 !   the module_data_mosaic_asecthp diameters are such that
  89 !       dlo/hi_sect*fact_apdiam give (cm)
  90 !   the module_data_mosaic_asecthp volumes are such that
  91 !       volumlo/hi_sect**fact_apdiam**3) give (cm^3)
  92 !   aerosol mass mixing ratios in rbox(massptr_aer(ll,n,itype,iphase))
  93 !       have units for which rbox*fact_apmassmr gives (g-AP/g-air)
  94 !   aerosol number mixing ratios in rbox(numptr_aer(n,itype,iphase))
  95 !       have units for which rbox*fact_apnumbmr gives (#/g-air)
  96 !   trace gas mixing ratios in rbox(...)
  97 !       have units for which rbox*fact_gasmr gives (mol/mol-air)
  98
  99 ! rbox0 holds mixrat values before gas-aerosol mass-transfer
 100 ! rbox  holds mixrat values after  gas-aerosol mass-transfer and movesect
 101         real(r8), intent(in)    :: rbox0(ntot_used)
 102         real(r8), intent(inout) :: rbox(ntot_used)
 103
 104 ! adrydens_box = aerosol dry density (units same as dens_aer)
 105 ! awetdens_box = aerosol wet density (units same as dens_aer)
 106         real(r8), intent(inout) :: adrydens_box(maxd_asize,maxd_atype)
 107         real(r8), intent(inout) :: awetdens_box(maxd_asize,maxd_atype)
 108 ! adrydpav_box = aerosol mean dry diameter (units same as dlo_sect)
 109 ! awetdpav_box = aerosol mean wet diameter (units same as dlo_sect)
 110         real(r8), intent(inout) :: adrydpav_box(maxd_asize,maxd_atype)
 111         real(r8), intent(inout) :: awetdpav_box(maxd_asize,maxd_atype)
 112 ! adryqmas_box = aerosol total-dry-mass mixing ratio (units same as rbox)
 113         real(r8), intent(inout) :: adryqmas_box(maxd_asize,maxd_atype)
 114
 115
 116 !   local variables
 117         integer, parameter :: p1st = 1
 118
 119         integer :: it=1, jt=1, k=1
 120         integer :: l, ll
 121         integer :: isize, itype, iphase
 122         integer :: iconform_numb
 123         integer :: idiagbb
 124         integer :: ldiagaa, lundiagaa
 125         integer :: nsize
 126         integer, save :: ncount(10)
 127
 128 ! min h2so4 vapor for nuc calcs = 4.0e-16 mol/mol-air ~= 1.0e4 molecules/cm3,
 129         real(r8), parameter :: qh2so4_cutoff = 4.0e-16_r8
 130
 131 ! nh4hso4 values for a_zsr and b_zsr
 132         real(r8), parameter :: a_zsr_xx1 =  1.15510
 133         real(r8), parameter :: a_zsr_xx2 = -3.20815
 134         real(r8), parameter :: a_zsr_xx3 =  2.71141
 135         real(r8), parameter :: a_zsr_xx4 =  2.01155
 136         real(r8), parameter :: a_zsr_xx5 = -4.71014
 137         real(r8), parameter :: a_zsr_xx6 =  2.04616
 138         real(r8), parameter :: b_zsr_xx  = 29.4779
 139
 140         real(r8) :: aw
 141         real(r8) :: densdefault
 142         real(r8) :: dens_nh4so4a, dens_nh4so4a_kgm3
 143         real(r8) :: dtnuc
 144         real(r8) :: duma, dumb, dumc
 145         real(r8) :: fact_apvolu
 146         real(r8) :: h2so4_uptkrate
 147         real(r8) :: press_pa
 148         real(r8) :: qh2so4_avg, qh2so4_cur, qh2so4_del  ! [mol/mol-air]
 149         real(r8) :: qnh3_avg, qnh3_cur, qnh3_del        ! [mol/mol-air]
 150         real(r8) :: qso4a_del, qnh4a_del                ! [mol/mol-air]
 151         real(r8) :: qnuma_del                           ! [#/mol-air]
 152         real(r8) :: tmpa, tmpb, tmpc, tmp_q2, tmp_q3
 153         real(r8) :: wwdens, wwdpav, wwmass, wwvolu
 154         real(r8) :: xxdens, xxdpav, xxmass, xxnumb, xxvolu
 155
 156         real(r8),save :: tmpveca(10), tmpvecb(10), tmpvecc(10), tmpvecd(10), tmpvece(10), tmpvecf(10) !BSINGH - 05/28/2013(RCE updates - changed 'dum' to 'tmp')
 157         real(r8) :: dplom_nuc(maxd_asize), dphim_nuc(maxd_asize)
 158         real(r8) :: volumlo_nuc(maxd_asize), volumhi_nuc(maxd_asize)
 159
 160         character(len=100) :: msg
 161
 162
 163         if ( kh2so4 < 1 .or. kh2so4 > ntot_used .or. &
 164              knh3   < 1 .or. knh3   > ntot_used .or. &
 165              kso2   < 1 .or. kso2   > ntot_used ) then
 166            write(msg,'(a,3i15)') &
 167               'module_mosaic_newnucb bad kh2so4/knh3/kso2', &
 168               kh2so4, knh3, kso2
 169            call peg_message( lunerr, msg )
 170            call peg_error_fatal( lunerr, msg )
 171            call wrf_error_fatal(trim(adjustl(msg)))
 172         end if
 173
 174 !   check mnewnuc_flag1
 175         istat_newnuc = 0
 176         if (mnewnuc_flag1 <= 0) then !BSINGH - 05/28/2013(RCE updates)
 177            !BSINGH - 05/28/2013(RCE updates)
 178            return
 179         else if ( mnewnuc_flag1 ==  1 .or. &
 180              mnewnuc_flag1 ==  2 .or. &
 181              mnewnuc_flag1 ==  3 .or. &
 182              mnewnuc_flag1 == 11 .or. &
 183              mnewnuc_flag1 == 12 ) then
 184            continue
 185         else
 186            !BSINGH - 05/28/2013(RCE updates ENDS)
 187 !           if ((it .eq. its) .and. (jt .eq. jts))   &
 188                 call peg_message( lunerr,   &
 189                 '*** mosaic_newnuc_1box -- illegal mnewnuc_flag1' )
 190             istat_newnuc = -1
 191             return
 192             !BSINGH - 05/28/2013(RCE updates- got rid of an 'elseif' here)
 193         end if
 194
 195
 196 !   set variables that do not change
 197         dtnuc = dtchem
 198         idiagbb = idiagbb_in
 199         lundiagaa = 93!BSINGH - 05/28/2013(RCE updates)
 200
 201         itype = itype_newnuc
 202         iphase = ai_phase
 203         nsize = nsize_aer(itype)
 204         fact_apvolu = fact_apdiam*fact_apdiam*fact_apdiam
 205         densdefault = dens_aer(1,itype)*fact_apdens
 206         volumlo_nuc(1:nsize) = volumlo_sect(1:nsize,itype)*fact_apvolu
 207         volumhi_nuc(1:nsize) = volumhi_sect(1:nsize,itype)*fact_apvolu
 208         dplom_nuc(1:nsize) = dlo_sect(1:nsize,itype)*0.01*fact_apdiam   ! convert cm to m
 209         dphim_nuc(1:nsize) = dhi_sect(1:nsize,itype)*0.01*fact_apdiam   ! convert cm to m
 210
 211
 212 !   loop over subareas (currently only 1) and vertical levels
 213 !       do 2900 m = 1, nsubareas
 214 !       do 2800 k = kclm_calcbgn, kclm_calcend
 215
 216
 217 !   initialize diagnostics
 218 !       if ((it .eq. its) .and.   &
 219 !           (jt .eq. jts) .and. (k .eq. kclm_calcbgn)) then
 220             tmpveca(:) = 0.0         ! current grid param values
 221             tmpvecb(:) = +1.0e35     ! param minimums
 222             tmpvecc(:) = -1.0e35     ! param maximums
 223             tmpvecd(:) = 0.0         ! param averages
 224             tmpvece(:) = 0.0         ! param values for highest qnuma_del
 225             ncount(:) = 0
 226 !       end if
 227
 228
 229         ncount(1) = ncount(1) + 1
 230
 231         qh2so4_cur = max(0.0_r8,rbox(kh2so4))*fact_gasmr
 232 !   skip if h2so4 vapor < qh2so4_cutoff
 233         if (qh2so4_cur <= qh2so4_cutoff) goto 2700
 234
 235         qnh3_cur   = max(0.0_r8,rbox(knh3))*fact_gasmr
 236         qh2so4_avg = 0.5*( qh2so4_cur + max(0.0_r8,rbox0(kh2so4))*fact_gasmr )
 237         qnh3_avg   = 0.5*( qnh3_cur   + max(0.0_r8,rbox0(knh3))*fact_gasmr )
 238
 239
 240 ! calc h2so4_uptkrate = h2so4 first-order uptake rate to existing aerosol
 241 !                     = -d[ln(qh2so4)]/dt  (units = 1/s)
 242         tmp_q3 = qh2so4_cur
 243 !   tmp_q2 = qh2so4 before aerosol uptake
 244         tmp_q2 = max(0.0_r8,rbox0(kh2so4))*fact_gasmr
 245 !   following gets tmpb = log( max( 1, tmp_q2/tmp_q3 ) )
 246 !   but with some checks added to avoid floating point exceptions
 247         if (tmp_q2 > tmp_q3) then
 248            tmpc = tmp_q2 * exp( -20.0_r8 )
 249            if (tmp_q3 <= tmpc) then
 250               tmpb = 20.0_r8
 251            else
 252               tmpb = log( tmp_q2/tmp_q3 )
 253            end if
 254         else
 255            tmpb = 0.0
 256         end if
 257         h2so4_uptkrate = tmpb/dtchem
 258
 259
 260         qh2so4_del = 0.0
 261         qnh3_del = 0.0
 262         qnuma_del = 0.0
 263         qso4a_del = 0.0
 264         qnh4a_del = 0.0
 265
 266 !       dens_nh4so4a = dens_so4_aer
 267 !       dens_nh4so4a = dens_aer_mac(iso4_a)
 268         dens_nh4so4a = dens_nh4so4a_in*fact_apdens
 269
 270         isize = 0
 271
 272 !   make call to nucleation routine
 273         if (mnewnuc_flag1 /= 3) then!BSINGH - 05/28/2013(RCE updates)
 274
 275 !  new ternary (from cam3 code)
 276             ldiagaa = idiagbb!BSINGH - 05/28/2013(RCE updates)
 277             press_pa = patm_box*1.01325e5   ! convert atmos to pa
 278             dens_nh4so4a_kgm3 = dens_nh4so4a*1.0e3
 279
 280             call mer07_veh02_nuc_mosaic_1box(   &
 281                mnewnuc_flag1, dtnuc, temp_box, rh_box, press_pa,   &
 282                qh2so4_cur, qh2so4_avg, qnh3_cur, qnh3_avg, h2so4_uptkrate,   &
 283                nsize, maxd_asize, dplom_nuc, dphim_nuc,   &
 284                isize, qnuma_del, qso4a_del, qnh4a_del,   &
 285                qh2so4_del, qnh3_del, dens_nh4so4a_kgm3, &
 286                itstep, ldiagaa, lundiagaa )!BSINGH - 05/28/2013(RCE updates)
 287 !           subr mer07_veh02_nuc_mosaic_1box(   &
 288 !              newnuc_method_flagaa, dtnuc, temp_in, rh_in, press_in,   &
 289 !              qh2so4_cur, qh2so4_avg, qnh3_cur, h2so4_uptkrate,   &
 290 !              nsize, maxd_asize, dplom_sect, dphim_sect,   &
 291 !              isize_nuc, qnuma_del, qso4a_del, qnh4a_del,   &
 292 !              qh2so4_del, qnh3_del, dens_nh4so4a, &
 293 !              itstep, ldiagaa, lundiagaa )!BSINGH - 05/28/2013(RCE updates)
 294
 295             dens_nh4so4a = dens_nh4so4a_kgm3/1.0e3
 296
 297         else! (mnewnuc_flag1 == 3)!BSINGH - 05/28/2013(RCE updates)
 298
 299             call wexler_nuc_mosaic_1box(   &
 300                dtnuc, temp_box, rh_box, cair_box,   &
 301                qh2so4_avg, qh2so4_cur, qnh3_avg, qnh3_cur,   &
 302                nsize, maxd_asize, volumlo_nuc, volumhi_nuc,   &
 303                isize, qnuma_del, qso4a_del, qnh4a_del,   &
 304                qh2so4_del, qnh3_del, dens_nh4so4a )
 305
 306         !BSINGH - 05/28/2013(RCE updates- got rid of an 'else' here)
 307         end if
 308
 309
 310 !   temporary diagnostics
 311         tmpveca(1) = temp_box
 312         tmpveca(2) = rh_box
 313         tmpveca(3) = rbox(kso2)*fact_gasmr
 314         tmpveca(4) = qh2so4_avg
 315         tmpveca(5) = qnh3_avg
 316         tmpveca(6) = qnuma_del
 317         do l = 1, 6
 318             tmpvecb(l) = min( tmpvecb(l), tmpveca(l) )
 319             tmpvecc(l) = max( tmpvecc(l), tmpveca(l) )
 320             tmpvecd(l) = tmpvecd(l) + tmpveca(l)
 321             if (qnuma_del .gt. tmpvece(6)) tmpvece(l) = tmpveca(l)
 322         end do
 323
 324
 325 !   check for zero new particles
 326         if (qnuma_del .le. 0.0) goto 2700
 327
 328 !   check for valid isize
 329         if (isize .ne. 1) ncount(3) = ncount(3) + 1
 330         if ((isize .lt. 1) .or. (isize .gt. nsize)) then
 331             write(msg,93010) 'newnucxx bad isize_nuc' , it, jt, k,   &
 332                 isize, nsize
 333             call peg_message( lunerr, msg )
 334             goto 2700
 335         end if
 336 93010   format( a, 3i3, 1p, 9e10.2 )
 337
 338
 339         ncount(2) = ncount(2) + 1
 340
 341 !   update gas and aerosol so4 and nh3/nh4 mixing ratios
 342         rbox(kh2so4) = max( 0.0_r8, rbox(kh2so4) + qh2so4_del/fact_gasmr )
 343         rbox(knh3  ) = max( 0.0_r8, rbox(knh3  ) + qnh3_del/fact_gasmr )
 344
 345         l = lptr_so4_aer(isize,itype,iphase)
 346         if (l .ge. p1st) then
 347 !           rbox(l) = rbox(l) + qso4a_del
 348             tmpa = qso4a_del*mw_so4a/mw_air   ! g/g-air
 349             tmpvecf(1) = rbox(l)!BSINGH - 05/28/2013(RCE updates)
 350             rbox(l) = rbox(l) + tmpa/fact_apmassmr
 351             !BSINGH - 05/28/2013(RCE updates)
 352             tmpvecf(2) = rbox(l)
 353             !              (convert whatever to g/g-air) * (then to ug/g-air) * (then to ug/m3-air)
 354             tmpvecf(1:2) = tmpvecf(1:2) * fact_apmassmr  * 1.0e6              * (1.0e6_r8*cair_box*mw_air)
 355             tmpvecf(3) = tmpvecf(2) - tmpvecf(1)
 356             if (idiagbb >= 10) &
 357                  write(lundiagaa,'(a,1p,3e12.4)') 'so4 (ug/m3) o/n/d', tmpvecf(1:3)
 358             !BSINGH - 05/28/2013(RCE updates ENDS)
 359         end if
 360         l = lptr_nh4_aer(isize,itype,iphase)
 361         if (l .ge. p1st) then
 362 !           rbox(l) = rbox(l) + qnh4a_del
 363             tmpa = qnh4a_del*mw_nh4a/mw_air   ! g/g-air
 364             tmpvecf(1) = rbox(l)!BSINGH - 05/28/2013(RCE updates)
 365             rbox(l) = rbox(l) + tmpa/fact_apmassmr
 366             !BSINGH - 05/28/2013(RCE updates)
 367             tmpvecf(2) = rbox(l)
 368             !              (convert whatever to g/g-air) * (then to ug/g-air) * (then to ug/m3-air)
 369             tmpvecf(1:2) = tmpvecf(1:2) * fact_apmassmr  * 1.0e6              * (1.0e6_r8*cair_box*mw_air)
 370             tmpvecf(3) = tmpvecf(2) - tmpvecf(1)
 371             if (idiagbb >= 10) &
 372             write(lundiagaa,'(a,1p,3e12.4)') 'nh4 (ug/m3) o/n/d', tmpvecf(1:3)
 373             !BSINGH - 05/28/2013(RCE updates ENDS)
 374         end if
 375         l = numptr_aer(isize,itype,iphase)
 376 !       rbox(l) = rbox(l) + qnuma_del
 377         tmpa = qnuma_del/mw_air   ! #/g-air
 378         tmpvecf(1) = rbox(l)!BSINGH - 05/28/2013(RCE updates)
 379         rbox(l) = rbox(l) + tmpa/fact_apnumbmr
 380         !BSINGH - 05/28/2013(RCE updates)
 381         tmpvecf(2) = rbox(l)
 382         !              (convert whatever to #/g-air) * (then to #/cm3-air)
 383         tmpvecf(1:2) = tmpvecf(1:2) * fact_apnumbmr  * (cair_box*mw_air)
 384         tmpvecf(3) = tmpvecf(2) - tmpvecf(1)
 385         if (idiagbb >= 10) &
 386         write(lundiagaa,'(a,1p,3e12.4)') 'num (#/cm3) o/n/d', tmpvecf(1:3)
 387         !BSINGH - 05/28/2013(RCE updates ENDS)
 388         xxnumb = rbox(l)*fact_apnumbmr
 389
 390 !   update aerosol water, using mosaic parameterizations for nh4hso4
 391 !     duma = (mole-salt)/(mole-salt+water)
 392 !     dumb = (mole-salt)/(kg-water)
 393 !     dumc = (mole-water)/(mole-salt)
 394         l = waterptr_aer(isize,itype)
 395         if ((rh_box .gt. 0.10) .and. (l .ge. p1st)) then
 396             aw = min( rh_box, 0.98_r8 )
 397             if (aw .lt. 0.97) then
 398                 duma =       a_zsr_xx1 +   &
 399                         aw*( a_zsr_xx2 +   &
 400                         aw*( a_zsr_xx3 +   &
 401                         aw*( a_zsr_xx4 +   &
 402                         aw*( a_zsr_xx5 +   &
 403                         aw*  a_zsr_xx6 ))))
 404             else
 405                 dumb = -b_zsr_xx*log(aw)
 406                 dumb = max( dumb, 0.5_r8 )
 407                 duma = 1.0/(1.0 + 55.509/dumb)
 408             end if
 409             duma = max( duma, 0.01_r8 )
 410             dumc = (1.0 - duma)/duma
 411 !           rbox(l) = rbox(l) + qso4a_del*dumc
 412             tmpa = qso4a_del*dumc*mw_h2o/mw_air
 413             rbox(l) = rbox(l) + tmpa/fact_apmassmr
 414         end if
 415
 416
 417 !
 418 !   update dry mass, density, and volume,
 419 !   and check for mean dry-size within bounds
 420 !
 421         xxmass = adryqmas_box(isize,itype)*fact_apmassmr
 422         xxdens = adrydens_box( isize,itype)*fact_apdens
 423         iconform_numb = 1
 424
 425         if ((xxdens .lt. 0.1) .or. (xxdens .gt. 20.0)) then
 426 !   (exception) case of drydensity not valid
 427             continue
 428         else
 429 !   (normal) case of drydensity valid (which means drymass is valid also)
 430 !   so increment mass and volume with the so4 & nh4 deltas, then calc density
 431             xxvolu = xxmass/xxdens
 432 !           duma = qso4a_del*mw_aer_mac(iso4_a)   &
 433 !                + qnh4a_del*mw_aer_mac(inh4_a)
 434             duma = (qso4a_del*mw_so4a + qnh4a_del*mw_nh4a)/mw_air   ! g/g-air
 435             xxmass = xxmass + duma
 436             xxvolu  = xxvolu  + duma/dens_nh4so4a
 437             if (xxmass .le. smallmassbb) then
 438 !   do this to force calc of dry mass, volume from rbox
 439                 xxdens = 0.001
 440             else if (xxmass .gt. 1000.0*xxvolu) then
 441 !   in this case, density is too large.  setting density=1000 forces
 442 !   next IF block while avoiding potential divide by zero or overflow
 443                 xxdens = 1000.0
 444             else
 445                 xxdens = xxmass/xxvolu
 446             end if
 447         end if
 448
 449         if ((xxdens .lt. 0.1) .or. (xxdens .gt. 20.0)) then
 450 !   (exception) case of drydensity not valid (or drymass extremely small),
 451 !   so compute from dry mass, volume from rbox
 452             ncount(4) = ncount(4) + 1
 453             xxmass = 0.0
 454             xxvolu  = 0.0
 455             do ll = 1, ncomp_aer(itype)
 456                 l = massptr_aer(ll,isize,itype,iphase)
 457                 if (l .ge. p1st) then
 458 !                   duma = max( 0.0_r8, rbox(l) )*mw_aer(ll,itype)
 459 !                   xxmass = xxmass + duma
 460 !                   xxvolu = xxvolu + duma/dens_aer(ll,itype)
 461                     duma = max( 0.0_r8, rbox(l) )*fact_apmassmr
 462                     xxmass = xxmass + duma
 463                     xxvolu = xxvolu + duma/(dens_aer(ll,itype)*fact_apdens)
 464                 end if
 465             end do
 466         end if
 467
 468         if (xxmass .le. smallmassbb) then
 469 !   when drymass extremely small, use default density and bin center size,
 470 !   and zero out water
 471             ncount(5) = ncount(5) + 1
 472             xxdens = densdefault
 473             xxvolu = xxmass/xxdens
 474             xxnumb = xxmass/(volumcen_sect(isize,itype)*fact_apvolu*xxdens)
 475             xxdpav = dcen_sect(isize,itype)*fact_apdiam
 476             wwdens = xxdens
 477             wwdpav = xxdpav
 478             iconform_numb = 0
 479             l = waterptr_aer(isize,itype)
 480             if (l .ge. p1st) rbox(l) = 0.0
 481             l = hyswptr_aer(isize,itype)
 482             if (l .ge. p1st) rbox(l) = 0.0
 483         else
 484             xxdens = xxmass/xxvolu
 485         end if
 486
 487         if (iconform_numb .gt. 0) then
 488 !   check for mean dry-size within bounds, and conform number if not
 489             if ( xxnumb .gt. xxvolu/volumlo_nuc(isize) ) then
 490                 ncount(6) = ncount(6) + 1
 491                 xxnumb = xxvolu/volumlo_nuc(isize)
 492                 xxdpav = dlo_sect(isize,itype)*fact_apdiam
 493             else if ( xxnumb .lt. xxvolu/volumhi_nuc(isize) ) then
 494                 ncount(7) = ncount(7) + 1
 495                 xxnumb = xxvolu/volumhi_nuc(isize)
 496                 xxdpav = dhi_sect(isize,itype)*fact_apdiam
 497             else
 498                 xxdpav = (xxvolu/(xxnumb*piover6))**third
 499             end if
 500
 501             tmpb = 0.0
 502             l = waterptr_aer(isize,itype)
 503 !           if (l .ge. p1st) tmpb = max(0.0_r8,rbox(l))*mw_water_aer
 504             if (l .ge. p1st) tmpb = max(0.0_r8,rbox(l))*fact_apmassmr
 505             wwmass = xxmass + tmpb
 506             wwvolu = xxvolu + tmpb/(dens_water_aer*fact_apdens)
 507             wwdens = wwmass/wwvolu
 508             wwdpav = xxdpav*((wwvolu/xxvolu)**third)
 509         end if
 510
 511 !   load dry mass, density, volume, and (possibly conformed) number
 512         l = numptr_aer(isize,itype,iphase)
 513         rbox(l) = xxnumb/fact_apnumbmr
 514         adrydens_box(isize,itype) = xxdens/fact_apdens
 515         adrydpav_box(isize,itype) = xxdpav/fact_apdiam
 516         adryqmas_box(isize,itype) = xxmass/fact_apmassmr
 517         awetdens_box(isize,itype) = wwdens/fact_apdens
 518         awetdpav_box(isize,itype) = wwdpav/fact_apdiam
 519
 520
 521 2700    continue
 522
 523 !   temporary diagnostics
 524         if (idiagbb .ge. 100) then
 525 !       if ((idiagbb .ge. 100) .and.   &
 526 !           (it .eq. ite) .and.    &
 527 !           (jt .eq. jte) .and. (k .eq. kclm_calcend)) then
 528             !BSINGH - 05/28/2013(RCE updates - got rid of lundiagaa)
 529             write(lundiagaa,'(2a)')   &
 530                 'newnucbb names  temp      rh        ',   &
 531                 'so2       h2so4_avg nh3_avg   numa_del'
 532             if (idiagbb .ge. 110) then
 533               write(lundiagaa,93020) 'newnucbb mins ', tmpvecb(1:6)
 534               write(lundiagaa,93020) 'newnucbb maxs ', tmpvecc(1:6)
 535               duma = max( 1, ncount(1) )
 536               write(lundiagaa,93020) 'newnucbb avgs ', tmpvecd(1:6)/duma
 537               write(lundiagaa,93020) 'newnucbb hinuc', tmpvece(1:6)
 538               write(lundiagaa,93020) 'newnucbb dtnuc', dtnuc
 539             end if
 540             write(lundiagaa,93030) 'newnucbb ncnt ', ncount(1:7)
 541         end if
 542 93020   format( a, 1p, 10e10.2 )
 543 93030   format( a, 1p, 10i10 )
 544
 545
 546 !2800   continue        ! k levels
 547 !2900   continue        ! subareas
 548
 549
 550         return
 551         end subroutine mosaic_newnuc_1box
 552
 553
 554
 555 !-----------------------------------------------------------------------
 556 !-----------------------------------------------------------------------
 557         subroutine mer07_veh02_nuc_mosaic_1box(   &
 558            newnuc_method_flagaa, dtnuc, temp_in, rh_in, press_in,   &
 559            qh2so4_cur, qh2so4_avg, qnh3_cur, qnh3_avg, h2so4_uptkrate,   &
 560            nsize, maxd_asize, dplom_sect, dphim_sect,   &
 561            isize_nuc, qnuma_del, qso4a_del, qnh4a_del,   &
 562            qh2so4_del, qnh3_del, dens_nh4so4a, &
 563            itstep, ldiagaa, lundiagaa ) !BSINGH - 05/28/2013(RCE updates)
 564 !.......................................................................
 565 !
 566 ! calculates new particle production from homogeneous nucleation
 567 !    over timestep dtnuc, using nucleation rates from either
 568 !    merikanto et al. (2007) h2so4-nh3-h2o ternary parameterization
 569 !    vehkamaki et al. (2002) h2so4-h2o binary parameterization
 570 !
 571 ! the new particles are "grown" to the lower-bound size of the host code's
 572 !    smallest size bin.  (this "growth" is somewhat ad hoc, and would not be
 573 !    necessary if the host code's size bins extend down to ~1 nm.)
 574 !
 575 !    if the h2so4 and nh3 mass mixing ratios (mixrats) of the grown new
 576 !    particles exceed the current gas mixrats, the new particle production
 577 !    is reduced so that the new particle mass mixrats match the gas mixrats.
 578 !
 579 !    the correction of kerminen and kulmala (2002) is applied to account
 580 !    for loss of the new particles by coagulation as they are
 581 !    growing to the "host code mininum size"
 582 !
 583 ! revision history
 584 !    coded by rc easter, pnnl, xx-apr-2007
 585 !
 586 ! key routines called: subr ternary_nuc_napari
 587 !
 588 ! references:
 589 !    merikanto, j., i. napari, h. vehkamaki, t. anttila,
 590 !     and m. kulmala, 2007, new parameterization of
 591 !     sulfuric acid-ammonia-water ternary nucleation
 592 !     rates at tropospheric conditions,
 593 !       j. geophys. res., 112, d15207, doi:10.1029/2006jd0027977
 594 !
 595 !    vehkamaki, h., m. kulmala, i. napari, k.e.j. lehtinen,
 596 !       c. timmreck, m. noppel and a. laaksonen, 2002,
 597 !       an improved parameterization for sulfuric acid-water nucleation
 598 !       rates for tropospheric and stratospheric conditions,
 599 !       j. geophys. res., 107, 4622, doi:10.1029/2002jd002184
 600 !
 601 !    kerminen, v., and m. kulmala, 2002,
 602 !       analytical formulae connecting the "real" and the "apparent"
 603 !       nucleation rate and the nuclei number concentration
 604 !       for atmospheric nucleation events
 605 !
 606 !.......................................................................
 607       implicit none
 608
 609 ! subr arguments (in)
 610         real(r8), intent(in) :: dtnuc             ! nucleation time step (s)
 611         real(r8), intent(in) :: temp_in           ! temperature, in k
 612         real(r8), intent(in) :: rh_in             ! relative humidity, as fraction
 613         real(r8), intent(in) :: press_in          ! air pressure (pa)
 614
 615         real(r8), intent(in) :: qh2so4_cur, qh2so4_avg
 616         real(r8), intent(in) :: qnh3_cur, qnh3_avg
 617              ! above 4 variables are gas h2so4 and gas nh3 mixing ratios (mol/mol-air)
 618              ! qxxx_cur = current value (after gas chem and condensation)
 619              ! qxxx_avg = estimated average value (for simultaneous source/sink calcs)
 620         real(r8), intent(in) :: h2so4_uptkrate    ! h2so4 uptake rate to aerosol (1/s)
 621
 622         integer,  intent(in) :: newnuc_method_flagaa    ! 1=merikanto et al (2007) ternary
 623                                                         ! 2=vehkamaki et al (2002) binary
 624         integer,  intent(in) :: nsize                   ! number of aerosol size bins
 625         integer,  intent(in) :: maxd_asize              ! dimension for dplom_sect, ...
 626         real(r8), intent(in) :: dplom_sect(maxd_asize)  ! dry diameter at lower bnd of bin (m)
 627         real(r8), intent(in) :: dphim_sect(maxd_asize)  ! dry diameter at upper bnd of bin (m)
 628         integer,  intent(in) :: itstep !BSINGH - 05/28/2013(RCE updates)
 629         integer,  intent(in) :: ldiagaa, lundiagaa
 630
 631 ! subr arguments (inout & out)
 632         integer,  intent(out) :: isize_nuc        ! size bin into which new particles go
 633         real(r8), intent(out) :: qnuma_del        ! change to aerosol number mixing ratio (#/mol-air)
 634         real(r8), intent(out) :: qso4a_del        ! change to aerosol so4 mixing ratio (mol/mol-air)
 635         real(r8), intent(out) :: qnh4a_del        ! change to aerosol nh4 mixing ratio (mol/mol-air)
 636         real(r8), intent(out) :: qh2so4_del       ! change to gas h2so4 mixing ratio (mol/mol-air)
 637         real(r8), intent(out) :: qnh3_del         ! change to gas nh3 mixing ratio (mol/mol-air)
 638                                                   ! aerosol changes are > 0; gas changes are < 0
 639         real(r8), intent(inout) :: dens_nh4so4a   ! dry-density of the new nh4-so4 aerosol mass (kg/m3)
 640                                                   ! use 'in' value only if it is between 1500-2000 kg/m3
 641
 642 ! subr arguments (out) passed via common block
 643 !    these are used to duplicate the outputs of yang zhang's original test driver
 644 !    they are not really needed in wrf-chem
 645         real(r8) :: ratenuclt        ! j = ternary nucleation rate from napari param. (cm-3 s-1)
 646         real(r8) :: rateloge         ! ln (j)
 647         real(r8) :: cnum_h2so4       ! number of h2so4 molecules in the critical nucleus
 648         real(r8) :: cnum_nh3         ! number of nh3   molecules in the critical nucleus
 649         real(r8) :: cnum_tot         ! total number of molecules in the critical nucleus
 650         real(r8) :: radius_cluster   ! the radius of cluster (nm)
 651
 652
 653 ! local variables
 654         integer :: i
 655         integer :: igrow
 656         integer, save :: icase = 0, icase_reldiffmax = 0
 657 !       integer, parameter :: ldiagaa = -1
 658         integer :: lun
 659         integer :: newnuc_method_flagaa2
 660
 661         real(r8), parameter :: onethird = 1.0/3.0
 662         real(r8), parameter :: avogad = 6.022e23   ! avogadro number (molecules/mol)
 663         real(r8), parameter :: mw_air = 28.966     ! dry-air mean molecular weight (g/mol)
 664
 665         real(r8), parameter :: accom_coef_h2so4 = 0.65   ! accomodation coef for h2so4 conden
 666
 667 ! dry densities (kg/m3) molecular weights of aerosol
 668 ! ammsulf, ammbisulf, and sulfacid (from mosaic  dens_electrolyte values)
 669 !       real(r8), parameter :: dens_ammsulf   = 1.769e3
 670 !       real(r8), parameter :: dens_ammbisulf = 1.78e3
 671 !       real(r8), parameter :: dens_sulfacid  = 1.841e3
 672 ! use following to match cam3 modal_aero densities
 673         real(r8), parameter :: dens_ammsulf   = 1.770e3
 674         real(r8), parameter :: dens_ammbisulf = 1.770e3
 675         real(r8), parameter :: dens_sulfacid  = 1.770e3
 676         real(r8), parameter :: dens_water     = 1.0e3
 677
 678 ! molecular weights (g/mol) of aerosol ammsulf, ammbisulf, and sulfacid
 679 !    for ammbisulf and sulfacid, use 114 & 96 here rather than 115 & 98
 680 !    because we don't keep track of aerosol hion mass
 681         real(r8), parameter :: mw_ammsulf   = 132.0
 682         real(r8), parameter :: mw_ammbisulf = 114.0
 683         real(r8), parameter :: mw_sulfacid  =  96.0
 684 ! molecular weights of aerosol sulfate and ammonium
 685         real(r8), parameter :: mw_so4a      =  96.0
 686         real(r8), parameter :: mw_nh4a      =  18.0
 687         real(r8), parameter :: mw_water     =  18.0
 688
 689         real(r8), save :: reldiffmax = 0.0
 690
 691         real(r8) cair                     ! dry-air molar density (mol/m3)
 692         real(r8) cs_prime_kk              ! kk2002 "cs_prime" parameter (1/m2)
 693         real(r8) cs_kk                    ! kk2002 "cs" parameter (1/s)
 694         real(r8) dens_part                ! "grown" single-particle dry density (kg/m3)
 695         real(r8) dfin_kk, dnuc_kk         ! kk2002 final/initial new particle wet diameter (nm)
 696         real(r8) dpdry_clus               ! critical cluster diameter (m)
 697         real(r8) dpdry_part               ! "grown" single-particle dry diameter (m)
 698         real(r8) tmpa, tmpb, tmpc, tmpe, tmpq
 699         real(r8) tmpa1, tmpb1
 700         real(r8) tmp_m1, tmp_m2, tmp_m3, tmp_n1, tmp_n2, tmp_n3
 701         real(r8) tmp_spd                  ! h2so4 vapor molecular speed (m/s)
 702         real(r8) factor_kk
 703         real(r8) fogas, foso4a, fonh4a, fonuma
 704         real(r8) freduce                  ! reduction factor applied to nucleation rate
 705                                           ! due to limited availability of h2so4 & nh3 gases
 706         real(r8) freducea, freduceb
 707         real(r8) gamma_kk                 ! kk2002 "gamma" parameter (nm2*m2/h)
 708         real(r8) gr_kk                    ! kk2002 "gr" parameter (nm/h)
 709         real(r8) kgaero_per_moleso4a      ! (kg dry aerosol)/(mol aerosol so4)
 710         real(r8) mass_part                ! "grown" single-particle dry mass (kg)
 711         real(r8) molenh4a_per_moleso4a    ! (mol aerosol nh4)/(mol aerosol so4)
 712         real(r8) nh3ppt, nh3ppt_bb        ! actual and bounded nh3 (ppt)
 713         real(r8) nu_kk                    ! kk2002 "nu" parameter (nm)
 714         real(r8) qmolnh4a_del_max         ! max production of aerosol nh4 over dtnuc (mol/mol-air)
 715         real(r8) qmolso4a_del_max         ! max production of aerosol so4 over dtnuc (mol/mol-air)
 716         real(r8) ratenuclt_bb             ! nucleation rate (#/m3/s)
 717         real(r8) ratenuclt_kk             ! nucleation rate after kk2002 adjustment (#/m3/s)
 718         real(r8) rh_bb                    ! bounded value of rh_in
 719         real(r8) so4vol_in                ! concentration of h2so4 for nucl. calc., molecules cm-3
 720         real(r8) so4vol_bb                ! bounded value of so4vol_in
 721         real(r8) temp_bb                  ! bounded value of temp_in
 722         real(r8) voldry_clus              ! critical-cluster dry volume (m3)
 723         real(r8) voldry_part              ! "grown" single-particle dry volume (m3)
 724         real(r8) wetvol_dryvol            ! grown particle (wet-volume)/(dry-volume)
 725         real(r8) wet_volfrac_so4a         ! grown particle (dry-volume-from-so4)/(wet-volume)
 726
 727
 728
 729 !
 730 ! if h2so4 vapor < qh2so4_cutoff
 731 ! exit with new particle formation = 0
 732 !
 733         isize_nuc = 1
 734         qnuma_del = 0.0
 735         qso4a_del = 0.0
 736         qnh4a_del = 0.0
 737         qh2so4_del = 0.0
 738         qnh3_del = 0.0
 739 !       if (qh2so4_avg .le. qh2so4_cutoff) return   ! this no longer needed
 740 !       if (qh2so4_cur .le. qh2so4_cutoff) return   ! this no longer needed
 741
 742
 743 !
 744 ! make call to vehkamaki parameterization routine
 745 !
 746
 747 ! calc h2so4 in molecules/cm3 and nh3 in ppt
 748         cair = press_in/(temp_in*8.3144)
 749         so4vol_in = qh2so4_avg * cair * avogad * 1.0e-6
 750         nh3ppt    = qnh3_avg * 1.0e12
 751         newnuc_method_flagaa2 = 0!BSINGH - 05/28/2013(RCE updates)
 752
 753         if ((newnuc_method_flagaa == 1) .and. (nh3ppt >= 0.1)) then
 754             if (so4vol_in < 5.0e4) return
 755 ! make call to merikanto ternary parameterization routine
 756 ! (when nh3ppt < 0.1, use binary param instead)
 757             temp_bb = max( 235.0_r8, min( 295.0_r8, temp_in ) )
 758             rh_bb = max( 0.05_r8, min( 0.95_r8, rh_in ) )
 759             so4vol_bb = max( 5.0e4_r8, min( 1.0e9_r8, so4vol_in ) )
 760             nh3ppt_bb = max( 0.1_r8, min( 1.0e3_r8, nh3ppt ) )
 761             call ternary_nuc_merik2007(   &
 762                 temp_bb, rh_bb, so4vol_bb, nh3ppt_bb,   &
 763                 rateloge,   &
 764                 cnum_tot, cnum_h2so4, cnum_nh3, radius_cluster )
 765             newnuc_method_flagaa2 = 1
 766
 767         else if ((newnuc_method_flagaa == 1) .or. &
 768                  (newnuc_method_flagaa == 2)) then
 769
 770             if (so4vol_in < 1.0e4) return
 771 ! make call to vehkamaki binary parameterization routine
 772             temp_bb = max( 230.15_r8, min( 305.15_r8, temp_in ) )
 773             rh_bb = max( 1.0e-4_r8, min( 1.0_r8, rh_in ) )
 774             so4vol_bb = max( 1.0e4_r8, min( 1.0e11_r8, so4vol_in ) )
 775             call binary_nuc_vehk2002(   &
 776                 temp_bb, rh_bb, so4vol_bb,   &
 777                 ratenuclt, rateloge,   &
 778                 cnum_h2so4, cnum_tot, radius_cluster )
 779             cnum_nh3 = 0.0
 780             newnuc_method_flagaa2 = 2
 781             !BSINGH - 05/28/2013(RCE updates)
 782          else if ((newnuc_method_flagaa == 11) .or. &
 783               (newnuc_method_flagaa == 12)) then
 784             call pbl_nuc_wang2008( so4vol_in,   &
 785                  newnuc_method_flagaa, newnuc_method_flagaa2,   &
 786                  ratenuclt, rateloge,   &
 787                  cnum_tot, cnum_h2so4, cnum_nh3, radius_cluster )
 788             !BSINGH - 05/28/2013(RCE updates ENDS)
 789         else
 790             return
 791         end if
 792         if (ldiagaa > 0) &
 793         write(lundiagaa,'(a,i8,1p,e10.2)') 'newnuc - itstep, ratenucl =', itstep, exp( rateloge )!BSINGH - 05/28/2013(RCE updates)
 794 ! if nucleation rate is less than 1e-6 #/m3/s ~= 0.1 #/cm3/day,
 795 ! exit with new particle formation = 0
 796         if (rateloge  .le. -13.82_r8) return
 797 !       if (ratenuclt .le. 1.0e-6) return
 798         ratenuclt = exp( rateloge )
 799         ratenuclt_bb = ratenuclt*1.0e6_r8
 800
 801
 802 ! wet/dry volume ratio - use simple kohler approx for ammsulf/ammbisulf
 803         tmpa = max( 0.10_r8, min( 0.95_r8, rh_in ) )
 804         wetvol_dryvol = 1.0 - 0.56/log(tmpa)
 805
 806
 807 ! determine size bin into which the new particles go
 808 ! (probably it will always be bin #1, but ...)
 809         voldry_clus = ( max(cnum_h2so4,1.0_r8)*mw_so4a + cnum_nh3*mw_nh4a ) /   &
 810                       (1.0e3*dens_sulfacid*avogad)
 811         dpdry_clus = (voldry_clus*6.0/pi)**onethird
 812
 813         isize_nuc = 1
 814         dpdry_part = dplom_sect(1)
 815         if (dpdry_clus <= dplom_sect(1)) then
 816            igrow = 1   ! need to clusters to larger size
 817         else if (dpdry_clus >= dphim_sect(nsize)) then
 818            igrow = 0
 819            isize_nuc = nsize
 820            dpdry_part = dphim_sect(nsize)
 821         else
 822            igrow = 0
 823            do i = 1, nsize
 824               if (dpdry_clus < dphim_sect(i)) then
 825                  isize_nuc = i
 826                  dpdry_part = dpdry_clus
 827                  dpdry_part = min( dpdry_part, dphim_sect(i) )
 828                  dpdry_part = max( dpdry_part, dplom_sect(i) )
 829                  exit
 830               end if
 831            end do
 832         end if
 833         voldry_part = (pi/6.0)*(dpdry_part**3)
 834         if (ldiagaa > 0) &
 835              write(lundiagaa,'(a,i8,i4)') 'newnuc - itstep, igrow    =', itstep, igrow!BSINGH - 05/28/2013(RCE updates)
 836
 837
 838 !
 839 ! determine composition and density of the "grown particles"
 840 ! the grown particles are assumed to be liquid
 841 !    (since critical clusters contain water)
 842 !    so any (nh4/so4) molar ratio between 0 and 2 is allowed
 843 ! assume that the grown particles will have
 844 !    (nh4/so4 molar ratio) = min( 2, (nh3/h2so4 gas molar ratio) )
 845 !
 846         if (igrow .le. 0) then
 847 ! no "growing" so pure sulfuric acid
 848            tmp_n1 = 0.0
 849            tmp_n2 = 0.0
 850            tmp_n3 = 1.0
 851         else if (qnh3_cur .ge. qh2so4_cur) then
 852 ! combination of ammonium sulfate and ammonium bisulfate
 853 ! tmp_n1 & tmp_n2 = mole fractions of the ammsulf & ammbisulf
 854            tmp_n1 = (qnh3_cur/qh2so4_cur) - 1.0
 855            tmp_n1 = max( 0.0_r8, min( 1.0_r8, tmp_n1 ) )
 856            tmp_n2 = 1.0 - tmp_n1
 857            tmp_n3 = 0.0
 858         else
 859 ! combination of ammonium bisulfate and sulfuric acid
 860 ! tmp_n2 & tmp_n3 = mole fractions of the ammbisulf & sulfacid
 861            tmp_n1 = 0.0
 862            tmp_n2 = (qnh3_cur/qh2so4_cur)
 863            tmp_n2 = max( 0.0_r8, min( 1.0_r8, tmp_n2 ) )
 864            tmp_n3 = 1.0 - tmp_n2
 865         end if
 866
 867         tmp_m1 = tmp_n1*mw_ammsulf
 868         tmp_m2 = tmp_n2*mw_ammbisulf
 869         tmp_m3 = tmp_n3*mw_sulfacid
 870         dens_part = (tmp_m1 + tmp_m2 + tmp_m3)/   &
 871            ((tmp_m1/dens_ammsulf) + (tmp_m2/dens_ammbisulf)   &
 872                                   + (tmp_m3/dens_sulfacid))
 873 ! use 'in' value only if it is between 1500-2000 kg/m3
 874         if (abs(dens_nh4so4a-1750.0) .le. 250.0) then
 875             dens_part = dens_nh4so4a
 876         else
 877             dens_nh4so4a = dens_part
 878         end if
 879
 880         mass_part  = voldry_part*dens_part
 881 ! (mol aerosol nh4)/(mol aerosol so4)
 882         molenh4a_per_moleso4a = 2.0*tmp_n1 + tmp_n2
 883 ! (kg dry aerosol)/(mol aerosol so4)
 884         kgaero_per_moleso4a = 1.0e-3*(tmp_m1 + tmp_m2 + tmp_m3)
 885
 886 ! fraction of wet volume due to so4a
 887         tmpb = 1.0 + molenh4a_per_moleso4a*17.0/98.0
 888         wet_volfrac_so4a = 1.0 / ( wetvol_dryvol * tmpb )
 889
 890
 891 !
 892 ! calc kerminen & kulmala (2002) correction
 893 !
 894         if (igrow <=  0) then
 895             factor_kk = 1.0
 896
 897         else
 898 ! "gr" parameter (nm/h) = condensation growth rate of new particles
 899 ! use kk2002 eqn 21 for h2so4 uptake, and correct for nh3 & h2o uptake
 900             tmp_spd = 14.7*sqrt(temp_in)   ! h2so4 molecular speed (m/s)
 901             gr_kk = 3.0e-9*tmp_spd*mw_sulfacid*so4vol_in/   &
 902                     (dens_part*wet_volfrac_so4a)
 903
 904 ! "gamma" parameter (nm2/m2/h)
 905 ! use kk2002 eqn 22
 906 !
 907 ! dfin_kk = wet diam (nm) of grown particle having dry dia = dpdry_part (m)
 908             dfin_kk = 1.0e9 * dpdry_part * (wetvol_dryvol**onethird)
 909 ! dnuc_kk = wet diam (nm) of cluster
 910             dnuc_kk = 2.0*radius_cluster
 911             dnuc_kk = max( dnuc_kk, 1.0_r8 )
 912 ! neglect (dmean/150)**0.048 factor,
 913 ! which should be very close to 1.0 because of small exponent
 914             gamma_kk = 0.23 * (dnuc_kk)**0.2   &
 915                      * (dfin_kk/3.0)**0.075   &
 916                      * (dens_part*1.0e-3)**(-0.33)   &
 917                      * (temp_in/293.0)**(-0.75)
 918
 919 ! "cs_prime parameter" (1/m2)
 920 ! instead kk2002 eqn 3, use
 921 !     cs_prime ~= tmpa / (4*pi*tmpb * h2so4_accom_coef)
 922 ! where
 923 !     tmpa = -d(ln(h2so4))/dt by conden to particles   (1/h units)
 924 !     tmpb = h2so4 vapor diffusivity (m2/h units)
 925 ! this approx is generally within a few percent of the cs_prime
 926 !     calculated directly from eqn 2,
 927 !     which is acceptable, given overall uncertainties
 928 ! tmpa = -d(ln(h2so4))/dt by conden to particles   (1/h units)
 929             tmpa = h2so4_uptkrate * 3600.0
 930             tmpa1 = tmpa
 931             tmpa = max( tmpa, 0.0_r8 )
 932 ! tmpb = h2so4 gas diffusivity (m2/s, then m2/h)
 933             tmpb = 6.7037e-6 * (temp_in**0.75) / cair
 934             tmpb1 = tmpb         ! m2/s
 935             tmpb = tmpb*3600.0   ! m2/h
 936             cs_prime_kk = tmpa/(4.0*pi*tmpb*accom_coef_h2so4)
 937             cs_kk = cs_prime_kk*4.0*pi*tmpb1
 938
 939 ! "nu" parameter (nm) -- kk2002 eqn 11
 940             nu_kk = gamma_kk*cs_prime_kk/gr_kk
 941 ! nucleation rate adjustment factor (--) -- kk2002 eqn 13
 942             factor_kk = exp( (nu_kk/dfin_kk) - (nu_kk/dnuc_kk) )
 943
 944         end if
 945         ratenuclt_kk = ratenuclt_bb*factor_kk
 946
 947
 948 ! max production of aerosol dry mass (kg-aero/m3-air)
 949         tmpa = max( 0.0_r8, (ratenuclt_kk*dtnuc*mass_part) )
 950 ! max production of aerosol so4 (mol-so4a/mol-air)
 951         tmpe = tmpa/(kgaero_per_moleso4a*cair)
 952 ! max production of aerosol so4 (mol/mol-air)
 953 ! based on ratenuclt_kk and mass_part
 954         qmolso4a_del_max = tmpe
 955
 956 ! check if max production exceeds available h2so4 vapor
 957         freducea = 1.0
 958         if (qmolso4a_del_max .gt. qh2so4_cur) then
 959            freducea = qh2so4_cur/qmolso4a_del_max
 960         end if
 961
 962 ! check if max production exceeds available nh3 vapor
 963         freduceb = 1.0
 964         if (molenh4a_per_moleso4a .ge. 1.0e-10) then
 965 ! max production of aerosol nh4 (ppm) based on ratenuclt_kk and mass_part
 966            qmolnh4a_del_max = qmolso4a_del_max*molenh4a_per_moleso4a
 967            if (qmolnh4a_del_max .gt. qnh3_cur) then
 968               freduceb = qnh3_cur/qmolnh4a_del_max
 969            end if
 970         end if
 971         freduce = min( freducea, freduceb )
 972
 973 ! if adjusted nucleation rate is less than 1e-12 #/m3/s ~= 0.1 #/cm3/day,
 974 ! exit with new particle formation = 0
 975         if (freduce*ratenuclt_kk .le. 1.0e-12) return
 976
 977
 978 ! note:  suppose that at this point, freduce < 1.0 (no gas-available
 979 !    constraints) and molenh4a_per_moleso4a < 2.0
 980 ! if the gas-available constraints is do to h2so4 availability,
 981 !    then it would be possible to condense "additional" nh3 and have
 982 !    (nh3/h2so4 gas molar ratio) < (nh4/so4 aerosol molar ratio) <= 2
 983 ! one could do some additional calculations of
 984 !    dens_part & molenh4a_per_moleso4a to realize this
 985 ! however, the particle "growing" is a crude approximate way to get
 986 !    the new particles to the host code's minimum particle size,
 987 ! are such refinements worth the effort?
 988
 989
 990 ! changes to h2so4 & nh3 gas (in mol/mol-air), limited by amounts available
 991         tmpa = 0.9999
 992         qh2so4_del = min( tmpa*qh2so4_cur, freduce*qmolso4a_del_max )
 993         qnh3_del   = min( tmpa*qnh3_cur, qh2so4_del*molenh4a_per_moleso4a )
 994         qh2so4_del = -qh2so4_del
 995         qnh3_del   = -qnh3_del
 996
 997 ! changes to so4 & nh4 aerosol (in mol/mol-air)
 998         qso4a_del = -qh2so4_del
 999         qnh4a_del =   -qnh3_del
1000 ! change to aerosol number (in #/mol-air)
1001         qnuma_del = 1.0e-3*(qso4a_del*mw_so4a + qnh4a_del*mw_nh4a)/mass_part
1002
1003 ! do the following (tmpa, tmpb, tmpc) calculations as a check
1004 ! max production of aerosol number (#/mol-air)
1005         tmpa = max( 0.0_r8, (ratenuclt_kk*dtnuc/cair) )
1006 ! adjusted production of aerosol number (#/mol-air)
1007         tmpb = tmpa*freduce
1008 ! relative difference from qnuma_del
1009         tmpc = (tmpb - qnuma_del)/max(tmpb, qnuma_del, 1.0e-35_r8)
1010
1011
1012 !
1013 ! diagnostic output to fort.41
1014 ! (this should be commented-out or deleted in the wrf-chem version)
1015 !
1016         if (ldiagaa <= 99) return!BSINGH - 05/28/2013(RCE updates)
1017         lun = lundiagaa
1018
1019         icase = icase + 1
1020         if (abs(tmpc) .gt. abs(reldiffmax)) then
1021            reldiffmax = tmpc
1022            icase_reldiffmax = icase
1023         end if
1024 !       do lun = 41, 51, 10
1025 !       do lun = 6, 6
1026 !          write(lun,'(/)')
1027            write(lun,'(a,2i9,1p,e10.2)')   &
1028                'vehkam bin-nuc icase, icase_rdmax =',   &
1029                icase, icase_reldiffmax, reldiffmax
1030            if (freduceb .lt. freducea) then
1031               if (abs(freducea-freduceb) .gt.   &
1032                    3.0e-7*max(freduceb,freducea)) write(lun,'(a,1p,2e15.7)')   &
1033                  'freducea, b =', freducea, freduceb
1034            end if
1035 !       end do
1036
1037 ! output factors so that output matches that of ternucl03
1038 !       fogas  = 1.0e6                     ! convert mol/mol-air to ppm
1039 !       foso4a = 1.0e9*mw_so4a/mw_air      ! convert mol-so4a/mol-air to ug/kg-air
1040 !       fonh4a = 1.0e9*mw_nh4a/mw_air      ! convert mol-nh4a/mol-air to ug/kg-air
1041 !       fonuma = 1.0e3/mw_air              ! convert #/mol-air to #/kg-air
1042         fogas  = 1.0
1043         foso4a = 1.0
1044         fonh4a = 1.0
1045         fonuma = 1.0
1046
1047 !       do lun = 41, 51, 10
1048 !       do lun = 6, 6
1049
1050         write(lun,'(a,2i5)') 'newnuc_method_flagaa/aa2',   &
1051            newnuc_method_flagaa, newnuc_method_flagaa2
1052
1053         write(lun,9210)
1054         write(lun,9201) temp_in, rh_in,   &
1055            ratenuclt, 2.0*radius_cluster*1.0e-7, dpdry_part*1.0e2,   &
1056            voldry_part*1.0e6, float(igrow)
1057         write(lun,9215)
1058         write(lun,9201)   &
1059            qh2so4_avg*fogas, qnh3_avg*fogas,  &
1060            qh2so4_cur*fogas, qnh3_cur*fogas,  &
1061            qh2so4_del*fogas, qnh3_del*fogas,  &
1062            qso4a_del*foso4a, qnh4a_del*fonh4a
1063
1064         write(lun,9220)
1065         write(lun,9201)   &
1066            dtnuc, dens_nh4so4a*1.0e-3,   &
1067            (qnh3_cur/qh2so4_cur), molenh4a_per_moleso4a,   &
1068            qnuma_del*fonuma, tmpb*fonuma, tmpc, freduce
1069
1070 !       end do
1071
1072 !       lun = 51
1073 !       lun = 6
1074         write(lun,9230)
1075         write(lun,9201)   &
1076            press_in, cair*1.0e-6, so4vol_in,   &
1077            wet_volfrac_so4a, wetvol_dryvol, dens_part*1.0e-3
1078
1079         if (igrow > 0) then
1080         write(lun,9240)
1081         write(lun,9201)   &
1082            tmp_spd, gr_kk, dnuc_kk, dfin_kk,   &
1083            gamma_kk, tmpa1, tmpb1, cs_kk
1084
1085         write(lun,9250)
1086         write(lun,9201)   &
1087            cs_prime_kk, nu_kk, factor_kk, ratenuclt,   &
1088            ratenuclt_kk*1.0e-6
1089         end if
1090
1091 9201    format ( 1p, 40e10.2  )
1092 9210    format (   &
1093         '      temp        rh',   &
1094         '   ratenuc  dia_clus ddry_part',   &
1095         ' vdry_part     igrow' )
1096 9215    format (   &
1097         '  h2so4avg  h2so4pre',   &
1098         '  h2so4cur   nh3_cur',   &
1099         '  h2so4del   nh3_del',   &
1100         '  so4a_del  nh4a_del' )
1101 9220    format (    &
1102         '     dtnuc    dens_a   nh/so g   nh/so a',   &
1103         '  numa_del  numa_dl2   reldiff   freduce' )
1104 9230    format (   &
1105         '  press_in      cair so4_volin',   &
1106         ' wet_volfr wetv_dryv dens_part' )
1107 9240    format (   &
1108         '   tmp_spd     gr_kk   dnuc_kk   dfin_kk',   &
1109         '  gamma_kk     tmpa1     tmpb1     cs_kk' )
1110 9250    format (   &
1111         ' cs_pri_kk     nu_kk factor_kk ratenuclt',   &
1112         ' ratenu_kk' )
1113
1114
1115         return
1116         end subroutine mer07_veh02_nuc_mosaic_1box
1117
1118
1119 !BSINGH - 05/28/2013(RCE updates)
1120
1121 !-----------------------------------------------------------------------
1122 !-----------------------------------------------------------------------
1123         subroutine pbl_nuc_wang2008( so4vol,   &
1124             newnuc_method_flagaa, newnuc_method_flagaa2,   &
1125             ratenucl, rateloge,   &
1126             cnum_tot, cnum_h2so4, cnum_nh3, radius_cluster )
1127 !
1128 ! calculates boundary nucleation nucleation rate
1129 ! using the first or second-order parameterization in
1130 !     wang, m., and j.e. penner, 2008,
1131 !        aerosol indirect forcing in a global model with particle nucleation,
1132 !        atmos. chem. phys. discuss., 8, 13943-13998
1133 !
1134         implicit none
1135
1136 ! subr arguments (in)
1137         real(r8), intent(in) :: so4vol            ! concentration of h2so4 (molecules cm-3)
1138         integer, intent(in)  :: newnuc_method_flagaa
1139                                 ! [11,12] value selects [first,second]-order parameterization
1140
1141 ! subr arguments (inout)
1142         integer, intent(inout)  :: newnuc_method_flagaa2
1143         real(r8), intent(inout) :: ratenucl         ! binary nucleation rate, j (# cm-3 s-1)
1144         real(r8), intent(inout) :: rateloge         ! log( ratenucl )
1145
1146         real(r8), intent(inout) :: cnum_tot         ! total number of molecules
1147                                                     ! in the critical nucleus
1148         real(r8), intent(inout) :: cnum_h2so4       ! number of h2so4 molecules
1149         real(r8), intent(inout) :: cnum_nh3         ! number of nh3 molecules
1150         real(r8), intent(inout) :: radius_cluster   ! the radius of cluster (nm)
1151
1152
1153 ! local variables
1154         real(r8) :: tmp_diam, tmp_mass, tmp_volu
1155         real(r8) :: tmp_rateloge, tmp_ratenucl
1156
1157 ! executable
1158
1159
1160 ! nucleation rate
1161         if (newnuc_method_flagaa == 11) then
1162            tmp_ratenucl = 1.0e-6_r8 * so4vol
1163         else if (newnuc_method_flagaa == 12) then
1164            tmp_ratenucl = 1.0e-12_r8 * (so4vol**2)
1165         else
1166            return
1167         end if
1168         tmp_ratenucl = tmp_ratenucl * adjust_factor_pbl_ratenucl
1169         tmp_rateloge = log( tmp_ratenucl )
1170
1171 ! exit if pbl nuc rate is lower than (incoming) ternary/binary rate
1172 !       if (tmp_rateloge <= rateloge) return
1173
1174         rateloge = tmp_rateloge
1175         ratenucl = tmp_ratenucl
1176         newnuc_method_flagaa2 = newnuc_method_flagaa
1177
1178 ! following wang 2002, assume fresh nuclei are 1 nm diameter
1179 !    subsequent code will "grow" them to aitken mode size
1180         radius_cluster = 0.5_r8
1181
1182 ! assume fresh nuclei are pure h2so4
1183 !    since aitken size >> initial size, the initial composition
1184 !    has very little impact on the results
1185         tmp_diam = radius_cluster * 2.0e-7_r8   ! diameter in cm
1186         tmp_volu = (tmp_diam**3) * (pi/6.0_r8)  ! volume in cm^3
1187         tmp_mass = tmp_volu * 1.8_r8            ! mass in g
1188         cnum_h2so4 = (tmp_mass / 98.0_r8) * 6.023e23_r8   ! no. of h2so4 molec assuming pure h2so4
1189         cnum_tot = cnum_h2so4
1190         cnum_nh3 = 0.0_r8
1191
1192
1193         return
1194         end subroutine pbl_nuc_wang2008
1195
1196 !BSINGH - 05/28/2013(RCE updates ENDS)
1197
1198
1199 !-----------------------------------------------------------------------
1200 !-----------------------------------------------------------------------
1201         subroutine binary_nuc_vehk2002( temp, rh, so4vol,   &
1202             ratenucl, rateloge,   &
1203             cnum_h2so4, cnum_tot, radius_cluster )
1204 !
1205 ! calculates binary nucleation rate and critical cluster size
1206 ! using the parameterization in
1207 !     vehkamaki, h., m. kulmala, i. napari, k.e.j. lehtinen,
1208 !        c. timmreck, m. noppel and a. laaksonen, 2002,
1209 !        an improved parameterization for sulfuric acid-water nucleation
1210 !        rates for tropospheric and stratospheric conditions,
1211 !        j. geophys. res., 107, 4622, doi:10.1029/2002jd002184
1212 !
1213         implicit none
1214
1215 ! subr arguments (in)
1216         real(r8), intent(in) :: temp              ! temperature (k)
1217         real(r8), intent(in) :: rh                ! relative humidity (0-1)
1218         real(r8), intent(in) :: so4vol            ! concentration of h2so4 (molecules cm-3)
1219
1220 ! subr arguments (out)
1221         real(r8), intent(out) :: ratenucl         ! binary nucleation rate, j (# cm-3 s-1)
1222         real(r8), intent(out) :: rateloge         ! log( ratenucl )
1223
1224         real(r8), intent(out) :: cnum_h2so4       ! number of h2so4 molecules
1225                                                   ! in the critical nucleus
1226         real(r8), intent(out) :: cnum_tot         ! total number of molecules
1227                                                   ! in the critical nucleus
1228         real(r8), intent(out) :: radius_cluster   ! the radius of cluster (nm)
1229
1230
1231 ! local variables
1232         real(r8) :: crit_x
1233         real(r8) :: acoe, bcoe, ccoe, dcoe, ecoe, fcoe, gcoe, hcoe, icoe, jcoe
1234         real(r8) :: tmpa, tmpb
1235
1236 ! executable
1237
1238
1239 ! calc sulfuric acid mole fraction in critical cluster
1240         crit_x = 0.740997 - 0.00266379 * temp   &
1241                - 0.00349998 * log (so4vol)   &
1242                + 0.0000504022 * temp * log (so4vol)   &
1243                + 0.00201048 * log (rh)   &
1244                - 0.000183289 * temp * log (rh)   &
1245                + 0.00157407 * (log (rh)) ** 2.0   &
1246                - 0.0000179059 * temp * (log (rh)) ** 2.0   &
1247                + 0.000184403 * (log (rh)) ** 3.0   &
1248                - 1.50345e-6 * temp * (log (rh)) ** 3.0
1249
1250
1251 ! calc nucleation rate
1252         acoe    = 0.14309+2.21956*temp   &
1253                 - 0.0273911 * temp**2.0   &
1254                 + 0.0000722811 * temp**3.0 + 5.91822/crit_x
1255
1256         bcoe    = 0.117489 + 0.462532 *temp   &
1257                 - 0.0118059 * temp**2.0   &
1258                 + 0.0000404196 * temp**3.0 + 15.7963/crit_x
1259
1260         ccoe    = -0.215554-0.0810269 * temp   &
1261                 + 0.00143581 * temp**2.0   &
1262                 - 4.7758e-6 * temp**3.0   &
1263                 - 2.91297/crit_x
1264
1265         dcoe    = -3.58856+0.049508 * temp   &
1266                 - 0.00021382 * temp**2.0   &
1267                 + 3.10801e-7 * temp**3.0   &
1268                 - 0.0293333/crit_x
1269
1270         ecoe    = 1.14598 - 0.600796 * temp   &
1271                 + 0.00864245 * temp**2.0   &
1272                 - 0.0000228947 * temp**3.0   &
1273                 - 8.44985/crit_x
1274
1275         fcoe    = 2.15855 + 0.0808121 * temp   &
1276                 -0.000407382 * temp**2.0   &
1277                 -4.01957e-7 * temp**3.0   &
1278                 + 0.721326/crit_x
1279
1280         gcoe    = 1.6241 - 0.0160106 * temp   &
1281                 + 0.0000377124 * temp**2.0   &
1282                 + 3.21794e-8 * temp**3.0   &
1283                 - 0.0113255/crit_x
1284
1285         hcoe    = 9.71682 - 0.115048 * temp   &
1286                 + 0.000157098 * temp**2.0   &
1287                 + 4.00914e-7 * temp**3.0   &
1288                 + 0.71186/crit_x
1289
1290         icoe    = -1.05611 + 0.00903378 * temp   &
1291                 - 0.0000198417 * temp**2.0   &
1292                 + 2.46048e-8  * temp**3.0   &
1293                 - 0.0579087/crit_x
1294
1295         jcoe    = -0.148712 + 0.00283508 * temp   &
1296                 - 9.24619e-6  * temp**2.0   &
1297                 + 5.00427e-9 * temp**3.0   &
1298                 - 0.0127081/crit_x
1299
1300         tmpa     =     (   &
1301                   acoe   &
1302                 + bcoe * log (rh)   &
1303                 + ccoe * ( log (rh))**2.0   &
1304                 + dcoe * ( log (rh))**3.0   &
1305                 + ecoe * log (so4vol)   &
1306                 + fcoe * (log (rh)) * (log (so4vol))   &
1307                 + gcoe * ((log (rh) ) **2.0)   &
1308                        * (log (so4vol))   &
1309                 + hcoe * (log (so4vol)) **2.0   &
1310                 + icoe * log (rh)   &
1311                        * ((log (so4vol)) **2.0)   &
1312                 + jcoe * (log (so4vol)) **3.0   &
1313                 )
1314         rateloge = tmpa
1315         tmpa = min( tmpa, log(1.0e38_r8) )
1316         ratenucl = exp ( tmpa )
1317 !       write(*,*) 'tmpa, ratenucl =', tmpa, ratenucl
1318
1319
1320
1321 ! calc number of molecules in critical cluster
1322         acoe    = -0.00295413 - 0.0976834*temp   &
1323                 + 0.00102485 * temp**2.0   &
1324                 - 2.18646e-6 * temp**3.0 - 0.101717/crit_x
1325
1326         bcoe    = -0.00205064 - 0.00758504*temp   &
1327                 + 0.000192654 * temp**2.0   &
1328                 - 6.7043e-7 * temp**3.0 - 0.255774/crit_x
1329
1330         ccoe    = +0.00322308 + 0.000852637 * temp   &
1331                 - 0.0000154757 * temp**2.0   &
1332                 + 5.66661e-8 * temp**3.0   &
1333                 + 0.0338444/crit_x
1334
1335         dcoe    = +0.0474323 - 0.000625104 * temp   &
1336                 + 2.65066e-6 * temp**2.0   &
1337                 - 3.67471e-9 * temp**3.0   &
1338                 - 0.000267251/crit_x
1339
1340         ecoe    = -0.0125211 + 0.00580655 * temp   &
1341                 - 0.000101674 * temp**2.0   &
1342                 + 2.88195e-7 * temp**3.0   &
1343                 + 0.0942243/crit_x
1344
1345         fcoe    = -0.038546 - 0.000672316 * temp   &
1346                 + 2.60288e-6 * temp**2.0   &
1347                 + 1.19416e-8 * temp**3.0   &
1348                 - 0.00851515/crit_x
1349
1350         gcoe    = -0.0183749 + 0.000172072 * temp   &
1351                 - 3.71766e-7 * temp**2.0   &
1352                 - 5.14875e-10 * temp**3.0   &
1353                 + 0.00026866/crit_x
1354
1355         hcoe    = -0.0619974 + 0.000906958 * temp   &
1356                 - 9.11728e-7 * temp**2.0   &
1357                 - 5.36796e-9 * temp**3.0   &
1358                 - 0.00774234/crit_x
1359
1360         icoe    = +0.0121827 - 0.00010665 * temp   &
1361                 + 2.5346e-7 * temp**2.0   &
1362                 - 3.63519e-10 * temp**3.0   &
1363                 + 0.000610065/crit_x
1364
1365         jcoe    = +0.000320184 - 0.0000174762 * temp   &
1366                 + 6.06504e-8 * temp**2.0   &
1367                 - 1.4177e-11 * temp**3.0   &
1368                 + 0.000135751/crit_x
1369
1370         cnum_tot = exp (   &
1371                   acoe   &
1372                 + bcoe * log (rh)   &
1373                 + ccoe * ( log (rh))**2.0   &
1374                 + dcoe * ( log (rh))**3.0   &
1375                 + ecoe * log (so4vol)   &
1376                 + fcoe * (log (rh)) * (log (so4vol))   &
1377                 + gcoe * ((log (rh) ) **2.0)   &
1378                        * (log (so4vol))   &
1379                 + hcoe * (log (so4vol)) **2.0   &
1380                 + icoe * log (rh)   &
1381                        * ((log (so4vol)) **2.0)   &
1382                 + jcoe * (log (so4vol)) **3.0   &
1383                 )
1384
1385         cnum_h2so4 = cnum_tot * crit_x
1386
1387 !   calc radius (nm) of critical cluster
1388         radius_cluster = exp( -1.6524245 + 0.42316402*crit_x   &
1389                               + 0.3346648*log(cnum_tot) )
1390
1391
1392       return
1393       end subroutine binary_nuc_vehk2002
1394
1395
1396
1397 !----------------------------------------------------------------------
1398 !----------------------------------------------------------------------
1399 subroutine ternary_nuc_merik2007( t, rh, c2, c3, j_log, ntot, nacid, namm, r )
1400 !subroutine ternary_fit(          t, rh, c2, c3, j_log, ntot, nacid, namm, r )
1401 ! *************************** ternary_fit.f90 ********************************
1402 ! joonas merikanto, 2006
1403 !
1404 ! fortran 90 routine that calculates the parameterized composition
1405 ! and nucleation rate of critical clusters in h2o-h2so4-nh3 vapor
1406 !
1407 ! warning: the fit should not be used outside its limits of validity
1408 ! (limits indicated below)
1409 !
1410 ! in:
1411 ! t:     temperature (k), limits 235-295 k
1412 ! rh:    relative humidity as fraction (eg. 0.5=50%) limits 0.05-0.95
1413 ! c2:    sulfuric acid concentration (molecules/cm3) limits 5x10^4 - 10^9 molecules/cm3
1414 ! c3:    ammonia mixing ratio (ppt) limits 0.1 - 1000 ppt
1415 !
1416 ! out:
1417 ! j_log: logarithm of nucleation rate (1/(s cm3))
1418 ! ntot:  total number of molecules in the critical cluster
1419 ! nacid: number of sulfuric acid molecules in the critical cluster
1420 ! namm:  number of ammonia molecules in the critical cluster
1421 ! r:     radius of the critical cluster (nm)
1422 !  ****************************************************************************
1423 implicit none
1424
1425 real(r8), intent(in) :: t, rh, c2, c3
1426 real(r8), intent(out) :: j_log, ntot, nacid, namm, r
1427 real(r8) :: j, t_onset
1428
1429 t_onset=143.6002929064716 + 1.0178856665693992*rh + &
1430    10.196398812974294*log(c2) - &
1431    0.1849879416839113*log(c2)**2 - 17.161783213150173*log(c3) + &
1432    (109.92469248546053*log(c3))/log(c2) + &
1433    0.7734119613144357*log(c2)*log(c3) - 0.15576469879527022*log(c3)**2
1434
1435 if(t_onset.gt.t) then
1436
1437    j_log=-12.861848898625231 + 4.905527742256349*c3 - 358.2337705052991*rh -&
1438    0.05463019231872484*c3*t + 4.8630382337426985*rh*t + &
1439    0.00020258394697064567*c3*t**2 - 0.02175548069741675*rh*t**2 - &
1440    2.502406532869512e-7*c3*t**3 + 0.00003212869941055865*rh*t**3 - &
1441    4.39129415725234e6/log(c2)**2 + (56383.93843154586*t)/log(c2)**2 -&
1442    (239.835990963361*t**2)/log(c2)**2 + &
1443    (0.33765136625580167*t**3)/log(c2)**2 - &
1444    (629.7882041830943*rh)/(c3**3*log(c2)) + &
1445    (7.772806552631709*rh*t)/(c3**3*log(c2)) - &
1446    (0.031974053936299256*rh*t**2)/(c3**3*log(c2)) + &
1447    (0.00004383764128775082*rh*t**3)/(c3**3*log(c2)) + &
1448    1200.472096232311*log(c2) - 17.37107890065621*t*log(c2) + &
1449    0.08170681335921742*t**2*log(c2) - 0.00012534476159729881*t**3*log(c2) - &
1450    14.833042158178936*log(c2)**2 + 0.2932631303555295*t*log(c2)**2 - &
1451    0.0016497524241142845*t**2*log(c2)**2 + &
1452    2.844074805239367e-6*t**3*log(c2)**2 - 231375.56676032578*log(c3) - &
1453    100.21645273730675*rh*log(c3) + 2919.2852552424706*t*log(c3) + &
1454    0.977886555834732*rh*t*log(c3) - 12.286497122264588*t**2*log(c3) - &
1455    0.0030511783284506377*rh*t**2*log(c3) + &
1456    0.017249301826661612*t**3*log(c3) + 2.967320346100855e-6*rh*t**3*log(c3) + &
1457    (2.360931724951942e6*log(c3))/log(c2) - &
1458    (29752.130254319443*t*log(c3))/log(c2) + &
1459    (125.04965118142027*t**2*log(c3))/log(c2) - &
1460    (0.1752996881934318*t**3*log(c3))/log(c2) + &
1461    5599.912337254629*log(c2)*log(c3) - 70.70896612937771*t*log(c2)*log(c3) + &
1462    0.2978801613269466*t**2*log(c2)*log(c3) - &
1463    0.00041866525019504*t**3*log(c2)*log(c3) + 75061.15281456841*log(c3)**2 - &
1464    931.8802278173565*t*log(c3)**2 + 3.863266220840964*t**2*log(c3)**2 - &
1465    0.005349472062284983*t**3*log(c3)**2 - &
1466    (732006.8180571689*log(c3)**2)/log(c2) + &
1467    (9100.06398573816*t*log(c3)**2)/log(c2) - &
1468    (37.771091915932004*t**2*log(c3)**2)/log(c2) + &
1469    (0.05235455395566905*t**3*log(c3)**2)/log(c2) - &
1470    1911.0303773001353*log(c2)*log(c3)**2 + &
1471    23.6903969622286*t*log(c2)*log(c3)**2 - &
1472    0.09807872005428583*t**2*log(c2)*log(c3)**2 + &
1473    0.00013564560238552576*t**3*log(c2)*log(c3)**2 - &
1474    3180.5610833308*log(c3)**3 + 39.08268568672095*t*log(c3)**3 - &
1475    0.16048521066690752*t**2*log(c3)**3 + &
1476    0.00022031380023793877*t**3*log(c3)**3 + &
1477    (40751.075322248245*log(c3)**3)/log(c2) - &
1478    (501.66977622013934*t*log(c3)**3)/log(c2) + &
1479    (2.063469732254135*t**2*log(c3)**3)/log(c2) - &
1480    (0.002836873785758324*t**3*log(c3)**3)/log(c2) + &
1481    2.792313345723013*log(c2)**2*log(c3)**3 - &
1482    0.03422552111802899*t*log(c2)**2*log(c3)**3 + &
1483    0.00014019195277521142*t**2*log(c2)**2*log(c3)**3 - &
1484    1.9201227328396297e-7*t**3*log(c2)**2*log(c3)**3 - &
1485    980.923146020468*log(rh) + 10.054155220444462*t*log(rh) - &
1486    0.03306644502023841*t**2*log(rh) + 0.000034274041225891804*t**3*log(rh) + &
1487    (16597.75554295064*log(rh))/log(c2) - &
1488    (175.2365504237746*t*log(rh))/log(c2) + &
1489    (0.6033215603167458*t**2*log(rh))/log(c2) - &
1490    (0.0006731787599587544*t**3*log(rh))/log(c2) - &
1491    89.38961120336789*log(c3)*log(rh) + 1.153344219304926*t*log(c3)*log(rh) - &
1492    0.004954549700267233*t**2*log(c3)*log(rh) + &
1493    7.096309866238719e-6*t**3*log(c3)*log(rh) + &
1494    3.1712136610383244*log(c3)**3*log(rh) - &
1495    0.037822330602328806*t*log(c3)**3*log(rh) + &
1496    0.0001500555743561457*t**2*log(c3)**3*log(rh) - &
1497    1.9828365865570703e-7*t**3*log(c3)**3*log(rh)
1498
1499    j=exp(j_log)
1500
1501    ntot=57.40091052369212 - 0.2996341884645408*t + &
1502    0.0007395477768531926*t**2 - &
1503    5.090604835032423*log(c2) + 0.011016634044531128*t*log(c2) + &
1504    0.06750032251225707*log(c2)**2 - 0.8102831333223962*log(c3) + &
1505    0.015905081275952426*t*log(c3) - 0.2044174683159531*log(c2)*log(c3) + &
1506    0.08918159167625832*log(c3)**2 - 0.0004969033586666147*t*log(c3)**2 + &
1507    0.005704394549007816*log(c3)**3 + 3.4098703903474368*log(j) - &
1508    0.014916956508210809*t*log(j) + 0.08459090011666293*log(c3)*log(j) - &
1509    0.00014800625143907616*t*log(c3)*log(j) + 0.00503804694656905*log(j)**2
1510
1511    r=3.2888553966535506e-10 - 3.374171768439839e-12*t + &
1512    1.8347359507774313e-14*t**2 + 2.5419844298881856e-12*log(c2) - &
1513    9.498107643050827e-14*t*log(c2) + 7.446266520834559e-13*log(c2)**2 + &
1514    2.4303397746137294e-11*log(c3) + 1.589324325956633e-14*t*log(c3) - &
1515    2.034596219775266e-12*log(c2)*log(c3) - 5.59303954457172e-13*log(c3)**2 - &
1516    4.889507104645867e-16*t*log(c3)**2 + 1.3847024107506764e-13*log(c3)**3 + &
1517    4.141077193427042e-15*log(j) - 2.6813110884009767e-14*t*log(j) + &
1518    1.2879071621313094e-12*log(c3)*log(j) - &
1519    3.80352446061867e-15*t*log(c3)*log(j) - 1.8790172502456827e-14*log(j)**2
1520
1521    nacid=-4.7154180661803595 + 0.13436423483953885*t - &
1522    0.00047184686478816176*t**2 - &
1523    2.564010713640308*log(c2) + 0.011353312899114723*t*log(c2) + &
1524    0.0010801941974317014*log(c2)**2 + 0.5171368624197119*log(c3) - &
1525    0.0027882479896204665*t*log(c3) + 0.8066971907026886*log(c3)**2 - &
1526    0.0031849094214409335*t*log(c3)**2 - 0.09951184152927882*log(c3)**3 + &
1527    0.00040072788891745513*t*log(c3)**3 + 1.3276469271073974*log(j) - &
1528    0.006167654171986281*t*log(j) - 0.11061390967822708*log(c3)*log(j) + &
1529    0.0004367575329273496*t*log(c3)*log(j) + 0.000916366357266258*log(j)**2
1530
1531    namm=71.20073903979772 - 0.8409600103431923*t + &
1532    0.0024803006590334922*t**2 + &
1533    2.7798606841602607*log(c2) - 0.01475023348171676*t*log(c2) + &
1534    0.012264508212031405*log(c2)**2 - 2.009926050440182*log(c3) + &
1535    0.008689123511431527*t*log(c3) - 0.009141180198955415*log(c2)*log(c3) + &
1536    0.1374122553905617*log(c3)**2 - 0.0006253227821679215*t*log(c3)**2 + &
1537    0.00009377332742098946*log(c3)**3 + 0.5202974341687757*log(j) - &
1538    0.002419872323052805*t*log(j) + 0.07916392322884074*log(c3)*log(j) - &
1539    0.0003021586030317366*t*log(c3)*log(j) + 0.0046977006608603395*log(j)**2
1540
1541 else
1542 ! nucleation rate less that 5e-6, setting j_log arbitrary small
1543    j_log=-300.
1544 end if
1545
1546 return
1547
1548 end  subroutine ternary_nuc_merik2007
1549
1550
1551
1552 !-----------------------------------------------------------------------
1553 !-----------------------------------------------------------------------
1554         subroutine wexler_nuc_mosaic_1box(   &
1555            dtnuc, temp_in, rh_in, cair,   &
1556            qh2so4_avg, qh2so4_cur, qnh3_avg, qnh3_cur,   &
1557            nsize, maxd_asize, volumlo_sect, volumhi_sect,   &
1558            isize_nuc, qnuma_del, qso4a_del, qnh4a_del,   &
1559            qh2so4_del, qnh3_del, dens_nh4so4a )
1560 !.......................................................................
1561 !
1562 ! calculates new particle production from h2so4-h2o binary nucleation
1563 !    over timestep dtnuc, using the wexler et al. (1994) parameterization
1564 !
1565 ! the size of new particles is the lower-bound size of the host code's
1566 !    smallest size bin.  their composition is so4 and nh4, since the nuclei
1567 !    would incorporate nh3 as they grow from ~1 nm to the lower-bound size.
1568 !    (the new particle composition can be forced to pure so4 by setting
1569 !    the qnh3_avg & qnh3_cur input arguments to 0.0).
1570 !
1571 ! revision history
1572 !    coded by rc easter, pnnl, 20-mar-2006
1573 !
1574 ! key routines called: none
1575 !
1576 ! references:
1577 !    wexler, a. s., f. w. lurmann, and j. h. seinfeld,
1578 !       modelling urban and regional aerosols -- i.  model development,
1579 !       atmos. environ., 28, 531-546, 1994.
1580 !
1581 !.......................................................................
1582         implicit none
1583
1584 ! subr arguments (in)
1585         real(r8), intent(in) :: dtnuc             ! nucleation time step (s)
1586         real(r8), intent(in) :: temp_in           ! temperature, in k
1587         real(r8), intent(in) :: rh_in             ! relative humidity, as fraction
1588         real(r8), intent(in) :: cair              ! dry-air molar density (mole-air/cm3)
1589
1590         real(r8), intent(in) :: qh2so4_avg, qh2so4_cur   ! gas h2so4 mixing ratios (ppm)
1591         real(r8), intent(in) :: qnh3_avg, qnh3_cur       ! gas nh3 mixing ratios (ppm)
1592              ! qxxx_cur = current value (at end of condensation)
1593              ! qxxx_avg = average value (from start to end of condensation)
1594
1595         integer, intent(in) :: nsize                    ! number of aerosol size bins
1596         integer, intent(in) :: maxd_asize               ! dimension for volumlo_sect, ...
1597         real(r8), intent(in) :: volumlo_sect(maxd_asize)    ! dry volume at lower bnd of bin (cm3)
1598         real(r8), intent(in) :: volumhi_sect(maxd_asize)    ! dry volume at upper bnd of bin (cm3)
1599
1600 ! subr arguments (out)
1601         integer, intent(out) :: isize_nuc     ! size bin into which new particles go
1602         real(r8), intent(out) :: qnuma_del        ! change to aerosol number mixing ratio (#/kg)
1603         real(r8), intent(out) :: qso4a_del        ! change to aerosol so4 mixing ratio (ug/kg)
1604         real(r8), intent(out) :: qnh4a_del        ! change to aerosol nh4 mixing ratio (ug/kg)
1605         real(r8), intent(out) :: qh2so4_del       ! change to gas h2so4 mixing ratio (ppm)
1606         real(r8), intent(out) :: qnh3_del         ! change to gas nh3 mixing ratio (ppm)
1607                                               ! aerosol changes are > 0; gas changes are < 0
1608
1609 ! subr arguments (inout)
1610         real(r8), intent(inout) :: dens_nh4so4a   ! dry-density of the new nh4-so4 aerosol mass (g/cm3)
1611                                               ! use 'in' value only if it is between 1.5-2.0 g/cm3
1612
1613 ! local variables
1614         integer i
1615         integer, save :: icase = 0, icase_reldiffmax = 0
1616
1617         real(r8), parameter :: pi = 3.1415926536
1618         real(r8), parameter :: avogad = 6.022e23   ! avogadro number (molecules/mole)
1619         real(r8), parameter :: mw_air = 28.966     ! dry-air mean molecular weight (g/mole)
1620
1621 ! dry densities (g/cm3) molecular weights of aerosol
1622 ! ammsulf, ammbisulf, and sulfacid (from mosaic  dens_electrolyte values)
1623         real(r8), parameter :: dens_ammsulf   = 1.769
1624         real(r8), parameter :: dens_ammbisulf = 1.78
1625         real(r8), parameter :: dens_sulfacid  = 1.841
1626
1627 ! molecular weights (g/mole) of aerosol ammsulf, ammbisulf, and sulfacid
1628 !    for ammbisulf and sulfacid, use 114 & 96 here rather than 115 & 98
1629 !    because we don't keep track of aerosol hion mass
1630         real(r8), parameter :: mw_ammsulf   = 132.0
1631         real(r8), parameter :: mw_ammbisulf = 114.0
1632         real(r8), parameter :: mw_sulfacid  =  96.0
1633 ! molecular weights of aerosol sulfate and ammonium
1634         real(r8), parameter :: mw_so4a      =  96.0
1635         real(r8), parameter :: mw_nh4a      =  18.0
1636
1637         real(r8), save :: reldiffmax = 0.0
1638
1639         real(r8) ch2so4_crit              ! critical h2so4 conc (ug/m3)
1640         real(r8) dens_part                ! "grown" single-particle dry density (g/cm3)
1641         real(r8) duma, dumb, dumc, dume
1642         real(r8) dum_m1, dum_m2, dum_m3, dum_n1, dum_n2, dum_n3
1643         real(r8) fogas, foso4a, fonh4a, fonuma
1644         real(r8) mass_part                ! "grown" single-particle mass (g)
1645         real(r8) molenh4a_per_moleso4a    ! (mole aerosol nh4)/(mole aerosol so4)
1646         real(r8) qh2so4_crit              ! critical h2so4 mixrat (ppm)
1647         real(r8) qh2so4_avail             ! amount of h2so4 available for new particles (ppm)
1648         real(r8) vol_part                 ! "grown" single-particle volume (cm3)
1649
1650
1651 !
1652 ! initialization output arguments with "zero nucleation" values
1653 !
1654         isize_nuc = 1
1655         qnuma_del = 0.0
1656         qso4a_del = 0.0
1657         qnh4a_del = 0.0
1658         qh2so4_del = 0.0
1659         qnh3_del = 0.0
1660
1661 !
1662 ! calculate critical h2so4 concentration (ug/m3) and mixing ratio (mole/mole-air)
1663 !
1664         ch2so4_crit = 0.16 * exp( 0.1*temp_in - 3.5*rh_in - 27.7 )
1665 ! ch2so4 = (ug-h2so4/m3-air)
1666 ! ch2so4*1.0e-12/mwh2so4 = (mole-h2so4/cm3-air)
1667         qh2so4_crit = (ch2so4_crit*1.0e-12/98.0)/cair
1668         qh2so4_avail = qh2so4_cur - qh2so4_crit
1669
1670 ! if "available" h2so4 vapor < 4.0e-18 mole/mole-air ~= 1.0e2 molecules/cm3,
1671 ! exit with new particle formation = 0
1672         if (qh2so4_avail .le. 4.0e-18) then
1673            return
1674         end if
1675
1676 ! determine size bin into which the new particles go
1677         isize_nuc = 1
1678         vol_part = volumlo_sect(1)
1679
1680 !
1681 ! determine composition and density of the "grown particles"
1682 ! the grown particles are assumed to be liquid
1683 !    (since critical clusters contain water)
1684 !    so any (nh4/so4) molar ratio between 0 and 2 is allowed
1685 ! assume that the grown particles will have
1686 !    (nh4/so4 molar ratio) = min( 2, (nh3/h2so4 gas molar ratio) )
1687 !
1688         if (qnh3_cur .ge. qh2so4_avail) then
1689 ! combination of ammonium sulfate and ammonium bisulfate
1690 ! dum_n1 & dum_n2 = mole fractions of the ammsulf & ammbisulf
1691            dum_n1 = (qnh3_cur/qh2so4_avail) - 1.0
1692            dum_n1 = max( 0.0_r8, min( 1.0_r8, dum_n1 ) )
1693            dum_n2 = 1.0 - dum_n1
1694            dum_n3 = 0.0
1695         else
1696 ! combination of ammonium bisulfate and sulfuric acid
1697 ! dum_n2 & dum_n3 = mole fractions of the ammbisulf & sulfacid
1698            dum_n1 = 0.0
1699            dum_n2 = (qnh3_cur/qh2so4_avail)
1700            dum_n2 = max( 0.0_r8, min( 1.0_r8, dum_n2 ) )
1701            dum_n3 = 1.0 - dum_n2
1702         end if
1703
1704         dum_m1 = dum_n1*mw_ammsulf
1705         dum_m2 = dum_n2*mw_ammbisulf
1706         dum_m3 = dum_n3*mw_sulfacid
1707         dens_part = (dum_m1 + dum_m2 + dum_m3)/   &
1708            ((dum_m1/dens_ammsulf) + (dum_m2/dens_ammbisulf)   &
1709                                   + (dum_m3/dens_sulfacid))
1710 ! 25-jul-2006 - use 'in' value only if it is between 1.5-2.0 g/cm3
1711         if (abs(dens_nh4so4a-1.75) .le. 0.25) then
1712             dens_part = dens_nh4so4a
1713         else
1714             dens_nh4so4a = dens_part
1715         end if
1716         mass_part  = vol_part*dens_part
1717         molenh4a_per_moleso4a = 2.0*dum_n1 + dum_n2
1718
1719
1720 ! changes to h2so4 & nh3 gas (in mole/mole-air), limited by amounts available
1721         duma = 0.9999
1722         qh2so4_del = min( duma*qh2so4_cur, qh2so4_avail )
1723         qnh3_del   = min( duma*qnh3_cur, qh2so4_del*molenh4a_per_moleso4a )
1724         qh2so4_del = -qh2so4_del
1725         qnh3_del   = -qnh3_del
1726
1727 ! changes to so4 & nh4 aerosol (in mole/mole-air)
1728         qso4a_del = -qh2so4_del
1729         qnh4a_del =   -qnh3_del
1730 ! change to aerosol number (in #/mole-air)
1731         qnuma_del = (qso4a_del*mw_so4a + qnh4a_del*mw_nh4a)/mass_part
1732
1733
1734         return
1735         end subroutine wexler_nuc_mosaic_1box
1736
1737
1738
1739
1740 !-----------------------------------------------------------------------
1741
1742
1743
1744         end module module_mosaic_newnucb