!**********************************************************************************
! This computer software was prepared by Battelle Memorial Institute, hereinafter
! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
!
! MOSAIC module: see module_mosaic_driver.F for references and terms of use
!**********************************************************************************
module module_mosaic_newnuc
use module_peg_util
implicit none
contains
!-----------------------------------------------------------------------
subroutine mosaic_newnuc_1clm( istat_newnuc, &
it, jt, kclm_calcbgn, kclm_calcend, &
idiagbb_in, &
dtchem, dtnuc_in, rsub0, &
id, ktau, ktauc, its, ite, jts, jte, kts, kte )
!
! calculates new particle nucleation for grid points in
! the i=it, j=jt vertical column over timestep dtnuc_in
! works with mosaic sectional aerosol packages
!
! when newnuc_method = 1 (internal parameter), uses
! h2so4-nh3-h2o ternary nucleation from napari et al., 2002
! *** this option is currently disabled because the ternary
! parameterization was recently shown to be invalid
! when newnuc_method = 2 (internal parameter), uses
! h2so4-h2o binary nucleation from wexler et al., 1994
!
use module_data_mosaic_asect
use module_data_mosaic_other
use module_state_description, only: param_first_scalar
! subr arguments
integer, intent(inout) :: istat_newnuc ! =0 if no problems
integer, intent(in) :: &
it, jt, kclm_calcbgn, kclm_calcend, &
idiagbb_in, &
id, ktau, ktauc, its, ite, jts, jte, kts, kte
real, intent(in) :: dtchem, dtnuc_in
real, intent(in) :: rsub0(l2maxd,kmaxd,nsubareamaxd)
! rsub0 holds mixrat values before the aerchemistry calcs
! NOTE - much information is passed via arrays in
! module_data_mosaic_asect and module_data_mosaic_other
!
! rsub (inout) - trace gas and aerosol mixing ratios
! aqmassdry_sub, aqvoldry_sub (inout) -
! aerosol dry-mass and dry-volume mixing ratios
! adrydens_sub (inout) - aerosol dry density
! rsub(ktemp,:,:), cairclm, relhumclm (in) -
! air temperature, molar density, and relative humidity
! local variables
integer, parameter :: newnuc_method = 2
integer :: k, l, ll, m
integer :: isize, itype, iphase
integer :: iconform_numb
integer :: idiagbb
integer :: nsize, ntau_nuc
integer, save :: ncount(10)
integer :: p1st
real, parameter :: densdefault = 2.0
real, parameter :: smallmassbb = 1.0e-30
! nh4hso4 values for a_zsr and b_zsr
real, parameter :: a_zsr_xx1 = 1.15510
real, parameter :: a_zsr_xx2 = -3.20815
real, parameter :: a_zsr_xx3 = 2.71141
real, parameter :: a_zsr_xx4 = 2.01155
real, parameter :: a_zsr_xx5 = -4.71014
real, parameter :: a_zsr_xx6 = 2.04616
real, parameter :: b_zsr_xx = 29.4779
real :: aw
real :: cair_box
real :: dens_nh4so4a, dtnuc
real :: duma, dumb, dumc
real :: rh_box
real :: qh2so4_avg, qh2so4_cur, qh2so4_del
real :: qnh3_avg, qnh3_cur, qnh3_del
real :: qnuma_del, qso4a_del, qnh4a_del
real :: temp_box
real :: xxdens, xxmass, xxnumb, xxvolu
real,save :: dumveca(10), dumvecb(10), dumvecc(10), dumvecd(10), dumvece(10)
real :: volumlo_nuc(maxd_asize), volumhi_nuc(maxd_asize)
character(len=100) :: msg
p1st = PARAM_FIRST_SCALAR
! check newnuc_method
istat_newnuc = 0
if (newnuc_method .ne. 2) then
if ((it .eq. its) .and. (jt .eq. jts)) &
call peg_message( lunerr, &
'*** mosaic_newnuc_1clm -- illegal newnuc_method' )
istat_newnuc = -1
return
end if
! when dtnuc_in > dtchem, do not perform nucleation calcs
! on every chemistry time step
ntau_nuc = nint( dtnuc_in/dtchem )
ntau_nuc = max( 1, ntau_nuc )
if (mod(ktau,ntau_nuc) .ne. 0) return
dtnuc = dtchem*ntau_nuc
! set variables that do not change
idiagbb = idiagbb_in
itype = 1
iphase = ai_phase
nsize = nsize_aer(itype)
volumlo_nuc(1:nsize) = volumlo_sect(1:nsize,itype)
volumhi_nuc(1:nsize) = volumhi_sect(1:nsize,itype)
! loop over subareas (currently only 1) and vertical levels
do 2900 m = 1, nsubareas
do 2800 k = kclm_calcbgn, kclm_calcend
! initialize diagnostics
if ((it .eq. its) .and. &
(jt .eq. jts) .and. (k .eq. kclm_calcbgn)) then
dumveca(:) = 0.0 ! current grid param values
dumvecb(:) = +1.0e35 ! param minimums
dumvecc(:) = -1.0e35 ! param maximums
dumvecd(:) = 0.0 ! param averages
dumvece(:) = 0.0 ! param values for highest qnuma_del
ncount(:) = 0
end if
ncount(1) = ncount(1) + 1
if (afracsubarea(k,m) .lt. 1.e-4) goto 2700
cair_box = cairclm(k)
temp_box = rsub(ktemp,k,m)
rh_box = relhumclm(k)
qh2so4_cur = max(0.0,rsub(kh2so4,k,m))
qnh3_cur = max(0.0,rsub(knh3,k,m))
qh2so4_avg = 0.5*( qh2so4_cur + max(0.0,rsub0(kh2so4,k,m)) )
qnh3_avg = 0.5*( qnh3_cur + max(0.0,rsub0(knh3,k,m)) )
qh2so4_del = 0.0
qnh3_del = 0.0
qnuma_del = 0.0
qso4a_del = 0.0
qnh4a_del = 0.0
dens_nh4so4a = dens_so4_aer
isize = 0
! make call to nucleation routine
if (newnuc_method .eq. 1) then
call ternary_nuc_mosaic_1box( &
dtnuc, temp_box, rh_box, cair_box, &
qh2so4_avg, qh2so4_cur, qnh3_avg, qnh3_cur, &
nsize, maxd_asize, volumlo_nuc, volumhi_nuc, &
isize, qnuma_del, qso4a_del, qnh4a_del, &
qh2so4_del, qnh3_del, dens_nh4so4a )
else if (newnuc_method .eq. 2) then
call wexler_nuc_mosaic_1box( &
dtnuc, temp_box, rh_box, cair_box, &
qh2so4_avg, qh2so4_cur, qnh3_avg, qnh3_cur, &
nsize, maxd_asize, volumlo_nuc, volumhi_nuc, &
isize, qnuma_del, qso4a_del, qnh4a_del, &
qh2so4_del, qnh3_del, dens_nh4so4a )
else
istat_newnuc = -1
return
end if
! temporary diagnostics
dumveca(1) = temp_box
dumveca(2) = rh_box
dumveca(3) = rsub(kso2,k,m)
dumveca(4) = qh2so4_avg
dumveca(5) = qnh3_avg
dumveca(6) = qnuma_del
do l = 1, 6
dumvecb(l) = min( dumvecb(l), dumveca(l) )
dumvecc(l) = max( dumvecc(l), dumveca(l) )
dumvecd(l) = dumvecd(l) + dumveca(l)
if (qnuma_del .gt. dumvece(6)) dumvece(l) = dumveca(l)
end do
! check for zero new particles
if (qnuma_del .le. 0.0) goto 2700
! check for valid isize
if (isize .ne. 1) ncount(3) = ncount(3) + 1
if ((isize .lt. 1) .or. (isize .gt. nsize)) then
write(msg,93010) 'newnucxx bad isize_nuc' , it, jt, k, &
isize, nsize
call peg_message( lunerr, msg )
goto 2700
end if
93010 format( a, 3i3, 1p, 9e10.2 )
ncount(2) = ncount(2) + 1
! update gas and aerosol so4 and nh3/nh4 mixing ratios
rsub(kh2so4,k,m) = max( 0.0, rsub(kh2so4,k,m) + qh2so4_del )
rsub(knh3, k,m) = max( 0.0, rsub(knh3, k,m) + qnh3_del )
l = lptr_so4_aer(isize,itype,iphase)
if (l .ge. p1st) then
rsub(l,k,m) = rsub(l,k,m) + qso4a_del
end if
l = lptr_nh4_aer(isize,itype,iphase)
if (l .ge. p1st) then
rsub(l,k,m) = rsub(l,k,m) + qnh4a_del
end if
l = numptr_aer(isize,itype,iphase)
rsub(l,k,m) = rsub(l,k,m) + qnuma_del
xxnumb = rsub(l,k,m)
! update aerosol water, using mosaic parameterizations for nh4hso4
! duma = (mole-salt)/(mole-salt+water)
! dumb = (mole-salt)/(kg-water)
! dumc = (mole-water)/(mole-salt)
l = waterptr_aer(isize,itype)
if ((rh_box .gt. 0.10) .and. (l .ge. p1st)) then
aw = min( rh_box, 0.98 )
if (aw .lt. 0.97) then
duma = a_zsr_xx1 + &
aw*( a_zsr_xx2 + &
aw*( a_zsr_xx3 + &
aw*( a_zsr_xx4 + &
aw*( a_zsr_xx5 + &
aw* a_zsr_xx6 ))))
else
dumb = -b_zsr_xx*log(aw)
dumb = max( dumb, 0.5 )
duma = 1.0/(1.0 + 55.509/dumb)
end if
duma = max( duma, 0.01 )
dumc = (1.0 - duma)/duma
rsub(l,k,m) = rsub(l,k,m) + qso4a_del*dumc
end if
!
! update dry mass, density, and volume,
! and check for mean dry-size within bounds
!
xxmass = aqmassdry_sub(isize,itype,k,m)
xxdens = adrydens_sub( isize,itype,k,m)
iconform_numb = 1
if ((xxdens .lt. 0.1) .or. (xxdens .gt. 20.0)) then
! (exception) case of drydensity not valid
continue
else
! (normal) case of drydensity valid (which means drymass is valid also)
! so increment mass and volume with the so4 & nh4 deltas, then calc density
xxvolu = xxmass/xxdens
duma = qso4a_del*mw_so4_aer + qnh4a_del*mw_nh4_aer
xxmass = xxmass + duma
xxvolu = xxvolu + duma/dens_nh4so4a
if (xxmass .le. smallmassbb) then
! do this to force calc of dry mass, volume from rsub
xxdens = 0.001
else if (xxmass .gt. 1000.0*xxvolu) then
! in this case, density is too large. setting density=1000 forces
! next IF block while avoiding potential divide by zero or overflow
xxdens = 1000.0
else
xxdens = xxmass/xxvolu
end if
end if
if ((xxdens .lt. 0.1) .or. (xxdens .gt. 20.0)) then
! (exception) case of drydensity not valid (or drymass extremely small),
! so compute from dry mass, volume from rsub
ncount(4) = ncount(4) + 1
xxmass = 0.0
xxvolu = 0.0
do ll = 1, ncomp_aer(itype)
l = massptr_aer(ll,isize,itype,iphase)
if (l .ge. p1st) then
duma = max( 0.0, rsub(l,k,m) )*mw_aer(ll,itype)
xxmass = xxmass + duma
xxvolu = xxvolu + duma/dens_aer(ll,itype)
end if
end do
end if
if (xxmass .le. smallmassbb) then
! when drymass extremely small, use default density and bin center size,
! and zero out water
ncount(5) = ncount(5) + 1
xxdens = densdefault
xxvolu = xxmass/xxdens
xxnumb = xxmass/(volumcen_sect(isize,itype)*xxdens)
iconform_numb = 0
l = waterptr_aer(isize,itype)
if (l .ge. p1st) rsub(l,k,m) = 0.0
l = hyswptr_aer(isize,itype)
if (l .ge. p1st) rsub(l,k,m) = 0.0
else
xxdens = xxmass/xxvolu
end if
if (iconform_numb .gt. 0) then
! check for mean dry-size within bounds, and conform number if not
if (xxnumb .gt. xxvolu/volumlo_sect(isize,itype)) then
ncount(6) = ncount(6) + 1
xxnumb = xxvolu/volumlo_sect(isize,itype)
else if (xxnumb .lt. xxvolu/volumhi_sect(isize,itype)) then
ncount(7) = ncount(7) + 1
xxnumb = xxvolu/volumhi_sect(isize,itype)
end if
end if
! load dry mass, density, volume, and (possibly conformed) number
l = numptr_aer(isize,itype,iphase)
rsub(l,k,m) = xxnumb
adrydens_sub( isize,itype,k,m) = xxdens
aqmassdry_sub(isize,itype,k,m) = xxmass
aqvoldry_sub( isize,itype,k,m) = xxvolu
2700 continue
! temporary diagnostics
if ((idiagbb .ge. 100) .and. &
(it .eq. ite) .and. &
(jt .eq. jte) .and. (k .eq. kclm_calcend)) then
if (idiagbb .ge. 110) then
write(msg,93020) 'newnucbb mins ', dumvecb(1:6)
call peg_message( lunerr, msg )
write(msg,93020) 'newnucbb maxs ', dumvecc(1:6)
call peg_message( lunerr, msg )
duma = max( 1, ncount(1) )
write(msg,93020) 'newnucbb avgs ', dumvecd(1:6)/duma
call peg_message( lunerr, msg )
write(msg,93020) 'newnucbb hinuc', dumvece(1:6)
call peg_message( lunerr, msg )
write(msg,93020) 'newnucbb dtnuc', dtnuc
call peg_message( lunerr, msg )
end if
write(msg,93030) 'newnucbb ncnt ', ncount(1:7)
call peg_message( lunerr, msg )
end if
93020 format( a, 1p, 10e10.2 )
93030 format( a, 1p, 10i10 )
2800 continue ! k levels
2900 continue ! subareas
return
end subroutine mosaic_newnuc_1clm
!-----------------------------------------------------------------------
! note - subrs ternary_nuc_mosaic_1box, ternary_nuc_napari, wexler_nuc_mosaic_1box
! taken from ternucl04.f90 on 22-jul-2006
!-----------------------------------------------------------------------
subroutine ternary_nuc_mosaic_1box( &
dtnuc, temp_in, rh_in, cair, &
qh2so4_avg, qh2so4_cur, qnh3_avg, qnh3_cur, &
nsize, maxd_asize, volumlo_sect, volumhi_sect, &
isize_nuc, qnuma_del, qso4a_del, qnh4a_del, &
qh2so4_del, qnh3_del, dens_nh4so4a )
!.......................................................................
!
! calculates new particle production from h2so4-nh3-h2o ternary nucleation
! over timestep dtnuc, using nucleation rates from the
! napari et al. (2002) parameterization
!
! the new particles are "grown" to the lower-bound size of the host code's
! smallest size bin. (this "growth" is purely ad hoc, and would not be
! necessary if the host code's size bins extend down to ~1 nm.)
! if the h2so4 and nh3 mass mixing ratios (mixrats) of the grown new
! particles exceed the current gas mixrats, the new particle production
! is reduced so that the new particle mass mixrats match the gas mixrats.
!
! revision history
! coded by rc easter, pnnl, 20-mar-2006
!
! key routines called: subr ternary_nuc_napari
!
! references:
! napari, i., m. noppel, h. vehkamaki, and m. kulmala,
! parameterization of ternary nucleation rates for
! h2so4-nh3-h2o vapors. j. geophys. res., 107, 4381, 2002.
!
!.......................................................................
implicit none
! subr arguments (in)
real, intent(in) :: dtnuc ! nucleation time step (s)
real, intent(in) :: temp_in ! temperature, in k
real, intent(in) :: rh_in ! relative humidity, as fraction
real, intent(in) :: cair ! dry-air molar density (mole/cm3)
real, intent(in) :: qh2so4_avg, qh2so4_cur ! gas h2so4 mixing ratios (mole/mole-air)
real, intent(in) :: qnh3_avg, qnh3_cur ! gas nh3 mixing ratios (mole/mole-air)
! qxxx_cur = current value (at end of condensation)
! qxxx_avg = average value (from start to end of condensation)
integer, intent(in) :: nsize ! number of aerosol size bins
integer, intent(in) :: maxd_asize ! dimension for volumlo_sect, ...
real, intent(in) :: volumlo_sect(maxd_asize) ! dry volume at lower bnd of bin (cm3)
real, intent(in) :: volumhi_sect(maxd_asize) ! dry volume at upper bnd of bin (cm3)
! subr arguments (out)
integer, intent(out) :: isize_nuc ! size bin into which new particles go
real, intent(out) :: qnuma_del ! change to aerosol number mixing ratio (#/mole-air)
real, intent(out) :: qso4a_del ! change to aerosol so4 mixing ratio (mole/mole-air)
real, intent(out) :: qnh4a_del ! change to aerosol nh4 mixing ratio (mole/mole-air)
real, intent(out) :: qh2so4_del ! change to gas h2so4 mixing ratio (mole/mole-air)
real, intent(out) :: qnh3_del ! change to gas nh3 mixing ratio (mole/mole-air)
! aerosol changes are > 0; gas changes are < 0
! subr arguments (inout)
real, intent(inout) :: dens_nh4so4a ! dry-density of the new nh4-so4 aerosol mass (g/cm3)
! use 'in' value only if it is between 1.6-2.0 g/cm3
! subr arguments (out) passed via common block
! these are used to duplicate the outputs of yang zhang's original test driver
! they are not really needed in wrf-chem, so just make them local
real :: ratenuclt ! j = ternary nucleation rate from napari param. (cm-3 s-1)
real :: rateloge ! ln (j)
real :: cnum_h2so4 ! number of h2so4 molecules in the critical nucleus
real :: cnum_nh3 ! number of nh3 molecules in the critical nucleus
real :: cnum_tot ! total number of molecules in the critical nucleus
real :: radius_cluster ! the radius of cluster (nm)
! common / ternary_nuc_yzhang_cmn01 / &
! ratenuclt, rateloge, &
! cnum_h2so4, cnum_nh3, cnum_tot, radius_cluster
! local variables
integer i
integer, save :: icase = 0, icase_reldiffmax = 0
real, parameter :: pi = 3.1415926536
real, parameter :: avogad = 6.022e23 ! avogadro number (molecules/mole)
real, parameter :: mw_air = 28.966 ! dry-air mean molecular weight (g/mole)
! dry densities (g/cm3) molecular weights of aerosol
! ammsulf, ammbisulf, and sulfacid (from mosaic dens_electrolyte values)
real, parameter :: dens_ammsulf = 1.769
real, parameter :: dens_ammbisulf = 1.78
real, parameter :: dens_sulfacid = 1.841
! molecular weights (g/mole) of aerosol ammsulf, ammbisulf, and sulfacid
! for ammbisulf and sulfacid, use 114 & 96 here rather than 115 & 98
! because we don't keep track of aerosol hion mass
real, parameter :: mw_ammsulf = 132.0
real, parameter :: mw_ammbisulf = 114.0
real, parameter :: mw_sulfacid = 96.0
! molecular weights of aerosol sulfate and ammonium
real, parameter :: mw_so4a = 96.0
real, parameter :: mw_nh4a = 18.0
real, save :: reldiffmax = 0.0
real dens_part ! "grown" single-particle dry density (g/cm3)
real duma, dumb, dumc, dume
real dum_m1, dum_m2, dum_m3, dum_n1, dum_n2, dum_n3
real fogas, foso4a, fonh4a, fonuma
real freduce ! reduction factor applied to nucleation rate
! due to limited availability of h2so4 & nh3 gases
real freducea, freduceb
real gramaero_per_moleso4a ! (g dry aerosol)/(mole aerosol so4)
real mass_part ! "grown" single-particle mass (g)
real molenh4a_per_moleso4a ! (mole aerosol nh4)/(mole aerosol so4)
real nh3conc_in ! mixing ratio of nh3 for nucl. calc., pptv
real so4vol_in ! concentration of h2so4 for nucl. calc., molecules cm-3
real qmolnh4a_del_max ! max production of aerosol nh4 over dtnuc (mole/mole-air)
real qmolso4a_del_max ! max production of aerosol so4 over dtnuc (mole/mole-air)
real vol_cluster ! critical-cluster volume (cm3)
real vol_part ! "grown" single-particle volume (cm3)
!
! if h2so4 vapor < 4.0e-16 mole/moleair ~= 1.0e4 molecules/cm3,
! exit with new particle formation = 0
!
isize_nuc = 1
qnuma_del = 0.0
qso4a_del = 0.0
qnh4a_del = 0.0
qh2so4_del = 0.0
qnh3_del = 0.0
if (qh2so4_avg .le. 4.0e-16) return
if (qh2so4_cur .le. 4.0e-16) return
!
! make call to napari parameterization routine
!
! convert nh3 from mole/mole-air to ppt
! qnh3_cur = nh3conc(m)*1.0e-12
nh3conc_in = qnh3_avg/1.0e-12
! convert h2so4 from mole/mole-air to molecules/cm3
! qh2so4_cur = (so4vol(k)/avogad) / cair
so4vol_in = (qh2so4_avg) * cair * avogad
call ternary_nuc_napari( &
temp_in, rh_in, nh3conc_in, so4vol_in, &
ratenuclt, rateloge, &
cnum_h2so4, cnum_nh3, cnum_tot, radius_cluster )
! if nucleation rate is less than 0.1 particle/cm3/day,
! exit with new particle formation = 0
if (ratenuclt .le. 1.0e-6) return
! determine size bin into which the new particles go
! (probably it will always be bin #1, but ...)
vol_cluster = (pi*4.0/3.0)* (radius_cluster**3) * 1.0e-21
isize_nuc = 1
vol_part = volumlo_sect(1)
if (vol_cluster .le. volumlo_sect(1)) then
continue
else if (vol_cluster .ge. volumhi_sect(nsize)) then
isize_nuc = nsize
vol_part = volumhi_sect(nsize)
else
do i = 1, nsize
if (vol_cluster .lt. volumhi_sect(i)) then
isize_nuc = i
vol_part = vol_cluster
vol_part = min( vol_part, volumhi_sect(i) )
vol_part = max( vol_part, volumlo_sect(i) )
exit
end if
end do
end if
!
! determine composition and density of the "grown particles"
! the grown particles are assumed to be liquid
! (since critical clusters contain water)
! so any (nh4/so4) molar ratio between 0 and 2 is allowed
! assume that the grown particles will have
! (nh4/so4 molar ratio) = min( 2, (nh3/h2so4 gas molar ratio) )
!
if (qnh3_cur .ge. qh2so4_cur) then
! combination of ammonium sulfate and ammonium bisulfate
! dum_n1 & dum_n2 = mole fractions of the ammsulf & ammbisulf
dum_n1 = (qnh3_cur/qh2so4_cur) - 1.0
dum_n1 = max( 0.0, min( 1.0, dum_n1 ) )
dum_n2 = 1.0 - dum_n1
dum_n3 = 0.0
else
! combination of ammonium bisulfate and sulfuric acid
! dum_n2 & dum_n3 = mole fractions of the ammbisulf & sulfacid
dum_n1 = 0.0
dum_n2 = (qnh3_cur/qh2so4_cur)
dum_n2 = max( 0.0, min( 1.0, dum_n2 ) )
dum_n3 = 1.0 - dum_n2
end if
dum_m1 = dum_n1*mw_ammsulf
dum_m2 = dum_n2*mw_ammbisulf
dum_m3 = dum_n3*mw_sulfacid
dens_part = (dum_m1 + dum_m2 + dum_m3)/ &
((dum_m1/dens_ammsulf) + (dum_m2/dens_ammbisulf) &
+ (dum_m3/dens_sulfacid))
! 25-jul-2006 - use 'in' value only if it is between 1.6-2.0 g/cm3
if (abs(dens_nh4so4a-1.8) .le. 0.2) then
dens_part = dens_nh4so4a
else
dens_nh4so4a = dens_part
end if
mass_part = vol_part*dens_part
molenh4a_per_moleso4a = 2.0*dum_n1 + dum_n2 ! (mole aerosol nh4)/(mole aerosol so4)
gramaero_per_moleso4a = dum_m1 + dum_m2 + dum_m3 ! (g dry aerosol)/(mole aerosol so4)
! max production of aerosol dry mass (g-aero/cm3-air)
duma = max( 0.0, (ratenuclt*dtnuc*mass_part) )
! max production of aerosol so4 (mole-so4a/cm3-air)
dumc = duma/gramaero_per_moleso4a
! max production of aerosol so4 (mole-so4a/mole-air)
dume = dumc/cair
! max production of aerosol so4 (mole/mole-air)
! based on ratenuclt and mass_part
qmolso4a_del_max = dume
! check if max production exceeds available h2so4 vapor
freducea = 1.0
if (qmolso4a_del_max .gt. qh2so4_cur) then
freducea = qh2so4_cur/qmolso4a_del_max
end if
! check if max production exceeds available nh3 vapor
freduceb = 1.0
if (molenh4a_per_moleso4a .ge. 1.0e-10) then
! max production of aerosol nh4 (ppm) based on ratenuclt and mass_part
qmolnh4a_del_max = qmolso4a_del_max*molenh4a_per_moleso4a
if (qmolnh4a_del_max .gt. qnh3_cur) then
freduceb = qnh3_cur/qmolnh4a_del_max
end if
end if
freduce = min( freducea, freduceb )
! if adjusted nucleation rate is less than 0.1 particle/cm3/day,
! exit with new particle formation = 0
if (freduce*ratenuclt .le. 1.0e-6) return
! note: suppose that at this point, freduce = 1.0 (no gas-available
! constraints) and molenh4a_per_moleso4a < 2.0
! then it would be possible to condense "additional" nh3 and have
! (nh3/h2so4 gas molar ratio) < (nh4/so4 aerosol molar ratio) <= 2
! one could do some additional calculations of
! dens_part & molenh4a_per_moleso4a to realize this
! however, the particle "growing" is a crude approximate way to get
! the new particles to the host code's minimum particle size,
! so are such refinements worth the effort?
! changes to h2so4 & nh3 gas (in mole/mole-air), limited by amounts available
duma = 0.9999
qh2so4_del = min( duma*qh2so4_cur, freduce*qmolso4a_del_max )
qnh3_del = min( duma*qnh3_cur, qh2so4_del*molenh4a_per_moleso4a )
qh2so4_del = -qh2so4_del
qnh3_del = -qnh3_del
! changes to so4 & nh4 aerosol (in mole/mole-air)
qso4a_del = -qh2so4_del
qnh4a_del = -qnh3_del
! change to aerosol number (in #/mole-air)
qnuma_del = (qso4a_del*mw_so4a + qnh4a_del*mw_nh4a)/mass_part
return
end subroutine ternary_nuc_mosaic_1box
!-----------------------------------------------------------------------
subroutine ternary_nuc_napari( &
temp_in, rh_in, nh3conc_in, so4vol_in, &
ratenuclt, rateloge, &
cnum_h2so4, cnum_nh3, cnum_tot, radius_cluster )
!*********************************************************************
! purpose: calculate ternary nucleation rates based on *
! napari et al. (2002) parameterization *
! *
! revision history: *
! coded by yang zhang, ncsu, nov. 25, 2004 *
! *
! 14-mar-2006 rce - yang sent this on 10-mar-2005; *
! converted to lowercase; *
! moved the nucleation rate calcs *
! from the main program unit to this subr; *
! added temporary variables log_rh, log_nh3conc, log_so4vol; *
! in the nuc subr, apply limits to the input parameters; *
! converted to f90; *
! *
! reference: *
! napari, i., m, noppel, h. vehkamaki, . and . m. kulmala, *
! parameterization of ternary nucleation rates for *
! h2so4-nh3-h2o vapors. j. geophys. res., 107, 4381, 2002b. *
! *
! note: *
! advantage of this parameterization *
! this parameterization reproduces the ternary nucleation *
! rates obtained from the full model within the range of *
! one order of magnitude, with a cpu saving by a factor *
! of 10e5. *
! *
! limitations / assumptions of this parameterization *
! 1. the limits of validity are *
! t = 240-300 k, rh=0.05-0.95 *
! [h2so4]=10e4-10e9 cm-3 (4e-4 - 40 ppt at stp) *
! [nh3]=0.1-100 ppt, *
! j=10e-5 - 10e6 cm-3 s-1 *
! 2. it cannot be used to obtain binary h2o-h2so4 or h2o-nh3 *
! limit, due to the logarithmic dependencies on rh, *
! [h2so4], and [nh3] *
! 3. the fit is the worst at high temperatures ( > 298 k) *
! and low nucleation rates ( < 0.01 cm-3 s-1), but it is *
! more accurate at significant nucleation rates (0.01-0.1 *
! cm-3 s-1), thus adequate for simulating ternary *
! nucleation in the atmosphere *
! *
!*********************************************************************
implicit none
! subr arguments (in)
real, intent(in) :: temp_in ! temperature, in k
real, intent(in) :: rh_in ! relative humidity, as fraction
real, intent(in) :: nh3conc_in ! mixing ratio of nh3, pptv
real, intent(in) :: so4vol_in ! concentration of h2so4, cm-3
! subr arguments (out)
real, intent(out) :: ratenuclt ! ternary nucleation rate, j, cm-3 s-1
real, intent(out) :: rateloge ! ln (j)
real, intent(out) :: cnum_h2so4 ! number of h2so4 molecules
! in the critical nucleus
real, intent(out) :: cnum_nh3 ! number of nh3 molecules
! in the critical nucleus
real, intent(out) :: cnum_tot ! total number of molecules
! in the critical nucleus
real, intent(out) :: radius_cluster ! the radius of cluster, nm
! local variables
integer ncoeff
parameter ( ncoeff = 4 ) ! total fitting coefficients at each t
! corresponds to 3rd order polynomial
integer npoly
parameter ( npoly = 20 ) ! total number of polynomial functions
integer n ! loop index for functions of polynomial
! polynomial functions
real f ( npoly )
! coefficients of polynomials
real a ( ncoeff, npoly )
real temp, rh, nh3conc, so4vol ! bounded values of the input args
real log_rh, log_nh3conc, log_so4vol
! coefficients of polynomials fi(t), i = 1, 20
data a / -0.355297, -33.8449, 0.34536, -8.24007e-4, &
3.13735, -0.772861, 5.61204e-3, -9.74576e-6, &
19.0359, -0.170957, 4.79808e-4, -4.14699e-7, &
1.07605, 1.48932, -7.96052e-3, 7.61229e-6, &
6.0916, -1.25378, 9.39836e-3, -1.74927e-5, &
0.31176, 1.64009, -3.43852e-3, -1.09753e-5, &
-2.00738e-2, -0.752115, 5.25813e-3, -8.98038e-6, &
0.165536, 3.26623, -4.89703e-2, 1.46967e-4, &
6.52645, -0.258002, 1.43456e-3, -2.02036e-6, &
3.68024, -0.204098, 1.06259e-3, -1.2656e-6 , &
-6.6514e-2, -7.82382, 1.22938e-2, 6.18554e-5, &
0.65874, 0.190542, -1.65718e-3, 3.41744e-6, &
5.99321e-2, 5.96475, -3.62432e-2, 4.93337e-5, &
-0.732731, -1.84179e-2, 1.47186e-4, -2.37711e-7, &
0.728429, 3.64736, -2.7422e-2, 4.93478e-5, &
41.3016, -0.35752, 9.04383e-4, -5.73788e-7, &
-0.160336, 8.89881e-3, -5.39514e-5, 8.39522e-8, &
8.57868, -0.112358, 4.72626e-4, -6.48365e-7, &
5.30167e-2, -1.98815, 1.57827e-2, -2.93564e-5, &
-2.32736, 2.34646e-2, -7.6519e-5, 8.0459e-8 /
!
! apply limits to input parameters
!
temp = max( 240.15, min (300.15, temp_in ) )
rh = max( 0.05, min (0.95, rh_in ) )
so4vol = max( 1.0e4, min (1.0e9, so4vol_in ) )
nh3conc = max( 0.1, min (100.0, nh3conc_in ) )
!
! calculate the functions of polynomials, fi(t), i = 1, npoly
! at temperature j
!
do n = 1, npoly
f ( n ) = a ( 1, n ) + a ( 2, n ) * temp &
+ a ( 3, n ) * ( temp ) ** 2.0 &
+ a ( 4, n ) * ( temp ) ** 3.0
end do
!
! calculate ln (j)
!
log_rh = log ( rh )
log_nh3conc = log ( nh3conc )
log_so4vol = log ( so4vol )
rateloge = -84.7551 &
+ f ( 1 ) / log_so4vol &
+ f ( 2 ) * ( log_so4vol ) &
+ f ( 3 ) * ( log_so4vol ) **2.0 &
+ f ( 4 ) * ( log_nh3conc ) &
+ f ( 5 ) * ( log_nh3conc ) **2.0 &
+ f ( 6 ) * rh &
+ f ( 7 ) * ( log_rh ) &
+ f ( 8 ) * ( log_nh3conc / &
log_so4vol ) &
+ f ( 9 ) * ( log_nh3conc * &
log_so4vol ) &
+ f ( 10 ) * rh * &
( log_so4vol ) &
+ f ( 11 ) * ( rh / &
log_so4vol ) &
+ f ( 12 ) * ( rh * &
log_nh3conc ) &
+ f ( 13 ) * ( log_rh / &
log_so4vol ) &
+ f ( 14 ) * ( log_rh * &
log_nh3conc ) &
+ f ( 15 ) * (( log_nh3conc ) ** 2.0 &
/ log_so4vol ) &
+ f ( 16 ) * ( log_so4vol * &
( log_nh3conc ) ** 2.0 ) &
+ f ( 17 ) * (( log_so4vol ) ** 2.0 * &
log_nh3conc ) &
+ f ( 18 ) * ( rh * &
( log_nh3conc ) ** 2.0 ) &
+ f ( 19 ) * ( rh * log_nh3conc &
/ log_so4vol ) &
+ f ( 20 ) * (( log_so4vol ) ** 2.0 * &
( log_nh3conc ) ** 2.0 )
ratenuclt = exp ( rateloge )
!
! calculate number of molecules and critical radius of the cluster
!
cnum_h2so4 = 38.1645 + 0.774106 * rateloge &
+ 0.00298879 * ( rateloge ) ** 2.0 &
- 0.357605 * temp &
- 0.00366358 * temp * rateloge &
+ 0.0008553 * ( temp ) ** 2.0
cnum_nh3 = 26.8982 + 0.682905 * rateloge &
+ 0.00357521 * ( rateloge ) ** 2.0 &
- 0.265748 * temp &
- 0.00341895 * temp * rateloge &
+ 0.000673454 * ( temp ) ** 2.0
cnum_tot = 79.3484 + 1.7384 * rateloge &
+ 0.00711403 * ( rateloge ) ** 2.0 &
- 0.744993 * temp &
- 0.00820608 * temp * rateloge &
+ 0.0017855 * ( temp ) ** 2.0
radius_cluster = 0.141027 - 0.00122625 * rateloge &
- 7.82211e-6 * ( rateloge ) ** 2.0 &
- 0.00156727 * temp &
- 0.00003076 * temp * rateloge &
+ 0.0000108375 * ( temp ) ** 2.0
return
end subroutine ternary_nuc_napari
!-----------------------------------------------------------------------
subroutine wexler_nuc_mosaic_1box( &
dtnuc, temp_in, rh_in, cair, &
qh2so4_avg, qh2so4_cur, qnh3_avg, qnh3_cur, &
nsize, maxd_asize, volumlo_sect, volumhi_sect, &
isize_nuc, qnuma_del, qso4a_del, qnh4a_del, &
qh2so4_del, qnh3_del, dens_nh4so4a )
!.......................................................................
!
! calculates new particle production from h2so4-h2o binary nucleation
! over timestep dtnuc, using the wexler et al. (1994) parameterization
!
! the size of new particles is the lower-bound size of the host code's
! smallest size bin. their composition is so4 and nh4, since the nuclei
! would incorporate nh3 as they grow from ~1 nm to the lower-bound size.
! (the new particle composition can be forced to pure so4 by setting
! the qnh3_avg & qnh3_cur input arguments to 0.0).
!
! revision history
! coded by rc easter, pnnl, 20-mar-2006
!
! key routines called: none
!
! references:
! wexler, a. s., f. w. lurmann, and j. h. seinfeld,
! modelling urban and regional aerosols -- i. model development,
! atmos. environ., 28, 531-546, 1994.
!
!.......................................................................
implicit none
! subr arguments (in)
real, intent(in) :: dtnuc ! nucleation time step (s)
real, intent(in) :: temp_in ! temperature, in k
real, intent(in) :: rh_in ! relative humidity, as fraction
real, intent(in) :: cair ! dry-air molar density (mole-air/cm3)
real, intent(in) :: qh2so4_avg, qh2so4_cur ! gas h2so4 mixing ratios (ppm)
real, intent(in) :: qnh3_avg, qnh3_cur ! gas nh3 mixing ratios (ppm)
! qxxx_cur = current value (at end of condensation)
! qxxx_avg = average value (from start to end of condensation)
integer, intent(in) :: nsize ! number of aerosol size bins
integer, intent(in) :: maxd_asize ! dimension for volumlo_sect, ...
real, intent(in) :: volumlo_sect(maxd_asize) ! dry volume at lower bnd of bin (cm3)
real, intent(in) :: volumhi_sect(maxd_asize) ! dry volume at upper bnd of bin (cm3)
! subr arguments (out)
integer, intent(out) :: isize_nuc ! size bin into which new particles go
real, intent(out) :: qnuma_del ! change to aerosol number mixing ratio (#/kg)
real, intent(out) :: qso4a_del ! change to aerosol so4 mixing ratio (ug/kg)
real, intent(out) :: qnh4a_del ! change to aerosol nh4 mixing ratio (ug/kg)
real, intent(out) :: qh2so4_del ! change to gas h2so4 mixing ratio (ppm)
real, intent(out) :: qnh3_del ! change to gas nh3 mixing ratio (ppm)
! aerosol changes are > 0; gas changes are < 0
! subr arguments (inout)
real, intent(inout) :: dens_nh4so4a ! dry-density of the new nh4-so4 aerosol mass (g/cm3)
! use 'in' value only if it is between 1.6-2.0 g/cm3
! local variables
integer i
integer, save :: icase = 0, icase_reldiffmax = 0
real, parameter :: pi = 3.1415926536
real, parameter :: avogad = 6.022e23 ! avogadro number (molecules/mole)
real, parameter :: mw_air = 28.966 ! dry-air mean molecular weight (g/mole)
! dry densities (g/cm3) molecular weights of aerosol
! ammsulf, ammbisulf, and sulfacid (from mosaic dens_electrolyte values)
real, parameter :: dens_ammsulf = 1.769
real, parameter :: dens_ammbisulf = 1.78
real, parameter :: dens_sulfacid = 1.841
! molecular weights (g/mole) of aerosol ammsulf, ammbisulf, and sulfacid
! for ammbisulf and sulfacid, use 114 & 96 here rather than 115 & 98
! because we don't keep track of aerosol hion mass
real, parameter :: mw_ammsulf = 132.0
real, parameter :: mw_ammbisulf = 114.0
real, parameter :: mw_sulfacid = 96.0
! molecular weights of aerosol sulfate and ammonium
real, parameter :: mw_so4a = 96.0
real, parameter :: mw_nh4a = 18.0
real, save :: reldiffmax = 0.0
real ch2so4_crit ! critical h2so4 conc (ug/m3)
real dens_part ! "grown" single-particle dry density (g/cm3)
real duma, dumb, dumc, dume
real dum_m1, dum_m2, dum_m3, dum_n1, dum_n2, dum_n3
real fogas, foso4a, fonh4a, fonuma
real mass_part ! "grown" single-particle mass (g)
real molenh4a_per_moleso4a ! (mole aerosol nh4)/(mole aerosol so4)
real qh2so4_crit ! critical h2so4 mixrat (ppm)
real qh2so4_avail ! amount of h2so4 available for new particles (ppm)
real vol_part ! "grown" single-particle volume (cm3)
!
! initialization output arguments with "zero nucleation" values
!
isize_nuc = 1
qnuma_del = 0.0
qso4a_del = 0.0
qnh4a_del = 0.0
qh2so4_del = 0.0
qnh3_del = 0.0
!
! calculate critical h2so4 concentration (ug/m3) and mixing ratio (mole/mole-air)
!
ch2so4_crit = 0.16 * exp( 0.1*temp_in - 3.5*rh_in - 27.7 )
! ch2so4 = (ug-h2so4/m3-air)
! ch2so4*1.0e-12/mwh2so4 = (mole-h2so4/cm3-air)
qh2so4_crit = (ch2so4_crit*1.0e-12/98.0)/cair
qh2so4_avail = (qh2so4_cur - qh2so4_crit)*dtnuc
! if "available" h2so4 vapor < 4.0e-18 mole/mole-air ~= 1.0e2 molecules/cm3,
! exit with new particle formation = 0
if (qh2so4_avail .le. 4.0e-18) then
return
end if
! determine size bin into which the new particles go
isize_nuc = 1
vol_part = volumlo_sect(1)
!
! determine composition and density of the "grown particles"
! the grown particles are assumed to be liquid
! (since critical clusters contain water)
! so any (nh4/so4) molar ratio between 0 and 2 is allowed
! assume that the grown particles will have
! (nh4/so4 molar ratio) = min( 2, (nh3/h2so4 gas molar ratio) )
!
if (qnh3_cur .ge. qh2so4_avail) then
! combination of ammonium sulfate and ammonium bisulfate
! dum_n1 & dum_n2 = mole fractions of the ammsulf & ammbisulf
dum_n1 = (qnh3_cur/qh2so4_avail) - 1.0
dum_n1 = max( 0.0, min( 1.0, dum_n1 ) )
dum_n2 = 1.0 - dum_n1
dum_n3 = 0.0
else
! combination of ammonium bisulfate and sulfuric acid
! dum_n2 & dum_n3 = mole fractions of the ammbisulf & sulfacid
dum_n1 = 0.0
dum_n2 = (qnh3_cur/qh2so4_avail)
dum_n2 = max( 0.0, min( 1.0, dum_n2 ) )
dum_n3 = 1.0 - dum_n2
end if
dum_m1 = dum_n1*mw_ammsulf
dum_m2 = dum_n2*mw_ammbisulf
dum_m3 = dum_n3*mw_sulfacid
dens_part = (dum_m1 + dum_m2 + dum_m3)/ &
((dum_m1/dens_ammsulf) + (dum_m2/dens_ammbisulf) &
+ (dum_m3/dens_sulfacid))
! 25-jul-2006 - use 'in' value only if it is between 1.6-2.0 g/cm3
if (abs(dens_nh4so4a-1.8) .le. 0.2) then
dens_part = dens_nh4so4a
else
dens_nh4so4a = dens_part
end if
mass_part = vol_part*dens_part
molenh4a_per_moleso4a = 2.0*dum_n1 + dum_n2
! changes to h2so4 & nh3 gas (in mole/mole-air), limited by amounts available
duma = 0.9999
qh2so4_del = min( duma*qh2so4_cur, qh2so4_avail )
qnh3_del = min( duma*qnh3_cur, qh2so4_del*molenh4a_per_moleso4a )
qh2so4_del = -qh2so4_del
qnh3_del = -qnh3_del
! changes to so4 & nh4 aerosol (in mole/mole-air)
qso4a_del = -qh2so4_del
qnh4a_del = -qnh3_del
! change to aerosol number (in #/mole-air)
qnuma_del = (qso4a_del*mw_so4a + qnh4a_del*mw_nh4a)/mass_part
return
end subroutine wexler_nuc_mosaic_1box
!-----------------------------------------------------------------------
end module module_mosaic_newnuc