      subroutine aerochem_cf(h2o,tempk,press,cwc,cph,aph,ajno2,
     &                       cNV,con,convfac,
     &                       dtaq,dtaer,delcon,scNV,sddm,nfam,nsen,
     &                       ldoipr,ipa_cel,aero_flag,lpig,icig,ldark)
      use camxcom
      use chmstry
      use filunit
      use tracer
      use procan
c
c----CAMx v7Beta6 190902
c
c     AEROCHEM_CF drives the CF 2-section aerosol model in CAMx.
c     It calculates the following chemical transformations:
c        1. Condensible organic gasses to organic aerosol (SOAP or VBS)
c        2. Gaseous sulfate to aerosol sulfate (from gas-phase chem)
c        3. Inorganic aqueous reactions and aqueous formation of SOA
c           (CMAQ approach)
c        4. Inorganic gas-aerosol equilibrium partitioning for the
c           ammonium/nitrate/sulfate/sodium/chloride system (ISORROPIA)
c     The gas species are in ppm
c     The aerosol species are in (ug/m3)
c
c     Copyright 1996 - 2018
c     Ramboll 
c
c     Modifications: 
c        1/9/02         Minor bug fixes (units conversions)
c        11/19/02       Incorporated RADM aqueous chemistry and ISORROPIA
c        12/9/02        Incorporated SOAP routine
c        10/24/05       Cloud pH is now passed in/out; Henry's Law values
c                       for SO2 and O3 no longer passed out
c        8/23/06        Removed metastable option in ISOROPIA
c        12/26/06       AQ-CHEM now called each chemistry timestep
c        12/29/06       Revised for the updated SOA scheme
c        01/10/07       Added species mapping for CB05/SAPRC99
c        05/04/07       Added array (wt_save) to keep original total
c                       mass because ISOROPIA overwrites with lower
c                       bound value.
c        05/30/07       Skip call to RAQCHEM under high acid conditions.
c        05/13/08       Modified to skip ISOROPIA for certain PiG calls.
c        01/21/10       Added formation of mineral dust nitrate
c        07/31/13       Modified for new GREASD PiG chemistry 
c        08/23/13       Revised for DDM-PM; replaced hard-wired MW's with wtmol
c        03/18/14       Revised for benzene SOA
c        10/09/14       Added subgrid convective model
c        12/07/14       Revised for VBS
c        01/08/16       Updated for Revised VBS
c        08/25/16       Updated for new SOAP
c        09/16/16       Aqueous chemistry updates
c        09/14/17       Reestablished metastable option in ISOROPIA
c        10/13/17       Modified to get aerosol pH from ISOROPIA
c        10/16/17       Added particle-phase SOA photolysis to SOAP
c        11/22/17       Added PFE/PMN/PK/PCA/PMG
c                       Updated for revised WTMOL
c                       Added EQSAM option
c        10/26/18       Updated DDM for SOA photolysis & revised RADM
c
c     Input arguments:
c        h2o                 cell water vapor (ppm)
c        tempk               cell temperature (K)
c        press               cell pressure (mb)
c        cwc                 cloud water content (g/m3)
c        cph                 cloud water pH 
c        aph                 aerosol pH
c        ajno2               time-weighted accumulated NO2 photolysis rate ([1/hr]*[hr])
c        cNV                 NV product concentrations (ppm)
c        con                 species concentrations (ppm, ug/m3)
c        convfac             conversion factor: umol/m3 = ppm * convfac
c        dtaq                aqueous chemistry time step (hours)
c        dtaer               general PM chemistry time step (hours)
c        delcon              delta concentrations for source apportionment
c        scNV                NV product sensitivities
c        sddm                DDM-PM sensitivity array
c        nfam                number of DDM parameters
c        nsen                number of model species used in DDM-PM
c        ldoipr              flag to determine whether to do IPR calculation
c        ipa_cel             pointer into PA sub-domain array
c        aero_flag           1 - do only RAQCHEM (with ISOROPIA);
c                            2 - do SOAP, RAQCHEM, ISOROPIA
c        lpig                T - PiG call (do only RAQCHEM)
c        icig                0 - CiG not activated (default to aero_flag)
c                            1 - CiG call (do only RAQCHEM)
c                            2 - Grid call, CiG called prior (skip RAQCHEM)
c        ldark               darkness flag
c
c     Output arguments:
c        con                 species concentrations (ppm, ug/m3)
c        cph                 cloud water pH 
c        aph                 aerosol pH
c        delcon              delta concentrations for source apportionment
c        sddm                DDM-PM sensitivity array
c
c     Routines called: 
c        SOAP
c        RAQCHEM
c        ISOROPIA
c 
c     Called by: 
c        IRONCHEM
c        CIGDRIVE
c        CHEMDRIV 
c
      implicit none
      include 'camx.prm' 
      include 'camx_aero.inc'
      include 'soap.inc'
      include 'vbs.inc'
      include 'chmdbg.inc'
c
c========================== Process Analysis Begin =======================
c
      logical ldoipr
      integer ipa_cel
      real csulf,chno3,cnh3,cpso4,cpno3,cpnh4
c
c========================== Process Analysis End =======================
c
c========================= DDM Begin ===================================
c
      integer nfam, nsen, ierr, ip
      real scNV(nfam,NNV)
      real sddm(nfam,nsen)
c
c========================== DDM End ====================================
c
      real    h2o
      real    tempk
      real    press
      real    cwc
      real    cph
      real    aph
      real    ajno2
      real    cNV(NNV)
      real    con(*)
      real    convfac
      real    dtaq
      real    dtaer
      real    delcon(7,*)
      real    so2
      real    so4
      real    no3

      integer aero_flag
      integer icig

      logical lpig
      logical lnvoc
      logical lsvoc
      logical lradm
      logical lisor
      logical ldark ! darkness flag added for day-night difference in Fe_III
      logical ldosa

      integer ispc,is
      real qwatr,ev,es,rh
      real pres_pa
      real dt_sec
      real cold(MXSPEC), deln2o5
      real cw_kgm3
c
c-----Arrays for SOAP
c
      real soacon(NCG),cg(NCG),csatT(NCG)
      real preOM ! pre-existing OA mass (mole)
      real apoly(NSOAP), aphot(NSOAP+NNVOL)
      real fpoly, fphot
c
c-----Arrays for RADM aqueous chemistry
c
      real r_gas(15),r_aer(10),r_old(25)
      real r_1,r_2
      real csoac
c
c-----Variables for EQSAM4clim
c
      integer, parameter :: dp = SELECTED_REAL_KIND(12,307)
      real(dp) :: xTT,xAW,xWH2O,xYPa(5),xYMa(4)
      real(dp) :: xsPM,xaPM,xPMS,xPMT,xRHO,xVOL,xPH,xHp,xGF,
     &            xYPs(5),xYMs(4),xYG(4)
      real(dp), parameter :: Dd = 1.e-6_dp ! dry droplet diameter for the Kelvin effect (not used)
      integer, parameter :: jmeq = 1
c
c-----Variables for ISOROPIA need double precision
c
      integer imask(1,1)
      real*8 wi(5),wt(5),wt_save(5)
      real*8 rhi,tempi,cntrl(2)
      real*8 gasis(3),aerliq(12),aersld(9),other(6)
      character*15 scasi
c
c-----Variables for calcium nitrate
c
      real*8 dno3 ! calcium nitrate (mol/m3)
      logical lno3_isr ! flag to indicate non-zero nitrate for isorropia
c
c-----Variables for iron and manganese fractions in dust and other primary PM
c
      real fe_iii, mn_ii, fe_sol, mn_sol, fe_mmr, mn_mmr
c
c-----Speciation factors for elemental species (negative to fall back to background values)
c
      real, parameter :: f_FE_PRM = -1.0   , f_FE_CRS = -1.0    ! mass fraction of Fe in PRM or CRS
      real, parameter :: f_MN_PRM = -1.0   , f_MN_CRS = -1.0    ! mass fraction of Mn in PRM or CRS
      real, parameter :: f_CA_PRM = 0.0260 , f_CA_CRS = 0.0260  ! mass fraction of Ca in PRM or CRS
      real, parameter :: f_MG_PRM = -1.0   , f_MG_CRS = -1.0    ! mass fraction of Mn in PRM or CRS
      real, parameter :: f_K_PRM  = -1.0   , f_K_CRS  = -1.0    ! mass fraction of K in PRM or CRS
c 
c-----Entry point 
c
      lnvoc = .false.
      lsvoc = .false.
      lradm = .false.
      lisor = .false.
      ldosa = .false.

      if (.not.lpig .and. icig.ne.1) then
        lnvoc = .true.
        if (aero_flag.eq.2) then
          lsvoc = .true.
          lisor = .true.
        endif
        ldosa = ltrace
      endif
      if (cwc.ge.cwmin .and. tempk.gt.tamin) then
        so2 = con(kso2)*1.e-6
        so4 = (con(kpso4)/wtmol(1)/convfac)*1.e-6
        no3 = (con(kpno3)/wtmol(3)/convfac)*1.e-6
        if (so4.lt.5.0e-7 .and. no3.lt.1.5e-6 .and. so2.lt.1.0e-6) then
          lradm = .true.
          lisor = .true.
        endif
        if (lpig) lisor = .false.
        if (icig.eq.1) lisor = .false.
        if (icig.eq.2) lradm = .false.
      endif

      do ispc = 1,nspec
         cold(ispc) = con(ispc)
      enddo
      con(nspec+1) = 0.
      cph = 5.0
c
c-----Calculate relative humidity
c
      qwatr = 1.e-6*h2o*18./28.8 
      ev = qwatr*press/(qwatr + eps) 
      es = e0*exp((lv/rv)*(1./273. - 1./tempk)) 
      rh = 100.0*amin1(0.99,ev/es)
c
c-----Partitioning of non-volatile organic compounds.
c     Unit conversion is done here: [ppm] -> [ug/m3]
c
      if ( lnvoc ) then

        if ( lvbs ) then

          con(kPAS0) = con(kPAS0) + cNV(1) * convfac * mwPAS0
          con(kPBS0) = con(kPBS0) + cNV(2) * convfac * mwPBS0
          con(kPAP0) = con(kPAP0) + cNV(3) * convfac * mwPAP0
          con(kPCP0) = con(kPCP0) + cNV(4) * convfac * mwPCP0
          con(kPFP0) = con(kPFP0) + cNV(5) * convfac * mwPFP0
          if ( ldoipr ) then
            cipr(IPR_OAERO,ipa_cel,kPAS0) = cipr(IPR_OAERO,ipa_cel,kPAS0)
     &                                           + con(kPAS0)-cold(kPAS0)
            cipr(IPR_OAERO,ipa_cel,kPBS0) = cipr(IPR_OAERO,ipa_cel,kPBS0)
     &                                           + con(kPBS0)-cold(kPBS0)
            cipr(IPR_OAERO,ipa_cel,kPAP0) = cipr(IPR_OAERO,ipa_cel,kPAP0)
     &                                           + con(kPAP0)-cold(kPAP0)
            cipr(IPR_OAERO,ipa_cel,kPCP0) = cipr(IPR_OAERO,ipa_cel,kPCP0)
     &                                           + con(kPCP0)-cold(kPCP0)
            cipr(IPR_OAERO,ipa_cel,kPFP0) = cipr(IPR_OAERO,ipa_cel,kPFP0)
     &                                           + con(kPFP0)-cold(kPFP0)
          endif

        else ! SOAP

          con(ksopa) = con(ksopa) + cNV(1) * convfac * mwsopa
          con(ksopb) = con(ksopb) + cNV(2) * convfac * mwsopb
          if ( ldoipr ) then
            cipr(IPR_OAERO,ipa_cel,ksopa) = cipr(IPR_OAERO,ipa_cel,ksopa)
     &                                         + con(ksopa) - cold(ksopa)
            cipr(IPR_OAERO,ipa_cel,ksopb) = cipr(IPR_OAERO,ipa_cel,ksopb)
     &                                         + con(ksopb) - cold(ksopb)
          endif
c
c========================= DDM Begin ===================================
c
          if (lddm) then
            do ip = 1, nfam
              sddm(ip,jdpSOPA) = sddm(ip,jdpSOPA) + scNV(ip,1) * convfac * mwsopa
              sddm(ip,jdpSOPB) = sddm(ip,jdpSOPB) + scNV(ip,2) * convfac * mwsopb
            enddo
          endif
c
c========================== DDM End ====================================
c
        endif ! VBS or SOAP?

      endif ! lnvoc?
c
c-----Partitioning of semi-volatile organic compounds between the gas and
c     aerosol phases using the SOAP or VBS scheme.
c     CGs are in ppm; SOAs are in ug/m3; the unit conversion is done here.
c
      if ( lsvoc ) then

        if ( lvbs ) then

          cg(1)  = con(kVAS1) * convfac * mwVBS(1)
          cg(2)  = con(kVAS2) * convfac * mwVBS(2)
          cg(3)  = con(kVAS3) * convfac * mwVBS(3)
          cg(4)  = con(kVAS4) * convfac * mwVBS(4)
          cg(5)  = con(kVBS1) * convfac * mwVBS(5)
          cg(6)  = con(kVBS2) * convfac * mwVBS(6)
          cg(7)  = con(kVBS3) * convfac * mwVBS(7)
          cg(8)  = con(kVBS4) * convfac * mwVBS(8)
          cg(9)  = con(kVAP1) * convfac * mwVBS(9)
          cg(10) = con(kVAP2) * convfac * mwVBS(10)
          cg(11) = con(kVAP3) * convfac * mwVBS(11)
          cg(12) = con(kVAP4) * convfac * mwVBS(12)
          cg(13) = con(kVCP1) * convfac * mwVBS(13)
          cg(14) = con(kVCP2) * convfac * mwVBS(14)
          cg(15) = con(kVCP3) * convfac * mwVBS(15)
          cg(16) = con(kVCP4) * convfac * mwVBS(16)
          cg(17) = con(kVFP1) * convfac * mwVBS(17)
          cg(18) = con(kVFP2) * convfac * mwVBS(18)
          cg(19) = con(kVFP3) * convfac * mwVBS(19)
          cg(20) = con(kVFP4) * convfac * mwVBS(20)
          soacon(1)  = con(kPAS1)
          soacon(2)  = con(kPAS2)
          soacon(3)  = con(kPAS3)
          soacon(4)  = con(kPAS4)
          soacon(5)  = con(kPBS1)
          soacon(6)  = con(kPBS2)
          soacon(7)  = con(kPBS3)
          soacon(8)  = con(kPBS4)
          soacon(9)  = con(kPAP1)
          soacon(10) = con(kPAP2)
          soacon(11) = con(kPAP3)
          soacon(12) = con(kPAP4)
          soacon(13) = con(kPCP1)
          soacon(14) = con(kPCP2)
          soacon(15) = con(kPCP3)
          soacon(16) = con(kPCP4)
          soacon(17) = con(kPFP1)
          soacon(18) = con(kPFP2)
          soacon(19) = con(kPFP3)
          soacon(20) = con(kPFP4)
          preOM = con(kPAS0) / mwPAS0
     &          + con(kPBS0) / mwPBS0
     &          + con(kPAP0) / mwPAP0
     &          + con(kPCP0) / mwPCP0
     &          + con(kPFP0) / mwPFP0

          do ispc = 1, NVBS*NSET
            csatT(ispc) = csatV(ispc) * (cstempV(ispc)/tempk)
     &                                * exp((deltahV(ispc)/8.314)
     &                                * (1/cstempV(ispc)-1/tempk))
          enddo

          call soap(NSET*NVBS,mwVBS,soacon(:NSET*NVBS),cg(:NSET*NVBS),
     &              csatT(:NSET*NVBS),preOM,
     &              iout,igrdchm,ichm,jchm,kchm)

          con(kVAS1) = amax1(cg(1) /convfac/mwVBS(1) ,bdnl(kVAS1))
          con(kVAS2) = amax1(cg(2) /convfac/mwVBS(2) ,bdnl(kVAS2))
          con(kVAS3) = amax1(cg(3) /convfac/mwVBS(3) ,bdnl(kVAS3))
          con(kVAS4) = amax1(cg(4) /convfac/mwVBS(4) ,bdnl(kVAS4))
          con(kVBS1) = amax1(cg(5) /convfac/mwVBS(5) ,bdnl(kVBS1))
          con(kVBS2) = amax1(cg(6) /convfac/mwVBS(6) ,bdnl(kVBS2))
          con(kVBS3) = amax1(cg(7) /convfac/mwVBS(7) ,bdnl(kVBS3))
          con(kVBS4) = amax1(cg(8) /convfac/mwVBS(8) ,bdnl(kVBS4))
          con(kVAP1) = amax1(cg(9) /convfac/mwVBS(9) ,bdnl(kVAP1))
          con(kVAP2) = amax1(cg(10)/convfac/mwVBS(10),bdnl(kVAP2))
          con(kVAP3) = amax1(cg(11)/convfac/mwVBS(11),bdnl(kVAP3))
          con(kVAP4) = amax1(cg(12)/convfac/mwVBS(12),bdnl(kVAP4))
          con(kVCP1) = amax1(cg(13)/convfac/mwVBS(13),bdnl(kVCP1))
          con(kVCP2) = amax1(cg(14)/convfac/mwVBS(14),bdnl(kVCP2))
          con(kVCP3) = amax1(cg(15)/convfac/mwVBS(15),bdnl(kVCP3))
          con(kVCP4) = amax1(cg(16)/convfac/mwVBS(16),bdnl(kVCP4))
          con(kVFP1) = amax1(cg(17)/convfac/mwVBS(17),bdnl(kVFP1))
          con(kVFP2) = amax1(cg(18)/convfac/mwVBS(18),bdnl(kVFP2))
          con(kVFP3) = amax1(cg(19)/convfac/mwVBS(19),bdnl(kVFP3))
          con(kVFP4) = amax1(cg(20)/convfac/mwVBS(20),bdnl(kVFP4))
          con(kPAS1) = amax1(soacon(1) ,bdnl(kPAS1))
          con(kPAS2) = amax1(soacon(2) ,bdnl(kPAS2))
          con(kPAS3) = amax1(soacon(3) ,bdnl(kPAS3))
          con(kPAS4) = amax1(soacon(4) ,bdnl(kPAS4))
          con(kPBS1) = amax1(soacon(5) ,bdnl(kPBS1))
          con(kPBS2) = amax1(soacon(6) ,bdnl(kPBS2))
          con(kPBS3) = amax1(soacon(7) ,bdnl(kPBS3))
          con(kPBS4) = amax1(soacon(8) ,bdnl(kPBS4))
          con(kPAP1) = amax1(soacon(9) ,bdnl(kPAP1))
          con(kPAP2) = amax1(soacon(10),bdnl(kPAP2))
          con(kPAP3) = amax1(soacon(11),bdnl(kPAP3))
          con(kPAP4) = amax1(soacon(12),bdnl(kPAP4))
          con(kPCP1) = amax1(soacon(13),bdnl(kPCP1))
          con(kPCP2) = amax1(soacon(14),bdnl(kPCP2))
          con(kPCP3) = amax1(soacon(15),bdnl(kPCP3))
          con(kPCP4) = amax1(soacon(16),bdnl(kPCP4))
          con(kPFP1) = amax1(soacon(17),bdnl(kPFP1))
          con(kPFP2) = amax1(soacon(18),bdnl(kPFP2))
          con(kPFP3) = amax1(soacon(19),bdnl(kPFP3))
          con(kPFP4) = amax1(soacon(20),bdnl(kPFP4))
c
c========================= Process Analysis Begin ======================
c
          if ( ldoipr ) then
            cipr(IPR_OAERO,ipa_cel,kVAS1)=cipr(IPR_OAERO,ipa_cel,kVAS1)
     &                               + (con(kVAS1)-cold(kVAS1))*convfac
            cipr(IPR_OAERO,ipa_cel,kVAS2)=cipr(IPR_OAERO,ipa_cel,kVAS2)
     &                               + (con(kVAS2)-cold(kVAS2))*convfac
            cipr(IPR_OAERO,ipa_cel,kVAS3)=cipr(IPR_OAERO,ipa_cel,kVAS3)
     &                               + (con(kVAS3)-cold(kVAS3))*convfac
            cipr(IPR_OAERO,ipa_cel,kVAS4)=cipr(IPR_OAERO,ipa_cel,kVAS4)
     &                               + (con(kVAS4)-cold(kVAS4))*convfac
            cipr(IPR_OAERO,ipa_cel,kVBS1)=cipr(IPR_OAERO,ipa_cel,kVBS1)
     &                               + (con(kVBS1)-cold(kVBS1))*convfac
            cipr(IPR_OAERO,ipa_cel,kVBS2)=cipr(IPR_OAERO,ipa_cel,kVBS2)
     &                               + (con(kVBS2)-cold(kVBS2))*convfac
            cipr(IPR_OAERO,ipa_cel,kVBS3)=cipr(IPR_OAERO,ipa_cel,kVBS3)
     &                               + (con(kVBS3)-cold(kVBS3))*convfac
            cipr(IPR_OAERO,ipa_cel,kVBS4)=cipr(IPR_OAERO,ipa_cel,kVBS4)
     &                               + (con(kVBS4)-cold(kVBS4))*convfac
            cipr(IPR_OAERO,ipa_cel,kVAP1)=cipr(IPR_OAERO,ipa_cel,kVAP1)
     &                               + (con(kVAP1)-cold(kVAP1))*convfac
            cipr(IPR_OAERO,ipa_cel,kVAP2)=cipr(IPR_OAERO,ipa_cel,kVAP2)
     &                               + (con(kVAP2)-cold(kVAP2))*convfac
            cipr(IPR_OAERO,ipa_cel,kVAP3)=cipr(IPR_OAERO,ipa_cel,kVAP3)
     &                               + (con(kVAP3)-cold(kVAP3))*convfac
            cipr(IPR_OAERO,ipa_cel,kVAP4)=cipr(IPR_OAERO,ipa_cel,kVAP4)
     &                               + (con(kVAP4)-cold(kVAP4))*convfac
            cipr(IPR_OAERO,ipa_cel,kVCP1)=cipr(IPR_OAERO,ipa_cel,kVCP1)
     &                               + (con(kVCP1)-cold(kVCP1))*convfac
            cipr(IPR_OAERO,ipa_cel,kVCP2)=cipr(IPR_OAERO,ipa_cel,kVCP2)
     &                               + (con(kVCP2)-cold(kVCP2))*convfac
            cipr(IPR_OAERO,ipa_cel,kVCP3)=cipr(IPR_OAERO,ipa_cel,kVCP3)
     &                               + (con(kVCP3)-cold(kVCP3))*convfac
            cipr(IPR_OAERO,ipa_cel,kVCP4)=cipr(IPR_OAERO,ipa_cel,kVCP4)
     &                               + (con(kVCP4)-cold(kVCP4))*convfac
            cipr(IPR_OAERO,ipa_cel,kVFP1)=cipr(IPR_OAERO,ipa_cel,kVFP1)
     &                               + (con(kVFP1)-cold(kVFP1))*convfac
            cipr(IPR_OAERO,ipa_cel,kVFP2)=cipr(IPR_OAERO,ipa_cel,kVFP2)
     &                               + (con(kVFP2)-cold(kVFP2))*convfac
            cipr(IPR_OAERO,ipa_cel,kVFP3)=cipr(IPR_OAERO,ipa_cel,kVFP3)
     &                               + (con(kVFP3)-cold(kVFP3))*convfac
            cipr(IPR_OAERO,ipa_cel,kVFP4)=cipr(IPR_OAERO,ipa_cel,kVFP4)
     &                               + (con(kVFP4)-cold(kVFP4))*convfac
            cipr(IPR_OAERO,ipa_cel,kPAS1)=cipr(IPR_OAERO,ipa_cel,kPAS1)
     &                                  + con(kPAS1)-cold(kPAS1)
            cipr(IPR_OAERO,ipa_cel,kPAS2)=cipr(IPR_OAERO,ipa_cel,kPAS2)
     &                                  + con(kPAS2)-cold(kPAS2)
            cipr(IPR_OAERO,ipa_cel,kPAS3)=cipr(IPR_OAERO,ipa_cel,kPAS3)
     &                                  + con(kPAS3)-cold(kPAS3)
            cipr(IPR_OAERO,ipa_cel,kPAS4)=cipr(IPR_OAERO,ipa_cel,kPAS4)
     &                                  + con(kPAS4)-cold(kPAS4)
            cipr(IPR_OAERO,ipa_cel,kPBS1)=cipr(IPR_OAERO,ipa_cel,kPBS1)
     &                                  + con(kPBS1)-cold(kPBS1)
            cipr(IPR_OAERO,ipa_cel,kPBS2)=cipr(IPR_OAERO,ipa_cel,kPBS2)
     &                                  + con(kPBS2)-cold(kPBS2)
            cipr(IPR_OAERO,ipa_cel,kPBS3)=cipr(IPR_OAERO,ipa_cel,kPBS3)
     &                                  + con(kPBS3)-cold(kPBS3)
            cipr(IPR_OAERO,ipa_cel,kPBS4)=cipr(IPR_OAERO,ipa_cel,kPBS4)
     &                                  + con(kPBS4)-cold(kPBS4)
            cipr(IPR_OAERO,ipa_cel,kPAP1)=cipr(IPR_OAERO,ipa_cel,kPAP1)
     &                                  + con(kPAP1)-cold(kPAP1)
            cipr(IPR_OAERO,ipa_cel,kPAP2)=cipr(IPR_OAERO,ipa_cel,kPAP2)
     &                                  + con(kPAP2)-cold(kPAP2)
            cipr(IPR_OAERO,ipa_cel,kPAP3)=cipr(IPR_OAERO,ipa_cel,kPAP3)
     &                                  + con(kPAP3)-cold(kPAP3)
            cipr(IPR_OAERO,ipa_cel,kPAP4)=cipr(IPR_OAERO,ipa_cel,kPAP4)
     &                                  + con(kPAP4)-cold(kPAP4)
            cipr(IPR_OAERO,ipa_cel,kPCP1)=cipr(IPR_OAERO,ipa_cel,kPCP1)
     &                                  + con(kPCP1)-cold(kPCP1)
            cipr(IPR_OAERO,ipa_cel,kPCP2)=cipr(IPR_OAERO,ipa_cel,kPCP2)
     &                                  + con(kPCP2)-cold(kPCP2)
            cipr(IPR_OAERO,ipa_cel,kPCP3)=cipr(IPR_OAERO,ipa_cel,kPCP3)
     &                                  + con(kPCP3)-cold(kPCP3)
            cipr(IPR_OAERO,ipa_cel,kPCP4)=cipr(IPR_OAERO,ipa_cel,kPCP4)
     &                                  + con(kPCP4)-cold(kPCP4)
            cipr(IPR_OAERO,ipa_cel,kPFP1)=cipr(IPR_OAERO,ipa_cel,kPFP1)
     &                                  + con(kPFP1)-cold(kPFP1)
            cipr(IPR_OAERO,ipa_cel,kPFP2)=cipr(IPR_OAERO,ipa_cel,kPFP2)
     &                                  + con(kPFP2)-cold(kPFP2)
            cipr(IPR_OAERO,ipa_cel,kPFP3)=cipr(IPR_OAERO,ipa_cel,kPFP3)
     &                                  + con(kPFP3)-cold(kPFP3)
            cipr(IPR_OAERO,ipa_cel,kPFP4)=cipr(IPR_OAERO,ipa_cel,kPFP4)
     &                                  + con(kPFP4)-cold(kPFP4)
          endif
c
c========================== Process Analysis End =======================
c

        else ! SOAP

          cg(1)  = con(kcg1) * convfac * mwsoap(1)
          cg(2)  = con(kcg2) * convfac * mwsoap(2)
          cg(3)  = con(kcg3) * convfac * mwsoap(3)
          cg(4)  = con(kcg4) * convfac * mwsoap(4)
          soacon(1) = con(ksoa1)
          soacon(2) = con(ksoa2)
          soacon(3) = con(ksoa3)
          soacon(4) = con(ksoa4)
          preOM = con(kpoa) / mwpoa
     &          + con(ksopa) / mwsopa
     &          + con(ksopb) / mwsopb
c
c========================= DDM Begin ===================================
c
          if (lddm) then
            call DDMSOAINI (convfac,sddm,nfam,nsen)
          endif
c
c========================== DDM End ====================================
c
          do ispc = 1, NSOAP
            csatT(ispc) = csatS(ispc) * (cstempS(ispc)/tempk)
     &                                * exp((deltahS(ispc)/8.314)
     &                                * (1/cstempS(ispc)-1/tempk))
          enddo

          call soap(NSOAP,mwsoap,soacon(:NSOAP),cg(:NSOAP),
     &              csatT(:NSOAP),preOM,
     &              iout,igrdchm,ichm,jchm,kchm)
c
c-----Calculate SOA removal by particle-phase photolysis
c
          fphot = 1. - exp( -sfacJSOA * ajno2 ) ! fraction of removed SOA by photolysis
          do ispc = 1, NSOAP
            aphot(ispc) = soacon(ispc) * fphot
            soacon(ispc) = soacon(ispc) - aphot(ispc)
          enddo
          aphot(NSOAP+1) = con(ksopa) * fphot
          con(ksopa) = con(ksopa) - aphot(NSOAP+1)
          aphot(NSOAP+2) = con(ksopb) * fphot
          con(ksopb) = con(ksopb) - aphot(NSOAP+2)
c
c-----Calculate fraction of polymerized SOA
c
          fpoly = 1. - exp( -kpolya * dtaer ) ! polymerized fraction for anthro SOA
          do ispc = 1, NSOAA
            apoly(ispc) = soacon(ispc) * fpoly
            soacon(ispc) = soacon(ispc) - apoly(ispc)
            con(ksopa) = con(ksopa) + apoly(ispc)
          enddo
          fpoly = 1. - exp( -kpolyb * dtaer ) ! polymerized fraction for bio SOA
          do ispc = NSOAA + 1, NSOAP
            apoly(ispc) = soacon(ispc) * fpoly
            soacon(ispc) = soacon(ispc) - apoly(ispc)
            con(ksopb) = con(ksopb) + apoly(ispc)
          enddo
c
c========================= DDM Begin ===================================
c
          if (lddm) then
            call DDMSOAP (dtaer,ajno2,convfac,sddm,nfam,nsen,ierr,iout)
            if (ierr.lt.0) then
              write(iout,'(a,4i4)') ' igrd,i,j,k = ',
     &                                          igrdchm,ichm,jchm,kchm
              write(iout,*) ' T,DT = ',tempk,dtaer
              write(iout,*) ' CONVFAC = ',convfac
              write(iout,*) ' CG = ',con(kcg1),con(kcg2),
     &                               con(kcg3),con(kcg4)
              write(iout,*) ' SOA = ',con(ksoa1),con(ksoa2),
     &                                con(ksoa3),con(ksoa4),
     &                                con(ksopa),con(ksopb)
              write(iout,*) ' APHOT = ',(aphot(ispc),ispc=1,NSOAP+2)
              write(iout,*) ' APOLY = ',(apoly(ispc),ispc=1,NSOAP)
              write(iout,*) ' POA = ',con(kpoa)
              do ip = 1, nfam
                write(iout,*) ' SI(',ip,') = ',
     &                               (sddm(ip,jdpCG1+is-1),is=1,NSOAP),
     &                              (sddm(ip,jdpSOA1+is-1),is=1,NSOAP),
     &                               sddm(ip,jdpSOPA),sddm(ip,jdpSOPB),
     &                                                 sddm(ip,jdpPOA)
              enddo
              call camxerr()
            endif
          endif
c
c========================== DDM End ====================================
c
          con(kcg1)  = amax1(cg(1)/convfac/mwsoap(1),bdnl(kcg1))
          con(kcg2)  = amax1(cg(2)/convfac/mwsoap(2),bdnl(kcg2))
          con(kcg3)  = amax1(cg(3)/convfac/mwsoap(3),bdnl(kcg3))
          con(kcg4)  = amax1(cg(4)/convfac/mwsoap(4),bdnl(kcg4))
          con(ksoa1) = amax1(soacon(1),bdnl(ksoa1))
          con(ksoa2) = amax1(soacon(2),bdnl(ksoa2))
          con(ksoa3) = amax1(soacon(3),bdnl(ksoa3))
          con(ksoa4) = amax1(soacon(4),bdnl(ksoa4))
          con(ksopa) = amax1(con(ksopa),bdnl(ksopa))
          con(ksopb) = amax1(con(ksopb),bdnl(ksopb))
c
c========================= Process Analysis Begin ======================
c
          if ( ldoipr ) then
            cipr(IPR_OAERO,ipa_cel,kcg1) = cipr(IPR_OAERO,ipa_cel,kcg1)
     &                                 + (con(kcg1)-cold(kcg1))*convfac
            cipr(IPR_OAERO,ipa_cel,kcg2) = cipr(IPR_OAERO,ipa_cel,kcg2)
     &                                 + (con(kcg2)-cold(kcg2))*convfac
            cipr(IPR_OAERO,ipa_cel,kcg3) = cipr(IPR_OAERO,ipa_cel,kcg3)
     &                                 + (con(kcg3)-cold(kcg3))*convfac
            cipr(IPR_OAERO,ipa_cel,kcg4) = cipr(IPR_OAERO,ipa_cel,kcg4)
     &                                 + (con(kcg4)-cold(kcg4))*convfac
            cipr(IPR_OAERO,ipa_cel,ksoa1)=cipr(IPR_OAERO,ipa_cel,ksoa1)
     &                                 + con(ksoa1)-cold(ksoa1)
            cipr(IPR_OAERO,ipa_cel,ksoa2)=cipr(IPR_OAERO,ipa_cel,ksoa2)
     &                                 + con(ksoa2)-cold(ksoa2)
            cipr(IPR_OAERO,ipa_cel,ksoa3)=cipr(IPR_OAERO,ipa_cel,ksoa3)
     &                                 + con(ksoa3)-cold(ksoa3)
            cipr(IPR_OAERO,ipa_cel,ksoa4)=cipr(IPR_OAERO,ipa_cel,ksoa4)
     &                                 + con(ksoa4)-cold(ksoa4)
            cipr(IPR_OAERO,ipa_cel,ksopa)=cipr(IPR_OAERO,ipa_cel,ksopa)
     &                                 + con(ksopa)-cold(ksopa)
            cipr(IPR_OAERO,ipa_cel,ksopb)=cipr(IPR_OAERO,ipa_cel,ksopb)
     &                                 + con(ksopb)-cold(ksopb)
          endif
c
c========================== Process Analysis End =======================
c

        endif ! VBS or SOAP?

      endif ! lsvoc?
c
c-----Do RADM aqueous chemistry if CWC is above threshold
c     All conc units must be mol/mol (mixing ratio)
c
      if (lradm) then
        pres_pa = 100.*press
        dt_sec = dtaq*3600.
        cw_kgm3 = cwc/1000.
c
        r_gas(1)  = con(kso2)*1.e-6
        r_gas(2)  = con(khno3)*1.e-6
        r_gas(3)  = con(kn2o5)*1.e-6
        r_gas(4)  = co2*1.e-6
        r_gas(5)  = con(knh3)*1.e-6
        r_gas(6)  = con(khpo_c)*1.e-6
        r_gas(7)  = con(ko3)*1.e-6
        if (kfoa_c.lt.nspec+1) then
          r_gas(8) = con(kfoa_c)*1.e-6
        else
          r_gas(8) = foa*1.e-6
        endif
        if (kmhp_c.lt.nspec+1 .or. kohp_c.lt.nspec+1) then
          r_1 = con(kmhp_c) + con(kohp_c)
          r_gas(9) = r_1*1.e-6
          if( r_1 .NE. 0. ) r_1 = con(kmhp_c) / r_1
        else
          r_gas(9) = mhp*1.e-6
          r_1 = 0.0
        endif
        if (kpaa_c.lt.nspec+1 .or. kopa_c.lt.nspec+1) then
          r_2 = con(kpaa_c) + con(kopa_c)
          r_gas(10) = r_2*1.e-6
          if( r_2 .NE. 0. ) r_2 = con(kpaa_c) / r_2
        else
          r_gas(10) = paa*1.e-6
          r_2 = 0.0
        endif
        r_gas(11) = con(ksulf)*1.e-6
c
        if (kgly.lt.nspec+1) then
          r_gas(12) = con(kgly)*1.e-6
        else
          r_gas(12) = gly*1.e-6
        endif
        if (kmgly.lt.nspec+1) then
          r_gas(13) = con(kmgly)*1.e-6
        else
          r_gas(13) = mgly*1.e-6
        endif
        if (kglyd.lt.nspec+1) then
          r_gas(14) = con(kglyd)*1.e-6
        else
          r_gas(14) = glyd*1.e-6
        endif
        r_gas(15)  = con(koh)*1.e-6
c
        r_aer(1)  = (con(kpso4)/wtmol(1)/convfac)*1.e-6
        r_aer(2)  = (con(kpnh4)/wtmol(2)/convfac)*1.e-6
        r_aer(3)  = (con(kpno3)/wtmol(3)/convfac)*1.e-6

        if (f_CA_PRM.ge.0. .and. f_CA_CRS.ge.0.) then
          r_aer(4) = ((f_CA_PRM*con(kFPRM) + f_CA_CRS*con(kFCRS))/
     &                wtmol(8)/convfac)*1.e-6
        else
          r_aer(4) = (caco3/(wtmol(8)+wtmol(7))/convfac)*1.e-6
        endif
        if (kPCA.lt.nspec+1) then
          r_aer(4) = r_aer(4) + (con(kPCA)/wtmol(8)/convfac)*1.e-6
        endif

        if (f_MG_PRM.ge.0. .and. f_MG_CRS.ge.0.) then
          r_aer(5) = ((f_MG_PRM*con(kFPRM) + f_MG_CRS*con(kFCRS))/
     &                wtmol(9)/convfac)*1.e-6
        else
          r_aer(5) = (mgco3/(wtmol(9)+wtmol(7))/convfac)*1.e-6
        endif
        if (kPMG.lt.nspec+1) then
          r_aer(5) = r_aer(5) + (con(kPMG)/wtmol(9)/convfac)*1.e-6
        endif

        if (kNA.lt.nspec+1 .and. kPCL.lt.nspec+1) then
          if (con(kNA).gt.con(kPCL)) then
            r_aer(6) = (con(kPCL)/wtmol(4)/convfac)*1.e-6
          else
            r_aer(6) = (con(kNA)/wtmol(5)/convfac)*1.e-6
          endif
        else
          r_aer(6) = (nacl/(wtmol(5)+wtmol(4))/convfac)*1.e-6
        endif

        if ( ldark ) then
          fe_iii = 0.9  ! nighttime fraction partitioning of FE(III), GS 01july2011
        else
          fe_iii = 0.1  ! daytime fraction partitioning of FE(III), GS 01july2011
        endif
        fe_sol = 0.1    ! solubility of Fe, GS 01July2011
        if (f_FE_PRM.ge.0. .and. f_FE_CRS.ge.0.) then
          fe_mmr = ((f_FE_PRM*con(kFPRM) + f_FE_CRS*con(kFCRS))/
     &               wtmol(10)/convfac)*1.e-6
          r_aer(7) = fe_iii*fe_sol*fe_mmr ! GS 01July2011
        else
          r_aer(7) = (a3fe/wtmol(10)/convfac)*1.e-6
        endif
        if (kPFE.lt.nspec+1) then
          r_aer(7) = max(r_aer(7),
     &               fe_iii*fe_sol*(con(kPFE)/wtmol(10)/convfac)*1.e-6)
        endif

        mn_ii  = 1.0    ! fraction partitioning of MN(II), GS 01July2011
        mn_sol = 0.5    ! solubility of Mn, GS 28July2011
        if (f_MN_PRM.ge.0. .and. f_MN_CRS.ge.0.) then
          mn_mmr = ((f_MN_PRM*con(kFPRM) + f_MN_CRS*con(kFCRS))/ 
     &               wtmol(11)/convfac)*1.e-6
          r_aer(8) = mn_ii*mn_sol*mn_mmr ! GS 01July2011
        else
          r_aer(8) = (b2mn/wtmol(11)/convfac)*1.e-6
        endif
        if (kPMN.lt.nspec+1) then
          r_aer(8) = max(r_aer(8),
     &               mn_ii*mn_sol*(con(kPMN)/wtmol(11)/convfac)*1.e-6)
        endif

        if (f_K_PRM.ge.0. .and. f_K_CRS.ge.0.) then
          r_aer(9) = ((f_K_PRM*con(kFPRM) + f_K_CRS*con(kFCRS))/
     &                wtmol(6)/convfac)*1.e-6
        else
          r_aer(9) = (potcl/(wtmol(6)+wtmol(4))/convfac)*1.e-6
        endif
        if (kPK.lt.nspec+1) then
          r_aer(9) = r_aer(9) + (con(kPK)/wtmol(6)/convfac)*1.e-6
        endif
c
c-----Biogenic SOA tracer
c
        csoac = con(ksoac_c)
        r_aer(10) = (csoac/wtmol(12)/convfac)*1.e-6
c
c========================= DDM Begin ===================================
c
        if (lddm) then
          do is = 1, 15
            r_old(is) = r_gas(is)
          enddo
          do is = 1, 10
            r_old(is+15) = r_aer(is)
          enddo
          call DDMRADINI (con(kso2),con(kna),con(kpcl),
     &                    f_CA_PRM,f_CA_CRS,
     &                    convfac,wtmol,sddm,nfam,nsen,NWTMOL)
        endif
c
c========================== DDM End ====================================
c
        call raqchem(tempk,pres_pa,dt_sec,cw_kgm3,r_gas,r_aer,cph,
     &               idiag,iout,igrdchm,ichm,jchm,kchm,wtmol(12))
c
c========================= DDM Begin ===================================
c
        if (lddm) then
          call DDMRADM (r_1,r_gas(9)/r_old(9),r_2,r_gas(10)/r_old(10),
     &                  convfac,wtmol,sddm,nfam,nsen,NWTMOL,ierr,iout)
          if (ierr.lt.0) then
            write(iout,'(a,4i4)') ' igrd,i,j,k = ',
     &                                          igrdchm,ichm,jchm,kchm
            write(iout,*) ' T,P,CW,DT = ',tempk,pres_pa,cw_kgm3,dt_sec
            write(iout,*) ' CONVFAC = ',convfac
            write(iout,*) ' r1,r11,r2,r22 = ',r_1,r_gas(9)/r_old(9)
     &                                       ,r_2,r_gas(10)/r_old(10)
            write(iout,*) ' R_GAS = ',(r_old(is),is=1,15)
            write(iout,*) ' R_AER = ',(r_old(is+15),is=1,10)
            do ip = 1, nfam
              if (con(kna).gt.con(kpcl)) then
                r_1 = sddm(ip,jdpPCL )/wtmol(4)/convfac *1.e-6
              else
                r_1 = sddm(ip,jdpNA  )/wtmol(5)/convfac *1.e-6
              endif
              write(iout,*) ' SI(',ip,') = ',
     &      (sddm(ip,jdpNH3 )+sddm(ip,jdpPNH4)/wtmol(2)/convfac)*1.e-6,
     &      (sddm(ip,jdpHNO3)+sddm(ip,jdpN2O5)*2.0
     &                       +sddm(ip,jdpPNO3)/wtmol(3)/convfac)*1.e-6,
     &             sddm(ip,jdpFOA )*1.e-6, 0.0, sddm(ip,jdpSO2 )*1.e-6,
     &      (sddm(ip,jdpSULF)+sddm(ip,jdpPSO4)/wtmol(1)/convfac)*1.e-6,
     &                  sddm(ip,jdpH2O2)*1.e-6, sddm(ip,jdpO3  )*1.e-6,
     &                   ( sddm(ip,jdpMHP ) + sddm(ip,jdpOHP ) )*1.e-6,
     &                   ( sddm(ip,jdpPAA ) + sddm(ip,jdpOPA ) )*1.e-6,
     &          r_1,(f_CA_PRM*sddm(ip,jdpFPRM)
     &              +f_CA_CRS*sddm(ip,jdpFCRS))/wtmol(8)/convfac*1.e-6,
     &                    sddm(ip,jdpGLY)*1.e-6,sddm(ip,jdpMGLY)*1.e-6,
     &                     sddm(ip,jdpGLYD)*1.e-6,sddm(ip,jdpHO)*1.e-6,
     &                        sddm(ip,jdpSOPB)/wtmol(12)/convfac*1.e-6
            enddo
            call camxerr()
          endif
        endif
c
c========================== DDM End ====================================
c
        con(kso2)  = amax1(r_gas(1) *1.e6,bdnl(kso2))     ! SO2 (ppm)
        con(khno3) = amax1(r_gas(2) *1.e6,bdnl(khno3))    ! HNO3 (ppm)
        con(kn2o5) = amax1(r_gas(3) *1.e6,bdnl(kn2o5))    ! N2O5 gas (ppm)
        con(knh3)  = amax1(r_gas(5) *1.e6,bdnl(knh3))     ! NH3 (ppm)
        con(khpo_c)= amax1(r_gas(6) *1.e6,bdnl(khpo_c))   ! H2O2 (ppm)
        con(ko3)   = amax1(r_gas(7) *1.e6,bdnl(ko3))      ! O3 (ppm)
!bk        con(kfoa_c)= amax1(r_gas(8) *1.e6,bdnl(kfoa_c)) ! not changed by RADM
        con(kmhp_c)= amax1(r_gas(9) *1.e6*r_1,bdnl(kmhp_c))
        con(kohp_c)= amax1(r_gas(9) *1.e6*(1.-r_1),bdnl(kohp_c))
        con(kpaa_c)= amax1(r_gas(10)*1.e6*r_2,bdnl(kpaa_c))
        con(kopa_c)= amax1(r_gas(10)*1.e6*(1.-r_2),bdnl(kopa_c))
        con(ksulf) = amax1(r_gas(11)*1.e6,bdnl(ksulf))    ! H2SO4 (ppm)
c
        con(kgly)  = amax1(r_gas(12)*1.e6,bdnl(kgly))     ! GLY (ppm)
        con(kmgly) = amax1(r_gas(13)*1.e6,bdnl(kmgly))    ! MGLY (ppm)
        con(kglyd) = amax1(r_gas(14)*1.e6,bdnl(kglyd))    ! GLYD (ppm)
c
        con(kpso4) = amax1(r_aer(1)*convfac*wtmol(1)*1.e6,bdnl(kpso4)) ! PSO4 (ug/m3)
        con(kpnh4) = amax1(r_aer(2)*convfac*wtmol(2)*1.e6,bdnl(kpnh4)) ! PNH4 (ug/m3)
        con(kpno3) = amax1(r_aer(3)*convfac*wtmol(3)*1.e6,bdnl(kpno3)) ! PNO3 (ug/m3)
c
        con(ksoac_c) = amax1(r_aer(10)*convfac*wtmol(12)*1.e6,
     &                                                  bdnl(ksoac_c)) ! SOAC (ug/m3)
c
        con(nspec+1) = 0.
c
c========================= Process Analysis Begin ======================
c
        if ( ldoipr ) then
          cipr(IPR_AQCHEM,ipa_cel,kso2 )=cipr(IPR_AQCHEM,ipa_cel,kso2 )
     &                                +(con(kso2 )-cold(kso2 ))*convfac
          cipr(IPR_AQCHEM,ipa_cel,khno3)=cipr(IPR_AQCHEM,ipa_cel,khno3)
     &                                +(con(khno3)-cold(khno3))*convfac
          cipr(IPR_AQCHEM,ipa_cel,kn2o5)=cipr(IPR_AQCHEM,ipa_cel,kn2o5)
     &                                +(con(kn2o5)-cold(kn2o5))*convfac
          cipr(IPR_AQCHEM,ipa_cel,knh3 )=cipr(IPR_AQCHEM,ipa_cel,knh3 )
     &                                +(con(knh3 )-cold(knh3 ))*convfac
          cipr(IPR_AQCHEM,ipa_cel,khpo_c)
     &                                 =cipr(IPR_AQCHEM,ipa_cel,khpo_c)
     &                              +(con(khpo_c)-cold(khpo_c))*convfac
          cipr(IPR_AQCHEM,ipa_cel,ko3  )=cipr(IPR_AQCHEM,ipa_cel,ko3  )
     &                                +(con(ko3  )-cold(ko3  ))*convfac
!bk          if (kfoa_c.lt.nspec+1) cipr(IPR_AQCHEM,ipa_cel,kfoa_c) ! not changed by RADM
!bk     &                                 =cipr(IPR_AQCHEM,ipa_cel,kfoa_c)
!bk     &                              +(con(kfoa_c)-cold(kfoa_c))*convfac
          if (kmhp_c.lt.nspec+1) cipr(IPR_AQCHEM,ipa_cel,kmhp_c)
     &                                 =cipr(IPR_AQCHEM,ipa_cel,kmhp_c)
     &                              +(con(kmhp_c)-cold(kmhp_c))*convfac
          if (kohp_c.lt.nspec+1) cipr(IPR_AQCHEM,ipa_cel,kohp_c)
     &                                 =cipr(IPR_AQCHEM,ipa_cel,kohp_c)
     &                              +(con(kohp_c)-cold(kohp_c))*convfac
          if (kpaa_c.lt.nspec+1) cipr(IPR_AQCHEM,ipa_cel,kpaa_c)
     &                                 =cipr(IPR_AQCHEM,ipa_cel,kpaa_c)
     &                              +(con(kpaa_c)-cold(kpaa_c))*convfac
          if (kopa_c.lt.nspec+1) cipr(IPR_AQCHEM,ipa_cel,kopa_c)
     &                                 =cipr(IPR_AQCHEM,ipa_cel,kopa_c)
     &                              +(con(kopa_c)-cold(kopa_c))*convfac
          cipr(IPR_AQCHEM,ipa_cel,ksulf)=cipr(IPR_AQCHEM,ipa_cel,ksulf)
     &                                +(con(ksulf)-cold(ksulf))*convfac
          if (kgly.lt.nspec+1) cipr(IPR_AQCHEM,ipa_cel,kgly)
     &                                   =cipr(IPR_AQCHEM,ipa_cel,kgly)
     &                                  +(con(kgly)-cold(kgly))*convfac
          if (kmgly.lt.nspec+1) cipr(IPR_AQCHEM,ipa_cel,kmgly)
     &                                  =cipr(IPR_AQCHEM,ipa_cel,kmgly)
     &                                +(con(kmgly)-cold(kmgly))*convfac
          if (kglyd.lt.nspec+1) cipr(IPR_AQCHEM,ipa_cel,kglyd)
     &                                  =cipr(IPR_AQCHEM,ipa_cel,kglyd)
     &                                +(con(kglyd)-cold(kglyd))*convfac
          cipr(IPR_AQCHEM,ipa_cel,kpso4)=cipr(IPR_AQCHEM,ipa_cel,kpso4)
     &                                  +con(kpso4)-cold(kpso4)
          cipr(IPR_AQCHEM,ipa_cel,kpnh4)=cipr(IPR_AQCHEM,ipa_cel,kpnh4)
     &                                  +con(kpnh4)-cold(kpnh4)
          cipr(IPR_AQCHEM,ipa_cel,kpno3)=cipr(IPR_AQCHEM,ipa_cel,kpno3)
     &                                  +con(kpno3)-cold(kpno3)
          cipr(IPR_AQCHEM,ipa_cel,ksoac_c)
     &                                =cipr(IPR_AQCHEM,ipa_cel,ksoac_c)
     &                                  +con(ksoac_c)-csoac
        endif
c
c========================== Process Analysis End =======================
c
      endif ! lradm?
c
c======================== Source Apportion Begin =======================
c
      if( ldosa ) then
         do ispc=1,ngas
           delcon(2,ispc) = delcon(2,ispc) +
     &                     AMAX1(bdnl(ispc),con(ispc))*convfac
         enddo
         do ispc=ngas+1,nspec
           delcon(2,ispc) = delcon(2,ispc) +
     &                               AMAX1(bdnl(ispc),con(ispc))
         enddo
         deln2o5 = delcon(2,kn2o5) - delcon(1,kn2o5)
         delcon(4,khno3) = delcon(4,khno3) - deln2o5 * 2.0
         delcon(4,kn2o5) = delcon(4,kn2o5) + deln2o5
         if ( lsvoc .and. .not.lvbs ) then
           delcon(5,ksoa1) = -apoly(1)
           delcon(5,ksoa2) = -apoly(2)
           delcon(5,ksoa3) = -apoly(3)
           delcon(5,ksoa4) = -apoly(4)
           do ispc = 1, NSOAA
             delcon(5,ksopa) = delcon(5,ksopa) + apoly(ispc)
           enddo
           do ispc = NSOAA+1, NSOAP
             delcon(5,ksopb) = delcon(5,ksopb) + apoly(ispc)
           enddo
           delcon(6,ksoa1) = -aphot(1)
           delcon(6,ksoa2) = -aphot(2)
           delcon(6,ksoa3) = -aphot(3)
           delcon(6,ksoa4) = -aphot(4)
           delcon(6,ksopa) = -aphot(5)
           delcon(6,ksopb) = -aphot(6)
         endif
         if ( lradm ) then
           delcon(7,ksoac_c) = con(ksoac_c) - csoac
         endif
      endif
c
c======================== Source Apportion End =======================
c
c
c-----Inorganic aerosol equilibrium chemistry with ISOROPIA
c     Convert conc units to mol/m3 (double precision)
c
      if (lisor) then
c
c========================= Process Analysis Begin ======================
c
        if ( ldoipr ) then
          csulf = con(ksulf)
          chno3 = con(khno3)
          cnh3  = con(knh3 )
          cpso4 = con(kpso4)
          cpno3 = con(kpno3)
          cpnh4 = con(kpnh4)
        endif
c
c========================== Process Analysis End =======================
c
        if ( leqsam ) then ! EQSAM

        xTT = DBLE(tempk)
        xAW = DBLE(rh/100.)

        xYPa(1) = DBLE((con(knh3) *convfac + con(kpnh4)/wtmol(2))*1.e-6) ! NH4+  [mol/m3]
        if (kNA.lt.nspec+1 .and. kPCL.lt.nspec+1) then
          xYPa(2) = DBLE(                     con(kna )/wtmol(5) *1.e-6) ! Na+   [mol/m3]
          xYMa(4) = DBLE((con(khcl)*convfac + con(kpcl)/wtmol(4))*1.e-6) ! Cl-   [mol/m3]
        else
          xYPa(2) = DBLE(                nacl/(wtmol(5)+wtmol(4))*1.e-6)
          xYMa(4) = DBLE(con(khcl)*convfac*1.e-6) + xYPa(2)
        endif
        if (f_K_PRM.ge.0. .and. f_K_CRS.ge.0.) then
          xYPa(3) = DBLE((f_K_PRM*con(kFPRM) + f_K_CRS*con(kFCRS))
     &                                                 /wtmol(6) *1.e-6) ! K+    [mol/m3]
        else
          xYPa(3) = DBLE(               potcl/(wtmol(6)+wtmol(4))*1.e-6)
          xYMa(4) = xYMa(4) + xYPa(3)
        endif
        if (kPK.lt.nspec+1) then
          xYPa(3) = xYPa(3) + DBLE(           con(kpk )/wtmol(6) *1.e-6)
        endif
        if (f_CA_PRM.ge.0. .and. f_CA_CRS.ge.0.) then
          xYPa(4) = DBLE((f_CA_PRM*con(kFPRM) + f_CA_CRS*con(kFCRS))
     &                                                 /wtmol(8) *1.e-6) ! Ca++  [mol/m3]
        else
          xYPa(4) = DBLE(               caco3/(wtmol(8)+wtmol(7))*1.e-6)
        endif
        if (kPCA.lt.nspec+1) then
          xYPa(4) = xYPa(4) + DBLE(           con(kpca)/wtmol(8) *1.e-6)
        endif
        if (f_MG_PRM.ge.0. .and. f_MG_CRS.ge.0.) then
          xYPa(5) = DBLE((f_MG_PRM*con(kFPRM) + f_MG_CRS*con(kFCRS))
     &                                                 /wtmol(9) *1.e-6) ! Mg++  [mol/m3]
        else
          xYPa(5) = DBLE(               mgco3/(wtmol(9)+wtmol(7))*1.e-6)
        endif
        if (kPMG.lt.nspec+1) then
          xYPa(5) = xYPa(5) + DBLE(           con(kpmg)/wtmol(9) *1.e-6)
        endif
        xYMa(1) = DBLE((con(ksulf)*convfac + con(kpso4)/wtmol(1))*1.e-6) ! SO4-- [mol/m3]
        xYMa(2) = 0.0_dp                                                 ! HSO4- [mol/m3]
        xYMa(3) = DBLE((con(khno3)*convfac + con(kpno3)/wtmol(3))*1.e-6) ! NO3-  [mol/m3]

        xWH2O = 0.0_dp ! set to zero to make the EQSAM metastable flag effective

        wt_save(2) = xYMa(1)
        wt_save(3) = xYPa(1)
        wt_save(4) = xYMa(3)
        wt_save(5) = xYMa(4)

        call EQSAM4clim(xTT,xAW,xsPM,xaPM,xPMS,xPMT,xRHO,xVOL,xPH,xHp,
     &                  xGF,xWH2O,xYPa,xYPs,xYMa,xYMs,xYG,
     &                  imask,1,1,jmeq,Dd)

        con(ksulf) = bdnl(ksulf)
        con(kpso4) = wt_save(2) - DBLE(con(ksulf)*1.e-6*convfac)
        con(kpso4) = amax1(con(kpso4)*wtmol(1)*1.e6,bdnl(kpso4))

        con(khno3) = REAL(xYG(1))*1.e6/convfac
        con(khno3) = amax1(con(khno3),bdnl(khno3))
        con(kpno3) = wt_save(4) - DBLE(con(khno3)*1.e-6*convfac)
        con(kpno3) = amax1(con(kpno3)*wtmol(3)*1.e6,bdnl(kpno3))

        con(knh3 ) = REAL(xYG(3))*1.e6/convfac
        con(knh3 ) = amax1(con(knh3),bdnl(knh3))
        con(kpnh4) = wt_save(3) - DBLE(con(knh3)*1.e-6*convfac)
        con(kpnh4) = amax1(con(kpnh4)*wtmol(2)*1.e6,bdnl(kpnh4))

        if (khcl.lt.nspec+1) then
          con(khcl ) = REAL(xYG(2))*1.e6/convfac
          con(khcl ) = amax1(con(khcl),bdnl(khcl))
        endif
        if (kpcl.lt.nspec+1) then
          con(kpcl ) = wt_save(5) - DBLE(con(khcl)*1.e-6*convfac)
          con(kpcl ) = amax1(con(kpcl)*wtmol(4)*1.e6,bdnl(kpcl))
        endif

        con(kph2o) = amax1(REAL(xWH2O),bdnl(kph2o))

        aph = REAL(xPH)
        aph = MIN(MAX(aph,-2.0),10.0)

        else ! ISORROPIA

        rhi = DBLE(rh/100.)
        tempi = DBLE(tempk)
        cntrl(1) = 0.d0                   ! 0 = forward problem
        cntrl(2) = 1.d0                   ! 1 = metastable liquid state

        if (kNA.lt.nspec+1 .and. kPCL.lt.nspec+1) then
          wi(1) = DBLE(con(kna)/wtmol(5)*1.e-6)                        ! total sodium
          wi(5) = DBLE((con(khcl)*convfac + con(kpcl)/wtmol(4))*1.e-6) ! total chloride
        else
          wi(1) = DBLE(nacl/(wtmol(5)+wtmol(4))*1.e-6)                 ! total sodium
          wi(5) = DBLE(con(khcl)*convfac*1.e-6) + wi(1)                ! total chloride
        endif
        wi(2) = DBLE((con(ksulf)*convfac + con(kpso4)/wtmol(1))*1.e-6) ! total sulfate
        wi(3) = DBLE((con(knh3) *convfac + con(kpnh4)/wtmol(2))*1.e-6) ! total ammonium
        wi(4) = DBLE((con(khno3)*convfac + con(kpno3)/wtmol(3))*1.e-6) ! total nitrate

        if (f_CA_PRM.ge.0. .and. f_CA_CRS.ge.0.) then
          dno3 = DBLE((f_CA_PRM*con(kFPRM) + f_CA_CRS*con(kFCRS))
     &                                               /wtmol(8) *1.e-6)
        else
          dno3 = DBLE(caco3/(wtmol(8)+wtmol(7))*1.e-6)
        endif
        if (kPCA.lt.nspec+1) then
          dno3 = dno3 + DBLE(con(kPCA)/wtmol(8)*1.e-6) ! half of CaCO3 is assumed to be replaced by Ca(NO3)2          
        endif
        if (dno3.lt.wi(4)) then
          wi(4) = wi(4) - dno3 ! exclude calcium nitrate
          lno3_isr = .true.
        else ! all nitrate is calcium nitrate
          con(khno3) = bdnl(khno3)
          con(kpno3) = amax1( (REAL(wi(4))*1.e6 - con(khno3)*convfac)
     &                                  * wtmol(3), bdnl(kpno3) )
          wi(4) = 0.d0
          lno3_isr = .false.
        endif
c
        do is = 1, 5
          wt_save(is) = wi(is)
        enddo
c
c========================= DDM Begin ===================================
c
        if (lddm) then
          call DDMISRINI (lno3_isr,f_CA_PRM,f_CA_CRS,
     &                    convfac,wtmol,sddm,nfam,nsen,NWTMOL,iout)
        endif
c
c========================== DDM End ====================================
c
        call isoropia(wi,rhi,tempi,cntrl,wt,gasis,aerliq,aersld,
     &                scasi,other)
c
c========================= DDM Begin ===================================
c
        if (lddm) then
          call DDMISRPIA (lno3_isr,f_CA_PRM,f_CA_CRS,
     &                    convfac,wtmol,sddm,nfam,nsen,NWTMOL,ierr,iout)
          if (ierr.lt.0) then
            write(iout,'(a,4i4)') ' igrd,i,j,k = ',
     &                                          igrdchm,ichm,jchm,kchm
            write(iout,*) ' RH,T = ',rhi,tempi
            write(iout,*) ' CONVFAC = ',convfac
            write(iout,*) ' WI = ',(wt_save(is),is=1,5)
            write(iout,*) ' LNO3_ISR = ',lno3_isr,f_CA_PRM,f_CA_CRS,
     &                     (sddm(ip,jdpFPRM),sddm(ip,jdpFCRS),ip=1,nfam)
            do ip = 1, nfam
              if (lno3_isr) then
                r_1 = (          sddm(ip,jdpHNO3)*convfac
     &                +          sddm(ip,jdpPNO3)/wtmol(3)
     &                -(f_CA_PRM*sddm(ip,jdpFPRM)
     &                 +f_CA_CRS*sddm(ip,jdpFCRS))/wtmol(8))*1.e-6
              else
                r_1 = 0.0
              endif
              write(iout,*) ' SI(',ip,') = ',
     &                                sddm(ip,jdpNA  )/wtmol(5) *1.e-6,
     &      (sddm(ip,jdpSULF)*convfac+sddm(ip,jdpPSO4)/wtmol(1))*1.e-6,
     &      (sddm(ip,jdpNH3 )*convfac+sddm(ip,jdpPNH4)/wtmol(2))*1.e-6,
     &                                                             r_1,
     &      (sddm(ip,jdpHCL )*convfac+sddm(ip,jdpPCL )/wtmol(4))*1.e-6
            enddo
            call camxerr()
          endif
        endif
c
c========================== DDM End ====================================
c 
c-----Load results back to local CON array (ppm for gas, ug/m3 for aerosol)
c 
        con(ksulf) = bdnl(ksulf)                                  ! sulfuric acid gas
        con(kpso4) = wt_save(2) - DBLE(con(ksulf)*1.e-6*convfac)  ! sulfate
        con(kpso4) = amax1(con(kpso4)*wtmol(1)*1.e6,bdnl(kpso4))
c
        if (lno3_isr) then
          con(khno3) = REAL(gasis(2))*1.e6/convfac                ! nitric acid gas
          con(khno3) = amax1(con(khno3),bdnl(khno3))    
          con(kpno3) = wt_save(4) - DBLE(con(khno3)*1.e-6*convfac)! nitrate
          con(kpno3) = con(kpno3) + dno3                          ! add calcium nitrate back
          con(kpno3) = amax1(con(kpno3)*wtmol(3)*1.e6,bdnl(kpno3))
        endif
c
        con(knh3 ) = REAL(gasis(1))*1.e6/convfac                  ! ammonia gas
        con(knh3 ) = amax1(con(knh3),bdnl(knh3))    
        con(kpnh4) = wt_save(3) - DBLE(con(knh3)*1.e-6*convfac)   ! ammonium
        con(kpnh4) = amax1(con(kpnh4)*wtmol(2)*1.e6,bdnl(kpnh4))
c
        if (khcl.lt.nspec+1) then
          con(khcl) = REAL(gasis(3))*1.e6/convfac                 ! HCl gas
          con(khcl) = amax1(con(khcl),bdnl(khcl))
        endif
        if (kpcl.lt.nspec+1) then
          con(kpcl) = wt_save(5) - DBLE(con(khcl)*1.e-6*convfac)  ! chloride
          con(kpcl) = amax1(con(kpcl)*wtmol(4)*1.e6,bdnl(kpcl))
        endif
c
        con(kph2o) = REAL(aerliq(8))                              ! aerosol water
        con(kph2o) = amax1(con(kph2o)*18.*1.e6,bdnl(kph2o))
c
        aph = LOG10(REAL(aerliq(8)*18.D-3/MAX(aerliq(1),1.D-20))) ! aerosol pH
        aph = MIN(MAX(aph,-2.0),10.0)

        endif ! leqsam?
c
c========================= Process Analysis Begin ======================
c
        if ( ldoipr ) then
          cipr(IPR_IAERO,ipa_cel,ksulf) = cipr(IPR_IAERO,ipa_cel,ksulf)
     &                                  + (con(ksulf)-csulf)*convfac
          cipr(IPR_IAERO,ipa_cel,khno3) = cipr(IPR_IAERO,ipa_cel,khno3)
     &                                  + (con(khno3)-chno3)*convfac
          cipr(IPR_IAERO,ipa_cel,knh3 ) = cipr(IPR_IAERO,ipa_cel,knh3 )
     &                                  + (con(knh3 )-cnh3 )*convfac
          cipr(IPR_IAERO,ipa_cel,kpso4) = cipr(IPR_IAERO,ipa_cel,kpso4)
     &                                  + con(kpso4)-cpso4
          cipr(IPR_IAERO,ipa_cel,kpno3) = cipr(IPR_IAERO,ipa_cel,kpno3)
     &                                  + con(kpno3)-cpno3
          cipr(IPR_IAERO,ipa_cel,kpnh4) = cipr(IPR_IAERO,ipa_cel,kpnh4)
     &                                  + con(kpnh4)-cpnh4
          if (khcl.lt.nspec+1)
     &      cipr(IPR_IAERO,ipa_cel,khcl) = cipr(IPR_IAERO,ipa_cel,khcl)
     &                                 + (con(khcl)-cold(khcl))*convfac
          if (kpcl.lt.nspec+1)
     &      cipr(IPR_IAERO,ipa_cel,kpcl) = cipr(IPR_IAERO,ipa_cel,kpcl)
     &                                   + con(kpcl)-cold(kpcl)
          cipr(IPR_IAERO,ipa_cel,kph2o) = cipr(IPR_IAERO,ipa_cel,kph2o)
     &                                  + con(kph2o)-cold(kph2o)
        endif
c
c========================== Process Analysis End =======================
c
      endif ! lisor?
c
      return
c
      end
