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

#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