#include "tm5.inc" SUBROUTINE m7_nuck(kproma, kbdim, klev, pso4g, pelvoc, & ptp1, prelhum, papp1, panew, pa4delt, & ptime, paernl, pm6rp ) ! ! Authors: ! -------- ! J. Wilson, E. Vignati, JRC/EI, original source 09/2000 ! P. Stier, MPI f90-version, ! changes, ! comments, ! modularisation and ! implementation of ! Vehkamaeki (2002) 2001-2003 ! ! Purpose: ! -------- ! This routine calls the routines for the computation of the ! nucleation rate znucrate [molec. cm-3 s-1] and, for ! Vehkamaeki (2002), of the number of molecules in the critical ! cluster 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. Whilst this is done ! analytically for the old Kulmala (1998) parameterization, it has ! to be done numerically for the new Vehkamaeki (2002) scheme. ! ! Interface: ! ---------- ! *m7_nuck* is called from *m7_dnum* ! ! Method: ! ------- ! ! 1) Vehkamaeki et al. (2002): ! ! An analytical integration of the nucleation equation is not ! possible, therefor the nucleation rate is simply multiplied ! by dt, implying a fixed gas-phase H2SO4 concentration over ! the timestep. !@@@ The dependency of the results on the used timestep !@@@ should be checked carefully in a sensitivity study! ! The number of nucleated particles is then calculated by ! taking the calculated critical mass critn of newly nucleated ! particles and dividing the total mass of nucleated sulfate ! by this critical mass of one nucleated particle. ! ! 2) Kulmala et al. (1998): ! ! For the Kulmala et al. (1998) scheme the 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. ! ! ! 3) Paasonen et al. (2010): ! A semi-empirical parameterization for new particle formation ! mechanism involving low-volatility organic vapors. We use the ! mechanism from Eq. 15, which presumes homogeneous heteromolecular ! nucleation between sulfuric acid and organic molecules, or that ! one of the vapors activate the clusters the clusters of other ! compound parameterised as J=K[H2SO4][Org]. Diameter of produced particle ! is assumed to be 2 nm. ! ! Externals: ! ---------- ! nucl_kulmala ! nucl_vehkamaeki ! nucl_bl ! ! 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 ! Vignatti, E. (1999), Modelling Interactions between Aerosols and ! Gaseous Compounds in the Polluted Marine Atmosphere. PhD-Thesis, ! RISO National Laborartory Copenhagen, Riso-R-1163(EN) use GO, only : gol, goErr, goPr, goBug USE mo_time_control, ONLY: delta_time USE mo_control, ONLY: nrow USE mo_kind, ONLY: wp USE mo_aero_m7, ONLY: critn, naermod, nnucl, avo, nmod, isoans ! USE mo_aero_mem, ONLY: d_nuc_so4 USE mo_aero_nucl, ONLY: nucl_kulmala, nucl_vehkamaeki, nucl_bl USE chem_param, ONLY: xmelvoc USE mo_aero, ONLY: nsoa IMPLICIT NONE !--- Parameters: ! ! pso4g = mass of gas phase sulfate [molec. cm-3] ! pelvoc = mass of gas phase elvoc [molec cm-3] ! ptp1 = atmospheric temperature at time t+1 [K] ! prelhum = atmospheric relative humidity [%] ! pa4delt(:,:,1) = mass of H2SO4 added to the nucleation mode due ! to nucleation of H2SO4 over ztmst. ! Equilvalent to the integral of H2SO4 gas loss ! due to nucleation over timestep ztmst. [molec. cm-3] ! panew = number of nucleated particles over timestep ztmst ! panew=pa4delt/critn i.e. mass of formed sulfate ! divided by an assumed mass of a nucleus. [cm-3] ! ptime = TM5 timestep ! !--- Local variables: ! ! zso4g_new = temporay storage of gas phase sulfate [molec. cm-3] ! zncrit = number of molecules in the critical cluster [1] ! ! See comments! INTEGER :: kproma, kbdim, klev REAL :: ptime, time_step_len REAL :: pso4g(kbdim,klev), ptp1(kbdim,klev), & prelhum(kbdim,klev), panew(kbdim,klev), papp1(kbdim,klev), pelvoc(kbdim,klev) REAL :: pa4delt(kbdim,klev,naermod) REAL :: paernl(kbdim,klev,nmod), & pm6rp(kbdim,klev,nmod) ! Local variables: INTEGER :: jk, jl, jrow REAL :: ztmst, zqtmst, zf1, zeps REAL :: znucrate(kbdim,klev), & ! nucleation rate [m-3 s-1] !@@@ Make consistent with Kulmala zso4g_new(kbdim,klev), & ! new gas phase sulfate concentration [molec. cm-3] zalpha(kbdim,klev), & ! auxiliary coefficient for the analytical integration of Kulmala zbeta(kbdim,klev), & ! auxiliary coefficient for the analytical integration of Kulmala zncrit(kbdim,klev), & ! number of molecules in the critical cluster [1] zdiam(kbdim,klev) ! diameter of the critical cluster [nm] REAL :: bl_ppfr(kbdim,klev), bl_pns_sa(kbdim,klev), bl_pns_org(kbdim,klev), pmolecELVOC(kbdim,klev) !--- 0) Initialisations: ------------------------------------------------------ jrow=nrow(2) time_step_len = ptime ztmst=time_step_len zqtmst=1/ztmst zeps=EPSILON(1.0_wp) znucrate=0. zncrit=0. bl_ppfr=0. ! Convert ELVOC from ug/m3 to molec/cm3 !1e-12/xmbc*Avog DO jk=1, klev DO jl=1, kproma pmolecELVOC(jl,jk)=(pelvoc(jl,jk)*1.e-12*avo)/(xmelvoc) END DO END DO IF(nnucl==1) THEN CALL nucl_vehkamaeki(kproma, kbdim, klev, & ! ECHAM5 dimensions ptp1, prelhum, pso4g, & ! ECHAM5 temperature, relative humidity znucrate, zncrit, zdiam ) !--- Calculate updated gas phase concentration: ! ! N(t) = N(0) - znucrate * zncrit * dt ! [cm-3] = [cm-3] - [cm-3 s-1] * [1] * [s] DO jk=1, klev DO jl=1, kproma zso4g_new(jl,jk)=pso4g(jl,jk)-(znucrate(jl,jk)*zncrit(jl,jk)*ztmst) END DO END DO ELSE IF (nnucl==2) THEN zncrit(:,:)=critn CALL nucl_kulmala(kproma, kbdim, klev, & pso4g, ptp1, prelhum, & znucrate, zalpha, zbeta ) !--- 2) Analytical integration of the nucleation rate (Eq. 19) ---------------- ! over timestep ztmst assuming no new H2SO4(g) production: ! ! d(N_av/critn)/dt=exp(alpha + beta*ln(N_av)) => ... => ! ! N_av(t0+dt)= ! [N_av(t0)**(1-beta) + critn*(1-beta)exp(alpha)*dt]**(1/(1-beta) ! DO jk=1, klev DO jl=1, kproma IF (znucrate(jl,jk) .GT. 1e-10) THEN zf1 = pso4g(jl,jk)**(1.0-zbeta(jl,jk))-critn*EXP(zalpha(jl,jk))*(1.0-zbeta(jl,jk))*ztmst zso4g_new(jl,jk) = EXP(LOG(zf1)/(1.0 - zbeta(jl,jk))) ELSE zso4g_new(jl,jk) = pso4g(jl,jk) END IF END DO END DO ELSE IF (nsoa==2 .and. nnucl==3) THEN CALL nucl_vehkamaeki(kproma, kbdim, klev, & ! ECHAM5 dimensions ptp1, prelhum, pso4g, & ! ECHAM5 temperature, relative humidity znucrate, zncrit, zdiam ) CALL nucl_bl(kproma, kbdim, klev, & ! Dimensions pso4g, pmolecELVOC, ptp1, papp1, prelhum, & ! H2SO4 concentration, ELVOC concentration paernl, pm6rp, znucrate, zdiam, & ! bl_ppfr, bl_pns_sa, bl_pns_org, ztmst) ! Formation rate, number of H2SO4 molecules in cluster ELSE IF (nsoa==2 .and. nnucl==4) THEN CALL nucl_vehkamaeki(kproma, kbdim, klev, & ! ECHAM5 dimensions ptp1, prelhum, pso4g, & ! ECHAM5 temperature, relative humidity znucrate, zncrit, zdiam ) CALL nucl_bl(kproma, kbdim, klev, & ! Dimensions pso4g, pmolecELVOC, ptp1, papp1, prelhum, & ! H2SO4 concentration, ELVOC concentration paernl, pm6rp, znucrate, zdiam, & ! bl_ppfr, bl_pns_sa, bl_pns_org, ztmst) ! Formation rate, number of H2SO4 molecules in cluster ELSE write(gol,*) 'M7_nuck: Choice of nucleation scheme is not supported, nsoa=',nsoa,' nnucl= ',nnucl call goPr return END IF !--- 3) Calculate changes in gas-phase and aerosol mass and aerosol numbers: -- DO jk=1, klev DO jl=1, kproma IF( znucrate(jl,jk) > zeps .OR. bl_ppfr(jl,jk) > zeps ) THEN !--- 3.1) If BL nucleation is active, apply only BL nucleation: IF(nsoa==2 .and. (nnucl==3 .OR. nnucl==4)) THEN pa4delt(jl,jk,1) = bl_ppfr(jl,jk)*bl_pns_sa(jl,jk)*ztmst pa4delt(jl,jk,isoans) = bl_ppfr(jl,jk)*bl_pns_org(jl,jk)*ztmst*xmelvoc/(1.e-12*avo) pa4delt(jl,jk,1) = MAX(MIN(pso4g(jl,jk), pa4delt(jl,jk,1)),0.0) pa4delt(jl,jk,isoans) = MAX(MIN(pelvoc(jl,jk),pa4delt(jl,jk,isoans)),0.0) panew(jl,jk)=bl_ppfr(jl,jk)*ztmst !RM !IF(jk==klev .AND. jl==230) THEN ! WRITE(*,*) 'nuck, delt:',pa4delt(jl,jk,1),pa4delt(jl,jk,isoans),bl_ppfr(jl,jk) ! WRITE(*,*) 'nuck, orig :',pso4g(jl,jk),pelvoc(jl,jk),pmolecELVOC(jl,jk) !END IF ! !--- 3.4) Calculate changes in gas phase due to nucleation: ! pso4g(jl,jk)=pso4g(jl,jk)-pa4delt(jl,jk,1) pelvoc(jl,jk)=pelvoc(jl,jk)-pa4delt(jl,jk,isoans) !--- 4) Vertically integrate mass of nucleated sulfate for diagnostics ! Convert [molec. cm-3] to [kg(S) m-2]: ELSE ! BL nucleation not active !--- 3.2) Calculate mass of nucleated H2SO4 (equals the net ! gas phase H2SO4 loss): ! ! pa4delt(jl,jk,1) = (pso4g(jl,jk)-zso4g_new(jl,jk)) ! JadB: A problem occurs. We calculate gas_new, wich is gas_old - nucleation, ! and then calculate the nucleation again with gas_old - gas_new. ! Huge errors when the time step is small. ! That's why we take the nucleation as fresh as possible (from znucrate) and apply the limiters. !pa4delt(jl,jk,1) = (znucrate(jl,jk)*zncrit(jl,jk)*ztmst) pa4delt(jl,jk,1) = (znucrate(jl,jk)*zncrit(jl,jk)*ztmst) pa4delt(jl,jk,1) = Max(pa4delt(jl,jk,1),0.0) pa4delt(jl,jk,1) = Min(pa4delt(jl,jk,1),pso4g(jl,jk)) ! !--- 3.3) Calculate the number of nucleated particles (nucleated mass ! divided by the assumed mass of a critical cluster critn): panew(jl,jk)=pa4delt(jl,jk,1)/zncrit(jl,jk) ! !--- 3.4) Calculate changes in gas phase H2SO4 due to nucleation: ! pso4g(jl,jk)=pso4g(jl,jk)-pa4delt(jl,jk,1) !--- 4) Vertically integrate mass of nucleated sulfate for diagnostics ! Convert [molec. cm-3] to [kg(S) m-2]: END IF END IF END DO END DO END SUBROUTINE m7_nuck