ecearth3/sources/tm5mp/proj/cbm4/eqsam_param.F90

!
#include "tm5.inc"
!
module eqsam_param

  implicit none

  private

  public :: eqsam_v03d

contains


  subroutine eqsam_v03d(yi,yo,nca,nco,iopt,loop,imax,ipunit,in)
    !
    implicit none
    !___________________________________________________________________________________________________________________________________
    !      Written by Swen Metzger 3/11/99. Modified 2002, 2003.
    !
    !      Department of Atmospheric Chemistry, Max-Planck-Institute for Chemistry.
    !      email: metzger@mpch-mainz.mpg.de
    !      http://www.mpch-mainz.mpg.de/~metzger
    !
    !      COPYRIGHT 1999-2003
    !
    !      purpose
    !      -------
    !      EQSAM (EQuilibrium Simplified Aerosol Model) is a new and computationally efficient thermodynamic
    !      aerosol composition model that allows to calculate the gas/aerosol equilibrium partitioning,
    !      including aerosol water, sufficiently fast and accurate for global (or even regional) modeling.
    !      EQSAM is based on a number of parameterizations, including single solute molalities and activity 
    !      coefficients (AC). The thermodynamic framework (domains and subdomains, internally mixed aerosols) 
    !      is the same as of more sophisticated thermodynamic equilibrium models (EQMs), e.g. of ISORROPIA 
    !      (Nenes et al., 1998). Details are given in the references below (and the references therein).
    !
    !      The main assumption on which EQSAM/EQMs are based is thermodynamical and chemical equilibrium. 
    !      From this assumption it directly follows that the aerosol water activity (aw) equals the ambient 
    !      relative humidity (RH), if the water vapor pressure is sufficiently larger than the partial vapor
    !      pressure of the aerosol compounds. This is approximately true for tropospheric aerosols. Given the 
    !      large amount of water vapor present, water vapor and aerosol water equilibrate relatively faster 
    !      compared to all other aerosol compounds. This is subsequently also true for single aerosol compounds.
    !      The water activity of single solutes must also equal RH under this assumption. Therefore, the so 
    !      called ZSR-relation is (and can be) used to calculate the aerosol associated water mass (simply
    !      from the sum of all water mass fractions that are derived from measured single solute molalities). 
    !
    !      In contrast to other EQMs, EQSAM utilizes the fact that the RH fixes the water activity 
    !      (under the above assumptions) and the consequence that any changes in RH also causes changes in 
    !      the aerosol water mass and, hence, aerosol activity (including activity coefficients). Thus, an decrease
    !      (increase) in RH decrease (increases) the aerosol water mass (and water activity). This can change the
    !      aerosol composition, e.g. due to condensation (evaporation/crystallization), because the vapor pressure 
    !      above the aerosol reduces (increases). In turn, a vapor pressure reduction (increase) due to changes
    !      in the aerosol composition is compensated by an associated condensation (evaporation) of water vapor 
    !      to maintain the aerosol molality to remain constant (because aw=RH). Furthermore, the aerosol water 
    !      mainly depends on the aerosol mass and the type of solute, so that parameterizations of single solute 
    !      molalities and activity coefficients can be defined, only depending on the type of solute and RH.
    !      The advantage of using such parameterizations is that the entire aerosol equilibrium composition 
    !      can be solved analytically, i.e. non-iteratively, which considerably reduces the amount of CPU time 
    !      that is usually need for aerosol thermodynamic calculations (especially if an EQM is incorporated in
    !      an aerosol dynamical model that is in turn embedded in a high resolution regional or global model).
    !
    !      However, EQSAM should still be regarded as a starting point for further developments. There is still 
    !      room for improvements. For instance, this code is not yet numerically optimized (vectorized) and a 
    !      number of improvements with respect to an explicit treatment of additional equilibrium reactions,
    !      missing (or only implicit) dissociation, and a basic parameterization of the water uptake. 
    !      
    !      Note that EQSAM was originally developed to calculate the gas/aerosol equilibrium partitioning of the 
    !      ammonium-sulfate-nitrate-water system for climate models, excluding solid compounds. 
    !      This version (eqsam_v03d.f90) is extended with respect to sea salt. Solids/hysteresis are treated in a 
    !      simplified manner. Results of a box model comparison with ISORROPIA will be available from the web page.
    !      Please also note that the water uptake is based on additional (unpublished) parameterizations for single 
    !      solute molalities, which are derived from tabulated measurements used in ISORROPIA. Note further that 
    !      this extended version (eqsam_v03d.f90) is not yet published. A publication is in progress.
    !
    ! ToDo:
    !     Split ion-pairs into ions for water parameterizations (since info is actually available)
    !     Include uptake/dissociation of NH3, HNO3, HCl (mainly to get pH right at near neutral conditions)
    !     Extension to K+,Ca++,Mg++, CO2/(CO3)2--/HCO3-,SOA,etc.. (maybe not)
    !     Vectorization. Translation of hardcoded formulas in array syntax.
    !     I/O Interface and program structure clean up.
    !     EQSAM info webpage.
    !
    ! Version History:
    !
    !  eqsam_v03d.f90 (MPI-CH, June 2003): 
    !   - gama parameterizations now according to Metzger 2002 (JGR Appendix)
    !   - improved pH calculations (still restricted to strong acids)
    !   - removed bug that lead to too high nitrate formation at dry and cold regions (UT/LS) 
    !   - removed bug in solid/hysteresis calculations 
    !     (both bugs introduced in eqsam_v03b.f90 by cleaning up eqsam_v02a.f90)
    !   
    !  eqsam_v03c.f90 (MPI-CH, April 2003):
    !   - more accurate paramterizations of single solute molalities (Na, Cl species)
    !   - cleanded up RHD subdomain structure
    !   - improved water uptake (Na, Cl species)
    !
    !  eqsam_v03b.f90 (MPI-CH, March 2003): 
    !                 System extended to HCl,Cl-/Na+. 
    !                 Parameterization (fit) of additional HNO3 uptake removed.
    !                 Instead, complete analytical solution of equilibrium reactions, based on the AC-RH relationship.
    !  eqsam_v03.f90  (IMAU, October 1999): 
    !                 Test version (included in TM3).
    !  eqsam_v02a.f90 (IMAU, April 2000):
    !                 Box model version.
    !  eqsam_v02.f90  (IMAU, October 1999):
    !                 TM3 version.
    !                 Version including solids and additional HNO3 uptake on acidic aerosols (parameterized).
    !  eqsam_v01b.f90 (MPI-CH, January 2003):
    !                 Same as eqsam_v01a.f90 (additional lines though uncommented for test purposes only).
    !  eqsam_v01a.f90 (IMAU, April 2000):
    !                 Box model version.
    !  eqsam_v01.f90  (IMAU, October 1999):
    !                 TM3 version.
    !                 First and most basic version (without solids) for better vectorization (for global modeling).
    !                 System: NH3,NH4+/H2SO4+,HSO4-,SO4--/HNO3,NO3-, H2O 
    !                 based on equilibrium / internal mixture assumption / aw=rh / ZSR-relation
    !                 parameterization of activcity coefficients (AC), i.e. an AC-RH relationship
    !
    !      
    !      interface
    !      ---------
    !      call  eqsam_v03d(yi,yo,nca,nco,iopt,loop,imax,ipunit,in)
    !
    !      yi = input  array (imax, nca)
    !      yo = output array (imax, nco)
    !      imax = max loop (e.g. time steps)
    !      nca >= 11
    !      nco >= 35
    !      iopt = 1 metastable 
    !      iopt = 2 solids 
    !      iopt = 3 hysteresis (metastable/solids) for online calculations
    !      iopt = 31 hysteresis lower branch 
    !      iopt = 32 hysteresis upper branch 
    !      ipunit = I/O unit (can be skipped)
    !      in = array        (can be skipped)
    !         
    !      method
    !      ------
    !      equilibrium / internal mixture assumption / aw=rh
    !      System: NH3,NH4+/H2SO4+,HSO4-,SO4--/HNO3,NO3-, HCl,Cl-/Na+, H2O 
    !              (K+,Ca++,Mg++)
    !      external
    !      --------
    !      program    eqmd.f90    (driver only needed for the box model version)
    !      subroutine gribio.f90  (provides diagnostics output in grib/binary/ascii format)
    !      
    !      references
    !      ---------
    !      Swen Metzger Ph.D Thesis, University Utrecht, 2000.
    !         http://www.library.uu.nl/digiarchief/dip/diss/1930853/inhoud.htm
    !
    !      Metzger, S. M., F. J. Dentener, J. Lelieveld, and S. N. Pandis, 
    !         GAS/AEROSOL PARTITIONING I: A COMPUTATIONALLY EFFICIENT MODEL, 
    !         J Geophys. Res., 107, D16, 10.1029/2001JD001102, 2002
    !         http://www.agu.org/journals/jd/jd0216/2001JD001102/index.html
    !      Metzger, S. M., F. J. Dentener, A. Jeuken, and M. Krol, J. Lelieveld, 
    !         GAS/AEROSOL PARTITIONING II: GLOBAL MODELING RESULTS, 
    !         J Geophys. Res., 107, D16, 10.1029/2001JD001103, 2002.
    !         http://www.agu.org/journals/jd/jd0216/2001JD001103/index.html
    !___________________________________________________________________________________________________________________________________
    real,parameter                                :: RH_HIST_DW=1.50                                    ! mean value for mixture of wet (2) and dry (1) gridboxes (needed for HYSTERESIS)
    real,parameter                                :: T0=298.15,T1=298.0,AVO=6.03e23,R=82.0567e-6,     & ! in cu.m*atm/deg/mole
         r_kcal  = 1.986E-3                                 ! Ideal gas constant [kcal K-1.mole-1]
    real,parameter                                :: RHMAX=0.99,RHMIN=0.0001                            ! restrict to max / min RH
    real,parameter                                :: MWNH4=18.,MWSO4=96.,MWNO3=62.,MWCl=35.5            ! mole mass of species considered
    real,parameter                                :: MWNa=23.0,MWCa=40.1,MWN=14.0, MWS=32.1
    real,parameter                                :: MWH20=55.51*18.01,ZERO=0.0
    real,parameter                                :: GF1=0.25,GF2=0.50,GF3=0.40,GF4=1.00,K=2.           ! exponents of AC-RH functions
    !______________________________________________
    integer,parameter                             :: NPAIR=10 
    !
    integer                                       :: ii,il,IHYST
    integer,intent(in)                            :: nca,nco,imax,loop,ipunit
    integer,intent(inout)                         :: iopt
    !______________________________________________
    integer,dimension(6),intent(in)               :: in
    !______________________________________________
    real                                          :: T0T,TT,RH,PX,RHD,KAN,KAC,ZIONIC,RH_HIST,GAMA,GG,GF,GFN
    real                                          :: X00,X01,X02,X03,X04,X05,X08,X09,X10,X11
    real                                          :: X0,X1,X2,X3,X4,X5,X6,XK10,XK6
    real                                          :: ZFLAG,ZKAN,ZKAC,PH,COEF,HPLUS,AKW,XKW,MOLAL
    real                                          :: TNH4,TSO4,TNO3,TNa,TCl,TPo,TCa,TMg
    real                                          :: PNH4,PSO4,PNO3,PCl,PNa,GNO3,GNH3,GSO4,GHCl
    real                                          :: ASO4,ANO3,ANH4,ACl,ANa,SNH4,SSO4,SNO3,SCl,SNa
    real                                          :: WH2O,PM,PMs,PMt,RINC,DON,RATIONS,GR,NO3P,NH4P
    Real                                          :: H2O_NO3
    !_______________________________________________
    real,dimension(imax,nca),intent(in)           :: yi
    real,dimension(imax,nco),intent(out)          :: yo
    real,dimension(8)                             :: w1,w2
    real,dimension(8)                             :: RHDA,RHDE,RHDX,RHDZ    ! RHD / MRHD arrays for different aerosol types
    real,dimension(NPAIR)                         :: M0,MW,NW,ZW            ! arrays of ion pairs
    !
    ! salt solutes:
    !   1 = NACl,  2 = (NA)2SO4, 3 = NANO3,  4 = (NH4)2SO4,  5 = NH4NO3, 6 = NH4CL,   7 = 2H-SO4
    !   8 = NH4HSO4,   9 = NAHSO4, 10 = (NH4)3H(SO4)2
    !
    DATA MW(1:NPAIR)/ 58.5, 142.0, 88.0, 132.0, 80.0, 53.5, 98.0, 115.0, 120.0, 247.0/ ! mole mass of the salt solute
    DATA NW(1:NPAIR)/  2.0,   2.5,  2.5,   2.5,  3.5,  1.0,  4.5,   2.0,   2.0,   2.5/ ! square of max. dissocation number (not consistent)
    DATA ZW(1:NPAIR)/ 0.67,   1.0,  1.0,   1.0,  1.0,  1.0,  0.5,   1.0,   1.0,   1.0/ ! exponents of water activity functions
    !
    DATA RHDA(1:8)/0.32840, 0.4906, 0.6183, 0.7997, 0.67500, 0.5000, 0.4000, 0.0000/ ! RHD / MRHD values as of ISORROPIA / SCAPE (T=298.15K)
    DATA RHDE(1:8)/-1860.0, -431.0, 852.00, 80.000, 262.000, 3951.0, 384.00, 0.0000/ ! Temp. coeff.
    !___________________________________________________________________________________________________________________________________
    IHYST=2
    IF(IOPT.EQ.31) THEN      ! SOLID HYSTORY
       IHYST=1
       IOPT=3
    ELSEIF(IOPT.EQ.32) THEN  ! WET   HISTORY
       IHYST=2
       IOPT=3
    ENDIF
    !kt write(ipunit,*)'eqsam_v03d ...'
    !kt print*,'                                                          '
    !kt print*,'              EQuilibrium Simplified Aerosol Model (EQSAM)'
    !kt print*,'                        for global modeling               ' 
    !kt print*,'                                 by                       '
    !kt print*,'                         Swen Metzger, MPI-CH             '
    !kt print*,'                         Copyright 1999-2003              '
    !kt print*,'                    >> metzger@mpch-mainz.mpg.de <<       '
    !kt print*,'                     last change:  04. June, 2003         '
    !kt print*,'                           (version 3.0d)                 '
    !kt print*,'                   gas/aerosol calculations assuming      '
    !kt print*,'                  System: NH3,NH4+/H2SO4+,HSO4-,SO4--     '
    !kt print*,'                      HNO3,NO3-, HCl,Cl-/Na+, H2O         '
    !kt if(iopt.eq.1) then
    !kt print*,'                         metastable aeorsols              '
    !kt elseif(iopt.eq.2) then
    !kt print*,'                            solid aeorsols                '
    !kt elseif(iopt.eq.3) then
    !kt print*,'                             hysteresis                   '
    !kt print*,'                         (metastable/solids)              '
    !kt if(IHYST.eq.1) then
    !kt print*,'                            solid hystory                 '
    !kt elseif(IHYST.eq.2) then
    !kt print*,'                             wet hystory                  '
    !kt endif
    !kt endif
    !kt print*,'                                                          '
    !kt print*,'loop over ',loop,' data sets'
    !kt print*,'   '
    !___________________________________________________________________________________________________________________________________
    yo=0.;w1=0.;w2=0.                                                                                   ! init/reset
    !___________________________________________________________________________________________________________________________________
    do il=1,loop

       ! get old relative humidity to calculate aerosol hysteresis (online only)

       RH_HIST = 2.                                                                                     ! WET HISTORY (DEFAULT)
       IF(IHYST.EQ.1.OR.IOPT.EQ.2)  RH_HIST = 1.                                                        ! SET TO SOLIDS

       !  meteorology
       TT = yi(il,1)                    ! T                      [K]
       RH = yi(il,2)                    ! RH                     [0-1]
       PX = yi(il,11)                   ! p                      [hPa]
       !
       ! gas+aerosol:
       w1(1) = yi(il,6)                 ! Na+ (ss  + xsod) (a)   [umol/m^3 air]
       w1(2) = yi(il,4)                 ! H2SO4    + SO4-- (p)   [umol/m^3 air]
       w1(3) = yi(il,3)                 ! NH3  (g) + NH4+  (p)   [umol/m^3 air]
       w1(4) = yi(il,5)                 ! HNO3 (g) + NO3-  (p)   [umol/m^3 air]
       w1(5) = yi(il,7)                 ! HCl  (g) + Cl-   (p)   [umol/m^3 air]
       w1(6) = yi(il, 8)                ! K+   (p) from Dust     [umol/m^3 air]
       w1(7) = yi(il, 9)                ! Ca++ (p) from Dust     [umol/m^3 air]
       w1(8) = yi(il,10)                ! Mg++ (p) from Dust     [umol/m^3 air]
       !______________________________________________

       zflag=1.

       w1=w1*1.0e-6                     ! [mol/m^3 air]

       TNa   = w1(1)                    ! total input sodium   (g+p) 
       TSO4  = w1(2)                    ! total input sulfate  (g+p) 
       TNH4  = w1(3)                    ! total input ammonium (g+p)
       TNO3  = w1(4)                    ! total input nitrate  (g+p) 
       TCl   = w1(5)                    ! total input chloride (g+p) 
       TPo   = w1(6)                    ! total input potasium (g+p) 
       TCa   = w1(7)                    ! total input calcium  (g+p)
       TMg   = w1(8)                    ! total input magnesium(g+p)

       ! SULFATE RICH

       if((w1(1)+w1(3)+w1(6)+2.*(w1(7)+w1(8))).le.(2.*w1(2))) then
          zflag=3.
       endif

       ! SULFATE VERY RICH CASE if (NH4+Na+K+2(Ca+Mg))/SO4 < 1

       if((w1(1)+w1(3)+w1(6)+2.*(w1(7)+w1(8))).le.w1(2)) then
          zflag=4.
       endif

       ! SULFATE NEUTRAL CASE

       if((w1(1)+w1(3)+w1(6)+2.*(w1(7)+w1(8))).gt.(2.*w1(2))) then
          zflag=2.
       endif

       ! SULFATE POOR AND CATION POOR CASE

       if((w1(1)+w1(6)+2.*(w1(7)+w1(8))).gt.(2.*w1(2))) then       
          zflag=1.
       endif

       IF ( RH .LT. RHMIN ) RH=RHMIN
       IF ( RH .GT. RHMAX ) RH=RHMAX

       ! CALCULATE TEMPERATURE DEPENDENCY FOR SOME RHDs

       RHDX(:)=RHDA(:)*exp(RHDE(:)*(1./TT-1./T0))
       RHDZ(:)=RHDX(:)

       ! ACCOUNT FOR VARIOUS AMMOMIUM/SODIUM SULFATE SALTS ACCORDING TO MEAN VALUE AS OF ISORROPIA

       GG=2.0                          ! (Na)2SO4 / (NH4)2SO4 IS THE PREFFERED SPECIES FOR SULFATE DEFICIENT CASES
       IF(ZFLAG.EQ.3.) THEN
          IF(RH.LE.RHDZ(7)) THEN       ! ACCOUNT FOR MIXTURE OF (NH4)2SO4(s) & NH4HSO4(s) & (NH4)3H(SO4)2(s) 
             GG=1.677                  !                        (Na)2SO4 &  NaHSO4
             !           GG=1.5
          ELSEIF(RH.GT.RHDZ(7).AND.RH.LE.RHDZ(5)) THEN ! MAINLY (Na)2SO4 / (NH4)2SO4(s) & (NH4)3H(SO4)2(s)
             GG=1.75
             !           GG=1.5
          ELSEIF(RH.GE.RHDZ(5)) THEN   ! (NH4)2SO4(S) & NH4HSO4(S) & SO4-- & HSO4-
             GG=1.5                    !  (Na)2SO4 &  NaHSO4
          ENDIF
       ENDIF
       IF(ZFLAG.EQ.4.) GG=1.0          ! IF SO4 NEUTRALIZED, THEN ONLY AS NaHSO4 / NH4HSO4(S) OR  HSO4- / H2SO4

       RHD=RH
       IF(IOPT.EQ.2.OR.RH_HIST.LT.RH_HIST_DW) THEN   ! GET RHD FOR SOLIDS / HYSTERESIS
          !
          ! GET LOWEST DELIQUESCENCE RELATIVE HUMIDITIES ACCORDING TO THE CONCENTRATION DOMAIN (APROXIMATION) 
          ! BASED ON RHD / MRHD ISORROPIA/SCAPE
          !
          w2(:)=1.
          do ii=1,8
             if(w1(ii).le.1.e-12) w2(ii)=0.  ! skip compound in RHD calculation if value is concentration is zero or rather small
          enddo

          ! GET LOWEST RHD ACCORDING TO THE CONCENTRATION DOMAIN

          ! zflag=1. (cation rich)  ...
          ! 1. sea salt      aerosol          : RHDX(1)=MgCl2
          ! 2. mineral dust  aerosol          : RHDX(2)=Ca(NO3)2
          !
          ! zflag=2. (sulfate neutral) ...
          ! 3. ammonium + nitrate             : RHDX(3)= NH4NO3
          ! 4. ammonium + sulfate             : RHDX(4)=(NH4)2SO4        
          ! 5. ammonium + sulfate mixed salt  : RHDX(5)=(NH4)3H(SO4)2, (NH4)2SO4        
          ! 6. ammonium + nitrate  + sulfate  : RHDX(6)=(NH4)2SO4, NH4NO3, NA2SO4, NH4CL
          !
          ! zflag=3. (sulfate poor) ...
          ! 7. ammonium + sulfate  (1:1,1.5)  : RHDX(7)= NH4HSO4
          !
          ! zflag=4. (sulfate very poor) ...
          ! 8. sulfuric acid                  : RHDX(8)= H2SO4       

          IF(ZFLAG.EQ.1.)THEN

             RHD=W2(1)+W2(5)                     ! Na+  dependency
             IF(RHD.EQ.0.)  RHDX(1)=1. 
             RHD=W2(6)+W2(7)+W2(8)               ! K+/Ca++/Mg++ dependency (incl. ss)
             IF(RHD.EQ.0.)  RHDX(2)=1. 

             RHD=MINVAL(RHDX(1:2))

          ELSEIF(ZFLAG.EQ.2.)THEN

             RHD=W2(3)*W2(4)                     ! NH4+ & NO3-  dependency
             IF(RHD.EQ.0.)  RHDX(3)=1. 
             RHD=W2(2)+W2(3)                     ! NH4+ & SO4-- dependency
             IF(GG.NE.2.)   RHD=0.               ! account only for pure (NH4)2SO4
             IF(RHD.EQ.0.)  RHDX(4)=1. 
             RHD=W2(2)+W2(3)                     ! NH4+ & SO4-- dependency
             IF(RHD.EQ.0.)  RHDX(5)=1. 
             RHD=W2(2)+W2(3)+W2(4)+W2(5)         ! (NH4)2SO4, NH4NO3, NA2SO4, NH4CL dependency
             IF(RHD.EQ.0.)  RHDX(6)=1. 

             !     RHD=MINVAL(RHDX(3:4))
             RHD=MINVAL(RHDX(3:6))

          ELSEIF(ZFLAG.EQ.3.)THEN

             RHD=W2(2)+W2(3)                     ! NH4+ & SO4-- dependency
             IF(RHD.EQ.0.)  RHDX(7)=1. 
             RHD=RHDX(7)

          ELSEIF(ZFLAG.EQ.4.)THEN

             RHD=W2(2)                           ! H2SO4 dependency (assume no dry aerosol)
             IF(RHD.EQ.0.)  RHDX(8)=1. 

             RHD=RHDX(8)

          ENDIF ! ZFLAG
       ENDIF ! SOLIDS

       ! GET WATER ACTIVITIES ACCORDING TO METZGER, 2000.
       ! FUNCTION DERIVED FROM ZSR RELATIONSHIP DATA (AS USED IN ISORROPIA)

       M0(:) = ((NW(:)*MWH20/MW(:)*(1./RH-1.)))**ZW(:)

       ! CALCULATE TEMPERATURE DEPENDENT EQUILIBRIUM CONSTANTS

       T0T=T0/TT
       COEF=1.0+LOG(T0T)-T0T

       ! EQUILIBRIUM CONSTANT NH4NO3(s) <==> NH3(g) + HNO3(g) [atm^2] (ISORROPIA)

       XK10 = 5.746e-17
       XK10= XK10 * EXP(-74.38*(T0T-1.0) + 6.120*COEF)
       KAN = XK10/(R*TT)/(R*TT)

       ! EQUILIBRIUM CONSTANT  NH4CL(s) <==> NH3(g) + HCL(g) [atm^2] (ISORROPIA)

       XK6  = 1.086e-16
       XK6 = XK6 * EXP(-71.00*(T0T-1.0) + 2.400*COEF)
       KAC = XK6/(R*TT)/(R*TT)

       !
       ! CALCULATE AUTODISSOCIATION CONSTANT (KW) FOR WATER H2O <==> H(aq) + OH(aq) [mol^2/kg^2] (ISORROPIA)

       XKW  = 1.010e-14
       XKW = XKW *EXP(-22.52*(T0T-1.0) + 26.920*COEF)

       ! GET MEAN MOLAL IONIC ACTIVITY COEFF ACCORDING TO METZGER, 2002.

       GAMA=0.0
       IF(RH.GE.RHD)                GAMA=(RH**ZFLAG/(1000./ZFLAG*(1.-RH)+ZFLAG))
       GAMA = GAMA**GF1            ! ONLY GAMA TYPE OF NH4NO3, NaCl, etc. NEEDED SO FAR

       GAMA=0.0
       GFN=K*K                      ! K=2, i.e. condensation of 2 water molecules per 1 mole ion pair
       GF=GFN*GF1                   ! = GFN[=Nw=4] * GF1[=(1*1^1+1*1^1)/2/Nw=1/4] = 1
       ! ONLY GAMA TYPE OF NH4NO3, NH4Cl, etc. needed so far

       IF(RH.GE.RHD) GAMA=RH**GF/((GFN*MWH20*(1./RH-1.)))**GF1

       GAMA = min(GAMA,1.0)        ! FOCUS ON 0-1 SCALE
       GAMA = max(GAMA,0.0)
       GAMA = (1.-GAMA)**K          ! transplate into aqueous phase equillibrium and account for 
       ! enhanced uptake of aerosol precursor gases with increasing RH
       ! (to match the results of ISORROPIA)

       ! CALCULATE RHD DEPENDENT EQ: IF RH <  RHD => NH4NO3(s) <==> NH3 (g) + HNO3(g) (ISORROPIA)
       !                             IF RH >> RHD => HNO3  (g)   -> NO3 (aq)

       X00  = MAX(ZERO,MIN(TNa,GG*TSO4))       ! MAX SODIUM   SULFATE
       X0   = MAX(ZERO,MIN(TNH4,GG*TSO4-X00))  ! MAX AMMOMIUM SULFATE
       X01  = MAX(ZERO,MIN(TNa-X00, TNO3))     ! MAX SODIUM   NITRATE
       X1   = MAX(ZERO,MIN(TNH4-X0,TNO3-X01))  ! MAX AMMOMIUM NITRATE
       !
       X02  = MAX(ZERO,MIN(TNa-X01-X00,TCl))   ! MAX SODIUM   CHLORIDE
       X03  = MAX(ZERO,MIN(TNH4-X0-X1,TCl-X02))! MAX AMMOMIUM CHLORIDE

       X2   = MAX(TNH4-X1-X0-X03,ZERO)         ! INTERIM RESIDUAL NH3
       X3   = MAX(TNO3-X1-X01,ZERO)            ! INTERIM RESIDUAL HNO3
       X04  = MAX(TSO4-(X0+X00)/GG,ZERO)       ! INTERIM RESIDUAL H2SO4
       X05  = MAX(TCl-X03-X02,ZERO)            ! INTERIM RESIDUAL HCl
       !     X06  = MAX(TNa-X02-X01-X00,ZERO)        ! INTERIM RESIDUAL Na (should be zero for electro-neutrality in input data)
       !
       ZKAN=2.
       IF(RH.GE.RHD) ZKAN=ZKAN*GAMA

       X4   = X2 + X3
       X5   = SQRT(X4*X4+KAN*ZKAN*ZKAN)
       X6   = 0.5*(-X4+X5)
       X6   = MIN(X1,X6)

       GHCl = X05                              ! INTERIM RESIDUAl HCl
       GNH3 = X2 + X6                          ! INTERIM RESIDUAl NH3
       GNO3 = X3 + X6                          ! RESIDUAl HNO3
       GSO4 = X04                              ! RESIDUAl H2SO4
       PNa  = X02 + X01 + X00                  ! RESIDUAl Na (neutralized)

       ZKAC=2.
       IF(RH.GE.RHD) ZKAC=ZKAC*GAMA

       X08   = GNH3 + GHCl
       X09   = SQRT(X08*X08+KAC*ZKAC*ZKAC)
       X10   = 0.5*(-X08+X09)
       X11   = MIN(X03,X10)

       GHCl = GHCl + X11                       ! RESIDUAL HCl
       GNH3 = GNH3 + X11                       ! RESIDUAL NH3

       ! GO SAVE ...

       IF(GHCl.LT.0.)   GHCl=0.
       IF(GSO4.LT.0.)   GSO4=0.
       IF(GNH3.LT.0.)   GNH3=0.
       IF(GNO3.LT.0.)   GNO3=0.
       IF(PNa.LT.0.)    PNa=0.
       IF(GSO4.GT.TSO4) GSO4=TSO4
       IF(GNH3.GT.TNH4) GNH3=TNH4
       IF(GNO3.GT.TNO3) GNO3=TNO3
       IF(GHCl.GT.TCl)  GHCl=TCl
       IF(PNa.GT.TNa)   PNa=TNa
       !     IF(PNa.LT.TNa)   print*,il,' PNa.LT.TNa => no electro-neutrality in input data! ',PNa,TNa


       ! DEFINE AQUEOUSE PHASE (NO SOLID NH4NO3 IF NO3/SO4>1, TEN BRINK, ET AL., 1996, ATMOS ENV, 24, 4251-4261)

       !     IF(TSO4.EQ.ZERO.AND.TNO3.GT.ZERO.OR.TNO3/TSO4.GE.1.) RHD=RH

       !     IF(IOPT.EQ.2.AND.RH.LT.RHD.OR.IOPT.EQ.2.AND.RH_HIST.LT.RH_HIST_DW) THEN        ! SOLIDS / HYSTERESIS
       IF(RH_HIST.EQ.1.AND.RH.LT.RHD) THEN        ! SOLIDS / HYSTERESIS

          ! EVERYTHING DRY, ONLY H2SO4 (GSO4) REMAINS IN THE AQUEOUSE PHASE

          ANH4 = 0.
          ASO4 = 0.
          ANO3 = 0.
          ACl  = 0.
          ANa  = 0.

       ELSE  !  SUPERSATURATED SOLUTIONS NO SOLID FORMATION

          ASO4 = TSO4 - GSO4
          ANH4 = TNH4 - GNH3
          ANO3 = TNO3 - GNO3
          ACl  = TCl  - GHCl
          ANa  = PNa

       ENDIF ! SOLIDS / HYSTERESIS

       ! CALCULATE AEROSOL WATER [kg/m^3(air)]
       !
       ! salt solutes:
       !   1 = NACl,  2 = (NA)2SO4, 3 = NANO3,  4 = (NH4)2SO4,  5 = NH4NO3, 6 = NH4CL,   7 = 2H-SO4
       !   8 = NH4HSO4,   9 = NAHSO4, 10 = (NH4)3H(SO4)2
       !
       IF(ZFLAG.EQ.1.) WH2O = ASO4/M0( 2) + ANO3/M0(3) + ACl/M0(6)
       IF(ZFLAG.EQ.2.) WH2O = ASO4/M0( 9) + ANO3/M0(5) + ACl/M0(6)
       IF(ZFLAG.EQ.3.) WH2O = ASO4/M0( 8) + ANO3/M0(5) + ACl/M0(6)
       IF(ZFLAG.EQ.4.) WH2O = ASO4/M0( 8) + GSO4/M0(7)

       ! H2O_NO3
       IF(ZFLAG.EQ.1.) H2O_NO3 = ANO3/M0(3)
       IF(ZFLAG.EQ.2.) H2O_NO3 = ANO3/M0(5)
       IF(ZFLAG.EQ.3.) H2O_NO3 = ANO3/M0(5)
       IF(ZFLAG.EQ.4.) H2O_NO3 = 0.0

       ! CALCULATE AQUEOUS PHASE PROPERTIES

       !     PH    = 9999.
       PH    = 7.
       MOLAL = 0.
       HPLUS = 0.
       ZIONIC= 0.

       IF(WH2O.GT.0.) THEN            

          ! CALCULATE AUTODISSOCIATION CONSTANT (KW) FOR WATER

          AKW=XKW*RH*WH2O*WH2O                         ! H2O <==> H+ + OH- with kw [mol^2/kg^2]
          AKW=AKW**0.5                                 ! [OH-] = [H+] [mol]

          ! trap zero [OH-] due to round-off:
          if ( AKW > 0.0 ) then

            ! Calculate hydrogen molality [mol/kg], i.e. H+ of the ions:
            !                                   Na+, NH4+, NO3-, Cl-, SO4--, HH-SO4- [mol/kg(water)]
            !                                   with [OH-] = kw/[H+]

            HPLUS = (-ANa/WH2O-ANH4/WH2O+ANO3/WH2O+ACl/WH2O+GG*ASO4/WH2O+GG*GSO4/WH2O+ & 
                 SQRT(( ANa/WH2O+ANH4/WH2O-ANO3/WH2O-ACl/WH2O-GG*ASO4/WH2O-GG*GSO4/WH2O)**2+XKW/AKW*WH2O))/2.

            ! trip zero [H+] due to round-off:
            if ( HPLUS > 0.0 ) then

              ! Calculate pH
              PH = -ALOG10(HPLUS)                             ! aerosol pH 

              ! Calculate ionic strength [mol/kg]

              ZIONIC=0.5*(ANa+ANH4+ANO3+ACl+ASO4*GG*GG+GSO4*GG*GG+XKW/AKW*WH2O*WH2O)
              ZIONIC=ZIONIC/WH2O                                           ! ionic strength [mol/kg]
              !        ZIONIC=min(ZIONIC,200.0)                                     ! limit for output
              !        ZIONIC=max(ZIONIC,0.0)

            end if

         end if  ! [OH-] > 0.0

       ENDIF ! AQUEOUS PHASE 
       !
       !-------------------------------------------------------
       ! calculate diagnostic output consistent with other EQMs ...
       !
       ASO4 = ASO4 + GSO4                                     ! assuming H2SO4 remains aqueous

       TNa   = TNa  * 1.e6                                    ! total input sodium   (g+p)  [umol/m^3]
       TSO4  = TSO4 * 1.e6                                    ! total input sulfate  (g+p)  [umol/m^3]
       TNH4  = TNH4 * 1.e6                                    ! total input ammonium (g+p)  [umol/m^3]
       TNO3  = TNO3 * 1.e6                                    ! total input nitrate  (g+p)  [umol/m^3]
       TCl   = TCl  * 1.e6                                    ! total input chloride (g+p)  [umol/m^3]
       TPo   = TPo  * 1.e6                                    ! total input potasium (g+p)  [umol/m^3]
       TCa   = TCa  * 1.e6                                    ! total input calcium  (g+p)  [umol/m^3]
       TMg   = TMg  * 1.e6                                    ! total input magnesium(g+p)  [umol/m^3]
       !
       ! residual gas:
       GNH3 = GNH3 * 1.e6                                     ! residual NH3
       GSO4 = GSO4 * 1.e6                                     ! residual H2SO4
       GNO3 = GNO3 * 1.e6                                     ! residual HNO3
       GHCl = GHCl * 1.e6                                     ! residual HCl

       ! total particulate matter (neutralized)
       PNH4=TNH4-GNH3                                         ! particulate ammonium   [umol/m^3]
       !kt      PNO3=TNO3-GNO3                                         ! particulate nitrate    [umol/m^3]
       PNO3=max(0.,TNO3-GNO3)                                         ! particulate nitrate    [umol/m^3]
       PCl =TCl -GHCl                                         ! particulate chloride   [umol/m^3]
       PNa =TNa                                               ! particulate sodium     [umol/m^3]
       PSO4=TSO4                                              ! particulate sulfate    [umol/m^3]

       ! liquid matter
       ASO4 = ASO4 * 1.e6                                     ! aqueous phase sulfate  [umol/m^3] 
       ANH4 = ANH4 * 1.e6                                     ! aqueous phase ammonium [umol/m^3]
       ANO3 = ANO3 * 1.e6                                     ! aqueous phase nitrate  [umol/m^3]
       ACl  = ACl  * 1.e6                                     ! aqueous phase chloride [umol/m^3]
       ANa  = ANa  * 1.e6                                     ! aqueous phase sodium   [umol/m^3]

       ! solid matter
       SNH4=PNH4-ANH4                                         ! solid phase ammonium   [umol/m^3]
       SSO4=PSO4-ASO4                                         ! solid phase sulfate    [umol/m^3]
       SNO3=PNO3-ANO3                                         ! solid phase nitrate    [umol/m^3]
       SCl =PCl -ACl                                          ! solid phase chloride   [umol/m^3]
       SNa =PNa -ANa                                          ! solid phase sodium     [umol/m^3]

       ! GO SAVE ...

       IF(SNH4.LT.0.)   SNH4=0.
       IF(SSO4.LT.0.)   SSO4=0.
       IF(SNO3.LT.0.)   SNO3=0.
       IF(SCl.LT.0.)    SCl=0.
       IF(SNa.LT.0.)    SNa=0.

       PM=SNH4+SSO4+SNO3+SNH4+SCl+SNa+ANH4+ASO4+ANO3+ACl+ANa  ! total PM [umol/m^3]
       PMs=SNH4*MWNH4+SSO4*MWSO4+SNO3*MWNO3+SCl*MWCl+SNa*MWNa ! dry particulate matter (PM)     [ug/m^3]
       PMt=PMs+ANH4*MWNH4+ASO4*MWSO4+ANO3*MWNO3+ACl*MWCl+  &
            ANa*MWNa                                           ! total (dry + wet) PM, excl. H20 [ug/m^3]

       WH2O = WH2O * 1.e9                                     ! convert aerosol water from [kg/m^3] to [ug/m^3]
       H2O_NO3 = H2O_NO3 * 1.e9                                     ! convert NO3 water from [kg/m^3] to [ug/m^3]
       IF(WH2O.LT.1.e-3) WH2O=0.

       ! UPDATE HISTORY RH FOR HYSTERESIS (ONLINE CALCULATIONS ONLY)

       RH_HIST=2.                                             ! wet
       IF(WH2O.EQ.0.) RH_HIST=1.                              ! dry

       RINC = 1.
       IF(PMt.GT.0.)   RINC = (WH2O/PMt+1)**(1./3.)           ! approx. radius increase due to water uptake
       IF(RINC.EQ.0.)  RINC = 1.

       RATIONS = 0.
       IF(PSO4.GT.0.) RATIONS = PNO3/PSO4                     ! nitrate / sulfate mol ratio

       GR = 0.
       IF(GNO3.GT.0.) GR = GNH3/GNO3                          ! gas ratio = residual NH3 / residual HNO3   [-]

       DON = 0.
       IF((PNO3+2.*PSO4).GT.0.) DON = 100.*PNH4/(PNO3+2.*PSO4)! degree of neutralization by ammonia : ammonium / total nitrate + sulfate  [%]

       NO3P = 0.
       IF(TNO3.GT.0.) NO3P = 100.*PNO3/TNO3                   ! nitrate  partitioning = nitrate  / total nitrate    [%]

       NH4P = 0.
       IF(TNH4.GT.0.) NH4P = 100.*PNH4/TNH4                   ! ammonium partitioning = ammonium / total ammonium   [%]
       !
       ! store aerosol species for diagnostic output:
       !______________________________________________________________
       ! Input values:
       yo(il, 1) = TT   - 273.15                                       ! T                                    [degC]
       yo(il, 2) = RH   * 100.00                                       ! RH                                      [%]
       yo(il, 3) = TNH4                                                ! total input ammonium (g+p)       [umol/m^3]
       yo(il, 4) = TSO4                                                ! total input sulfate  (g+p)       [umol/m^3]
       yo(il, 5) = TNO3                                                ! total input nitrate  (g+p)       [umol/m^3]
       yo(il, 6) = TNa                                                 ! total input sodium   (p)         [umol/m^3]
       yo(il,33) = TCl                                                 ! total input chloride (g+p)       [umol/m^3]
       yo(il, 7) = TPo                                                 ! total input potasium (p)         [umol/m^3]
       yo(il,34) = TCa                                                 ! total input calcium  (p)         [umol/m^3]
       yo(il,35) = TMg                                                 ! total input magnesium(p)         [umol/m^3]
       yo(il,25) = PX                                                  ! atmospheric pressure                  [hPa]
       ! Output values:
       yo(il, 8) = GHCL                                                ! residual HCl   (g)               [umol/m^3]
       yo(il, 9) = GNO3                                                ! residual HNO3  (g)               [umol/m^3]
       yo(il,10) = GNH3                                                ! residual NH3   (g)               [umol/m^3]
       yo(il,11) = GSO4                                                ! residual H2SO4 (aq)              [umol/m^3]
       yo(il,12) = WH2O                                                ! aerosol Water  (aq)                [ug/m^3]
       yo(il,13) = PH                                                  ! aerosol pH                            [log]
       yo(il,14) = ZFLAG                                               ! concnetration domain [1=SP,2=SN,3=SR,4=SVR]
       yo(il,15) = PM                                                  ! total particulate matter         [umol/m^3]
       yo(il,16) = SNH4                                                ! solid ammonium (s)               [umol/m^3]
       yo(il,17) = SNO3                                                ! solid nitrate  (s)               [umol/m^3]
       yo(il,18) = SSO4                                                ! solid sulfate  (s)               [umol/m^3]
       yo(il,19) = PNH4                                                ! particulate ammonium (p=a+s)     [umol/m^3]
       yo(il,20) = PNO3                                                ! particulate nitrate  (p=a+s)     [umol/m^3]
       yo(il,21) = PSO4                                                ! particulate sulfate  (p=a+s)     [umol/m^3]
       yo(il,22) = RATIONS                                             ! mol ratio Nitrate/Sulfate (p)           [-]
       yo(il,23) = GAMA                                                ! activity coefficient (e.g. NH4NO3)           [-]
       yo(il,24) = ZIONIC                                              ! ionic strength (aq)                [mol/kg]
       yo(il,26) = PMt                                                 ! total PM (liquids & solids)        [ug/m^3]
       yo(il,27) = PMs                                                 ! total PM (solid)                   [ug/m^3]
       yo(il,28) = RINC                                                ! radius increase (H2O/PMt+1)**(1/3)      [-]
       yo(il,29) = SCl                                                 ! solid chloride (s)               [umol/m^3]
       yo(il,30) = SNa                                                 ! solid sodium (s)                 [umol/m^3]
       yo(il,31) = PCl                                                 ! particulate chloride (p=a+s)     [umol/m^3]
       yo(il,32) = PNa                                                 ! particulate sodium (p=a+s)       [umol/m^3]
       yo(il,36) = GG
       yo(il,37) = H2O_NO3                                             ! Water associated with nitrate      [ug/m^3]
    enddo
    !
  end subroutine eqsam_v03d


end module eqsam_param