123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448 |
- #include "tm5.inc"
- MODULE mo_aero_nucl
- ! *mo_aero_nucl* holds routines calculating the nucleation rate
- !
- ! Author:
- ! -------
- ! Philip Stier, MPI-Met, Hamburg, 01/2003
- !
- ! Purpose:
- ! --------
- ! This module holds the routines for the calculation of the nucleation
- ! rate for the binary nucleation in the sulfate water system.
- ! These routines are called and applied within HAM from the routine
- ! m7_nuck.
- ! The old parameterization of Kulmala (1998), that apparently contains
- ! some errors is kept for consistency. It is recommended to use the
- ! new scheme of Vehkamaeki et al. (2002). However, this parameterization
- ! is no longer analytically integrable, the effect of this is to be
- ! scrutinized.
- ! The modularized version of the Kulmala parameterization has been tested
- ! to give bit-identical results with the old version.
- !
- ! References:
- ! -----------
- ! Vehkamaeki et al. (2002), An improved parameterization for sulfuric
- ! acid/water nucleation rates for tropospheric and stratospheric
- ! conditions, J. Geophys. Res, 107, D22, 4622
- ! Kulmala et al. (1998), Parameterizations for sulfuric acid/water
- ! nucleation rates. J. Geophys. Res., 103, No D7, 8301-8307
- ! Vignati, E. (1999), Modelling Interactions between Aerosols and
- ! Gaseous Compounds in the Polluted Marine Atmosphere. PhD-Thesis,
- ! RISO National Laborartory Copenhagen, Riso-R-1163(EN)
- IMPLICIT NONE
- CONTAINS
- SUBROUTINE nucl_kulmala(kproma, kbdim, klev, & ! ECHAM5 dims
- pso4g, ptp1, prelhum, &
- pbnrate, palpha, pbeta )
-
- ! Authors:
- ! --------
- ! P. Stier, MPI-Met, Hamburg, from the original f77 code
- ! in the routine m7_nuck 2001-2003
- ! J. Wilson, E. Vignati, JRC/EI, original source 09/2000
- !
- ! Purpose:
- ! --------
- ! This routine calculates the instananeous nucleation rate
- ! znucrate [molec. cm-3 s-1] from a given gas phase H2SO4 concentration
- ! pso4g [molec. cm-3]. It also calculates the integrated change of
- ! H2SO4 gas phase mass over one timestep due to nucleation
- ! pa4delt(:,:,1) [molec. cm-3] as well as the number of nucleated
- ! particles panew [1] during the timestep.
- !
- ! Interface:
- ! ----------
- ! *nucl_kulmala* is called from *m7_nuck*
- !
- ! Method:
- ! -------
- ! Kulmala et al. (1998) 's formula for binary nucleation is
- ! rewritten to take the form znucrate = exp[zalpha+ln(pso4g)*beta].
- ! Equation numbers are taken from Kulmala et al. (1998).
- ! After the calculation of the nucleation rate znucrate, it is
- ! integrated in 2) analytically over one timestep, i.e.:
- !
- ! Integration of:
- !
- ! znucrate=d(critn*znav)/dt=exp[zalpha + zbeta*ln(znav)]
- !
- ! gives znav(t0+dt, znav(t0)) where znav(t0) is pso4g.
- ! znav is temporarily stored in zso4g_new and confined between
- ! 0 and pso4g.
- ! The number of nucleated particles is then calculated by
- ! assuming a fixed critical mass critn of newly nucleated
- ! particles and dividing the total mass of nucleated sulfate
- ! by this critical mass of one nucleated particle.
- !
- ! Externals:
- ! ----------
- ! None
- USE mo_aero_m7, ONLY: bk
- IMPLICIT NONE
- INTEGER :: kproma, kbdim, klev
- REAL :: pso4g(kbdim,klev), ptp1(kbdim,klev), &
- prelhum(kbdim,klev), pbnrate(kbdim,klev), &
- palpha(kbdim,klev), pbeta(kbdim,klev)
- INTEGER :: jl, jk
- REAL :: znwv, zln_nac, ztk, zsupsat, &
- zpeh2o, zpeh2so4, zra, zxal, &
- ztkn, zssn, zdelta
- !---1) Calculation of the nucleation rate: ----------------------------
-
- DO jk=1,klev
- DO jl=1,kproma
- IF (pso4g(jl,jk) .GT. 1e-5) THEN
- ztk=ptp1(jl,jk)
- zsupsat=prelhum(jl,jk)
- !
- !--- 1.1) Restrict t, and rh to limits where the parameterization ok:
- !
- ztkn = MAX(ztk, 220.0)
- zssn = MIN(zsupsat, 0.90)
- !
- !--- 1.2) Equlibrium vapour pressures (Jaeker-Mirabel (1995), JGR):
- !
- !--- H2O equlibrium vapour pressure (Tabata):
- !
- zpeh2o=0.750064*(10.**(8.42926609-1827.17843/ &
- ztkn-71208.271/ztkn/ztkn))*1333./bk/ztkn
- !
- !--- H2SO4 equlibrium vapour pressure at 360
- !@@@ Check source: which Ayers?
- !
- zpeh2so4=EXP(-10156./ztkn+16.259)*7.6e2*1333./bk/ztkn
- !
- !--- H2SO4 equlibrium vapour pressure - correction of ayers
- ! by kulmala - currently not used
- !
- ! payers=exp(-10156/360+16.259)*7.6e2
- ! zpeh2so4=exp(log(payers)+10156*(-1./ztkn+1./360.+0.38/(905-360) * &
- ! (1+log(360./ztkn)-360./ztkn)))*1333/bk/ztkn
- !
- !--- 1.3) Relative acidity (0.0 -1.0):
- !
- zra=pso4g(jl,jk)/zpeh2so4
- !
- !--- 1.4) Water vapour molecule concentration [cm-3]:
- !
- znwv=zsupsat*zpeh2o
- !
- !--- 1.5) Factor delta in Eq. 22:
- !
- zdelta=1.0+(ztkn-273.15)/273.15
- !
- !--- 1.6) Molefraction of H2SO4 in the critical cluster
- ! minus the H2SO4(g) term in Eq. 17:
- !
- zxal = 1.2233-0.0154*zra/(zra+zssn)-0.0415*LOG(znwv)+ 0.0016*ztkn
- !
- !--- 1.7) Exponent of the critical cluster (Eq. 18):
- !
- zln_nac = -14.5125+0.1335*ztkn-10.5462*zssn+1958.4*zssn/ztkn
- !
- !--- 1.8) Sum of all terms in Eq. 20 containing H2SO4(g):
- !
- pbeta(jl,jk) = 25.1289 - 4890.8/ztkn + 7643.4*0.0102/ztkn - &
- 2.2479*zdelta*zssn - 1.9712*0.0102*zdelta/zssn
- !
- !--- 1.9) Sum all terms in Eq. 20 not containing H2SO4(g):
- !
- palpha(jl,jk) = zln_nac*(-25.1289 + 4890.8/ztkn + 2.2479*zdelta*zssn) - &
- 1743.3/ztkn + zxal*(7643.4/ztkn - 1.9712*zdelta/zssn)
- !
- !--- 1.10) Nucleation rate [cm-3 s-1] (Kulmala et al., 1998):
- !
- pbnrate(jl,jk) = EXP(palpha(jl,jk)+LOG(pso4g(jl,jk))*pbeta(jl,jk))
- ELSE
- palpha(jl,jk) =0.
- pbeta(jl,jk) =0.
- pbnrate(jl,jk)=0.
- END IF ! pso4g(jl,jk) .GT. 1e-5
- END DO ! kproma
- END DO !klev
- END SUBROUTINE nucl_kulmala
- SUBROUTINE nucl_vehkamaeki(kproma, kbdim, klev, & ! ECHAM5 dimensions
- ptp1, prhd, pmolecH2SO4, & ! ECHAM5 temperature, relative humidity
- pxtrnucr, pntot ) ! nucleation rate, number of molecules in the
- ! critical cluster
-
- !
- ! Authors:
- ! ---------
- ! C. TIMMRECK, MPI HAMBURG 2002
- !
- ! Purpose
- ! ---------
- ! Calculation of classical nucleation rate
- !
- ! calculation of the nucleation rate after Vehkamaeki et al. (2002)
- ! The calculation of the nucrate ZKNH2SO4 is in cm^-3 s^-1
- ! and a coarse approxmation for the first class
- !
- ! Modifications:
- ! --------------
- ! R. Hommel; rewrite in f90, adopted to ECHAM5; MPI HAMBURG; Dec. 2002
- ! P. Stier; bugfixes, modularisation and optimization; MPI HAMBURG; 2003-2004
- !
- ! H2SO4 still fixed to xxx molc/cm3, no sulfur cycle coupling yet
- !
- ! References:
- ! -----------
- ! Vehkamaeki et al. (2002), An improved parameterization for sulfuric
- ! acid/water nucleation rates for tropospheric and stratospheric
- ! conditions, J. Geophys. Res, 107, D22, 4622
- !
- ! Parameters
- ! ----------
- ! prho = prhop_neu in *sam*
- !
- ! prhd = relative humidity in sam_aeroprop & sam_nucl
- !
- ! pxtrnucr = nucleation rate in [1/m3s]
- ! xrhoc = density of the critical nucleus in kg/m^3
- ! zrxc = ?
- !
-
- USE mo_control, ONLY: nrow
- !----------------------------------------------------
- IMPLICIT NONE
- !----------------------------------------------------
-
- INTEGER :: kproma, kbdim, klev
-
- INTEGER :: jk, jl,jrow
-
- !----------------------------------------------------
- !
-
- REAL :: ptp1(kbdim,klev), prhd(kbdim,klev), &
- pxtrnucr(kbdim,klev), &
- pmolecH2SO4(kbdim,klev), & ! revisited, ok
- pntot(kbdim,klev)
-
- !----------------------------------------------------
- ! Local Arrays
-
- REAL :: zrxc(kbdim)
- REAL :: zrhoa, zrh, zt, x, zjnuc, zrc, &
- zxmole, zntot
- !--- 0) Initializations:
- jrow=nrow(2)
- DO jk=1, klev
- DO jl=1,kproma
- !----1.) Parameterization of nucleation rate after Vehkamaeki et al. (2002)
- ! t: temperature in K (190.15-300.15K)
- ! zrh: saturatio ratio of water (0.0001-1)
- ! zrhoa: sulfuric acid concentration in 1/cm3 (10^4-10^11 1/cm3)
- ! jnuc: nucleation rate in 1/cm3s (10^-7-10^10 1/cm3s)
- ! ntot: total number of molecules in the critical cluster (ntot>4)
- ! x: molefraction of H2SO4 in the critical cluster
- ! rc: radius of the critical cluster in nm
- ! Calculate nucleation only for valid thermodynamic conditions:
- IF( (pmolecH2SO4(jl,jk)>=1.E+4) .AND. &
- (prhd(jl,jk) >=1.E-4) .AND. &
- (ptp1(jl,jk)>=190.15 .AND. ptp1(jl,jk)<=300.15) ) THEN
- zrhoa=MIN(pmolecH2SO4(jl,jk),1.e11)
- zrh=MIN(prhd(jl,jk),1.0)
- zt=ptp1(jl,jk)
- ! Equation (11) - molefraction of H2SO4 in the critical cluster
- x=0.7409967177282139 - 0.002663785665140117*zt &
- + 0.002010478847383187*LOG(zrh) &
- - 0.0001832894131464668*zt*LOG(zrh) &
- + 0.001574072538464286*LOG(zrh)**2 &
- - 0.00001790589121766952*zt*LOG(zrh)**2 &
- + 0.0001844027436573778*LOG(zrh)**3 &
- - 1.503452308794887e-6*zt*LOG(zrh)**3 &
- - 0.003499978417957668*LOG(zrhoa) &
- + 0.0000504021689382576*zt*LOG(zrhoa)
- zxmole=x
- ! Equation (12) - nucleation rate in 1/cm3s
-
- zjnuc=0.1430901615568665 + 2.219563673425199*zt - &
- 0.02739106114964264*zt**2 + &
- 0.00007228107239317088*zt**3 + 5.91822263375044/x + &
- 0.1174886643003278*LOG(zrh) + 0.4625315047693772*zt*LOG(zrh) - &
- 0.01180591129059253*zt**2*LOG(zrh) + &
- 0.0000404196487152575*zt**3*LOG(zrh) + &
- (15.79628615047088*LOG(zrh))/x - &
- 0.215553951893509*LOG(zrh)**2 - &
- 0.0810269192332194*zt*LOG(zrh)**2 + &
- 0.001435808434184642*zt**2*LOG(zrh)**2 - &
- 4.775796947178588e-6*zt**3*LOG(zrh)**2 - &
- (2.912974063702185*LOG(zrh)**2)/x - &
- 3.588557942822751*LOG(zrh)**3 + &
- 0.04950795302831703*zt*LOG(zrh)**3 - &
- 0.0002138195118737068*zt**2*LOG(zrh)**3 + &
- 3.108005107949533e-7*zt**3*LOG(zrh)**3 - &
- (0.02933332747098296*LOG(zrh)**3)/x + &
- 1.145983818561277*LOG(zrhoa) - &
- 0.6007956227856778*zt*LOG(zrhoa) + &
- 0.00864244733283759*zt**2*LOG(zrhoa) - &
- 0.00002289467254710888*zt**3*LOG(zrhoa) - &
- (8.44984513869014*LOG(zrhoa))/x + &
- 2.158548369286559*LOG(zrh)*LOG(zrhoa) + &
- 0.0808121412840917*zt*LOG(zrh)*LOG(zrhoa) - &
- 0.0004073815255395214*zt**2*LOG(zrh)*LOG(zrhoa) - &
- 4.019572560156515e-7*zt**3*LOG(zrh)*LOG(zrhoa) + &
- (0.7213255852557236*LOG(zrh)*LOG(zrhoa))/x + &
- 1.62409850488771*LOG(zrh)**2*LOG(zrhoa) - &
- 0.01601062035325362*zt*LOG(zrh)**2*LOG(zrhoa) + &
- 0.00003771238979714162*zt**2*LOG(zrh)**2*LOG(zrhoa) + &
- 3.217942606371182e-8*zt**3*LOG(zrh)**2*LOG(zrhoa) - &
- (0.01132550810022116*LOG(zrh)**2*LOG(zrhoa))/x + &
- 9.71681713056504*LOG(zrhoa)**2 - &
- 0.1150478558347306*zt*LOG(zrhoa)**2 + &
- 0.0001570982486038294*zt**2*LOG(zrhoa)**2 + &
- 4.009144680125015e-7*zt**3*LOG(zrhoa)**2 + &
- (0.7118597859976135*LOG(zrhoa)**2)/x - &
- 1.056105824379897*LOG(zrh)*LOG(zrhoa)**2 + &
- 0.00903377584628419*zt*LOG(zrh)*LOG(zrhoa)**2 - &
- 0.00001984167387090606*zt**2*LOG(zrh)*LOG(zrhoa)**2 + &
- 2.460478196482179e-8*zt**3*LOG(zrh)*LOG(zrhoa)**2 - &
- (0.05790872906645181*LOG(zrh)*LOG(zrhoa)**2)/x - &
- 0.1487119673397459*LOG(zrhoa)**3 + &
- 0.002835082097822667*zt*LOG(zrhoa)**3 - &
- 9.24618825471694e-6*zt**2*LOG(zrhoa)**3 + &
- 5.004267665960894e-9*zt**3*LOG(zrhoa)**3 - &
- (0.01270805101481648*LOG(zrhoa)**3)/x
-
- zjnuc=EXP(zjnuc) ! add. Eq. (12) [1/(cm^3s)]
- ! Equation (13) - total number of molecules in the critical cluster
- zntot=-0.002954125078716302 - 0.0976834264241286*zt + &
- 0.001024847927067835*zt**2 - 2.186459697726116e-6*zt**3 - &
- 0.1017165718716887/x - 0.002050640345231486*LOG(zrh) - &
- 0.007585041382707174*zt*LOG(zrh) + &
- 0.0001926539658089536*zt**2*LOG(zrh) - &
- 6.70429719683894e-7*zt**3*LOG(zrh) - &
- (0.2557744774673163*LOG(zrh))/x + &
- 0.003223076552477191*LOG(zrh)**2 + &
- 0.000852636632240633*zt*LOG(zrh)**2 - &
- 0.00001547571354871789*zt**2*LOG(zrh)**2 + &
- 5.666608424980593e-8*zt**3*LOG(zrh)**2 + &
- (0.03384437400744206*LOG(zrh)**2)/x + &
- 0.04743226764572505*LOG(zrh)**3 - &
- 0.0006251042204583412*zt*LOG(zrh)**3 + &
- 2.650663328519478e-6*zt**2*LOG(zrh)**3 - &
- 3.674710848763778e-9*zt**3*LOG(zrh)**3 - &
- (0.0002672510825259393*LOG(zrh)**3)/x - &
- 0.01252108546759328*LOG(zrhoa) + &
- 0.005806550506277202*zt*LOG(zrhoa) - &
- 0.0001016735312443444*zt**2*LOG(zrhoa) + &
- 2.881946187214505e-7*zt**3*LOG(zrhoa) + &
- (0.0942243379396279*LOG(zrhoa))/x - &
- 0.0385459592773097*LOG(zrh)*LOG(zrhoa) - &
- 0.0006723156277391984*zt*LOG(zrh)*LOG(zrhoa) + &
- 2.602884877659698e-6*zt**2*LOG(zrh)*LOG(zrhoa) + &
- 1.194163699688297e-8*zt**3*LOG(zrh)*LOG(zrhoa) - &
- (0.00851515345806281*LOG(zrh)*LOG(zrhoa))/x - &
- 0.01837488495738111*LOG(zrh)**2*LOG(zrhoa) + &
- 0.0001720723574407498*zt*LOG(zrh)**2*LOG(zrhoa) - &
- 3.717657974086814e-7*zt**2*LOG(zrh)**2*LOG(zrhoa) - &
- 5.148746022615196e-10*zt**3*LOG(zrh)**2*LOG(zrhoa) + &
- (0.0002686602132926594*LOG(zrh)**2*LOG(zrhoa))/x - &
- 0.06199739728812199*LOG(zrhoa)**2 + &
- 0.000906958053583576*zt*LOG(zrhoa)**2 - &
- 9.11727926129757e-7*zt**2*LOG(zrhoa)**2 - &
- 5.367963396508457e-9*zt**3*LOG(zrhoa)**2 - &
- (0.007742343393937707*LOG(zrhoa)**2)/x + &
- 0.0121827103101659*LOG(zrh)*LOG(zrhoa)**2 - &
- 0.0001066499571188091*zt*LOG(zrh)*LOG(zrhoa)**2 + &
- 2.534598655067518e-7*zt**2*LOG(zrh)*LOG(zrhoa)**2 - &
- 3.635186504599571e-10*zt**3*LOG(zrh)*LOG(zrhoa)**2 + &
- (0.0006100650851863252*LOG(zrh)*LOG(zrhoa)**2)/x + &
- 0.0003201836700403512*LOG(zrhoa)**3 - &
- 0.0000174761713262546*zt*LOG(zrhoa)**3 + &
- 6.065037668052182e-8*zt**2*LOG(zrhoa)**3 - &
- 1.421771723004557e-11*zt**3*LOG(zrhoa)**3 + &
- (0.0001357509859501723*LOG(zrhoa)**3)/x
-
- zntot=EXP(zntot) ! add. Eq. (13)
-
- ! Equation (14) - radius of the critical cluster in nm
- zrc=EXP(-1.6524245+0.42316402*x+0.33466487*LOG(zntot)) ! [nm]
- ! Conversion [nm -> m]
- zrxc(jl)=zrc*1e-9
-
- !----1.2) Limiter
- IF(zjnuc<1.e-7 .OR. zntot<4.0) zjnuc=0.0
- ! limitation to 1E+10 [1/cm3s]
-
- zjnuc=MIN(zjnuc,1.e10)
- pxtrnucr(jl,jk) = zjnuc
- ! convert total number of molecules in the critical cluster
- ! to number of sulfate molecules:
- pntot(jl,jk)=zntot*zxmole
- ELSE ! pmolecH2SO4, ptp1 , prhd out of range
- pntot(jl,jk) =0.0
- pxtrnucr(jl,jk)=0.0
- END IF
- END DO ! kproma
- END DO ! klev
- END SUBROUTINE nucl_vehkamaeki
- END MODULE mo_aero_nucl
|