ecearth3/sources/tm5mp/proj/cb05/m7_nuck.F90

#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