ecearth3/sources/tm5mp/proj/cbm4/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:
  ! -----------
  ! 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