!----------------------------------------------------------------------
! --- API_CHEM -- CALPUFF subroutines added for API (MCHEM=6,7)
!----------------------------------------------------------------------
!
! --- CALPUFF    Version: TNG-7.1.0    Level: 140521          API_CHEM
!
! --- PURPOSE: This collection of routines was prepared for API by
!              AER to add chemical transformation and aerosol options
!              to CALPUFF.  They are bundled as an include-file to
!              facilitate future modular upgrades.  Note that another
!              set of subroutines (isorropia.for) has not been placed
!              in api_chem.for and must be 'included' to complete
!              the code for the related options.
!
! --- IMPLEMENTATION:  Several modifications made to these new options
!              have been assigned either new user-selected modeling
!              option variables (changed via the control file) or
!              local logical operators that may be switched in an
!              individual subroutine to restore the original logic.
!              CALPUFF will need to be recompiled if a local logical
!              variable is changed.
!
!              The following variable settings will restore the original
!              features and logic of the API code (for testing):
!
!              Control File
!                     MLWC = 0 Do not read gridded cloud water from file
!                     MNH3 = 1 Read monthly ammonia vertical profiles
!                  MAVGNH3 = 0 Do not average ammonia vertical profiles
!                              across puff
!                   RNITE1 = 0.0 for no heterogeneous SO2 transformation
!              Local Logical in Subroutine
!                     L_KGPCM = .FALSE. in CHEMRIV6, CHEMRIV7
!                 L_TNO3FLOOR = .FALSE. in CHEMRIV6, CHEMRIV7
!                 L_RAINCLOUD = .FALSE. in CHEMRIV6, CHEMRIV7
!                     L_SCAV6 = .FALSE. in WET
!
! -----------------------------
! --- CONTENT:
! -----------------------------
! --- AER routines based on existing CALPUFF routines
!      subroutine chemriv6
!      subroutine chemriv7
!      subroutine chmriv6
!      subroutine chmriv7
!      subroutine setbckoc
!
! --- AER ISORROPIA interface routine
!     subroutine isodriver
!     (include 'isorropia.for')
!
! --- AER CALTECH SOA routines
!     subroutine soadriver
!     subroutine caltech_soa
!     subroutine soasub
!
! --- AQUEOUS-CHEMISTRY ROUTINES BASED ON RADM/CMAQ
!     subroutine aqradm
!     function hlconst
!     function index1
! -----------------------------
!
!
! --- UPDATE
!
! **********************************************************************
! --- Exponent, Inc. Updates:
! **********************************************************************
! --- V6.41-V6.42_x1.1 140521  : Call ISORROPIA with the METASTABLE
!                                control mode ON
!                                Modified: ISODRIVER
!                      140521  : Use iterative procedure in calling
!                                CALTECH_SOA to revise absorbing
!                                organic mass based on current aerosol
!                                mass as equilibrium is sought
!                                Modified: CHEMRIV7, SOADRIVER
!                      140521  : Use cloud fraction to set puff mass
!                                altered when using AUX-file LWC
!                                Use local temperature, pressure for AQ
!                                Adjust mass transformed by fraction of
!                                puff that overlays cloud water layers
!                                Modified: CHEMRIV6, CHEMRIV7
!                      140521  : Enforce NO3<=TNO3 to avoid negative
!                                HNO3 concentrations due to precision in
!                                HNO3=TNO3-NO3 operation
!                                Also guard against NO3<0
!                                Modified: CHEMRIV6, CHEMRIV7
!                      140521  : Add minimum RH and SO4 for ISORROPIA
!                                Modified: CHEMRIV6, CHEMRIV7
! **********************************************************************
!
! --- V6.4-V6.41   110301 (DGS): Treat aqueous-phase transformation
!                                case of precip without cloud water
!                                (from AUX file) by using a default
!                                lwc=0.5g/m3 with at least 10% cloud
!                                cover
!                                Modified: CHEMRIV6, CHEMRIV7
!
! --- V6.302-V6.4  101025 (DGS): Update ISORROPIA to V2.1
!                                Replaced: (include file isorropia.for)
!                                          (include file isrpia.inc)
!                                Modified: ISODRIVER
!
! --- V5.8-V6.302  100917 (DGS): Restructure value**-n to value**(-n)
!                                to satisfy LF95 compiler
!                                Modified: CHMRIV6, CHMRIV7
!                  100917 (DGS): Place NH3 profiles into /CHEMDAT/
!                                and add chembkz to chembk logic
!                  100917 (DGS): pass RSHET(fraction/hr) for setting
!                                heterogeneous SO2 reaction rate
!                                Modified: CHEMRIV6, CHEMRIV7
!                                          CHMRIV6, CHMRIV7
!                  100917 (DGS): Pass local cloud water mixing ratio
!                                when available from CALMET 3D file
!                                Modified: CHEMRIV6, CHEMRIV7
!                  100917 (DGS): Add local logical to convert LWC
!                                passed to AQRADM from g/m3 to kg/m3
!                                Modified: CHEMRIV6, CHEMRIV7
!                  100917 (DGS): Add local logical to use the minimum
!                                TNO3 as a floor rather than as a
!                                breakpoint for NO3=0.0
!                                Modified: CHEMRIV6, CHEMRIV7
!                  100917 (DGS): Add local logical to set a minimum
!                                cloud cover (10%) whenever there is
!                                liquid precipitation, and
!                                add but do not activate code to
!                                change effective scavenging
!                                coefficients to limit their action
!                                in any timestep to just the cloud
!                                fraction affecting the puff mass.
!                                Modified: CHEMRIV6, CHEMRIV7
!----------------------------------------------------------------------
      subroutine chemriv6(delt,qin,coz,ctnh3,ch2o2,maqchem,rshet,temp,         &
     &                    rhum,rhoair,pivol,zlen,zpuf,cldamr,fzcld,            &
     &                    cldt,cldp,cloud,prate,zcoef,nspec,ldb1,io6,          &
     &                    rh_isrp,so4_isrp,                                    &
     &                    rate,scav)
!----------------------------------------------------------------------
!
! --- CALPUFF    Version: TNG-7.1.0    Level: 140521          CHEMRIV6
!                P. Karamchandani, AER (Adapted from chemriv)
!
! --- PURPOSE:  This routine sets up call to RIVAD chemical
!               transformation subroutines, and passes results to
!               calling program
!
! --- UPDATES:
!
! --- V6.41-V6.42_x1.1 140521  : Use cloud fraction to set puff mass
!                                altered when using AUX-file LWC
!                                Use local temperature, pressure for AQ
!                                Adjust mass transformed by fraction of
!                                puff that overlays cloud water layers
!                                Condition NO3 to be no larger than TNO3
!                                and no smaller than zero
!                                Add minimum RH and SO4 for ISORROPIA
!
! --- V6.302-V6.41 110301 (DGS): Treat aqueous-phase transformation
!                                case of precip without cloud water
!                                (from AUX file) by using a default
!                                lwc=0.5g/m3 with at least 10% cloud
!                                cover
! --- V5.8-V6.302  100917 (DGS): Pass RSHET(fraction/hr) for setting
!                                heterogeneous SO2 reaction rate
!                         (DGS): Add local cloud water mixing ratio
!                                to input argument list (MLWC option)
!                         (DGS): Add local logical to convert LWC
!                                passed to AQRADM from g/m3 to kg/m3
!                         (DGS): Add local logical to use the minimum
!                                TNO3 as a floor rather than as a
!                                breakpoint for NO3=0.0
!                         (DGS): Add local logical to set a minimum
!                                cloud cover (10%) whenever there is
!                                liquid precipitation for use in the
!                                aqueous phase reaction/scavenging,
!                                and add but do not activate code to
!                                change the effective scavenging
!                                coefficients by limiting their action
!                                in any timestep to just the cloud
!                                fraction affecting the puff mass.
!
! --- INPUTS:
!         DELT - real    - integration time interval (hours)
!  QIN(mxspec) - real    - Pollutant mass (g) in the puff
!                            QIN(1) = SO2
!                            QIN(2) = SO4
!                            QIN(3) = NO
!                            QIN(4) = NO2
!                            QIN(5) = TNO3 (HNO3 + NO3)
!                            QIN(6) = NH4NO3
!                            QIN(+) = Not Used Here
!
!          COZ - real    - puff ozone concentration (ppb)
!        CTNH3 - real    - background ammonia concentration (ppb)
!        CH2O2 - real    - puff H2O2 concentration (ppb)
!      MAQCHEM - integer - Aqueous phase transformation flag
!                            0 = aqueous phase transformation
!                                not modeled
!                            1 = transformation rates adjusted
!                                for aqueous phase reactions
!        RSHET - real    - SO2 heterogeneous loss rate (fraction/hr)
!         TEMP - real    - temperature (deg. K)
!         RHUM - real    - relative humidity (percent)
!       RHOAIR - real    - surface air density (kg/m**3)
!        PIVOL - real    - Reciprocal of puff volume (1/m**3)
!         ZLEN - real    - Puff vertical length scale (m)
!         ZPUF - real    - Puff/Slug elevation (m MSL)
!       CLDAMR - real    - Average cloud water mixing ratio (g/kg)
!
! --- 6.42_x1.1
!        FZCLD - real    - Fraction of puff layer in cloud
!         CLDT - real    - Average cloud temperature (K)
!         CLDP - real    - Average cloud pressure (atm)
!      RH_ISRP - real    - Minimum relative humidity (%)
!                          for ISORROPIA
!     SO4_ISRP - real    - Minimum SO4 (g/m3) for ISORROPIA
!
!        CLOUD - real    - Cloud cover (tenths)
!        PRATE - real    - Precip. rate (mm/hr)
!        ZCOEF - real    - Cosine of solar zenith angle
!        NSPEC - real    - number of species
!         LDB1 - logical - Control variable for printing of debug
!                          information
!          IO6 - integer - Fortran unit number of printed output
!
!
! --- OUTPUT:
!  QIN(mxspec) - real    - Pollutant mass (g) in the puff
!                            QIN(1) = SO2
!                            QIN(2) = SO4
!                            QIN(3) = NO
!                            QIN(4) = NO2
!                            QIN(5) = TNO3 (HNO3 + NO3)
!                            QIN(6) = NH4NO3
!                            QIN(+) = Not Used Here
!      RATE(4) - real    - Transformation rates (percent/hour)
!                            R(1) -- SO2 loss rate
!                            R(2) -- NOX loss rate
!                            R(3) -- TNO3 formation rate
!                            R(4) -- NO  loss rate
!     SCAV(6) - real    -  Scavenging coefficients (1/s)
!                            SCAV(1) -- SO2
!                            SCAV(2) -- SO4
!                            SCAV(3) -- NO
!                            SCAV(4) -- NO2
!                            SCAV(5) -- HNO3
!                            SCAV(6) -- NO3
!
!
! --- CHEMRIV6 called by: CHEM
! --- CHEMRIV6 calls:     PHOT, CHMRIV6, ISODRIVER, AQRADM
!----------------------------------------------------------------------
!
      implicit none

! --- Arguments
! --- wangzhm save
      integer nspec
      integer maqchem, io6
      real delt, coz, ctnh3, ch2o2, temp, rhum, rhoair, pivol, zlen
      real zpuf, cloud, prate, zcoef, rshet
      real qin(nspec), scav(nspec)
      real rate(4)
      real cldamr

! --- 6.42_x1.1
      real cldt,cldp,fzcld
      real,save :: fcloud2
      real rh_isrp,so4_isrp
! --- wangzhm save
      real,save :: cldta,cldpa
      real,save :: confcta,rhoaira
      real,save :: test

      logical ldb1

! --- wangzhm save
! --- Locals
      real,save :: ppb(6),ppbi(6),q(6)
      real rmwt(6)
      real,save :: lwc    ! Liquid water content, g/m3
      real,save :: cozm, ch2o2m, ctnh3m  ! Concentrations in mols/mols air units
      real,save :: tno3floor
! --- 6.42_x1.1
      real,save :: tso4min
! --- Local controls
      logical l_kgpcm, l_tno3floor, l_raincloud

! --- Concentration array (moles/mole of air) for AQRADM
      real,save :: con(6)

! --- Note: TNO3 is weighted as NO3
      data rmwt/64.,96.,30.,46.,62.,62./

      real,save :: dt, vkgmol, confct, rk1, zcoefb, o3ppm
      real,save :: presur, f, ppbix,  ppbi4, ppbx
      real,save :: delno, tso4, tno3, tnh3, rhfrac, pno3
      real,save :: patm, taucld, rhoairm
      integer,save :: is, iss, j, ii, i

      real,save :: pcloud, fcloud

! ----------------------
! --- Set local controls
! ----------------------
! --- These enable 3 changes to the code to be reverted to their
! --- original condition.
! --- 1.  Selecting l_kgpcm=.TRUE. converts the LWC passed to AQRADM
! ---     from g/m3 to kg/m3, since AQRADM assumes it to be kg/m3.
! ---     Selecting l_kgpcm=.FALSE. retains the original code and
! ---     passes LWC to AQRADM in g/m3.
      data l_kgpcm/.TRUE./
! --- 2.  Selecting l_tno3floor=.TRUE. uses the cut-off TNO3
! ---     concentration as a floor, and computes NO3 corresponding
! ---     to this floor for all TNO3 less than this floor.  The
! ---     resulting ratio of NO3/TNO3(floor) is then multiplied by
! ---     the actual TNO3 to obtain the final NO3 concentration.
! ---     Selecting l_tno3floor=.FALSE. retains the original code
! ---     and NO3=0.0 for all TNO3<TNO3(floor).
      data l_tno3floor/.TRUE./
! --- 3.  Selecting l_raincloud=.TRUE. sets a minimum cloud cover
! ---     of 10% whenever liquid precipitation is non-zero.  This
! ---     only applies to the aqueous phase option and forces wet
! ---     removal whenever there is rain.  When cloud cover is zero,
! ---     aqueous conversion and wet removal are also zero.
! ---     (This also changes the effective scavenging coefficients
! ---     by limiting their action in any timestep to just the
! ---     cloud fraction affecting the puff mass.)--NA
      data l_raincloud/.TRUE./
! ----------------------

! --- wangzhm omp: set threadsprivate
      !$OMP THREADPRIVATE(cldta,cldpa,confcta,rhoaira,test,ppb,ppbi,q,&
      !$OMP& lwc,cozm,ch2o2m,ctnh3m,tno3floor,tso4min,&
      !$OMP& con,dt,vkgmol,confct,rk1,zcoefb,o3ppm,presur,f,ppbix,ppbi4,ppbx,delno,tso4,tno3,&
      !$OMP& tnh3,rhfrac,pno3,patm,taucld,rhoairm,is,iss,j,ii,i,pcloud,fcloud,fcloud2)

! --- Transfer mass data to local array
      iss = MIN(nspec,6)
      do is = 1, iss
         q(is) = qin(is)
      end do
      do is = (iss+1), 6
         q(is) = 0.0
      end do

! --- Compute initial concentrations as PPB
! --- kg-Molar volume (m^3) at ambient T,P computed from air density
      vkgmol = 28.97/rhoair
! --- Conversion factor to ppm in RIVAD
      confct = (1.0e-3)*vkgmol*pivol
      do is = 1, iss
         ppbi(is) = (1.0e09)*confct*q(is)/rmwt(is)
      end do

! --- Get NO2 photolysis rate
      call PHOT(zpuf,cloud,zcoef,rk1,zcoefb)

! --- Do conversions needed for RIVAD input arguments
! --- RIVAD weights TNO3 as HNO3: scale mass by 63/62=1.016129
      q(5) = q(5)*1.016129
! --- ozone in PPM rather than PPB
      o3ppm = 1.0e-03*coz
! --- Implied pressure (mb) -- use T0=273. P0=1013. as in RIVAD
      patm = (22.4/vkgmol)*(temp/273.)
      presur = 1013.*patm

      if(ldb1) then
         write(io6,*)
         write(io6,*)'CHEMRIV6:'
         write(io6,*)'    rk1,rhoair,O3ppb = ',rk1,rhoair,coz
         write(io6,*)'     TNH3ppb,H2O2ppb = ',ctnh3, ch2o2
         write(io6,*)'    presur,temp,rhum = ',presur,temp,rhum
         write(io6,*)'  cloud,prate,cldamr = ',cloud,prate,cldamr
         write(io6,*)'     zlen,zmsl,zcoef = ',zlen,zpuf,zcoef
         write(io6,*)' vkgmol,confct,pivol = ',vkgmol,confct,pivol

! --- 6.42_x1.1
         write(io6,*)'     fzcld,cldt,cldp = ',fzcld,cldt,cldp
         write(io6,*)'    rh_isrp,so4_isrp = ',rh_isrp,so4_isrp

         write(io6,*)
!         write(io6,101)
101      format('ppb =  SO2',7x,'SO4',7x,'NO',7x,'NO2',7x,'TNO3',              &
     &              6x,'NO3')
         write(io6,*)'Starting Concs -----'
         write(io6,101)
         write(io6,'(3x,6e10.2)') (ppbi(j),j=1,6)
         write(io6,*)'Starting Puff Mass (g) -----'
         write(io6,'(3x,6f10.2)') (qin(j),j=1,6)
      endif

! --- Call RIVAD module
      call CHMRIV6(q,temp,presur,rhum,o3ppm,rk1,zcoef,delt,confct,             &
     &             rshet)
! --- ozone in PPB
      coz = 1.0e03*o3ppm

! --- Re-weight TNO3 as NO3
      q(5) = q(5)*0.984127

      if(ldb1) then
         write(io6,*)'After CHMRIV6, Puff Mass (g) -----'
         write(io6,'(3x,6f10.2)') (q(j),j=1,6)
      endif

      if (maqchem.EQ.1) then
!
! --- initialize scavenging coefficients
         do ii = 1, 6
            scav(ii) = 0.
         end do

! --- Use local cloud water mixing ratio if valid, or
! --- Assign liquid water content if there is cloud cover and the
! --- temperature is above freezing
         fcloud = cloud * 0.1    ! fractional cloud cover
         pcloud = 0.0
         if(l_raincloud .AND. temp.GT.273.15) then
! ---       Force at least 10% cloud cover when there is rain
            if(prate.GT.0.0) fcloud=MAX(fcloud,0.1)
         endif

         if(fcloud.GT.0.0) pcloud = prate/fcloud

! --- 6.42_x1.1
         rhoaira=rhoair
         confcta=confct
         cldta=temp
         cldpa=patm
! ---    Fraction of puff within cloud horizontally and vertically
         fcloud2=fcloud
         if(cldamr.GE.0.0) then
! ---       Valid cloud water provided (use it)
! ---       Use average cloud temperature and pressure
            cldta=cldt
            cldpa=cldp
            rhoaira=rhoair*(cldpa/patm)*(temp/cldta)
            confcta=(1.0e-3)*(28.97/rhoaira)*pivol
! ---       Convert: g/kg(air) * rhoair(kg/m3) ==> g/m3(air)
            lwc=cldamr*rhoaira
! ---       Set minimum cloud fraction to 10%
            fcloud=MAX(fcloud,0.1)
            pcloud=prate/fcloud
! ---       Fraction of puff within cloud horizontally and vertically
            fcloud2=fcloud*fzcld

         elseif (fcloud.gt.0. .and. temp.gt.273.15) then
            if (prate .gt. 0.) then
               lwc = 0.5
            else
               lwc = 0.1
            end if
         else
            lwc = 0.
         end if

         if(l_kgpcm) then
! ---       Units for lwc in AQRADM are kg/m3
            lwc=0.001*lwc
         endif

         if (lwc .gt. 0.) then
!
! --- Calculate concs in moles/mole air units
            do ii = 1, 6

! --- 6.42_x1.1
               con(ii) = MAX(confcta*q(ii)/rmwt(ii),0.)

            end do
!
! --- Get HNO3 conc from total nitrate and PM nitrate
            con(5) = MAX(con(5) - con(6),0.)
!
            cozm = coz * 1.E-9
            ch2o2m = ch2o2 * 1.E-9
            ctnh3m = ctnh3 * 1.E-9
!
! --- Call aqueous-phase chemistry module
            taucld = delt * 3600.   ! timestep in seconds

! --- 6.42_x1.1
            rhoairm = rhoaira * 1.E3/28.97   ! air density in moles/m3
            if(ldb1) then
               write(io6,*)'AQRADM:cldta,cldpa,pcloud,lwc,cozm,'//             &
     &                     'ch2o2m,ctnh3m,con6= '
               write(io6,*)'called:   ',cldta,cldpa,pcloud,lwc,cozm,           &
     &                      ch2o2m,ctnh3m,con
            endif
            call AQRADM(cldta,cldpa,taucld,pcloud,lwc,cozm,ch2o2m,             &
     &                  ctnh3m,con,rhoairm,zlen,scav)
            if(ldb1) then
               write(io6,*)'returned: ',cldta,cldpa,pcloud,lwc,cozm,           &
     &                      ch2o2m,ctnh3m,con
               write(io6,*)'    scav: ',scav
            endif

! --- Assign adjusted SO2 and SO4 concs; adjust for cloud cover
            q(1) = q(1)*(1. - fcloud2) + fcloud2*con(1)*rmwt(1)/confcta
            q(2) = q(2)*(1. - fcloud2) + fcloud2*con(2)*rmwt(2)/confcta
!
! --- Adjusted oxidant concs
            coz = coz*(1. - fcloud2) + fcloud2*cozm*1.E9
            ch2o2 = ch2o2*(1. - fcloud2) + fcloud2*ch2o2m*1.E9

!
! --- Adjusted scavenging rates
            if (prate .gt. 0.) then
! ---          Possible alternate method (NA)
! ---           if(l_raincloud .AND. fcloud.LT.0.99) then
! ---              do ii = 1, 6
! ---                 scav(ii)=-ALOG(1.0 - fcloud*
! --- &                        (1.0 - EXP(-scav(ii)*taucld) ))/taucld
! ---              end do
! ---           else
                  do ii = 1, 6
                     scav(ii) = scav(ii) * fcloud
                  end do
! ---           endif
            end if
         end if
      end if

! --- Condition mass results
      do is = 1,iss
         q(is) = MAX(q(is),0.0)
      end do

! --- Compute ending concentrations as PPB
      do is = 1,iss
         ppb(is) = (1.0e09)*confct*q(is)/rmwt(is)
      end do

! --- Compute conversion rates (%/hr) for QA review
      do i = 1,4
         rate(i) = 0.0
      end do

! --- Compute these logs only if debug output is ON
      if(ldb1) then
         dt = .01*delt
! ---    SOX:
         if(ppbi(1).GT.0.0 .AND. ppb(1).GT.0.0) rate(1) = -ALOG(ppb(1)/        &
     &                                                   ppbi(1))/dt
! ---    NOX:
         ppbix = ppbi(4) + ppbi(3)
         ppbx = ppb(4) + ppb(3)
         if(ppbix.GT.0.0 .AND. ppbx.GT.0.0) rate(2) = -ALOG(ppbx/              &
     &                                                   ppbix)/dt
! ---    TNO3:  Conversion of NO2 after NO:NO2 process
         ppbi4 = ppbi(4) + ppbi(3) - ppb(3)
         if(ppbi4.GT.0.0 .AND. ppb(4).GT.0.0) rate(3) = -ALOG(ppb(4)/          &
     &                                                   ppbi4)/dt
! ---    NO:
         delno = ppb(3) - ppbi(3)
         if(delno.LT.0.0) then
            if(ppbi(3).GT.0.0 .AND. ppb(3).GT.0.0) rate(4) =                   &
     &                          ALOG(ppb(3)/ppbi(3))/dt
         elseif(delno.GT.0.0) then
            if(ppbi4.GT.0.0 .AND. ppb(4).GT.0.0) rate(4) = -ALOG(ppbi4/        &
     &                                                    ppbi(4))/dt
         endif
      endif

! --- Inorganic aerosol equilibrium with ISORROPIA
      if(q(5).gt.0.0)then
!
! --- concs in mols/m3

! --- 6.42_x1.1
! ---   Apply so4(g/m3) concentration constraint from control file
        tso4=q(2)*pivol
        tso4=MAX(tso4,so4_isrp)/rmwt(2)

        tso4 = MAX(tso4,1.E-12)
        tno3 = q(5) * pivol / rmwt(6)
        tnh3 = ctnh3 * rhoair / 28.97E6

! --- 6.42_x1.1
! ---   Apply RH(%) constraint from control file to compute RH fraction
        rhfrac=0.01*MAX(rhum,rh_isrp)

        if(tno3.gt.1.E-10)then
          if(ldb1) then

! --- 6.42_x1.1
             write(io6,*)'ISODRIVER:tso4,tno3,tnh3,rhfrac,temp     = '
             write(io6,*)'called:   ',tso4,tno3,tnh3,rhfrac,temp

          endif
          call isodriver(tso4,tno3,tnh3,rhfrac,temp,pno3)
          if(ldb1) then
             write(io6,*)'returned: ',tso4,tno3,tnh3,rhfrac,temp,pno3
          endif

        elseif(l_tno3floor) then
          tno3floor=1.E-10
          if(ldb1) then
             write(io6,*)                                                      &
     &             'ISODRIVER:tso4,tno3floor,tnh3,rhfrac,temp,pno3= '

! --- 6.42_x1.1
             write(io6,*)'called:   ',                                         &
     &                        tso4,tno3floor,tnh3,rhfrac,temp

          endif
          call isodriver(tso4,tno3floor,tnh3,rhfrac,temp,pno3)
          if(ldb1) then
             write(io6,*)'returned: ',                                         &
     &                        tso4,tno3floor,tnh3,rhfrac,temp,pno3
          endif
          pno3=pno3*(tno3/tno3floor)

        else
          pno3 = 0.
        endif
        q(6) = pno3 * rmwt(6) / pivol

! --- 6.42_x1.1
! ---   Condition NO3 (q(6)) to be no smaller than zero
        if(q(6).LT.0.0) then
           write(io6,*)'CHEMRIV6 Warning: reset NO3(g) from ',                 &
     &                 q(6),' to ZERO'
           q(6)=0.0
        endif
! ---   Condition NO3 to be no larger than TNO3 (q(5))
! ---   (q(5)>0.0 in this block)
        test=q(6)/q(5)-1.0
        if(test.GT.1.0e-06) then
           write(io6,*)'CHEMRIV6 Warning: reset NO3(g) from ',                 &
     &                 q(6),' to ',q(5)
        endif
        q(6)=MIN(q(5),q(6))

      endif

! --- Transfer revised mass data to original array
      do is = 1,iss
         qin(is) = q(is)
      end do

      if(ldb1) then
         write(io6,*)'After ISODRIVER, Puff Mass (g) -----'
         write(io6,'(3x,6f10.2)') (qin(j),j=1,6)
         write(io6,*)
      endif

      return
      end

!----------------------------------------------------------------------
      subroutine chemriv7(delt,qin,coz,ctnh3,ch2o2,maqchem,rshet,temp,         &
     &                    rhum,rhoair,pivol,zlen,zpuf,cldamr,fzcld,            &
     &                    cldt,cldp,cloud,prate,zcoef,nspec,ldb1,io6,          &
     &                    rh_isrp,so4_isrp,                                    &
     &                    rate,scav)
!----------------------------------------------------------------------
!
! --- CALPUFF    Version: TNG-7.1.0    Level: 140521          CHEMRIV7
!                P. Karamchandani, AER (Adapted from chemriv)
!
! --- PURPOSE:  This routine sets up call to RIVAD chemical
!               transformation subroutines, and passes results to
!               calling program
!
! --- UPDATES:
!
! --- V6.41-V6.42_x1.1 140521  : Use cloud fraction to set puff mass
!                                altered when using AUX-file LWC
!                                Use local temperature, pressure for AQ
!                                Adjust mass transformed by fraction of
!                                puff that overlays cloud water layers
!                                Condition NO3 to be no larger than TNO3
!                                and no smaller than zero
!                                Add list file unit to SOADRIVER
!                                Add minimum RH and SO4 for ISORROPIA
!
! --- V6.302-V6.41 110301 (DGS): Treat aqueous-phase transformation
!                                case of precip without cloud water
!                                (from AUX file) by using a default
!                                lwc=0.5g/m3 with at least 10% cloud
!                                cover
! --- V5.8-V6.302  100917 (DGS): Pass RSHET(fraction/hr) for setting
!                                heterogeneous SO2 reaction rate
!                         (DGS): Add local cloud water mixing ratio
!                                to input argument list (MLWC option)
!                         (DGS): Add local logical to convert LWC
!                                passed to AQRADM from g/m3 to kg/m3
!                         (DGS): Add local logical to use the minimum
!                                TNO3 as a floor rather than as a
!                                breakpoint for NO3=0.0
!                         (DGS): Add local logical to set a minimum
!                                cloud cover (10%) whenever there is
!                                liquid precipitation for use in the
!                                aqueous phase reaction/scavenging,
!                                and add but do not activate code to
!                                change the effective scavenging
!                                coefficients by limiting their action
!                                in any timestep to just the cloud
!                                fraction affecting the puff mass.
!
! --- INPUTS:
!         DELT - real    - integration time interval (hours)
!  QIN(mxspec) - real    - Pollutant mass (g) in the puff
!                            QIN(1) = SO2
!                            QIN(2) = SO4
!                            QIN(3) = NO
!                            QIN(4) = NO2
!                            QIN(5) = TNO3 (HNO3 + NO3)
!                            QIN(6) = NH4NO3
!                            QIN(7)  = Primary OC (POC)
!                            QIN(8)  = TOL (Toluene)
!                            QIN(9)  = TOLAER1 (Condensable product)
!                            QIN(10) = TOLAER2 (Condensable product)
!                            QIN(11) = ATOLA1 (SOA 1 from TOL)
!                            QIN(12) = ATOLA2 (SOA 2 from TOL)
!                            QIN(13) = XYL (Xylene)
!                            QIN(14) = XYLAER1 (Condensable product)
!                            QIN(15) = XYLAER2 (Condensable product)
!                            QIN(16) = AXYLA1 (SOA 1 from XYL)
!                            QIN(17) = AXYLA2 (SOA 2 from XYL)
!                            QIN(18) = ALKH (Higher alkanes)
!                            QIN(19) = ALKHAER (Condensable product)
!                            QIN(20) = ALKHA (SOA 1 from ALKH)
!                            QIN(21) = PAH
!                            QIN(22) = PAHAER1 (Condensable product)
!                            QIN(23) = PAHAER2 (Condensable product)
!                            QIN(24) = APAHA1 (SOA 1 from PAH)
!                            QIN(25) = APAHA2 (SOA 2 from PAH)
!
!          COZ - real    - puff ozone concentration (ppb)
!        CTNH3 - real    - background ammonia concentration (ppb)
!        CH2O2 - real    - puff H2O2 concentration (ppb)
!      MAQCHEM - integer - Aqueous phase transformation flag
!                            0 = aqueous phase transformation
!                                not modeled
!                            1 = transformation rates adjusted
!                                for aqueous phase reactions
!        RSHET - real    - SO2 heterogeneous loss rate (fraction/hr)
!         TEMP - real    - temperature (deg. K)
!         RHUM - real    - relative humidity (percent)
!       RHOAIR - real    - surface air density (kg/m**3)
!        PIVOL - real    - Reciprocal of puff volume (1/m**3)
!         ZLEN - real    - Puff vertical length scale (m)
!         ZPUF - real    - Puff/Slug elevation (m MSL)
!       CLDAMR - real    - Average cloud water mixing ratio (g/kg)
!
! --- 6.42_x1.1
!        FZCLD - real    - Fraction of puff layer in cloud
!         CLDT - real    - Average cloud temperature (K)
!         CLDP - real    - Average cloud pressure (atm)
!      RH_ISRP - real    - Minimum relative humidity (%)
!                          for ISORROPIA
!     SO4_ISRP - real    - Minimum SO4 (g/m3) for ISORROPIA
!
!        CLOUD - real    - Cloud cover (tenths)
!        PRATE - real    - Precip. rate (mm/hr)
!        ZCOEF - real    - Cosine of solar zenith angle
!        NSPEC - real    - number of species
!         LDB1 - logical - Control variable for printing of debug
!                          information
!          IO6 - integer - Fortran unit number of printed output
!
!
! --- OUTPUT:
!  QIN(mxspec) - real    - Pollutant mass (g) in the puff
!                            QIN(1)  = SO2
!                            QIN(2)  = SO4
!                            QIN(3)  = NO
!                            QIN(4)  = NO2
!                            QIN(5)  = TNO3 (HNO3 + NO3)
!                            QIN(6)  = NH4NO3
!                            QIN(7)  = Primary OC (POC)
!                            QIN(8)  = TOL (Toluene)
!                            QIN(9)  = TOLAER1 (Condensable product)
!                            QIN(10) = TOLAER2 (Condensable product)
!                            QIN(11) = ATOLA1 (SOA 1 from TOL)
!                            QIN(12) = ATOLA2 (SOA 2 from TOL)
!                            QIN(13) = XYL (Xylene)
!                            QIN(14) = XYLAER1 (Condensable product)
!                            QIN(15) = XYLAER2 (Condensable product)
!                            QIN(16) = AXYLA1 (SOA 1 from XYL)
!                            QIN(17) = AXYLA2 (SOA 2 from XYL)
!                            QIN(18) = ALKH (Higher alkanes)
!                            QIN(19) = ALKHAER (Condensable product)
!                            QIN(20) = ALKHA (SOA 1 from ALKH)
!                            QIN(21) = PAH
!                            QIN(22) = PAHAER1 (Condensable product)
!                            QIN(23) = PAHAER2 (Condensable product)
!                            QIN(24) = APAHA1 (SOA 1 from PAH)
!                            QIN(25) = APAHA2 (SOA 2 from PAH)
!      RATE(8) - real    - Transformation rates (percent/hour)
!                            R(1) -- SO2 loss rate
!                            R(2) -- NOX loss rate
!                            R(3) -- TNO3 formation rate
!                            R(4) -- NO  loss rate
!                            R(5) -- TOL loss rate
!                            R(6) -- XYL loss rate
!                            R(7) -- ALKH loss rate
!                            R(8) -- PAH loss rate
!     SCAV(6) - real    -  Scavenging coefficients (1/s)
!                            SCAV(1) -- SO2
!                            SCAV(2) -- SO4
!                            SCAV(3) -- NO
!                            SCAV(4) -- NO2
!                            SCAV(5) -- HNO3
!                            SCAV(6) -- NO3
!
! --- CHEMRIV7 called by: CHEM
! --- CHEMRIV7 calls:     PHOT, CHMRIV7, ISODRIVER, SOADRIVER, AQRADM
!----------------------------------------------------------------------
!
      implicit none

! --- wangzhm save

! --- Arguments
      integer nspec
      integer maqchem, io6
      real delt, coz, ctnh3, ch2o2, temp, rhum, rhoair, pivol, zlen
      real zpuf, cloud, prate, zcoef, rshet
      real qin(nspec), scav(nspec)
      real rate(8)
      real cldamr

! --- 6.42_x1.1
      real cldt,cldp,fzcld
      real,save :: fcloud2
      real rh_isrp,so4_isrp

      real,save :: cldta,cldpa
      real,save :: confcta,rhoaira
      real,save :: test

      logical ldb1

! --- Locals
      real,save :: ppb(25),ppbi(25),q(25)
      real rmwt(25)
      real,save :: lwc    ! Liquid water content, g/m3
      real,save :: cozm, ch2o2m, ctnh3m  ! Concentrations in mols/mols air units
      real,save :: tno3floor
! --- 6.42_x1.1
      real,save :: tso4min
! --- Local controls
      logical l_kgpcm, l_tno3floor, l_raincloud

! --- Concentration array (moles/mole of air) for AQRADM
      real,save :: con(6)

! --- Note: TNO3 is weighted as NO3
      data rmwt/64.,96.,30.,46.,62.,62.,180.,                                  &
     &          5*92.,5*106.,3*226.,5*156./

      real,save :: dt, vkgmol, confct, rk1, zcoefb, o3ppm
      real,save :: presur, f, ppbix,  ppbi4, ppbx
      real,save :: delno, tso4, tno3, tnh3, rhfrac, pno3
      real,save :: patm, taucld, rhoairm
      integer,save :: is, iss, j, ii, i
!
      real,save :: pcloud, fcloud

! ----------------------
! --- Set local controls
! ----------------------
! --- These enable 2 changes to the code to be reverted to their
! --- original condition.
! --- 1.  Selecting l_kgpcm=.TRUE. converts the LWC passed to AQRADM
! ---     from g/m3 to kg/m3, since AQRADM assumes it to be kg/m3.
! ---     Selecting l_kgpcm=.FALSE. retains the original code and
! ---     passes LWC to AQRADM in g/m3.
      data l_kgpcm/.TRUE./
! --- 2.  Selecting l_tno3floor=.TRUE. uses the cut-off TNO3
! ---     concentration as a floor, and computes NO3 corresponding
! ---     to this floor for all TNO3 less than this floor.  The
! ---     resulting ratio of NO3/TNO3(floor) is then multiplied by
! ---     the actual TNO3 to obtain the final NO3 concentration.
! ---     Selecting l_tno3floor=.FALSE. retains the original code
! ---     and NO3=0.0 for all TNO3<TNO3(floor).
      data l_tno3floor/.TRUE./
! --- 3.  Selecting l_raincloud=.TRUE. sets a minimum cloud cover
! ---     of 10% whenever liquid precipitation is non-zero.  This
! ---     only applies to the aqueous phase option and forces wet
! ---     removal whenever there is rain.  When cloud cover is zero,
! ---     aqueous conversion and wet removal are also zero.
! ---     (This also changes the effective scavenging coefficients
! ---     by limiting their action in any timestep to just the
! ---     cloud fraction affecting the puff mass.)--NA
      data l_raincloud/.TRUE./
! ----------------------
! --- wangzhm omp: set threadsprivate
      !$OMP THREADPRIVATE(cldta,cldpa,confcta,rhoaira,test,ppb,ppbi,q,&
      !$OMP& lwc,cozm,ch2o2m,ctnh3m,tno3floor,tso4min,&
      !$OMP& con,dt,vkgmol,confct,rk1,zcoefb,o3ppm,presur,f,ppbix,ppbi4,ppbx,delno,&
      !$OMP& tso4,tno3,tnh3,rhfrac,pno3,patm,taucld,rhoairm,is,iss,j,ii,i,pcloud,fcloud,fcloud2)

! --- Transfer mass data to local array
      iss = MIN(nspec,25)
      do is = 1, iss
         q(is) = qin(is)
      end do
      do is = (iss+1), 25
         q(is) = 0.0
      end do

! --- Compute initial concentrations as PPB
! --- kg-Molar volume (m^3) at ambient T,P computed from air density
      vkgmol = 28.97/rhoair
! --- Conversion factor to ppm in RIVAD
      confct = (1.0e-3)*vkgmol*pivol
      do is = 1, iss
         ppbi(is) = (1.0e09)*confct*q(is)/rmwt(is)
      end do

! --- Get NO2 photolysis rate
      call PHOT(zpuf,cloud,zcoef,rk1,zcoefb)

! --- Do conversions needed for RIVAD input arguments
! --- RIVAD weights TNO3 as HNO3: scale mass by 63/62=1.016129
      q(5) = q(5)*1.016129
! --- ozone in PPM rather than PPB
      o3ppm = 1.0e-03*coz
! --- Implied pressure (mb) -- use T0=273. P0=1013. as in RIVAD
      patm = (22.4/vkgmol)*(temp/273.)
      presur = 1013.*patm

      if(ldb1) then
         write(io6,*)'CHEMRIV7:'
         write(io6,*)'    rk1,rhoair,O3ppb = ',rk1,rhoair,coz
         write(io6,*)'     TNH3ppb,H2O2ppb = ',ctnh3, ch2o2
         write(io6,*)'    presur,temp,rhum = ',presur,temp,rhum
         write(io6,*)'  cloud,prate,cldamr = ',cloud,prate,cldamr
         write(io6,*)'     zlen,zmsl,zcoef = ',zlen,zpuf,zcoef
         write(io6,*)' vkgmol,confct,pivol = ',vkgmol,confct,pivol

! --- 6.42_x1.1
         write(io6,*)'     fzcld,cldt,cldp = ',fzcld,cldt,cldp
         write(io6,*)'    rh_isrp,so4_isrp = ',rh_isrp,so4_isrp

         write(io6,*)
!         write(io6,101)
101      format('ppb =  SO2',7x,'SO4',7x,'NO',7x,'NO2',7x,'TNO3',              &
     &              6x,'NO3')
         write(io6,*)'Starting Concs -----'
         write(io6,101)
         write(io6,'(3x,6e10.2)') (ppbi(j),j=1,6)
      endif

! --- Call RIVAD module
      call CHMRIV7(q,temp,presur,rhum,o3ppm,rk1,zcoef,delt,confct,             &
     &             rshet)
! --- ozone in PPB
      coz = 1.0e03*o3ppm

! --- Re-weight TNO3 as NO3
      q(5) = q(5)*0.984127

      if(ldb1) then
         write(io6,*)'After CHMRIV6, Puff Mass (g) -----'
         write(io6,'(3x,6f10.2)') (q(j),j=1,6)
      endif

      if (maqchem.EQ.1) then
!
! --- initialize scavenging coefficients
         do ii = 1, 6
            scav(ii) = 0.
         end do

! --- Use local cloud water mixing ratio if valid, or
! --- Assign liquid water content if there is cloud cover and the
! --- temperature is above freezing
         fcloud = cloud * 0.1    ! fractional cloud cover
         pcloud = 0.0
         if(l_raincloud .AND. temp.GT.273.15) then
! ---       Force at least 10% cloud cover when there is rain
            if(prate.GT.0.0) fcloud=MAX(fcloud,0.1)
         endif

         if(fcloud.GT.0.0) pcloud = prate/fcloud

! --- 6.42_x1.1
         rhoaira=rhoair
         confcta=confct
         cldta=temp
         cldpa=patm
! ---    Fraction of puff within cloud horizontally and vertically
         fcloud2=fcloud
         if(cldamr.GE.0.0) then
! ---       Valid cloud water provided (use it)
! ---       Use average cloud temperature and pressure
            cldta=cldt
            cldpa=cldp
            rhoaira=rhoair*(cldpa/patm)*(temp/cldta)
            confcta=(1.0e-3)*(28.97/rhoaira)*pivol
! ---       Convert: g/kg(air) * rhoair(kg/m3) ==> g/m3(air)
            lwc=cldamr*rhoaira
! ---       Set minimum cloud fraction to 10%
            fcloud=MAX(fcloud,0.1)
            pcloud=prate/fcloud
! ---       Fraction of puff within cloud horizontally and vertically
            fcloud2=fcloud*fzcld

         elseif (fcloud.gt.0. .and. temp.gt.273.15) then
            if (prate .gt. 0.) then
               lwc = 0.5
            else
               lwc = 0.1
            end if
         else
            lwc = 0.
         end if

         if(l_kgpcm) then
! ---       Units for lwc in AQRADM are kg/m3
            lwc=0.001*lwc
         endif

         if (lwc .gt. 0.) then
!
! --- Calculate concs in moles/mole air units
            do ii = 1, 6

! --- 6.42_x1.1
               con(ii) = MAX(confcta*q(ii)/rmwt(ii),0.)

            end do
!
! --- Get HNO3 conc from total nitrate and PM nitrate
            con(5) = MAX(con(5) - con(6),0.)
!
            cozm = coz * 1.E-9
            ch2o2m = ch2o2 * 1.E-9
            ctnh3m = ctnh3 * 1.E-9
!
! --- Call aqueous-phase chemistry module
            taucld = delt * 3600.   ! timestep in seconds

! --- 6.42_x1.1
            rhoairm = rhoaira * 1.E3/28.97   ! air density in moles/m3
            if(ldb1) then
               write(io6,*)'AQRADM:cldta,cldpa,pcloud,lwc,cozm,'//             &
     &                     'ch2o2m,ctnh3m,con6= '
               write(io6,*)'called:   ',cldta,cldpa,pcloud,lwc,cozm,           &
     &                      ch2o2m,ctnh3m,con
            endif
            call AQRADM(cldta,cldpa,taucld,pcloud,lwc,cozm,ch2o2m,             &
     &                  ctnh3m,con,rhoairm,zlen,scav)
            if(ldb1) then
               write(io6,*)'returned: ',cldta,cldpa,pcloud,lwc,cozm,           &
     &                      ch2o2m,ctnh3m,con
               write(io6,*)'    scav: ',scav
            endif

! --- Assign adjusted SO2 and SO4 concs; adjust for cloud cover
            q(1) = q(1)*(1. - fcloud2) + fcloud2*con(1)*rmwt(1)/confcta
            q(2) = q(2)*(1. - fcloud2) + fcloud2*con(2)*rmwt(2)/confcta
!
! --- Adjusted oxidant concs
            coz = coz*(1. - fcloud2) + fcloud2*cozm*1.E9
            ch2o2 = ch2o2*(1. - fcloud2) + fcloud2*ch2o2m*1.E9

!
! --- Adjusted scavenging rates
            if (prate .gt. 0.) then
! ---          Possible alternate method (NA)
! ---           if(l_raincloud .AND. fcloud.LT.0.99) then
! ---              do ii = 1, 6
! ---                 scav(ii)=-ALOG(1.0 - fcloud*
! --- &                        (1.0 - EXP(-scav(ii)*taucld) ))/taucld
! ---              end do
! ---           else
                  do ii = 1, 6
                     scav(ii) = scav(ii) * fcloud
                  end do
! ---           endif
            end if
         end if
      end if

! --- Condition mass results
      do is = 1,iss
         q(is) = MAX(q(is),0.0)
      end do

! --- Compute ending concentrations as PPB
      do is = 1,iss
         ppb(is) = (1.0e09)*confct*q(is)/rmwt(is)
      end do

! --- Compute conversion rates (%/hr) for QA review
      do i = 1,8
         rate(i) = 0.0
      end do

! --- Compute these logs only if debug output is ON
      if(ldb1) then
         dt = .01*delt
! ---    SOX:
         if(ppbi(1).GT.0.0 .AND. ppb(1).GT.0.0) rate(1) = -ALOG(ppb(1)/        &
     &                                                   ppbi(1))/dt
! ---    NOX:
         ppbix = ppbi(4) + ppbi(3)
         ppbx = ppb(4) + ppb(3)
         if(ppbix.GT.0.0 .AND. ppbx.GT.0.0) rate(2) = -ALOG(ppbx/              &
     &                                                   ppbix)/dt
! ---    TNO3:  Conversion of NO2 after NO:NO2 process
         ppbi4 = ppbi(4) + ppbi(3) - ppb(3)
         if(ppbi4.GT.0.0 .AND. ppb(4).GT.0.0) rate(3) = -ALOG(ppb(4)/          &
     &                                                   ppbi4)/dt
! ---    NO:
         delno = ppb(3) - ppbi(3)
         if(delno.LT.0.0) then
            if(ppbi(3).GT.0.0 .AND. ppb(3).GT.0.0) rate(4) =                   &
     &                          ALOG(ppb(3)/ppbi(3))/dt
         elseif(delno.GT.0.0) then
            if(ppbi4.GT.0.0 .AND. ppb(4).GT.0.0) rate(4) = -ALOG(ppbi4/        &
     &                                                    ppbi(4))/dt
! ---    TOL:
         if(ppbi(8).GT.0.0 .AND. ppb(8).GT.0.0) rate(5) = -ALOG(ppb(8)/        &
     &                                                    ppbi(8))/dt
! ---    XYL:
         if(ppbi(13).GT.0.0 .AND. ppb(13).GT.0.0) rate(6) =                    &
     &                                  -ALOG(ppb(13)/ppbi(13))/dt
! ---    ALKH:
         if(ppbi(18).GT.0.0 .AND. ppb(18).GT.0.0) rate(7) =                    &
     &                                  -ALOG(ppb(18)/ppbi(18))/dt
! ---    PAH:
         if(ppbi(21).GT.0.0 .AND. ppb(21).GT.0.0) rate(8) =                    &
     &                                  -ALOG(ppb(21)/ppbi(21))/dt
         endif
      endif

! --- Inorganic aerosol equilibrium with ISORROPIA
      if(q(5).gt.0.0)then
!
! --- concs in mols/m3

! --- 6.42_x1.1
! ---   Apply so4(g/m3) concentration constraint from control file
        tso4=q(2)*pivol
        tso4=MAX(tso4,so4_isrp)/rmwt(2)

        tso4 = MAX(tso4,1.E-12)
        tno3 = q(5) * pivol / rmwt(6)
        tnh3 = ctnh3 * rhoair / 28.97E6

! --- 6.42_x1.1
! ---   Apply RH(%) constraint from control file to compute RH fraction
        rhfrac=0.01*MAX(rhum,rh_isrp)

        if(tno3.gt.1.E-10)then
          if(ldb1) then

! --- 6.42_x1.1
             write(io6,*)'ISODRIVER:tso4,tno3,tnh3,rhfrac,temp     = '
             write(io6,*)'called:   ',tso4,tno3,tnh3,rhfrac,temp

          endif
          call isodriver(tso4,tno3,tnh3,rhfrac,temp,pno3)
          if(ldb1) then
             write(io6,*)'returned: ',tso4,tno3,tnh3,rhfrac,temp,pno3
          endif

        elseif(l_tno3floor) then
          tno3floor=1.E-10
          if(ldb1) then
             write(io6,*)                                                      &
     &             'ISODRIVER:tso4,tno3floor,tnh3,rhfrac,temp,pno3= '

! --- 6.42_x1.1
             write(io6,*)'called:   ',                                         &
     &                        tso4,tno3floor,tnh3,rhfrac,temp

          endif
          call isodriver(tso4,tno3floor,tnh3,rhfrac,temp,pno3)
          if(ldb1) then
             write(io6,*)'returned: ',                                         &
     &                        tso4,tno3floor,tnh3,rhfrac,temp,pno3
          endif
          pno3=pno3*(tno3/tno3floor)

        else
          pno3 = 0.
        endif
        q(6) = pno3 * rmwt(6) / pivol

! --- 6.42_x1.1
! ---   Condition NO3 (q(6)) to be no smaller than zero
        if(q(6).LT.0.0) then
           write(io6,*)'CHEMRIV6 Warning: reset NO3(g) from ',                 &
     &                 q(6),' to ZERO'
           q(6)=0.0
        endif
! ---   Condition NO3 to be no larger than TNO3 (q(5))
! ---   (q(5)>0.0 in this block)
        test=q(6)/q(5)-1.0
        if(test.GT.1.0e-06) then
           write(io6,*)'CHEMRIV6 Warning: reset NO3(g) from ',                 &
     &                 q(6),' to ',q(5)
        endif
        q(6)=MIN(q(5),q(6))

      endif
!
! --- SOA equilibrium
! --- 6.42_x1.1
      call soadriver(io6,q,temp,presur,pivol)
!
! --- Transfer revised mass data to original array
      do is = 1,iss
         qin(is) = q(is)
      end do

      return
      end

!----------------------------------------------------------------------
      subroutine chmriv6(PM,TEMPER,PRESUR,RH,O3PUFF,RK1,ZCOEF,                 &
     &                   TSTP,CONFCT,RSHET)
!----------------------------------------------------------------------
!
! --- CALPUFF    Version: TNG-7.1.0    Level: 100917           CHMRIV6
!                P. Karamchandani, AER (Adapted from chemriv)
!
! --- ADOPTED FROM ARM3 (see banner below)
!                Sulfate & nitrate conversion uses exponential fcn
!                Unused variables are deactivated
!
! --- UPDATE
! --- V5.8-V6.302  100917 (DGS): Restructure value**-n to value**(-n)
!                                to satisfy LF95 compiler
!                  100917 (DGS): Add heterogeneous reaction rate for
!                                SO4 and pass RSHET(fraction/hr)
! --- V5.8         071025  (PK): Original
!
!----------------------------------------------------------------------
      IMPLICIT NONE
      SAVE
!
!DECK.CHMRIV
!
!  DATE: NOVEMBER 1987
!  VERSION: ARM3-1.0
!
!     CALCULATE TRANSFORMATION OF SO2 TO SO4 AND NO2 TO HNO3-NO3
!     FOR ONE TIME STEP, TSTP
!     BASED ON RIVAD CHEMICAL MECHANISM
!     ADPATED FROM THE RIVAD MODEL
!
!
!  INPUT ARGUMENTS:
!                   PM     R  MASS OF SPECIES (I) IN PUFF (G)
!                             SO2, SO4, NO, NO2, HNO3-NO3, TSP
!                   TEMPER R  TEMPERATURE AT PLUME HEIGHT (K)
!                   PRESUR R  PRESSURE AT PLUME HEIGHT (MB)
!                   RH     R  RELATIVE HUMIDITY AT PLUME (%)
!                   O3PUFF R  PUFF OZONE CONCENTRATION (PPM)
!                   RK1    R  NO2 PHOTOLYSIS RATE CONSTANT (PPM/MIN)
!                   ZCOEF  R  COSINE OF SOLAR ZENITH ANGLE
!                   TSTP   R  TIME STEP (HOURS)
!                   CONFCT R  CONVERSION FACTOR FOR UG/M3 TO PPM-MWT
!                   RSHET  R  SO2 heterogeneous loss rate (fraction/hr)
!
!  OUTPUT ARGUMENTS:
!                   PM     R  NEW MASS OF PUFF DUE TO CHANGES BASED ON
!
!
!  SUBROUTINES CALLED:
!
!  CALLED BY:   CHEMRIV6
!

! --- Constants
      REAL COEF1                  ! Molec/cc to ppm conv factor coefficient
      PARAMETER ( COEF1 = 7.33981E+15 )

      REAL CONSTC                 ! Constant for falloff type reaction
      PARAMETER ( CONSTC = 0.6 )

      REAL TI300                  ! 1.0 / 300.
      PARAMETER ( TI300 = 1.0 / 300.0 )

! --- Arguments
      REAL PM(6)
      REAL TEMPER, PRESUR, RH, O3PUFF, RK1, ZCOEF, TSTP, CONFCT, RSHET

! --- wangzhm save

! --- Locals
      REAL,save :: TFACT, R26, ROHM, QJ, PHIKK, RCONST
      INTEGER,save :: L
      REAL,save :: CNO, CNO2, CNOX, TAMB, PAMB, RH1, H2O, SUM
      REAL,save :: XNO, XNOX, XO3, XOX, XNO2, RNO2X, XNO3, XN2O5
      REAL,save :: RSULF, RNITR, RKX, RNO3, RNITRN, RNITRD, SULFN
      REAL,save :: XOHMAX, XOH
      REAL,save :: A0, A, B, C

      REAL,save :: RK0               ! K0 in falloff rate expressions
      REAL,save :: RKINF             ! KINF in falloff rate expressions
      REAL,save :: XEND              ! Exponent in falloff rate expressions
!
      REAL,save :: KSO2OH            ! SO2 + OH rate constant
      REAL,save :: KNO2OH            ! NO2 + OH rate constant
!
      REAL*8,save :: TINV, CFACT, RFACT
!
      REAL,save :: PPM(5),SMWT(5)
      LOGICAL,save :: LFIRST
      DATA LFIRST/.TRUE./
      DATA SMWT/64.0,96.0,30.0,46.0,63.0/
!
!     FIRST TIME THROUGH SET UP SOME GLOBAL RATE CONSTANTS
!
! --- wangzhm omp: set threadsprivate
      !$OMP THREADPRIVATE(TFACT,R26,ROHM,QJ,PHIKK,RCONST,L,&
      !$OMP& CNO,CNO2,CNOX,TAMB,PAMB,RH1,H2O,SUM,&
      !$OMP& XNO,XNOX,XO3,XOX,XNO2,RNO2X,XNO3,XN2O5,&
      !$OMP& RSULF,RNITR,RKX,RNO3,RNITRN,RNITRD,SULFN,&
      !$OMP& XOHMAX,XOH,A0,A,B,C,RK0,RKINF,XEND,KSO2OH,KNO2OH,&
      !$OMP& TINV,CFACT,RFACT,ppm,smwt,lfirst)

      IF (LFIRST) THEN
         TFACT = 1./273.
         R26 = 1./26.4
         ROHM = 4.87E-7/1.32E-3
         LFIRST = .FALSE.
      END IF
!
      DO L = 1,5
         PPM(L) = 1.E6*CONFCT*PM(L)/SMWT(L)
      END  DO
!
      CNO = PPM(3)
      CNO2 = PPM(4)
      CNOX = CNO + CNO2
!
      TAMB = TEMPER
      PAMB = PRESUR/1013.0
      TINV = 1. / TAMB
      CFACT = COEF1 * PAMB * TINV
!
      IF (RK1.LE.0.0) THEN
         QJ = 0.0
         PHIKK = 0.0
      ELSE
         PHIKK = RK1*R26
         QJ = (1.338E-3)*ZCOEF**2.74
      ENDIF
      IF (TAMB.GE.273.) THEN
         RCONST = 18.02*(597.3-.566*(TAMB-273.))/1.9869
      ELSE
         RCONST = 6133.17
      END IF
      RH1 = MIN(RH,95.)
      RH1 = MAX(RH1,0.)
      H2O = (6030.*.01*RH1/PAMB)*EXP(RCONST*(TFACT-TINV))
!
!  DO SIMPLE CHEMISTRY
!
      XNOX = CNOX
      XOX = O3PUFF + CNO2
      SUM = XOX + XNOX + PHIKK
      XNO2 = 0.5*(SUM-SQRT(ABS(SUM*SUM-4.*XNOX*XOX)))
      XNO2 = MIN(XNO2,XNOX)
      XNO2 = MAX(XNO2,0.)
! dgs      XNO=XNOX-XNO2
! dgs      IF (XNO.LT.0.0) XNO=0.0
      RNO2X = 0.
      IF (XNOX.GT.0.) RNO2X = XNO2/XNOX

! dgs Use NO2/NOX ratio to test for zero NO (single precision)
      xno = xnox - xno2
      if(rno2x.GE.0.999999) XNO=0.0

      XO3 = XOX - XNO2
      XO3 = MAX(XO3,0.)

      IF (ZCOEF.LT.0.06975) THEN
!
!  NIGHTTIME CHEMISTRY
!
         RSULF = 0.
!
!  IMPLEMENT NEW NIGHTTIME CHEMISTRY - 1/86
!
         RNITR = 0.0
         XO3 = MAX(XO3,0.)
         IF (XO3.GT.0.) THEN
            RKX = 1780./(1.9E-6*H2O + 3.12)
            RNO3 = 0.086
            A0 = 0.59 + 1780. -3.12*RKX
            A = 2.*RKX*RNO3-A0
            B = RNO3 + 0.0474*XO3 +XNO2*A0
            C = -0.0474*XNO2*XO3
!
            XNO3 = -2.*C/(B+SQRT(B*B-4.*A*C))
            XN2O5 = XNO3*XNO2*RKX
            RNITRN = 2.*1.9E-6*H2O*XN2O5*60.*TSTP
            RNITRN = MIN(XNO2,RNITRN)
            PPM(5) = PPM(5) + RNITRN
            XNO2 = XNO2 - RNITRN
         END IF
      ELSE
!
!  DO DAYTIME CHEMISTRY
!
         XOHMAX = ROHM*QJ
         XOH = XOHMAX
         IF (PPM(1).NE.0..OR.XNO2.NE.0.)                                       &
     &   XOH = 2.*QJ*3.4E5*H2O*XO3/((4.45E10+3.4E5*H2O)*                       &
     &       (2000.*PPM(1) + 14000.*XNO2))
         XOH = MIN(XOH,XOHMAX)

! --- SO2 + OH rate constant (falloff expression)
         RK0 = 1.0E+06 * CFACT * 3.0E-31 * ( TAMB * TI300 )**(-3.3)
         RKINF = 1.5E-12
         XEND = 1.0 / ( 1.0 + ( LOG10( RK0 / RKINF ) )**2 )
         KSO2OH = ( RK0 / ( 1.0 + RK0 / RKINF ) ) * CONSTC**XEND

! --- NO2 + OH rate constant (falloff expression)
         RK0 = 1.0E+06 * CFACT * 2.6E-30 * ( TAMB * TI300 )**(-3.2)
         RKINF = 2.4E-11 * ( TAMB * TI300 )**(-1.3)
         XEND = 1.0 / ( 1.0 + ( LOG10( RK0 / RKINF ) )**2 )
         KNO2OH = ( RK0 / ( 1.0 + RK0 / RKINF ) ) * CONSTC**XEND
!
! --- Convert rate constants from molec-cc-1 s-1 to ppm-1 hr-1 units
         rfact = 2.64e19 * (pamb * tinv)
         kso2oh = kso2oh * rfact
         kno2oh = kno2oh * rfact
         rsulf = 1. - EXP(-kso2oh*xoh*tstp)
         rnitr = 1. - EXP(-kno2oh*xoh*tstp)
         RNITRD = RNITR*XNO2
         RNITRD = MIN(XNOX,RNITRD)
         PPM(5) = PPM(5)+RNITRD
         XNOX = XNOX-RNITRD
         XNOX =  MAX(XNOX,0.)
         XNO = XNOX*(1.-RNO2X)
         XNO2 = XNOX*RNO2X

      END IF
!
! --- Add heterogeneous rate wso4=rshet*ppm(1)*tstp
      SULFN = PPM(1)*RSULF + rshet*ppm(1)*tstp
      SULFN = MIN(PPM(1),SULFN)
      PPM(2) = PPM(2) + SULFN
      PPM(1) = PPM(1) - SULFN
!
!     UPDATE NEW STEADY-STATE NO AND NO2 VALUES
!
      PPM(3) = XNO
      PPM(4) = XNO2
!
!     UPDATE NEW O3 VALUE
      O3PUFF = XO3
!
!     CONVERT BACK TO GRAMS OF SPECIES IN PUFF
!
      DO L = 1,5
         PM(L) = 1.E-6*PPM(L)*SMWT(L)/CONFCT
      END DO
      RETURN
      END

!----------------------------------------------------------------------
      subroutine chmriv7(PM,TEMPER,PRESUR,RH,O3PUFF,RK1,ZCOEF,                 &
     &                   TSTP,CONFCT,RSHET)
!----------------------------------------------------------------------
!
! --- CALPUFF    Version: TNG-7.1.0    Level: 100917           CHMRIV7
!                P. Karamchandani, AER (Adapted from chemriv)
!
! --- ADOPTED FROM ARM3 (see banner below)
!                Sulfate & nitrate conversion uses exponential fcn
!                Unused variables are deactivated
!
! --- UPDATE
! --- V5.8-V6.302  100917 (DGS): Restructure value**-n to value**(-n)
!                                to satisfy LF95 compiler
!                         (DGS): Add heterogeneous reaction rate for
!                                SO4 and pass RSHET(fraction/hr)
! --- V5.8         071025  (PK): Original
!
!----------------------------------------------------------------------
      IMPLICIT NONE
      SAVE
!
!DECK.CHMRIV
!
!  DATE: NOVEMBER 1987
!  VERSION: ARM3-1.0
!
!     CALCULATE TRANSFORMATION OF SO2 TO SO4 AND NO2 TO HNO3-NO3
!     FOR ONE TIME STEP, TSTP
!     BASED ON RIVAD CHEMICAL MECHANISM
!     ADPATED FROM THE RIVAD MODEL
!
!
!  INPUT ARGUMENTS:
!                   PM     R  MASS OF SPECIES (I) IN PUFF (G)
!                             SO2, SO4, NO, NO2, HNO3-NO3, TSP
!                   TEMPER R  TEMPERATURE AT PLUME HEIGHT (K)
!                   PRESUR R  PRESSURE AT PLUME HEIGHT (MB)
!                   RH     R  RELATIVE HUMIDITY AT PLUME (%)
!                   O3PUFF R  PUFF OZONE CONCENTRATION (PPM)
!                   RK1    R  NO2 PHOTOLYSIS RATE CONSTANT (PPM/MIN)
!                   ZCOEF  R  COSINE OF SOLAR ZENITH ANGLE
!                   TSTP   R  TIME STEP (HOURS)
!                   CONFCT R  CONVERSION FACTOR FOR UG/M3 TO PPM-MWT
!                   RSHET  R  SO2 heterogeneous loss rate (fraction/hr)
!
!  OUTPUT ARGUMENTS:
!                   PM     R  NEW MASS OF PUFF DUE TO CHANGES BASED ON
!
!
!  SUBROUTINES CALLED:
!
!  CALLED BY:   CHEMRIV7
!

! --- Constants
      REAL COEF1                  ! Molec/cc to ppm conv factor coefficient
      PARAMETER ( COEF1 = 7.33981E+15 )

      REAL CONSTC                 ! Constant for falloff type reaction
      PARAMETER ( CONSTC = 0.6 )

      REAL TI300                  ! 1.0 / 300.
      PARAMETER ( TI300 = 1.0 / 300.0 )

! --- Arguments
      REAL PM(25)
      REAL TEMPER, PRESUR, RH, O3PUFF, RK1, ZCOEF, TSTP, CONFCT, RSHET

! --- Locals

! --- wangzhm save
      REAL,save :: TFACT, R26, ROHM, QJ, PHIKK, RCONST
      INTEGER,save :: L
      REAL,save :: CNO, CNO2, CNOX, TAMB, PAMB, RH1, H2O, SUM
      REAL,save :: XNO, XNOX, XO3, XOX, XNO2, RNO2X, XNO3, XN2O5
      REAL,save :: RSULF, RNITR, RKX, RNO3, RNITRN, RNITRD, SULFN
      REAL,save :: XOHMAX, XOH
      REAL,save :: A0, A, B, C

      REAL,save :: RK0               ! K0 in falloff rate expressions
      REAL,save :: RKINF             ! KINF in falloff rate expressions
      REAL,save :: XEND              ! Exponent in falloff rate expressions
!
      REAL,save :: KSO2OH            ! SO2 + OH rate constant
      REAL,save :: KNO2OH            ! NO2 + OH rate constant
      REAL,save :: KTOLOH            ! TOL + OH rate constant
      REAL,save :: KXYLOH            ! XYL + OH rate constant
      REAL,save :: KALKHOH           ! ALKH + OH rate constant
      REAL,save :: KPAHOH            ! PAH + OH rate constant
!
      REAL,save :: DTOL, DXYL, DALKH, DPAH
!
      REAL*8,save :: TINV, CFACT, RFACT
!
      REAL,save :: PPM(25),SMWT(25)
      LOGICAL,save :: LFIRST
      DATA LFIRST/.TRUE./
      DATA SMWT/64.0,96.0,30.0,46.0,63.0,62.0,180.0,                           &
     &          5*92.0,5*106.0,3*226.0,5*156.0/

! --- wangzhm omp: set threadsprivate
      !$OMP THREADPRIVATE(TFACT,R26,ROHM,QJ,PHIKK,RCONST,L,&
      !$OMP& CNO,CNO2,CNOX,TAMB,PAMB,RH1,H2O,SUM,XNO,XNOX,XO3,XOX,XNO2,RNO2X,XNO3,XN2O5,&
      !$OMP& RSULF,RNITR,RKX,RNO3,RNITRN,RNITRD,SULFN,XOHMAX,XOH,A0,A,B,C,RK0,RKINF,XEND,&
      !$OMP& KSO2OH,KNO2OH,KTOLOH,KXYLOH,KALKHOH,KPAHOH,DTOL,DXYL,DALKH,DPAH,TINV,CFACT,RFACT,&
      !$OMP& PPM,SMWT,LFIRST)
!
!     FIRST TIME THROUGH SET UP SOME GLOBAL RATE CONSTANTS
!
      IF (LFIRST) THEN
         TFACT = 1./273.
         R26 = 1./26.4
         ROHM = 4.87E-7/1.32E-3
         LFIRST = .FALSE.
      END IF
!
      DO L = 1,5
         PPM(L) = 1.E6*CONFCT*PM(L)/SMWT(L)
      END  DO
!
      DO L = 8,25
         PPM(L) = 1.E6*CONFCT*PM(L)/SMWT(L)
      END  DO
!
      CNO = PPM(3)
      CNO2 = PPM(4)
      CNOX = CNO + CNO2
!
      TAMB = TEMPER
      PAMB = PRESUR/1013.0
      TINV = 1. / TAMB
      CFACT = COEF1 * PAMB * TINV
!
      IF (RK1.LE.0.0) THEN
         QJ = 0.0
         PHIKK = 0.0
      ELSE
         PHIKK = RK1*R26
         QJ = (1.338E-3)*ZCOEF**2.74
      ENDIF
      IF (TAMB.GE.273.) THEN
         RCONST = 18.02*(597.3-.566*(TAMB-273.))/1.9869
      ELSE
         RCONST = 6133.17
      END IF
      RH1 = MIN(RH,95.)
      RH1 = MAX(RH1,0.)
      H2O = (6030.*.01*RH1/PAMB)*EXP(RCONST*(TFACT-TINV))
!
!  DO SIMPLE CHEMISTRY
!
      XNOX = CNOX
      XOX = O3PUFF + CNO2
      SUM = XOX + XNOX + PHIKK
      XNO2 = 0.5*(SUM-SQRT(ABS(SUM*SUM-4.*XNOX*XOX)))
      XNO2 = MIN(XNO2,XNOX)
      XNO2 = MAX(XNO2,0.)
! dgs      XNO=XNOX-XNO2
! dgs      IF (XNO.LT.0.0) XNO=0.0
      RNO2X = 0.
      IF (XNOX.GT.0.) RNO2X = XNO2/XNOX

! dgs Use NO2/NOX ratio to test for zero NO (single precision)
      xno = xnox - xno2
      if(rno2x.GE.0.999999) XNO = 0.0

      XO3 = XOX - XNO2
      XO3 = MAX(XO3,0.)
      IF (ZCOEF.LT.0.06975) THEN
!
!  NIGHTTIME CHEMISTRY
!
         RSULF = 0.
!
!  IMPLEMENT NEW NIGHTTIME CHEMISTRY - 1/86
!
         RNITR = 0.0
         XO3 = MAX(XO3,0.)
         IF (XO3.GT.0.) THEN
            RKX = 1780./(1.9E-6*H2O + 3.12)
            RNO3 = 0.086
            A0 = 0.59 + 1780. -3.12*RKX
            A = 2.*RKX*RNO3-A0
            B = RNO3 + 0.0474*XO3 +XNO2*A0
            C = -0.0474*XNO2*XO3
!
            XNO3 = -2.*C/(B+SQRT(B*B-4.*A*C))
            XN2O5 = XNO3*XNO2*RKX
            RNITRN = 2.*1.9E-6*H2O*XN2O5*60.*TSTP
            RNITRN = MIN(XNO2,RNITRN)
            PPM(5) = PPM(5) + RNITRN
            XNO2 = XNO2 - RNITRN
         END IF
      ELSE
!
!  DO DAYTIME CHEMISTRY
!
         XOHMAX = ROHM*QJ
         XOH = XOHMAX
         IF (PPM(1).NE.0..OR.XNO2.NE.0.)                                       &
     &   XOH = 2.*QJ*3.4E5*H2O*XO3/((4.45E10+3.4E5*H2O)*                       &
     &     (2000.*PPM(1) + 14000.*XNO2))
         XOH = MIN(XOH,XOHMAX)

! --- SO2 + OH rate constant (falloff expression)
         RK0 = 1.0E+06 * CFACT * 3.0E-31 * ( TAMB * TI300 )**(-3.3)
         RKINF = 1.5E-12
         XEND = 1.0 / ( 1.0 + ( LOG10( RK0 / RKINF ) )**2 )
         KSO2OH = ( RK0 / ( 1.0 + RK0 / RKINF ) ) * CONSTC**XEND
!
! --- NO2 + OH rate constant (falloff expression)
         RK0 = 1.0E+06 * CFACT * 2.6E-30 * ( TAMB * TI300 )**(-3.2)
         RKINF = 2.4E-11 * ( TAMB * TI300 )**(-1.3)
         XEND = 1.0 / ( 1.0 + ( LOG10( RK0 / RKINF ) )**2 )
         KNO2OH = ( RK0 / ( 1.0 + RK0 / RKINF ) ) * CONSTC**XEND
!
! --- TOL + OH rate constant
         ktoloh = 2.1e-12 * EXP(322.*tinv)
!
! --- XYL + OH rate constant
         kxyloh = 1.7e-11 * EXP(116.*tinv)
!
! --- ALKH + OH rate constant
         kalkhoh = 1.97e-11
!
! --- PAH + OH rate constant
         kpahoh = 7.7e-11
!
! --- Convert rate constants from molec-cc-1 s-1 to ppm-1 hr-1 units
         rfact = 2.64e19 * (pamb * tinv)
         kso2oh = kso2oh * rfact
         kno2oh = kno2oh * rfact
         rsulf = 1. - EXP(-kso2oh*xoh*tstp)
         rnitr = 1. - EXP(-kno2oh*xoh*tstp)
         RNITRD = RNITR*XNO2
         RNITRD = MIN(XNOX,RNITRD)
         PPM(5) = PPM(5) + RNITRD
         XNOX = XNOX - RNITRD
         XNOX = MAX(XNOX,0.)
         XNO = XNOX*(1.-RNO2X)
         XNO2 = XNOX*RNO2X

         ktoloh = ktoloh * rfact
         kxyloh = kxyloh * rfact
         kalkhoh = kalkhoh * rfact
         kpahoh = kpahoh * rfact
!
! --- Change in SOA precursor concentrations
         dtol = (1. - EXP(-ktoloh*xoh*tstp)) * ppm(8)
         dxyl = (1. - EXP(-kxyloh*xoh*tstp)) * ppm(13)
         dalkh = (1. - EXP(-kalkhoh*xoh*tstp)) * ppm(18)
         dpah = (1. - EXP(-kpahoh*xoh*tstp)) * ppm(21)

         ppm(8) = ppm(8) - dtol          ! TOL
         ppm(9) = ppm(9) + 0.071*dtol    ! TOLAER1
         ppm(10) = ppm(10) + 0.138*dtol  ! TOLAER2

         ppm(13) = ppm(13) - dxyl        ! XYL
         ppm(14) = ppm(14) + 0.038*dxyl  ! XYLAER1
         ppm(15) = ppm(15) + 0.167*dxyl  ! XYLAER2

         ppm(18) = ppm(18) - dalkh       ! ALKH
         ppm(19) = ppm(19) + 1.173*dalkh ! ALKHAER

         ppm(21) = ppm(21) - dpah        ! PAH
         ppm(22) = ppm(22) + 0.156*dpah  ! PAHAER1
         ppm(23) = ppm(23) + 0.777*dpah  ! PAHAER2
      END IF

! --- Add heterogeneous rate wso4=rshet*ppm(1)*tstp
      SULFN = PPM(1)*RSULF + rshet*ppm(1)*tstp
      SULFN = MIN(PPM(1),SULFN)
      PPM(2) = PPM(2) + SULFN
      PPM(1) = PPM(1) - SULFN
!
!     UPDATE NEW STEADY-STATE NO AND NO2 VALUES
!
      PPM(3) = XNO
      PPM(4) = XNO2
!
!     UPDATE NEW O3 VALUE
      O3PUFF = XO3
!
!     CONVERT BACK TO GRAMS OF SPECIES IN PUFF
!
      DO L = 1,5
         PM(L) = 1.E-6*PPM(L)*SMWT(L)/CONFCT
      END DO
!
      DO L = 8,25
         PM(L) = 1.E-6*PPM(L)*SMWT(L)/CONFCT
      END DO
!
      RETURN
      END

!---------------------------------------------------------------------
      subroutine setbckoc(ndathr)
!---------------------------------------------------------------------
!
! --- CALPUFF    Version: TNG-7.1.0    Level: 071025          SETBCKOC
!                P. Karamchandani, AER (Adapted from setsoa)
!
! --- PURPOSE:  Background OC concentration (ug/m3) for current hour
!
! --- INPUTS:
!            NDATHR - integer    - YYJJJHH date-time for hour
!
!     Common Block /CHEMDAT/ variables:
!        BCKPMF(12),OFRAC(12)
!
! --- OUTPUT:
!
!     Common Block /NEWSOA/ variables: bckoc
!
! --- SETBCKOC called by:  COMP
! --- SETBCKOC calls:      GRDAY
!----------------------------------------------------------------------
      include 'params.puf'
      include 'chemdat.puf'
      include 'newsoa.puf'

! --- Extract month from date-time
      iyr=ndathr/100000
      ijul=ndathr/100 - 1000*iyr
      call GRDAY(io6,iyr,ijul,imo,iday)

! --- Set background OC data for this month
      bckoc=bckpmf(imo)*ofrac(imo)

      return
      end

!----------------------------------------------------------------------
      subroutine isodriver(so4,no3,nh3,nrh,ntempk,pno3c)
!----------------------------------------------------------------------
!
! --- CALPUFF    Version: TNG-7.1.0    Level: 140521         ISODRIVER
!                P. Karamchandani, AER
!
! --- PURPOSE:  Driver routine to calculate gas/particle equilibrium
!               using ISORROPIA aerosol equilibrium module
!               This version uses double precision
!
! --- UPDATES:
!
! --- V6.4-V6.42_x1.1  140521  : Call ISOROPIA with METASTABLE ON
!
! --- V6.302-V6.4  101025 (DGS): Revise ISOROPIA arguments for V2.1
!
! --- INPUTS:
!            so4    - real - Total sulfate concentration (mole/m3)
!            no3    - real - Total nitrate concentration (mole/m3)
!            nh3    - real - Total ammonia concentration (mole/m3)
!            nrh    - real - Fractional relative humidity
!            ntempk - real - Temperature (K)
!
! --- OUTPUT:
!            pno3c  - real - Particle-phase equilibrium NO3
!                            concentration (mole/m3)
!
! --- ISODRIVER called by:  CHEMRIV6, CHEMRIV7
! --- ISODRIVER calls:      ISOROPIA
!----------------------------------------------------------------------

!***********************************************************************
!                                                                      *
!  REVISION HISTORY:                                                   *
!     Version 1.0 developed November 2006 by PK, AER, Inc. for         *
!     CALPUFF for API Contract No. 2006-102376                         *
!                                                                      *
!***********************************************************************
!
!...........  INCLUDES

      INCLUDE 'isrpia.inc'

!...........  ARGUMENTS and their descriptions

      REAL SO4             ! Total sulfate (moles/m3)
      REAL NO3             ! Total nitrate (moles/m3)
      REAL NH3             ! Total ammonia (moles/m3)
      REAL NRH             ! Relative humidity as fraction
      REAL NTEMPK          ! Temperature in Kelvin

      REAL PNO3C           ! Particle-phase NO3 concentration (mole/m**3)

! --- Locals

! --- wangzhm save
      REAL,save :: GNO3C           ! Gas-phase concentration of HNO3 in mole/m**3 air

      INTEGER NCTRL, NOTHER
!V1.7 PARAMETER(NCTRL = 2,NOTHER = 6)
!V2.1 PARAMETER(NCTRL = 2,NOTHER = 9)
      PARAMETER(NCTRL = 2,NOTHER = 9)

! --- Gas-phase concentration array in moles/m**3 air
      REAL*8,save ::  GAS(NGASAQ)

! --- Aqueous-phase concentration array in moles/m**3 air
      REAL*8,save ::  AERLIQ(NIONS+NGASAQ+2)

! --- Solid-phase concentration array in moles/m**3 air
      REAL*8,save ::  AERSLD(NSLDS)

! --- Flag for different types of problems solved and
! --- different state of aerosols (deliquescent or metastable)
      REAL*8,save ::  CNTRL(NCTRL)
                                                                               &
! Sol&tion information array (see ISOCOM.f for details)
      REAL*8,save ::  OTHER(NOTHER)
                                                                               &
! Tot&l species concentrations in moles/m**3 air
      REAL*8,save ::  WI(NCOMP)
                                                                               &
! Tem&erature and humidity
      REAL*8,save ::  TEMPI, RHI

! --- wangzhm omp: set threadsprivate
      !$OMP THREADPRIVATE(GNO3C,GAS,AERLIQ,AERSLD,CNTRL,OTHER,WI,TEMPI,RHI)

!***********************************************************************
!  begin body of subroutine
!
! *** Assign input total concentrations to WI

      WI(1) = 0.    ! Na concentration
      WI(2) = so4
      WI(3) = nh3
      WI(4) = no3
      WI(5) = 0.    ! Cl concentration
!V2.1 Additional 3 are not used here
      WI(6) = 0.    ! total calcium concentration
      WI(7) = 0.    ! total potassium concentration
      WI(8) = 0.    ! total magnesium concentration

      RHI   = NRH
      TEMPI = NTEMPK
!
! *** CALL ISORROPIA
!
      CNTRL(1) = 0.     ! 0 = FORWARD PROBLEM,  1 = REVERSE PROBLEM

! --- 6.42_x1.1
      CNTRL(2) = 1.     ! 0 = SOLID + LIQUID AEROSOL,  1 = METASTABLE
! --- wangzhm omp: critical
      !$OMP CRITICAL (call_isoropia)
      CALL ISOROPIA (WI,  RHI,  TEMPI,   CNTRL,                                &
     &               W,   GAS,  AERLIQ,  AERSLD,  SCASE,  OTHER)
      !$OMP END CRITICAL (call_isoropia)

!
! *** SAVE RESULTS
                                                                               &
! Gas&phase HNO3 (moles/m3)
      GNO3C  = GAS(2)
                                                                               &
! Par&icle-phase NO3 (moles/m3)
      PNO3C  = MAX(0.,NO3 - GNO3C)
!
      RETURN
      END

!----------------------------------------------------------------------
      subroutine soadriver(ilst,q,tempk,press,pivol)
!----------------------------------------------------------------------

!***********************************************************************
!  FUNCTION: Driver routine to calculate SOA gas-particle equilibrium  *
!                                                                      *
!
! --- UPDATES:
!
! --- V6.4-V6.42_x1.1  140521  : Iterate on CALTECH_SOA call to update
!                                absorbing organic mass
!
!  INPUTS:                                                             *
!     ilst   - List-file unit                                          *
!     q      - Mass of species in puff                                 *
!     tempk  - Temperature (K)                                         *
!     press  - Pressure at plume height (mb)                           *
!     pivol  - Reciprocal of puff volume (1/m**3)                      *
!                                                                      *
!  RETURN VALUES:                                                      *
!     q      - Adjusted mass of species in puff                        *
!                                                                      *
!  REVISION HISTORY:                                                   *
!     Version 1.0 developed November 2006 by PK, AER, Inc. for         *
!     CALPUFF for API Contract No. 2006-102376                         *
!***********************************************************************
      IMPLICIT NONE
!
!...........  INCLUDES

      INCLUDE 'newsoa.puf'
      INCLUDE 'soadat.puf'

!...........  ARGUMENTS and their descriptions

! --- 6.42_x1.1
      integer ilst         ! List-file unit

      REAL q(25)           ! Puff masses (g)
      REAL tempk           ! Temperature in Kelvin
      REAL press           ! Pressure in mb
      REAL pivol           ! reciprocal of puff volume (1/m3)

! --- wangzhm save
                                                                               &
! Loc&l variables
      integer,save :: iorg ! loop index                                                &
! Sec&ndary organic aerosol concentrations (ug/m3)
! --- Gas-phase compounds
      REAL,save :: gasorg(NORG)
! --- Particle-phase compounds
      REAL,save :: partorg(NORG)
! --- Total (gas-phase + particle-phase)
      REAL,save :: worg(NORG+1)  ! additional species at end for primary OC

! --- Conversion factor
      REAL,save :: cfact

! --- 6.42_x1.1
! --- OA mass iteration
      integer,save :: iter, niter
      real,save :: psum1, psum2, pdiff

! --- wangzhm omp: set threadsprivate
      !$OMP THREADPRIVATE(iorg,gasorg,partorg,worg,cfact,iter,niter,psum1,psum2,pdiff)

!***********************************************************************
!  begin body of subroutine
!
! *** Convert input puff masses (g) to puff concentrations (ug/m3)

      cfact = pivol * 1.E6
                                                                               &
! Gas&s
      gasorg(1) = q(9) * cfact    ! TOLAER1
      gasorg(2) = q(10) * cfact   ! TOLAER2
      gasorg(3) = q(14) * cfact   ! XYLAER1
      gasorg(4) = q(15) * cfact   ! XYLAER2
      gasorg(5) = q(19) * cfact   ! ALKHAER
      gasorg(6) = q(22) * cfact   ! PAHAER1
      gasorg(7) = q(23) * cfact   ! PAHAER2
                                                                               &
! Par&icles
      partorg(1) = q(11) * cfact   ! ATOLA1
      partorg(2) = q(12) * cfact   ! ATOLA2
      partorg(3) = q(16) * cfact   ! AXYLA1
      partorg(4) = q(17) * cfact   ! AXYLA2
      partorg(5) = q(20) * cfact   ! AALKHA
      partorg(6) = q(24) * cfact   ! APAHA1
      partorg(7) = q(25) * cfact   ! APAHA2
                                                                               &
! Tot&l
      do iorg = 1, NORG
        worg(iorg) = gasorg(iorg) + partorg(iorg)
      end do
                                                                               &
! Org&nic absorbing mass
      worg( NORG + 1 ) = MAX(q(7)*cfact + bckoc, 0.)
                                                                               &
! Cal&ulate organic aerosol formation using Caltech SOA Module

! --- 6.42_x1.1
! --- Iterate on this call until difference in particulate OC total <1%
! --- PDIFF is 1% criterion for difference
! --- NITER is set large but finite to trap endless loop
! --- ITER is current loop counter
      pdiff=0.01
      niter=10000
      iter=0
10    psum1=0.0
      do iorg = 1, NORG
        psum1=psum1+partorg(iorg)
      end do

      call caltech_soa(worg, gasorg, partorg, tempk)

! --- 6.42_x1.1
      iter=iter+1
      psum2=-psum1
      do iorg = 1, NORG
        psum2=psum2+partorg(iorg)
      end do
      if(psum1.GT.0.0) then
         if(iter.LE.niter) then
            if(ABS(psum2/psum1) .GT. pdiff) goto 10
         else
! ---       Report too many iterations warning and stop iterating
            write(ilst,*)'SOADRIVER:  Warning ...'
            write(ilst,*)'  Iterations exceed NITER = ',niter,                 &
     &                   ' when computing particulate SOA '
            write(ilst,*)'  Target fractional SOA change = ',pdiff
            write(ilst,*)'  Last fractional SOA change = ',psum2/psum1
         endif
      elseif(worg(norg+1).GT.0.0 .AND. iter.EQ.1) then
         goto 10
      endif

! --- Assign adjusted concentrations back to puff masses

      cfact = 1./cfact
! Gases
      q(9) = gasorg(1) * cfact    ! TOLAER1
      q(10) = gasorg(2) * cfact   ! TOLAER2
      q(14) = gasorg(3) * cfact   ! XYLAER1
      q(15) = gasorg(4) * cfact   ! XYLAER2
      q(19) = gasorg(5) * cfact   ! ALKHAER
      q(22) = gasorg(6) * cfact   ! PAHAER1
      q(23) = gasorg(7) * cfact   ! PAHAER2

! Particles
      q(11) = partorg(1) * cfact   ! ATOLA1
      q(12) = partorg(2) * cfact   ! ATOLA2
      q(16) = partorg(3) * cfact   ! AXYLA1
      q(17) = partorg(4) * cfact   ! AXYLA2
      q(20) = partorg(5) * cfact   ! AALKHA
      q(24) = partorg(6) * cfact   ! APAHA1
      q(25) = partorg(7) * cfact   ! APAHA2
!
      RETURN
      END

!----------------------------------------------------------------------
      SUBROUTINE CALTECH_SOA (WORG, GASORG, PARTORG, CURTEMP)
!----------------------------------------------------------------------

!***********************************************************************
!  FUNCTION: Program to simulate Secondary Organic Aerosols            *
!     Total no. of species = 4 + 34 = 38                               *
!     Partition equation                                               *
!     Kom, i = (Ai/Mo)/Gi                                              *
!     Ai = particle-phase concentration (ug/m**3 air)                  *
!     Gi = gas-phase concentration (ug/m**3 air)                       *
!     MSUM = sum Ai + primary organics                                 *
!     Kom, i has the units of m3/ug                                    *
!  PRECONDITION REQUIRED: called from subr. AEROEQ                     *
!  RETURN VALUES:                                                      *
!     PARTORG(I) - Particulate concentration of organic species i      *
!     GASORG(I)  - Gas-phase concentration of organic species i        *
!  REVISION HISTORY:                                                   *
!     Written by Betty K. Pun of AER, Inc. for EPRI's Aerosol          *
!         Module Implementation Project, May, 1999                     *
!     Revised by Yang Zhang of AER, Inc. for EPRI's Aerosol            *
!         Module Implementation Project based on Models3's             *
!         coding standard July, 1999                                   *
!     Revised by Betty K. Pun of AER, Inc. November, 99                *
!         for incorporation into 3-D model under EPRI.                 *
!         To reduce computational requirement and to                   *
!         facilitate numerical solution, simultaneous equations        *
!         are not solved in this version.  Instead, the partition of   *
!         each organic aerosol at each time depends on the amount of   *
!         material in the organic phase at the step before the partition*
!     Modified by Betty K. Pun, AER to include temperature dependence  *
!         based on generic Hvap, February/March 2002                   *
!     Revised by BKP March 2006, change reference temperature from     *
!             310 to 298K                                              *
!     Updated by PK December 2006, for implementation in CALPUFF       *
!             for API Contract No. 2006-102376                         *
!  REFERENCES:                                                         *
!    1. Odum et al.,  97.  EST 31:1890                                 *
!    2. Griffin et al.,  99.  JGR 104:3555                             *
!                                                                      *
!***********************************************************************

      IMPLICIT NONE

!...........  INCLUDES
      INCLUDE 'soadat.puf'

!...........  ARGUMENTS and their descriptions and some other variables

      REAL WORG(NORG+1)   ! Total organic species concentration
      REAL PARTORG(NORG)  ! Particulate concentration of organic species i
      REAL GASORG(NORG)   ! Gas-phase concentration of organic species i
      REAL CURTEMP        ! Puff temperature used in modifying parititon
                          ! constants

!..........  Local Variables

      REAL MP             ! Total organic species concentration
      REAL MSUM           ! Sum of particle-phase concentration (ug/m**3 air)
      INTEGER I           ! Organic species index

!***********************************************************************
!  begin body of subroutine

      MP = WORG(NORG + 1)              ! OC (primary+background)
      MSUM = MP

      DO I = 1,  NORG
         IF (WORG(I) .LT. 0.0) THEN
            WRITE (*,  95000) I
95000       FORMAT ('ERROR: Negative conc. read for species',  I2)
            WRITE (*,  *) 'Negative Conc. = ',WORG(I)
            call flush(6)
            STOP
         END IF
         MSUM = MSUM + PARTORG(I)      !add starting SOA to MSUM
      END DO

      IF (MSUM .GT. 0.0) THEN
         call SOASUB (MSUM, WORG, PARTORG, GASORG, CURTEMP)
      ENDIF

      DO I = 1,  NORG
         IF (PARTORG(I) .LT. 0.0) THEN
            PARTORG(I) = 0.0
            GASORG(I) = WORG(I)
         END IF
      END DO

      RETURN
      END

!----------------------------------------------------------------------
      SUBROUTINE SOASUB (MSUM, WORG, PARTORG, GASORG, CURTEMP)
!----------------------------------------------------------------------

!***********************************************************************
!  FUNCTION: Program to calculate the equilibria                       *
!  PRECONDITION REQUIRED: called from subr. CALTECHSOA                 *
!  RETURN VALUES:                                                      *
!     PARTORG(J) - Particule-phase concentrations (microgram/m**3)     *
!     GASORG(J) - Gas-phase concentrations (microgram/m**3)            *
!  REVISION HISTORY:                                                   *
!     Written by Betty K. Pun of AER, Inc. for EPRI's Aerosol          *
!         Module Implementation Project, May, 1999                     *
!     Revised by Yang Zhang of AER, Inc. for EPRI's Aerosol            *
!         Module Implementation Project based on Models3's             *
!         coding standard July, 1999                                   *
!     Modified by Betty K. Pun of AER, Inc, for Models-3 application   *
!         Analytical solution calculated November, 99                  *
!     Modified by Betty K. Pun, AER to include temperature dependence  *
!         based on generic Hvap, February/March 2002                   *
!***********************************************************************

      IMPLICIT NONE

!...........  INCLUDES

      INCLUDE 'soadat.puf'

!...........  ARGUMENTS and their descriptions

      REAL MSUM           ! Sum of particle-phase concentration (ug/m**3 air)
      REAL WORG(NORG+1)   ! Total organic species concentration
      REAL PARTORG(NORG)  ! Particulate concentration of organic species i
      REAL GASORG(NORG)   ! Gas-phase concentration of organic species i
      REAL CURTEMP        ! Grid cell temperature used in modifying parititon
                          ! constants
!...........  Local variables
      INTEGER J           ! Organic species index
      REAL KCORR(NORG)    ! Partition coefficients, corrected for temperature
                          ! [m**3/ug]
      REAL TEXPT          ! temperature where KOM (experimental values)
                          ! are obtained
      PARAMETER (TEXPT = 298.0)

! ..........  Other Variables

!***********************************************************************
!  begin body of subroutine

      IF (MSUM .EQ. 0.0) THEN
         WRITE(*,  *)    'Error in SOA: MSUM = 0.0'
         STOP
      END IF
!
! *** Calculate the equilibrium concentration
!     GASORG(J) = WORG(J) - PARTORG(J)
!     G         = C      -  A
!     (A/M)/G = K is equal to A/(C-A) = KM is equal to A = CKM - AKM
!     Therefore A = CKM / (1 + KM)

! add code to correct for CURTEMPerature, assuming KOM is experimental
! values obtained at 310 K (changed to a ref temp of 298K)
      DO J = 1, NORG
         KCORR(J) = EXP(HVAP(J)*(1/CURTEMP - 1/TEXPT)/RKJMOLK)                 &
     &                 * KOM(J) * CURTEMP/TEXPT
         PARTORG(J) = WORG(J) * KCORR(J)* MSUM /                               &
     &                    (1. + MSUM * KCORR(J))
         GASORG(J) = WORG (J) - PARTORG(J)
      END DO

      RETURN
      END

!----------------------------------------------------------------------
      subroutine aqradm ( temp, pres_atm, taucld, prcrate, wcavg, coz,         &
     &                    ch2o2, ctnh3, conc, rhoair, len, scav )
!----------------------------------------------------------------------
!
!  DESCRIPTION:
!    Compute concentration changes in cloud due to aqueous chemistry
!    Adapted from RADM Cloud implementation in CMAQ/SCICHEM for CALPUFF
!    by PK, AER, January 2007 for API Contract No. 2006-102376
!
!  Reference:
!     Walcek & Taylor, 1986, A theoretical Method for computing
!      vertical distributions of acidity and sulfate within cumulus
!      clouds, J. Atmos Sci.,  Vol. 43, no. 4 pp 339 - 355
!
!  Called by:  CHEMRIV6, CHEMRIV7
!
!  Calls the following functions:  HLCONST
!
!  ARGUMENTS     TYPE      I/O       DESCRIPTION
!  ---------     ----  ------------  --------------------------------
!   CONC(6)      real  input&output  Concentration for species i=1,11
!                                    (1) = SO2   conc (mol/mol of SO2)
!                                    (2) = SO4   conc (mol/mol of SO4)
!                                    (3) = NO    conc (mol/mol of NO)
!                                    (4) = NO2   conc (mol/mol of NO2)
!                                    (5) = HNO3  conc (mol/mol of HNO3)
!                                    (6) = NO3   conc (mol/mol of NO3)
!-----------------------------------------------------------------------

      IMPLICIT NONE

!...........PARAMETERS and their descriptions:
!      include 'params.puf'                                                    &
! tem&orarily set io6=2 here since params.puf variables are not declared,      &
! res&lting in compile errors.
      INTEGER IO6
      PARAMETER (IO6 = 2) ! from params.puf
                                                                               &
! num&er of oxidizing reactions
      INTEGER NUMOX
      PARAMETER (NUMOX =  3)   ! H2O2, O3, and Fe-Mn catalyzed
                                                                               &
! min&mum and maximum pH
      REAL PHMIN, PHMAX
      PARAMETER (PHMIN = 0.0001, PHMAX = 10.0)
                                                                               &
! min&mum concentration
      REAL CONCMIN
      PARAMETER (CONCMIN = 1.0E-30)
                                                                               &
! con&ert seconds to hours
      REAL SEC2HR
      PARAMETER (SEC2HR = 1.0 / 3600.0)
!                                                                              &
! Mol&r volume at STP [ L/mol ] Non MKS units
      REAL MOLVOL
      PARAMETER (MOLVOL = 22.41410)
                                                                               &
! Sta&dard Temperature [ K ]
      REAL STDTEMP
      PARAMETER (STDTEMP = 273.15)
                                                                               &
! den&ity of water at 20 C and 1 ATM (kg/m3)
      REAL H2ODENS
      PARAMETER (H2ODENS = 1000.0)
                                                                               &
! Mol&cular weights of Fe and Mn
      REAL MW_FE, MW_MN
      PARAMETER (MW_FE = 55.8, MW_MN = 54.9)

! CO2 background concentration (moles/mole air)
      REAL CCO2
      PARAMETER (CCO2 = 3.4E-04)   ! 340 ppm

!...........ARGUMENTS and their descriptions

      REAL         TEMP            ! temperature (K)
      REAL         PRES_ATM        ! pressure (atm)
      REAL         TAUCLD          ! timestep for cloud (s)
      REAL         PRCRATE         ! precip rate (mm/hr)
      REAL         WCAVG           ! liquid water content (kg/m3)
      REAL         COZ             ! ozone concentration (mol/mol)
      REAL         CH2O2           ! H2O2 concentration (mol/mol)
      REAL         CTNH3           ! total ammonium  concentration (mol/mol)
      REAL         CONC   ( 6 )    ! species concentrations (mol/molV)
      real         rhoair          ! air density, moles/m3

! --- puff length scale (for scavenging coefficient calculations) (m)
      real len

! --- Scavenging coefficients (1/s)
      real scav( 6 )

! --- wangzhm save

!...........LOCAL VARIABLES (scalars) and their descriptions:

      REAL,save ::       RT              ! gas const * temperature (liter atm/mol)
      REAL,save ::       RECIPAP1        ! one over pressure (/atm)
      REAL,save ::       ONE_OVER_TEMP   ! 1.0 / TEMP

! --- Gas liquid equilibria

      REAL,save ::         PH2O20          ! total H2O2 partial pressure (atm)
      REAL,save ::         PH2O2F          ! gas only H2O2 partial pressure (atm)
      REAL,save ::         H2O2H           ! Henry's Law Constant for H2O2
      REAL,save ::         H2O2L           ! H2O2 conc in cloudwater (mol/liter)

      REAL,save ::         PHNO30          ! total HNO3 partial pressure (atm)
      REAL,save ::         PHNO3F          ! gas only HNO3 partial pressure (atm)
      REAL,save ::         HNO3H           ! Henry's Law Constant for HNO3
      REAL,save ::         HNO31           ! First dissociation constant for HNO3
      REAL,save ::         HNO31H          ! HNO31*HNO3H
      REAL,save ::         HNO3L           ! HNO3 conc in cloudwater (mol/liter)

      REAL,save ::         PNH30           ! total NH3 partial pressure (atm)
      REAL,save ::         PNH3F           ! gas only NH3 partial pressure (atm)
      REAL,save ::         NH3H            ! Henry's Law Constant for NH3
      REAL,save ::         NH31            ! First dissociation constant for NH3
      REAL,save ::         NH3DH2O         !
      REAL,save ::         NH31HDH         !
      REAL,save ::         NH3L            ! NH3 conc in cloudwater (mol/liter)

      REAL,save ::         PO30            ! total O3 partial pressure (atm)
      REAL,save ::         PO3F            ! gas only O3 partial pressure (atm)
      REAL,save ::         O3H             ! Henry's Law Constant for O3
      REAL,save ::         O3L             ! O3 conc in cloudwater (mol/liter)

      REAL,save ::         PSO20           ! total SO2 partial pressure (atm)
      REAL,save ::         PSO2F           ! gas only SO2 partial pressure (atm)
      REAL,save ::         SO2H            ! Henry's Law Constant for SO2
      REAL,save ::         SO21            ! First dissociation constant for SO2
      REAL,save ::         SO22            ! Second dissociation constant for SO2
      REAL,save ::         SO212           ! SO21*SO22
      REAL,save ::         SO212H          ! SO21*SO22*SO2H
      REAL,save ::         SO21H           ! SO21*SO2H
      REAL,save ::         SO2L            ! SO2 conc in cloudwater (mol/liter)
      REAL,save ::         SO3             ! SO3= conc in cloudwater (mol/liter)

      REAL,save ::         PNO2            ! total NO2 partial pressure (atm)
      REAL,save ::         NO2H            ! Henry's Law Constant for NO2 (M/atm)
      REAL,save ::         NO2L            ! NO2 conc in cloudwater (mol/liter)

      REAL,save ::         PNO             ! total NO partial pressure (atm)
      REAL,save ::         NOH             ! Henry's Law Constant for NO (M/atm)
      REAL,save ::         NOL             ! NO conc in cloudwater (mol/liter)

      REAL,save ::         PCO20           ! total CO2 partial pressure (atm)
      REAL,save ::         PCO2F           ! gas only CO2 partial pressure (atm)
      REAL,save ::         CO2H            ! Henry's Law constant for CO2
      REAL,save ::         CO21            ! First dissociation constant for CO2
      REAL,save ::         CO22            ! Second dissociation constant for CO2
      REAL,save ::         CO212           ! CO21*CO22
      REAL,save ::         CO212H          ! CO2H*CO21*CO22
      REAL,save ::         CO21H           ! CO2H*CO21
      REAL,save ::         CO2L            ! CO2 conc in cloudwater (mol/liter)

      REAL,save ::         XL              ! conversion factor (liter-atm/mol)
      REAL,save ::         ONE_OVER_XL     ! 1.0 / XL
      REAL,save ::         PRES_ATM_OVER_XL     ! PRES_ATM / XL
      REAL,save ::         XLCO2           !
      REAL,save ::         XLH2O2          !
      REAL,save ::         XLHNO3          !
      REAL,save ::         XLNH3           !
      REAL,save ::         XLO3            !
      REAL,save ::         XLSO2           !

      REAL,save ::         HCO3            ! HCO3 conc in cloudwater (mol/liter)
      REAL,save ::         HSO3            ! HSO3 conc in cloudwater (mol/liter)
      REAL,save ::         HSO4            ! HSO4 concn in cloudwater (mol/liter)

      REAL,save ::         MN              ! Mn++ conc in cloudwater (mol/liter)
      REAL,save ::         MNA             ! initial Mn in cloudwater (mol/liter)
      REAL,save ::         NH4             ! NH4+ conc in cloudwater (mol/liter)
      REAL,save ::         NO3M            ! NO3- conc in cloudwater (mol/liter)
      REAL,save ::         OHM             ! OH- conc in cloudwater (mol/liter)
      REAL,save ::         SO4             ! SO4= conc in cloudwater (mol/liter)
      REAL,save ::         CO3             ! CO3= conc in cloudwater (mol/liter)

      REAL,save ::         A               ! iron's anion concentration
      REAL,save ::         B               ! manganese's anion concentration
      REAL,save ::         FE              ! Fe+++ conc in cloudwater (mol/liter)
      REAL,save ::         FEA             ! initial Fe in cloudwater (mol/liter)

      INTEGER,save ::      I20C            ! loop counter for do loop 20
      INTEGER,save ::      I30C            ! loop counter for do loop 30
      INTEGER,save ::      ITERAT          ! # iterations of aqueous chemistry solver
      INTEGER,save ::      I7777C          ! aqueous chem iteration counter
      INTEGER,save ::      ICNTAQ          ! aqueous chem iteration counter
      INTEGER,save ::      IOX             ! index over oxidation reactions

      REAL,save ::         ACT1            ! activity correction factor!single ions
      REAL,save ::         ACT2            ! activity factor correction!double ions
      REAL,save ::         ACTB            !

      REAL,save ::         FTST            !
      REAL,save ::         GM1             !
      REAL,save ::         GM1LOG          !
      REAL,save ::         GM2             ! activity correction factor
      REAL,save ::         GM2LOG          !

      REAL,save ::         FA              ! functional value ??
      REAL,save ::         FB              ! functional value ??
      REAL,save ::         AC              ! H+ concentration in cloudwater (mol/liter)
      REAL,save ::         AE              ! guess for H+ conc in cloudwater (mol/l)
      REAL,save ::         BB              ! lower limit guess of cloudwater pH
      REAL,save ::         HA              !
      REAL,save ::         HB              !
      REAL,save ::         STION           ! ionic strength

      REAL,save ::         H2OW            !
      REAL,save ::         HTST            !
      REAL,save ::         RATE            !
      REAL,save ::         RECIPA1         !
      REAL,save ::         RECIPA2         !

      REAL,save ::         RH2O2           !
      REAL,save ::         RMHP            !
      REAL,save ::         RPAA            !
      REAL,save ::         SIV             ! dissolved so2 in cloudwater (mol/liter)
      REAL,save ::         SK6             !
      REAL,save ::         SK6TS6          !
      REAL,save ::         DTS6            !
      REAL,save ::         TAC             !
      REAL,save ::         TEMP1           !
      REAL,save ::         TIMEW           ! cloud chemistry clock (sec)
      REAL,save ::         TOTOX           !
      REAL,save ::         TS6             ! SO4 conc in cloudwater (mol/liter)
      REAL,save ::         TSIV            !
      REAL,save ::         TST             !

      REAL,save ::         DSIVDT( 0:NUMOX ) ! rate of so2 oxid incloud (mol/liter/sec)
      REAL,save ::         DS4   ( 0:NUMOX ) ! S(IV) oxidized over timestep DTW(0)
      real,save ::         ds40old           ! total s(iv) oxidized
      REAL,save ::         DTW   ( 0:NUMOX ) ! cloud chemistry timestep (sec)

      real,save ::         depfactor

!
! --- Define background concs of iron and manganese (ug/m3)
      real cfe, cmn
      data cfe/0.01/, cmn/0.005/

! ... Specify fraction of activation for particles.
      real fracma
      data fracma/1.0/

!...........EXTERNAL FUNCTIONS and their descriptions:

      REAL HLCONST
      EXTERNAL HLCONST

!*********************************************************************
!     begin body of subroutine AQRADM

! --- wangzhm omp: critical
         !---$---OMP CRITICAL (call_aqradm)
! --- wangzhm omp: set threadsprivate
      !$OMP THREADPRIVATE(RT,RECIPAP1,ONE_OVER_TEMP,PH2O20,PH2O2F,H2O2H,H2O2L,&
      !$OMP& PHNO30,PHNO3F,HNO3H,HNO31,HNO31H,HNO3L,&
      !$OMP& PNH30,PNH3F,NH3H,NH31,NH3DH2O,NH31HDH,NH3L,&
      !$OMP& PO30,PO3F,O3H,O3L,PSO20,PSO2F,SO2H,SO21,SO22,SO212,SO212H,SO21H,SO2L,SO3,&
      !$OMP& PNO2,NO2H,NO2L,PNO,NOH,NOL,PCO20,PCO2F,CO2H,CO21,CO22,CO212,CO212H,CO21H,CO2L,&
      !$OMP& XL,ONE_OVER_XL,PRES_ATM_OVER_XL,XLCO2,XLH2O2,XLHNO3,XLNH3,XLO3,XLSO2,&
      !$OMP& HCO3,HSO3,HSO4,MN,MNA,NH4,NO3M,OHM,SO4,CO3,A,B,FE,FEA,&
      !$OMP& I20C,I30C,ITERAT,I7777C,ICNTAQ,IOX,ACT1,ACT2,ACTB,FTST,GM1,GM1LOG,GM2,GM2LOG,&
      !$OMP& FA,FB,AC,AE,BB,HA,HB,STION,H2OW,HTST,RATE,RECIPA1,RECIPA2,&
      !$OMP& RH2O2,RMHP,RPAA,SIV,SK6,SK6TS6,DTS6,TAC,TEMP1,TIMEW,TOTOX,TS6,TSIV,TST,&
      !$OMP& DSIVDT,DS4,ds40old,DTW,depfactor)

!...Check for bad temperature or pressure
      if ( temp .LE. 0.0 .or. pres_atm .LE. 0.0 ) then
         write(io6,*)'Invalid temp and/or pressure; T,P: ',temp,pres_atm
         stop 'Halted in AQRADM -- see list file.'
      end if

      one_over_temp = 1.0 / temp

!...Compute several conversion factors

      icntaq = 0
      iterat = 0
      rt = ( MOLVOL / STDTEMP ) * temp    ! r * t (liter atm / mol)
      xl   = wcavg * rt / H2ODENS         ! conversion factor (l-atm/mol)
      one_over_xl = 1.0 / xl
      pres_atm_over_xl = pres_atm / xl
      tst  = 0.999
      act1 = 1.0
      act2 = 1.0
      gm2  = 1.0
      timew = 0.0
      recipap1 = 1.0 / pres_atm

!...set equilibrium constants as a function of temperature
!...Henry's law constants

      so2h  = HLCONST( 'SO2', temp, .false., 0.0 )
      co2h  = HLCONST( 'CO2', temp, .false., 0.0 )
      nh3h  = HLCONST( 'NH3', temp, .false., 0.0 )
      h2o2h = HLCONST( 'H2O2', temp, .false., 0.0 )
      o3h   = HLCONST( 'O3', temp, .false., 0.0 )
      hno3h = HLCONST( 'HNO3', temp, .false., 0.0 )
      noh   = HLCONST( 'NO', temp, .false., 0.0 )
      no2h  = HLCONST( 'NO2', temp, .false., 0.0 )

!...Dissociation constants

      temp1 = one_over_temp - 1.0 / 298.0

      sk6   = 1.02e-02 * EXP(  2.72e+03 * temp1 )    ! Smith and Martell (1976)
      so21  = 1.30e-02 * EXP(  1.96e+03 * temp1 )    ! Smith and Martell (1976)
      so22  = 6.60e-08 * EXP(  1.50e+03 * temp1 )    ! Smith and Martell (1976)
      co21  = 4.30e-07 * EXP( -1.00e+03 * temp1 )    ! Smith and Martell (1976)
      co22  = 4.68e-11 * EXP( -1.76e+03 * temp1 )    ! Smith and Martell (1976)
      h2ow  = 1.00e-14 * EXP( -6.71e+03 * temp1 )    ! Smith and Martell (1976)
      nh31  = 1.70e-05 * EXP( -4.50e+02 * temp1 )    ! Smith and Martell (1976)
      hno31 = 1.54e+01 * EXP(  8.70e+03 * temp1 )    ! Schwartz (1984)

!...Kinetic oxidation rates
!...   From Chamedies (1982)

      rh2o2 = 8.0e+04 * EXP( -3650.0 * temp1 )

!...Make initializations

      do iox = 0, NUMOX
         dsivdt(iox) = 0.0
         dtw(iox)    = 0.0
         ds4(iox)    = 0.0
      end do

      fea = cfe * 1.E-6 * pres_atm_over_xl / (rhoair * MW_FE)
      mna = cmn * 1.E-6 * pres_atm_over_xl / (rhoair * MW_MN)

!...Set constant factors that will be used in later multiplications (moles/atm)

      xlh2o2  = h2o2h * xl
      xlo3    = o3h   * xl
      xlso2   = so2h  * xl
      xlnh3   = nh3h  * xl
      xlhno3  = hno3h * xl
      xlco2   = co2h  * xl

      so212   = so21  * so22
      so21h   = so21  * so2h
      so212h  = so212 * so2h
      co212   = co21  * co22
      co21h   = co21  * co2h
      co212h  = co22  * co21h
      nh3dh2o = nh31  / h2ow
      nh31hdh = nh3h  * nh3dh2o
      hno31h  = hno31 * hno3h

!...If kinetic calculations are made, return to this point

      i20c = 0
20    continue

      i20c = i20c + 1
      if ( i20c >= 1000 ) then
         write(io6,*) 'Excessive looping at I20C'
         stop 'Halted in AQRADM -- see list file.'
      end if

!...Initial gas phase partial pressures (atm)
      pnh30  = ctnh3 * pres_atm
      phno30 = ( conc(5)  + conc(6)*fracma ) * pres_atm
      pco20  = cco2 * pres_atm

! --- Reactive species

! --- H2O2
      ph2o20 = ch2o2 * pres_atm + xl * ds4( 1 )
! --- check if too much h2o2 has reacted (possible in plume with
! --- high SO2 concs)
      if ( ph2o20 .LT. 0. ) then
         ds4( 0 ) = ds4( 0 ) - ds4( 1 )
         ds4( 1 ) = -ch2o2 * pres_atm * one_over_xl
         ds4( 0 ) = ds4( 0 ) + ds4( 1 )
         ph2o20 = 0.
      end if

! --- O3
      po30   = coz * pres_atm + xl * ds4( 2 )
! --- check if too much o3 has reacted
      if ( po30 .LT. 0. ) then
         ds4( 0 ) = ds4( 0 ) - ds4( 2 )
         ds4( 2 ) = -coz * pres_atm * one_over_xl
         ds4( 0 ) = ds4( 0 ) + ds4( 2 )
         po30 = 0.
      end if

! --- SO2
      pso20  = conc( 1  ) * pres_atm + ds4( 0 ) * xl
! --- check if too much SO2 has reacted
      if ( pso20 .LT. 0. ) then
         ds40old = ds4( 0 )
         ds4( 0 ) = -conc( 1 ) * pres_atm * one_over_xl
         do iox = 1, NUMOX
            ds4( iox ) = ds4( iox ) * ds4( 0 ) / ds40old
         end do
         pso20 = 0.
         ph2o20 = ch2o2 * pres_atm + xl * ds4( 1 )
         po30   = coz * pres_atm + xl * ds4( 2 )
      end if

!...Don't allow gas concentrations to go below zero

!     pso20  = MAX( pso20,  0.0 )
!     ph2o20 = MAX( ph2o20, 0.0 )
!     po30   = MAX( po30,   0.0 )
      pnh30  = MAX( pnh30,  0.0 )
      pco20  = MAX( pco20,  0.0 )
      phno30 = MAX( phno30, 0.0 )

!...Molar concentrations of soluble aerosols
      ts6     = conc( 2 ) * fracma * pres_atm_over_xl                          &
     &        - ds4( 0 )   ! Sulfate

      fe      = fea
      mn      = mna
      a       = 3.0 * fe
      b       = 2.0 * mn

!...Don't allow aerosol concentrations to go below zero

      ts6     = MAX( ts6,     0.0 )

      sk6ts6 = sk6 * ts6

!...Find solution of the equation using a method of reiterative
!...bisections. Make initial guesses for pH between PHMIN  to  PHMAX.

      ha = PHMIN
      hb = PHMAX

      i7777c = 0
7777  continue

      i7777c = i7777c + 1
      if ( i7777c .GT. 1000 ) then
         write(io6,*)'Excessive looping at I7777C'
         stop 'Halted in AQRADM -- see list file.'
      end if

!     ha = MAX( ha - 0.8, 0.1 )
!     hb = MIN( hb + 0.8, 9.9 )
      ha = MAX( ha - 0.8, PHMIN )
      hb = MIN( hb + 0.8, PHMAX )
      ae = 10.0**( -ha )

      recipa1 = 1.0 / ( ae * act1 )
      recipa2 = 1.0 / ( ae * ae * act2 )

!...Calculate final gas phase partial pressure of SO2, NH3, HNO3 and
!...CO2 (atm)

      pso2f = pso20 / ( 1.0 + xlso2 * ( 1.0 + so21 * recipa1                   &
     &      + so212 * recipa2 ) )

      pnh3f = pnh30 / ( 1.0 + xlnh3 * ( 1.0 + nh3dh2o * ae ) )

      phno3f = phno30 / ( 1.0 + xlhno3 * ( 1.0 + hno31 * recipa1 ) )

      pco2f = pco20 / ( 1.0 + xlco2 * ( 1.0 + co21 * recipa1                   &
     &      + co212 * recipa2 ) )

!...Calculate liquid phase concentrations (moles/liter)

      so4  = sk6ts6 / ( ae * gm2 + sk6 )
      hso4 = ts6 - so4
      so3  = so212h  * pso2f  * recipa2
      hso3 = so21h   * pso2f  * recipa1
      co3  = co212h  * pco2f  * recipa2
      hco3 = co21h   * pco2f  * recipa1
      ohm  = h2ow    * recipa1
      nh4  = nh31hdh * pnh3f  * ae
      no3m = hno31h  * phno3f * recipa1

!...Compute functional value

      fa = ae + nh4 - 2.0 *  (co3 + so3 + so4 ) - ohm - hco3                   &
     &   - hso3 - no3m - hso4

!...Start iteration and bisection ****************<<<<<<<

      i30c = 0
30    continue

      i30c = i30c + 1
      if ( i30c .GT. 1000 ) then
         write(io6,*)'Excessive looping at I30C'
         stop 'Halted in AQRADM -- see list file.'
      end if

      bb = 0.5 * ( ha + hb )
      ae = 10.0**( -bb )

! --- don't solve for H+ if fa < 0 at first try
      if ( i7777c .EQ. 1 .and. fa .LT. 0. ) then

         bb = ha
         hb = ha
         ae = 10.0**( -bb )

      end if

      recipa1 = 1.0 / ( ae * act1 )
      recipa2 = 1.0 / ( ae * ae * act2 )

!...Calculate final gas phase partial pressure of SO2, NH3, HNO3 and
!...CO2 (atm)

      pso2f = pso20 / ( 1.0 + xlso2                                            &
     &	    * ( 1.0 + so21 * recipa1 + so212 * recipa2 ) )

      pnh3f = pnh30 / ( 1.0 + xlnh3 * ( 1.0 + nh3dh2o * ae ) )

      phno3f = phno30 / ( 1.0 + xlhno3 * ( 1.0 + hno31 * recipa1 ) )

      pco2f = pco20 / ( 1.0 + xlco2 * ( 1.0 + co21 * recipa1                   &
     &      + co212 * recipa2 ) )

!...Calculate liquid phase concentrations (moles/liter)

      so4  = sk6ts6 / ( ae * gm2 + sk6 )
      hso4 = ts6 - so4
      so3  = so212h  * pso2f  * recipa2
      hso3 = so21h   * pso2f  * recipa1
      co3  = co212h  * pco2f  * recipa2
      hco3 = co21h   * pco2f  * recipa1
      ohm  = h2ow    * recipa1
      nh4  = nh31hdh * pnh3f  * ae
      no3m = hno31h  * phno3f * recipa1

!...compute functional value

      fb = ae + nh4 - 2.0 * ( co3 + so3 + so4 ) - ohm - hco3                   &
     &   - hso3 - no3m - hso4

!...Calculate and check the sign of the product of the two functional values

      ftst = fa * fb
      if ( ftst .LE. 0.0 ) then
        hb = bb
      else
        ha = bb
        fa = fb
      end if

!...Check convergence of solutions

      htst = ha / hb
      if ( htst .LE. tst ) go to 30

!...end of zero-finding routine ****************<<<<<<<<<<<<

!...compute Ionic strength and activity coefficient by the Davies equation

      stion = 0.5 * (ae + nh4 + ohm + hco3 + hso3                              &
     &      + 4.0 * (so4 + co3 + so3)                                          &
     &      + no3m + hso4 + 9.0 * fe + a + b)
      gm1log = -0.509 * ( SQRT( stion )                                        &
     &       / ( 1.0 + SQRT( stion ) ) - 0.2 * stion )
      gm2log = gm1log * 4.0
      gm1  = 10.0**gm1log
      gm2  = MAX( 10.0**gm2log, 1.0e-30 )
      actb = act1
      act1 = MAX( gm1 * gm1, 1.0e-30 )
      act2 = MAX( gm1 * gm1 * gm2, 1.0e-30 )

!...check for convergence and possibly go to 7777, to recompute
!...  Gas and liquid phase concentrations

! --- don't solve for H+ if fa < 0 at first try
      if ( i7777c .EQ. 1 .and. fa .LT. 0. ) then
         actb = act1
      end if

      tac = ABS( actb - act1 ) / actb
      if ( tac >= 1.0e-2 ) then

         icntaq = icntaq + 1
         if ( icntaq .GT. 100 ) then
            write(io6,*)'Maximum iterations for pH calculation exceeded'
            write(io6,*)'Using last pH value'
         else
           go to 7777
         end if

      end if

!...return an error if the pH is not in range

!cc      if ( ( ha .lt. 0.02 ) .or. ( ha .gt. 9.49 ) ) then
      if ( ( ha .LT. PHMIN ) .or. ( ha .GT. PHMAX ) ) then
         write(io6,*)'pH value out of range: ',ha
         stop 'Halted in AQRADM -- see list file.'
      end if

!...Make those concentration calculations which can be made outside
!...  of the function.

      so2l = so2h * pso2f
      ac = 10.0**( -bb )
      siv = so3 + hso3 + so2l

!...Calculate final gas phase concentrations of oxidants (atm)

      ph2o2f = ph2o20 / ( 1.0 + xlh2o2 )
      po3f   = po30   / ( 1.0 + xlo3   )

      ph2o2f = MAX( ph2o2f, 0.0 )
      po3f   = MAX( po3f,   0.0 )

!...Calculate liquid phase concentrations (moles/liter)

      h2o2l = ph2o2f * h2o2h
      o3l   = po3f   * o3h
      nh3l  = pnh3f  * nh3h
      co2l  = pco2f  * co2h
      hno3l = phno3f * hno3h

!...if the maximum cloud lifetime has not been reached, then compute
!...the next timestep.

      if ( timew .LT. taucld ) then

!...make kinetics calculations
!...  note: DS4(i) and DSIV(I) are negative numbers!

        iterat = iterat + 1

!...Define the total S(iv) available for oxidation

        tsiv = pso20 * one_over_xl

!...Calculate sulfur iv oxidation rate due to H2O2

        dsivdt( 1 ) = -rh2o2 * h2o2l * so2l / ( 0.1 + ac )
        totox = ph2o20 * one_over_xl
        if ( ( dsivdt( 1 ) .EQ. 0.0 ) .or.                                     &
     &       ( tsiv  .LE. CONCMIN ) .or.                                       &
     &       ( totox .LE. CONCMIN ) ) then
          dtw(1) = taucld
        else
          dtw( 1 ) = -0.05 * MIN( totox, tsiv ) / dsivdt( 1 )
        end if

!...Calculate sulfur iv oxidation rate due to O3

        if ( bb .GE. 2.7 ) then
          dsivdt( 2 ) = -4.19e5 * ( 1.0 + 2.39e-4 / ac ) * o3l * siv
        else
          dsivdt( 2 ) = -1.9e4 * siv * o3l / SQRT( ac )
        end if
        totox = po30 * one_over_xl
        if ( ( dsivdt( 2 ) .EQ. 0.0 ) .or.                                     &
     &       ( tsiv  .LE. CONCMIN ) .or.                                       &
     &       ( totox .LE. CONCMIN ) ) then
          dtw( 2 ) = taucld
        else
          dtw( 2 ) = -0.01 * MIN( totox, tsiv ) / dsivdt( 2 )
        end if

!...Calculate sulfur iv oxidation rate due to O2 catalyzed by Mn++
!...  and Fe+++  See Table IV Walcek & Taylor ( 1986)

        if ( bb .GE. 4.0 )  then  ! 4.0  < ph

          if ( siv .LE. 1.0e-5 ) then
            dsivdt( 3 ) = -5000.0 * mn * hso3
          else
            dsivdt( 3 ) = -( 4.7 * mn * mn / ac                                &
     &                  + 1.0e7 * fe * siv * siv )
          end if  ! end of first pass through siv conc.

        else          ! ph < 4.0

          if ( siv .LE. 1.0e-5 ) then
            dsivdt( 3 ) = -3.0 * ( 5000.0 * mn * hso3                          &
     &                  + 0.82 * fe * siv / ac )
          else
            dsivdt( 3 ) = -( 4.7 * mn * mn / ac                                &
     &                  + ( 0.82 * fe * siv / ac )                             &
     &                  * ( 1.0 + 1.7e3 * mn**1.5 / ( 6.3e-6 + fe ) ) )
          end if ! end of second pass through siv conc.
        end if  ! end of pass through ph

        if ( ( dsivdt( 3 ) .EQ. 0.0 ) .or. ( tsiv .LE. CONCMIN ) ) then
          dtw( 3 ) = taucld
        else
          dtw( 3 ) = -0.1 * tsiv / dsivdt( 3 )
        end if

!...Calculate total sulfur iv oxidation rate

        dsivdt( 0 ) = 0.0
        do iox = 1, NUMOX
          dsivdt( 0 ) = dsivdt( 0 ) + dsivdt( iox )
        end do

!...Calculate a minimum time step required

        dtw( 0 ) = MIN( dtw( 1 ), dtw( 2 ), dtw( 3 ) )

!...check for large time step

        if ( dtw( 0 ) .GT. 8.0e+37 ) then
          write(io6,1001) dsivdt(0), ts6, dtw(0)
        else

!...calculate the change in sulfur iv for this time step

60        continue
          dts6 = ABS( dtw( 0 ) * ( -dsivdt( 0 ) ) )

!...If DSIV(0), sulfur iv oxidized during this time step would be
!...less than 5% of sulfur oxidized since time 0, then double DT

          if ( dtw( 0 ) .LE. taucld ) then
            if ( dts6 .LT. 0.05 * ts6 ) then
              dtw( 0 ) = dtw( 0 ) * 2.0
	      go to 60
            end if
          end if
        end if
        dtw( 0 ) = MIN( dtw( 0 ), taucld )

!...If the total time after this time increment will be greater than
!...  TAUCLD sec., then set DTW(0) so that total time will be TAUCLD

        if ( timew + dtw( 0 ) .GT. taucld ) dtw( 0 ) = taucld - timew
        if ( ts6 .LT. 1.0e-11 ) dtw( 0 ) = taucld - timew
        if ( iterat .GT. 100 ) dtw( 0 ) = taucld - timew

!...Set DSIV(I), I = 0,NUMOX, the amount of S(IV) oxidized by each
!... individual oxidizing agent, as well as the total.

        do iox = 0, NUMOX
           ds4( iox ) = ds4( iox ) + dtw( 0 ) * dsivdt( iox )
        end do

        timew = timew + dtw( 0 )

!...Return to make additional calculations

        go to 20
      end if

! --- Calculate liquid-phase concentrations of other species
      pno    = conc(3) * pres_atm
      pno2   = conc(4) * pres_atm

      nol    = pno   * noh   / ( 1.0 + noh    * xl )
      no2l   = pno2  * no2h  / ( 1.0 + no2h   * xl )

!...Compute the output concentrations

!...gas concentrations (mol/molV) (only for reactive species)

      conc(1) = (pso2f  + xl * siv)   * recipap1
      ch2o2   = (ph2o2f + xl * h2o2l) * recipap1
      coz     = (po3f   + xl * o3l)   * recipap1

! --- calculate scavenging coefficients for gases
      depfactor = prcrate * SEC2HR / ( rhoair * len )
      if ( conc( 1 ) > CONCMIN ) then
         scav( 1 ) = siv * depfactor / conc( 1 )
      end if
      if ( conc( 3 ) > CONCMIN ) then
         scav( 3 ) = nol * depfactor / conc( 3 )
      end if
      if ( conc( 4 ) > CONCMIN ) then
         scav( 4 ) = no2l * depfactor / conc( 4 )
      end if
      if ( conc( 5 ) > CONCMIN ) then
         scav( 5 ) = hno3l * depfactor / conc( 5 )
      end if

!...aerosol concentrations (mol/molV) (only reactive species-so4)

      conc( 2 ) = conc( 2 ) * ( 1.0 - fracma ) + ts6 * xl * recipap1

! --- calculate scavenging coefficients for aerosols

      if ( conc( 2 ) > CONCMIN ) then
         scav( 2 ) = ts6 * depfactor / conc( 2 )
      end if
      if ( conc( 6 ) > CONCMIN ) then
         scav( 6 ) = no3m * depfactor / conc( 6 )
      end if
! --- wangzhm omp: critical
         !---$---OMP END CRITICAL (call_aqradm)
      return

!...formats

1001  format (1X,'DSIVDT(0) =', F10.5,                                         &
     &       'TS6=', F10.5, 'DTW(0)=', F10.5)

      end

!----------------------------------------------------------------------
      REAL FUNCTION HLCONST ( NAME, TEMP, EFFECTIVE, HPLUS )
!----------------------------------------------------------------------

!
!  FUNCTION: return the Henry's law constant for the specified substance
!            at the given temperature
!  Adapted for CALPUFF from CMAQ version, PK, AER, Feb 2007

      IMPLICIT NONE

!...........PARAMETERS and their descriptions:
!      include 'params.puf'
! tem&orarily set io6=2 here since params.puf variables are not declared,
! res&lting in compile errors.
      INTEGER IO6
      PARAMETER (IO6 = 2) ! from params.puf

! Num&er of species
      INTEGER MXSPCS
      PARAMETER (MXSPCS = 8)
! Num&er of dissociating species
      INTEGER MXDSPCS
      PARAMETER (MXDSPCS = 8)

!...........ARGUMENTS and their descriptions

      CHARACTER*(*) NAME                ! name of substance
      REAL          TEMP                ! temperature (K)
      LOGICAL       EFFECTIVE           ! true=compute the effective henry's law
      REAL          HPLUS               ! hydrogen ion concentration (mol/l)

!...........SCRATCH LOCAL VARIABLES and their descriptions:
      CHARACTER*16 SUBNAME(MXSPCS)
      SAVE SUBNAME
! --- wangzhm save
      INTEGER,save ::       SPC                 ! species index
      INTEGER,save ::       LSO2                ! SO2 pointer
      INTEGER,save ::       LHSO3               ! HSO3 pointer
      INTEGER,save ::       LHNO3               ! HNO3 pointer
      INTEGER,save ::       LCO2                ! CO2 pointer
      INTEGER,save ::       LHCO3               ! HCO3 pointer
      INTEGER,save ::       LH2O2               ! H2O2 pointer
      INTEGER,save ::       LHO2                ! HO2 pointer
      INTEGER,save ::       LNH4OH              ! NH4OH pointer
      INTEGER,save ::       LH2O                ! H2O pointer

      REAL,save ::          HPLUSI              ! 1 / HPLUS
      REAL,save ::          HPLUS2I             ! 1 / HPLUS**2
      REAL,save ::          TFAC                ! (298-T)/(T*298)
      REAL,save ::          AKEQ1               ! temp var for dissociation constant
      REAL,save ::          AKEQ2               ! temp var for dissociation constant
      REAL,save ::          OHION               ! OH ion concentration
      REAL,save ::          KH                  ! temp var for henry's law constant
                                                                               &
! Hen&y's law constants at 298.15K (M/atm) (taken from Rolf Sanders'           &
! Com&ilation of Henry's Law Constants for Inorganic and Organic Species       &
! of &otential Importance in Environment Chemistry, 1999)
      REAL A(MXSPCS)
      SAVE A
                                                                               &
! Ent&alpy (like activation energy) (K) (taken from Rolf Sanders'              &
! Com&ilation of Henry's Law Constants for Inorganic and Organic Species       &
! of &otential Importance in Environment Chemistry, 1999)
      REAL E(MXSPCS)
      SAVE E
                                                                               &
! Dis&ociation constants at 298.15K (M or M2) (taken from Table 6.A.1,         &
! Sei&feld and Pandis, Atmospheric Chemistry and Physics, 1997)
      REAL B(MXDSPCS)
      SAVE B
                                                                               &
! -dH&R (K) (taken from Table 6.A.1,                                           &
! Sei&feld and Pandis, Atmospheric Chemistry and Physics, 1997)
      REAL D(MXDSPCS)
      SAVE D

      DATA SUBNAME(1), A(1), E(1) / 'O3', 1.2E-02, 2.7E+03 /  ! Chameides 1984
      DATA SUBNAME(2), A(2), E(2) / 'H2O2', 8.3E+04, 7.4E+03 /  ! O'Sullivan et
      DATA SUBNAME(3), A(3), E(3) / 'NH3', 6.1E+01, 4.2E+03 /  ! Clegg and Brimb
      DATA SUBNAME(4), A(4), E(4) / 'NO', 1.9E-03, 1.4E+03 /  ! Lide and Frederi
      DATA SUBNAME(5), A(5), E(5) / 'NO2', 1.2E-02, 2.5E+03 /  ! Chameides 1984
      DATA SUBNAME(6), A(6), E(6) / 'HNO3', 2.1E+05, 8.7E+03 /  ! Leieveld and C
      DATA SUBNAME(7), A(7), E(7) / 'SO2', 1.4E+00, 2.9E+03 /  ! Linde and Frede
      DATA SUBNAME(8), A(8), E(8) / 'CO2', 3.6E-02, 2.2E+03 /  ! Zheng et al. 19

      DATA LSO2,  B(1), D(1) /  1, 1.30E-02,  1.96E+03 /  ! SO2*H2O<=>HSO3+H
      DATA LHSO3, B(2), D(2) /  2, 6.60E-08,  1.50E+03 /  ! HSO3<=>SO3+H
      DATA LHNO3, B(3), D(3) /  3, 1.54E+01,  8.70E+03 /  ! HNO3(aq)<=>NO3+H
      DATA LCO2,  B(4), D(4) /  4, 4.30E-07, -1.00E+03 /  ! CO2*H2O<=>HCO3+H
      DATA LHCO3, B(5), D(5) /  5, 4.68E-11, -1.76E+03 /  ! HCO3<=>CO3+H
      DATA LH2O2, B( 6), D(6) /  6, 2.20E-12, -3.73E+03 /  ! H2O2(aq)<=>HO2+H
      DATA LNH4OH, B(7), D(7) /  7, 1.70E-05, -4.50E+02 /  ! NH4*OH<=>NH4+OH
      DATA LH2O,   B(8), D(8) /  8, 1.00E-14, -6.71E+03 /  ! H2O<=>H+OH

!...........EXTERNAL FUNCTIONS and their descriptions:
                                                                               &
! Fun&tion to look up name in table
      INTEGER INDEX1
      EXTERNAL INDEX1

! --- wangzhm omp: set threadsprivate
      !$OMP THREADPRIVATE(SUBNAME,SPC,LSO2,LHSO3,LHNO3,LCO2,LHCO3,LH2O2,LHO2,LNH4OH,LH2O,&
      !$OMP& HPLUSI,HPLUS2I,TFAC,AKEQ1,AKEQ2,OHION,KH,A,E,B,D)

!-----------------------------------------------------------------------
!  begin body of subroutine HLCONST

      SPC = INDEX1( NAME, MXSPCS, SUBNAME )

!...error if species not found in table

      IF ( SPC <= 0 ) THEN
         write(io6,*)TRIM(NAME) //                                             &
     &         ' not found in Henrys Law Constant table'
         stop 'Halted in HLCONST -- see list file.'
      END IF

!...compute the Henry's Law Constant
      TFAC = (298.0 - TEMP) / (298.0 * TEMP)
      KH = A(SPC) * EXP(E(SPC) * TFAC)
      HLCONST = KH

!...compute the effective Henry's law constants

      IF (EFFECTIVE) THEN

         IF ( HPLUS <= 0.0 ) THEN
            write(io6,*)'Negative or Zero [H+] concentration specified '
            stop 'Halted in HLCONST -- see list file.'
         END IF

         HPLUSI = 1.0 / HPLUS
         HPLUS2I = HPLUSI * HPLUSI

         IF (TRIM(NAME) .EQ. 'SO2') THEN

            AKEQ1 = B(LSO2) * EXP(D(LSO2) * TFAC)     !SO2H2O <=> HSO3- + H+
            AKEQ2 = B(LHSO3) * EXP(D(LHSO3) * TFAC)   !HSO3- <=> SO3= + H+
            HLCONST = KH * (1.0 + AKEQ1*HPLUSI + AKEQ1*AKEQ2*HPLUS2I)

         ELSE IF (TRIM(NAME) .EQ. 'HNO3') THEN

            AKEQ1 = B(LHNO3) * EXP(D(LHNO3) * TFAC)   !HNO3(aq) <=> NO3- + H+
            HLCONST = KH * (1.0 + AKEQ1*HPLUSI)

         ELSE IF (TRIM(NAME) .EQ. 'CO2') THEN

            AKEQ1 = B(LCO2) * EXP(D(LCO2)  * TFAC)    !CO2H2O <=> HCO3- + H+
            AKEQ2 = B(LHCO3) * EXP(D(LHCO3) * TFAC)   !HCO3- <=> CO3= + H+
            HLCONST = KH * (1.0 + AKEQ1*HPLUSI + AKEQ1*AKEQ2*HPLUS2I)

         ELSE IF (TRIM(NAME) .EQ. 'H2O2') THEN

            AKEQ1 = B(LH2O2) * EXP(D(LH2O2) * TFAC)   !H2O2(aq) <=> HO2- + H+
            HLCONST = KH * (1.0 + AKEQ1*HPLUSI)

         ELSE IF (TRIM(NAME) .EQ. 'NH3') THEN

            AKEQ1 = B(LNH4OH) * EXP(D(LNH4OH) * TFAC) !NH4OH <=> NH4+ + OH-
            AKEQ2 = B(LH2O) * EXP(D(LH2O) * TFAC)
            OHION = AKEQ2 * HPLUSI
            HLCONST = KH * (1.0 + AKEQ1/OHION)

         END IF

      END IF

      RETURN
      END

!----------------------------------------------------------------------
      INTEGER FUNCTION INDEX1 (NAME, N, NLIST)
!----------------------------------------------------------------------

!***********************************************************************
!  subroutine body starts at line 39
!
!  FUNCTION:
!
!    Searches for NAME in list NLIST and returns the subscript
!    (1...N) at which it is found, or returns 0 when NAME not
!    found in NLIST
!
!  PRECONDITIONS REQUIRED:  none
!
!  SUBROUTINES AND FUNCTIONS CALLED:  none
!
!  Based on index1 routine from Models-3 I/O Library
!
!***********************************************************************

      IMPLICIT NONE

!.......   Arguments and their descriptions:

      CHARACTER*(*) NAME        !  Character string being searched for
      INTEGER       N           !  Length of array to be searched
      CHARACTER*(*) NLIST(*)    !  array to be searched

!.......   Local variable:

      INTEGER       I   !  loop counter

!.....................................................................
!.......   begin body of INDEX1()

      DO 100 I = 1, N

          IF (TRIM(NAME) .EQ. TRIM(NLIST(I))) THEN    ! Found NAME in NLIST
              INDEX1 = I
              RETURN
          END IF

100   CONTINUE

      INDEX1 = 0        !  not found
      RETURN

      END