#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: ! ----------- ! Paasonen et al. (2010), On the roles of sulphuric acid and low-volatility ! organic vapours in the initial steps of atmospheric new particle ! formation, Atmos. Chem. Phys., 10, 11223-11242, 2010 ! 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 real, public ::d_form 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, pdiam ) ! 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), & pdiam(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 pdiam(jl,jk) = 2.0*zrc ELSE ! pmolecH2SO4, ptp1 , prhd out of range pntot(jl,jk) =0.0 pxtrnucr(jl,jk)=0.0 pdiam(jl,jk) =0.0 END IF END DO ! kproma END DO ! klev END SUBROUTINE nucl_vehkamaeki SUBROUTINE nucl_bl(kproma, kbdim, klev, & ! Dimensions pmolecH2SO4, pmolecELVOC, ptp1, papp1, prhd, & ! H2SO4 concentration, ELVOC concentration, temperature, pressure, RH paernl, pm6rp, binnucrate, bindiam, & ! ppfr, pns_sa, pns_org, ztmst) ! Number concs, radii, formation rate, number of H2SO4 molecules in cluster ! ! Authors: ! --------- ! R. MAKKONEN, UNIV. HELSINKI 2016-2017 ! ! Purpose ! --------- ! Calculation of 2 nm formation rate based on semi-empirical equation. ! ! References: ! ----------- ! Paasonen et al. (2010) ACP ! Kerminen and Kulmala (2002) Aer. Sci. ! USE mo_control, ONLY: nrow USE chem_param, ONLY: xmh2so4, xmelvoc USE mo_aero_m7, ONLY: pi, dh2so4,nnucl, avo, condensation_sink, nmod, doc !---------------------------------------------------- IMPLICIT NONE !---------------------------------------------------- INTEGER :: kproma, kbdim, klev INTEGER :: jk, jl,jrow !---------------------------------------------------- ! REAL :: pns_org(kbdim,klev), & pns_sa(kbdim,klev), & ppfr(kbdim,klev), & binnucrate(kbdim,klev),& bindiam(kbdim,klev) REAL :: ptp1(kbdim,klev), prhd(kbdim,klev), papp1(kbdim,klev), & pxtrnucr(kbdim,klev), & pmolecH2SO4(kbdim,klev), & ! Sulfuric acid concentration pntot(kbdim,klev), & pmolecELVOC(kbdim,klev) ! ELVOC concentration REAL :: paernl(kbdim,klev,nmod), & pm6rp(kbdim,klev,nmod) REAL :: zaernl(nmod),zm6rp(nmod) REAL :: ztmst !---------------------------------------------------- ! Local Arrays REAL :: gr_sa, gr_org, nmolec, d_nuc, lambda, factor, cs_prime, fracELVOCnucl, org_ss, binfactor REAL :: vtot, nuclrate, formrate, n_sa, n_org, delvoc, vfr_sa, vfr_org, diffcoef_org, cs(nmod) !--- 0) Initializations: jrow=nrow(2) ppfr =0. pns_sa =1. pns_org =1. d_nuc=2. ! Size at parameterized nucleation rate (nm) ! d_form= ! Size at simulated formation rate (nm) delvoc=doc ! Density of ELVOC g/cm3 diffcoef_org=1.5E-5 ! Diffusion coefficient for calculation of CS from CS', unit m2/s !--- Total volume of d_form -sized spherical particle (cm3) vtot = (4./3.) * pi * (0.5*d_form*1.E-7)**3 DO jk=1, klev DO jl=1,kproma ! Mean free path (m) lambda=0.066 *(1.01325E+5/papp1(jl,jk))*(ptp1(jl,jk)/293.15)*1.E-06 ! Condensation sink prime (m-2) !zaernl(:)=paernl(jl,jk,:) !zm6rp(:)=pm6rp(jl,jk,:) !CALL condensation_sink_prime(zaernl(:), 0.01*zm6rp(:), lambda, cs) CALL condensation_sink_prime(paernl(jl,jk,:), 0.01*pm6rp(jl,jk,:), lambda, cs) cs_prime=MAX(MIN(sum(cs),1.E10),1.e-20) ! Steady-state concentration of ELVOCs assuming above condensation sink org_ss=pmolecELVOC(jl,jk)/(4*pi*diffcoef_org*ztmst*cs_prime) org_ss=MAX(MIN(org_ss,pmolecELVOC(jl,jk)),0.0) ! Semi-empirical nucleation rate (d_nuc sized clusters). Only used for rate, not for composition information IF(nnucl==3) THEN nuclrate = 11.0 * 1.E-14 * pmolecH2SO4(jl,jk) * org_ss ! Paasonen et al. (2010) Eq. 15 d_nuc=2.0 ELSE IF(nnucl==4) THEN nuclrate = 3.27 * 1.E-21 * (pmolecH2SO4(jl,jk)**2) * org_ss ! Riccobono et al. (2014) d_nuc=1.7 END IF ! Calculate "Kerminen-Kulmala factor" for binary nucleation from critical cluster size to d_nuc ! Safety check for division by zero and floating overflow (when binnucrate==0, can bindiam>d_nuc) if (binnucrate(jl,jk)>1e-20) then CALL kkfactor(ptp1(jl,jk),pmolecH2SO4(jl,jk),org_ss,pm6rp(jl,jk,:),paernl(jl,jk,:),bindiam(jl,jk),d_nuc,cs_prime,gr_sa,gr_org,binfactor) else binfactor=0.0 end if ! Calculate "Kerminen-Kulmala factor", a proportionality factor to account for particles losses between d_nuc and d_form CALL kkfactor(ptp1(jl,jk),pmolecH2SO4(jl,jk),org_ss,pm6rp(jl,jk,:),paernl(jl,jk,:),d_nuc,d_form,cs_prime,gr_sa,gr_org,factor) ! Formation rate of d_form clusters: binary nucleation (at d_nuc) plus semi-empirical nucleation (at d_nuc) formrate = (binnucrate(jl,jk)*binfactor + nuclrate)*factor ! security check IF (formrate< 1e-20) then cycle end IF ! Assume that volume ratio v_org/v_sa in forming particle is proportional to GR_org^3 / GR_sa^3 IF(gr_sa+gr_org .GT. 1.E-15) THEN vfr_sa = gr_sa**3 / (gr_sa+gr_org)**3 vfr_org = gr_org**3 / (gr_sa+gr_org)**3 vfr_sa = MAX(MIN(vfr_sa / (vfr_sa+vfr_org), 1.0), 0.0) vfr_org = 1 - vfr_sa ELSE vfr_sa=0.5 vfr_org=0.5 END IF ! Number of sulphuric acid and organic molecules in d_form sized particle ! # = 1 cm3 mol-1 g cm-3 g-1 mol n_sa = vfr_sa * vtot * avo * dh2so4 / xmh2so4 n_org = vfr_org * vtot * avo * delvoc / xmelvoc ! Sanity checks and limiters ! Limiting formation rate based on both sulphuric acid and ELVOC concentration. ! Limit 1a: sulfuric acid deficit IF( n_sa * formrate * ztmst .GT. pmolecH2SO4(jl,jk)) THEN ! Try to maintain formation rate: calculate maximum VFR using all available sulfuric acid vfr_sa = pmolecH2SO4(jl,jk) * xmh2so4/ ( formrate * vtot * avo * dh2so4 * ztmst) vfr_org = 1 - vfr_sa ! Update number of molecules in d_form sized particle n_sa = vfr_sa * vtot * avo * dh2so4 / xmh2so4 ! (==[H2SO4]/(dT*J) ) n_org = vfr_org * vtot * avo * delvoc / xmelvoc ! Limit 1b: J can not be maintained by organics IF( n_org * formrate * ztmst .GT. org_ss) THEN ! Calculate J from available H2SO4 and ELVOC !cm-3 s-1= mol cm-3 s-1 ( cm-3 g mol-1 g-1 cm3 ) formrate = 1.0 / (avo * vtot * ztmst) * (pmolecH2SO4(jl,jk) * xmh2so4 / dh2so4 + org_ss * xmelvoc / delvoc) n_sa = pmolecH2SO4(jl,jk) / (formrate * ztmst) n_org = org_ss / (formrate * ztmst) !NOT-NEEDED-EXCEPT-4-DEBUG vfr_sa = vtot * avo * dh2so4 / (xmh2so4 * n_sa) ! (==[H2SO4]/(dT*J) ) !NOT-NEEDED-EXCEPT-4-DEBUG vfr_org = vtot * avo * delvoc / (xmelvoc * n_org) END IF ! Limit 1b ! Limit 2a: ELVOC deficit ELSE IF( n_org * formrate * ztmst .GT. org_ss) THEN ! Try to maintain formation rate: calculate maximum VFR using all available ELVOC ! 1 = cm-3 g mol-1 cm3 s cm-3 mol g-1 cm3 s-1 vfr_org = MAX(MIN(org_ss * xmelvoc/ ( formrate * vtot * avo * delvoc * ztmst), 1.0), 0.0) vfr_sa = 1 - vfr_org ! Update number of molecules in d_form sized particle n_sa = vfr_sa * vtot * avo * dh2so4 / xmh2so4 n_org = vfr_org * vtot * avo * delvoc / xmelvoc ! (==[ELVOC(SS)]/(dT*J) ) ! Limit 2b: J can not be maintained by sulfuric acid IF( n_sa * formrate * ztmst .GT. pmolecH2SO4(jl,jk) ) THEN ! Calculate J from available H2SO4 and ELVOC formrate = 1.0 / (avo * vtot * ztmst) * (pmolecH2SO4(jl,jk) * xmh2so4 / dh2so4 + org_ss * xmelvoc / delvoc) n_sa = pmolecH2SO4(jl,jk) / (formrate * ztmst) n_org = org_ss / (formrate * ztmst) !NOT-NEEDED-EXCEPT-4-DEBUG vfr_sa = vtot * avo * dh2so4 / (xmh2so4 * n_sa) ! (==[H2SO4]/(dT*J) ) !NOT-NEEDED-EXCEPT-4-DEBUG vfr_org = vtot * avo * delvoc / (xmelvoc * n_org) END IF ! Limit 2b END IF ! Limit 2a pns_sa(jl,jk) = n_sa ! * formrate ! kerrotaan jo m7_nuck? pns_org(jl,jk) = n_org ! * formrate ! kerrotaan jo m7_nuck? ppfr(jl,jk) = formrate !DBG IF(jk==klev .AND. jl==230) THEN !DBG! IF(3.GT.2) THEN ! .AND. cs_prime .GT. 1.E-4 .AND. ptp1(jl,jk) .GT. 293.0 .AND. nuclrate .GT. 1) THEN !DBG WRITE(*,*) 'blnuc: gr_org, n_org, pmolecELVOC(jl,jk),org_ss', gr_org, n_org, pmolecELVOC(jl,jk),org_ss !DBG WRITE(*,*) 'blnuc: gr_sa, n_sa, pmolecH2SO4(jl,jk)', gr_sa, n_sa, pmolecH2SO4(jl,jk) !DBG WRITE(*,*) 'blnuc: nuclrate,formrate,factor,cs',nuclrate,formrate,factor,cs_prime !DBG WRITE(*,*) 'blnuc: vfr ', vfr_sa, vfr_org !DBG WRITE(*,*) 'blnuc: vtot ', vtot, ztmst !DBG END IF END DO END DO END SUBROUTINE nucl_bl SUBROUTINE kkfactor(ptp1, ph2so4, porg, pm6rp, paernl, d_nuc, d_form, cs_prime, gr_sa, gr_org, factor) USE mo_aero_m7, ONLY: nmod, sigma, pi, dh2so4 USE chem_param, ONLY: xmh2so4, xmelvoc REAL,INTENT(IN) :: d_nuc, d_form, paernl(nmod), pm6rp(nmod), ptp1, ph2so4, porg, cs_prime REAL,INTENT(OUT) :: gr_sa, gr_org, factor REAL :: dorg, lambda, cs(nmod), gamma, vmol ! Gamma-parameter for K&K CALL kk_gamma_parameter(paernl(:), 0.01*pm6rp(:), d_nuc, d_form, dh2so4, ptp1, gamma) !--- Growth rates, calculated using Eq. 21 in Kerminen&Kulmala 2002, assuming 1 g/cm3 density for nuclei !--- Growth rate due to H2SO4 vmol=SQRT(3.0 * 8.314472 * ptp1 / (xmh2so4*0.001) ) ! RMS Molecular speed (m/s) gr_sa=((3.E-9)/(1000.))*(vmol*xmh2so4*ph2so4) ! Growth rate (nm/h) gr_sa=MAX(MIN(gr_sa,1.E2),0.) ! Limit GR to <100 nm/h ! Growth rate from ELVOCs vmol=SQRT(3.0 * 8.314472 * ptp1 / (xmelvoc*0.001) ) ! RMS Molecular speed (m/s) gr_org=((3.E-9)/(1000.))*(vmol*xmelvoc*porg) ! Growth rate (nm/h) gr_org=MAX(MIN(gr_org,1.E2),0.) ! Limit GR to <100 nm/h ! Safety check for division by zero IF((gr_sa+gr_org) .GT. 1.E-10)THEN factor=EXP((1./d_form-1./d_nuc)*(gamma*cs_prime/(gr_sa+gr_org))) else factor=0.0 end if END SUBROUTINE kkfactor SUBROUTINE kk_gamma_parameter(paernl, pm6rp, d_nuc, d_form, rho, t, gamma) USE mo_aero_m7, ONLY: nmod, sigma, wh2so4, rerg, dh2so4 REAL :: gamma0 REAL :: gamma REAL,intent(in) :: d_nuc, d_form REAL :: d_mean, rho, t REAL :: paernl(nmod), pm6rp(nmod) gamma0=0.23 !! [gamma]=(nm2*m2)/h ! Number mean diameter for pre-existing particle population (now mass mean diameter) ! [d_mean] = nm d_mean=(1.E9)*(paernl(2)*pm6rp(2)+paernl(3)*pm6rp(3)+paernl(4)*pm6rp(4)+ & paernl(5)*pm6rp(5)+paernl(6)*pm6rp(6)+paernl(7)*pm6rp(7))/6. gamma=gamma0*((d_nuc/1.)**0.2)*((d_form/3.)**0.075)*((d_mean/150.)**0.048)*(rho**(-0.33))*(t/293.) END SUBROUTINE kk_gamma_parameter SUBROUTINE condensation_sink_prime(N,D,LAMBDA,OUTPUT) ! Calculate 'condensation sink' using eq. 6 in Kerminen et. al 2004 ! eq. 3 in Kerminen et. al 2002 USE mo_aero_m7, ONLY: nmod REAL, INTENT(OUT) :: OUTPUT(nmod) !output coefficients REAL, INTENT(IN) :: N(nmod) !number concentrations [#/cm3] REAL, INTENT(IN) :: D(nmod) !current geom median diameters of the modes [m] REAL, INTENT(IN) :: LAMBDA ![m] REAL :: KN ! Knudsen number REAL :: CS(nmod) INTEGER :: jmod DO jmod=1,nmod IF (N(jmod)>1.e0) THEN KN=(2*LAMBDA)/D(jmod) CS(jmod)=((D(jmod))*(N(jmod)*1.E6)*(1+KN))/(1+0.377*KN+1.33*KN*(1+KN)) CS(jmod)=MAX(CS(jmod),0.) CS(jmod)=MIN(CS(jmod),1.E15) ELSE CS(jmod)=0. END IF END DO !Note: the 0.5 factor is applied already here OUTPUT=0.5*CS END SUBROUTINE condensation_sink_prime END MODULE mo_aero_nucl