      subroutine khetero(tcell,pcell,H2O,con,lddm,nfam,sddm)
      use chmstry
      use filunit
      implicit none
c
c----CAMx v7Beta6 190902
c
c     KHETERO estimates the rate constants for heterogeneous reactions of N2O5:
c       N2O5 + H2O -> 2 HNO3
c       N2O5 + HCL -> CLN2 + HNO3
c     KHETERO also sets the rate constant for hydrolysis of organic nitrate:
c       NTR2       -> HNO3
c       INTR       -> HNO3
c     KHETERO also sets the rate constant for heterogeneous reaction of SO2:
c       SO2        -> SULF
c
c     Parameterization of N2O5 on aqueous particles:
c     Bertram and Thornton (2009) Atmos. Chem. Phys., 9, 8351-8363.
c
c     Mean molecular speed formula - Seinfeld and Pandis (1998)
c     Atmospheric Chemistry and Physics, p453.
c
c     NTR2 hydrolysis rate:
c     Liu et al. (2012) Aerosol Sci. Technol., 46, 1359-1369.
c
c     INTR hydrolysis rate:
c     Fisher et al. (2016) Atmos. Chem. Phys., 16, 5969-5991.
c
c     Organic nitrate partitioning ratio:
c     Rollins et al. (2013) J. Geophys. Res., 118, 6651-6662.
c
c     Heterogeneous SO2 reaction rate:
c     Zheng et al. (2015) Atmos. Chem. Phys., 15, 2031-2049.
c
c     Molecular diffusion coefficient:
c     Tang et al. (2014) Atmos. Chem. Phys., 14, 9233-9247.
c     (https://sites.google.com/site/mingjintang/home/diffusion)
c
c     Copyright 1996 - 2018
c     Ramboll
c
c     Modifications:
c        01/30/2006  - bkoo -     Original development
c        07/04/2007  - bkoo -     Modified to use pointer to N2O5 hydrolysis
c                                 Extended for CMU_AERO
c        03/29/2011  -cemery-     Revised to allow for inert PM with gas-phase
c                                 chemistry
c        03/28/2012  -cemery-     Revised to set rate of 0.25/hr in the
c                                 absence of PM chemistry
c        09/21/2013  - bkoo -     Revised for DDM-PM
c        10/14/2013  - bkoo -     Added hydrolysis of organic nitrate
c        01/03/2014  - bkoo -     Updated DDM-PM for organic nitrate hydrolysis
c        03/18/2014  - bkoo -     Updated for benzene SOA
c        12/07/2014  - bkoo -     Modified for VBS
c                                 Cap khetN2O5 sens
c        01/08/2016  - bkoo -     Updated for Revised VBS
c        02/29/2016  - bkoo -     Updated heterogeneous rxn of N2O5
c        07/20/2016  - bkoo -     Added heterogeneous hydrolysis of INTR
c        08/25/2016  - bkoo -     Updated for new SOAP
c        09/02/2016  - bkoo -     Added het SO2 -> SULF; revised het rxn pointers
c        11/28/2016  - bkoo -     Revised N2O5 uptake rxn formula
c
c     Input arguments:
c        tcell          temperature (K)
c        pcell          pressure (mb)
c        H2O            water vapor concentration (ppm)
c        con            species concentrations (ppm or ug/m3)
c        lddm           logical flag for DDM sensitivities
c        nfam           number of sensitivity families
c        sddm           sensitivity matrix (ppm; ug/m3; 1/hr or 1/ppm-hr)
c                          sddm(:,nspec+1:nspec+2) [1/ppm-hr]
c                          sddm(:,nspec+3:nspec+5) [1/hr]
c
c     Output arguments:
c        sddm
c
c     Routines called:
c        NONE
c
c     Called by:
c        CHEMDRIV
c
      include 'camx.prm'
      include 'camx_aero.inc'
      include 'vbs.inc'
c
c-----Arguments
c
      real  tcell, pcell, H2O
      real  con(*)
c
c-----Local variables:
c
      real, parameter :: khon0  = 0.16667       ! [1/hr]; lifetime = 6 hr (particle-phase NTR2 hydrolysis)
      real, parameter :: khon20 = 1.0           ! [1/hr]; lifetime = 1 hr (particle-phase INTR hydrolysis)
      real, parameter :: alfa1  = 0.34,         ! partioning parameters for total ON
     &                   alfa2  = 0.66,
     &                   cstar1 = 0.73,
     &                   cstar2 = 1000.
      real, parameter :: fntr2  = 0.71          ! NTR2 / total ON
      real, parameter :: fintr  = 0.08          ! INTR / total ON

      real :: fcond                             ! condensable ON / total ON
      real :: khon                              ! psuedo gas-phase NTR2 hydrolysis rate constant [1/hr]
      real :: khon2                             ! psuedo gas-phase INTR hydrolysis rate constant [1/hr]
      real :: fp                                ! fraction of NTR2 or INTR in the particle phase
                                                ! (assuming the same fp for NTR2 and INTR)
      real :: ctoa                              ! total OA concentration [ug/m3]

      real :: qwatr, ev, es, rh                 ! relative humidity calculation

      real, parameter :: Rgc    = 8.314         ! gas constant [J/mol-K]
      real, parameter :: Pi     = 3.1415927     ! Pi
      real, parameter :: mwN2O5 = 0.108         ! MW of N2O5 [kg/mol]
      real, parameter :: mwSO2  = 0.064         ! MW of SO2 [kg/mol]
      real, parameter :: rtau   = 1.66667e-3    ! reciprocal of minimum lifetime of N2O5 (10 min = 600 sec) [1/s]

      real, parameter :: Afac   = 3.2e-8        ! empirical pre-factor for gamma [s]
      real, parameter :: beta   = 1.15e6        ! parameter for k2f' [1/s]
      real, parameter :: delta  = 1.3e-1        ! parameter for k2f' [1/M]
      real, parameter :: k3k2b  = 6.0e-2        ! k3/k2b
      real, parameter :: k4k2b  = 29.0          ! k4/k2b
      real, parameter :: molarPH2O = 55.        ! molar concentration of aerosol water [M]; assume diluted aqueous phase

      real, parameter :: diff_N2O5 = 65.        ! Diffusivity of N2O5 at 296 K [Torr-cm2/s]
      real, parameter :: diff_SO2 = 94.         ! Diffusivity of SO2 at 296 K [Torr-cm2/s]
      real, parameter :: g1_SO2 = 2.e-5         ! lower limit of SO2 uptake coefficient
      real, parameter :: g2_SO2 = 5.e-5         ! upper limit of SO2 uptake coefficient

      real, parameter :: h2s = 3600.0           ! hour to second

      real :: vt                                ! total dry volume of fine particles
      real :: diam                              ! diameter of a fine particle [um]
      real :: nden                              ! number density of fine particles [#/cm3_air]
      real :: saden                             ! surface area density [um2/cm3_air]
      real :: cbarN2O5                          ! mean molecular speed of N2O5 [m/s]
      real :: cbarSO2                           ! mean molecular speed of SO2 [m/s]
      real :: k2fprime                          ! k2f' = k2f * [H2O] [1/s]
      real :: dc_N2O5                           ! diffusion coefficient of N2O5 [m2/s]
      real :: dc_SO2                            ! diffusion coefficient of SO2 [m2/s]
      real :: gamma_N2O5                        ! uptake coefficient of N2O5
      real :: gamma_SO2                         ! uptake coefficient of SO2
      real :: y_clno2                           ! product yield for ClNO2

      real*8 :: r1, r2, rdenom
      real :: rscl
      real :: khet0, khet1, khet2, khet5

      integer :: l, isempty
      integer :: kwtr, isec, iaero
c
c======================== DDM Begin =======================
c
      logical :: lddm
      integer :: nfam
      integer :: ip
      real :: sddm(nfam,MXSPEC+NHETRXN)
      real :: stoa,sfp
      real :: const0,sad0
      real :: dvt,dsad,dr1r
c
c======================== DDM End =======================
c
c-----No PM chemistry:
c       N2O5 + H2O -> 2 HNO3      : reset to the IUPAC 2006 rate
c       N2O5 + HCL -> CLN2 + HNO3 : use the rate set in the CHEMPARAM input file
c       NTR2       -> HNO3        : use the rate set in the CHEMPARAM input file
c       INTR       -> HNO3        : use the rate set in the CHEMPARAM input file
c       SO2        -> SULF        : use the rate set in the CHEMPARAM input file
c
      if (aeropt.eq.'NONE' .or. aeropt.eq.'INERT') then
        rk(ihetrxn(1)) = 2.5e-22 * 8.87e16    ! [cm3/molecule-sec] -> [1/ppm-hr]
        return
      endif
c
c-----Calculate RH
c
      qwatr = 1.e-6*H2O*18./28.8
      ev = qwatr*pcell/(qwatr + eps)
      es = e0*exp((lv/rv)*(1./273. - 1./tcell))
      rh = amin1(1.,ev/es)
c
c-----Hydrolysis of organic nitrate in the organic solution phase
c
      if (ihetrxn(3).gt.0) then
        if (aeropt.eq.'CF') then
          ctoa = con(kpoa)  + con(ksoa1) + con(ksoa2)
     &                      + con(ksoa3) + con(ksoa4)
     &                      + con(ksopa) + con(ksopb)
          if (lvbs) then
            ctoa = con(kpas0)+con(kpas1)+con(kpas2)+con(kpas3)+con(kpas4)
     &            +con(kpbs0)+con(kpbs1)+con(kpbs2)+con(kpbs3)+con(kpbs4)
     &            +con(kpap0)+con(kpap1)+con(kpap2)+con(kpap3)+con(kpap4)
     &            +con(kpcp0)+con(kpcp1)+con(kpcp2)+con(kpcp3)+con(kpcp4)
     &            +con(kpfp0)+con(kpfp1)+con(kpfp2)+con(kpfp3)+con(kpfp4)
          endif
        elseif (aeropt.eq.'CMU') then
          ctoa = 0.0
          do isec = 1, nbin
            ctoa = ctoa + con(kpoa_1+isec-1)
     &                  + con(ksoa1_1+isec-1) + con(ksoa2_1+isec-1)
     &                  + con(ksoa3_1+isec-1) + con(ksoa4_1+isec-1)
     &                  + con(ksopa_1+isec-1) + con(ksopb_1+isec-1)
          enddo
        else
          write(iout,'(//,a)') 'ERROR in KHETERO: invalid AEROPT'
          call camxerr()
        endif
        fp = 0.0
        if (ctoa.ge.0.0001) then
          fp = alfa1 / (1. + cstar1/ctoa) + alfa2 / (1. + cstar2/ctoa)
          fcond = fntr2
          if (ihetrxn(4).gt.0) fcond = fcond + fintr
          fp = fp / fcond
        endif
        if (rh.lt.0.2) then
          rk(ihetrxn(3)) = 0.0
          if (ihetrxn(4).gt.0) rk(ihetrxn(4)) = 0.0
        elseif (rh.ge.0.4) then
          rk(ihetrxn(3)) = fp * khon0
          if (ihetrxn(4).gt.0) rk(ihetrxn(4)) = fp * khon20
        else
          rk(ihetrxn(3)) = fp * khon0 * (rh - 0.2) / 0.2
          if (ihetrxn(4).gt.0) rk(ihetrxn(4)) = fp * khon20
     &                                * (rh - 0.2) / 0.2
        endif
c
c======================== DDM Begin =======================
c
        if ( lddm .and. aeropt.eq.'CF' ) then
          do ip = 1, nfam
            stoa = sddm(ip,kpoa)
     &           + sddm(ip,ksoa1) + sddm(ip,ksoa2)
     &           + sddm(ip,ksoa3) + sddm(ip,ksoa4)
     &           + sddm(ip,ksopa) + sddm(ip,ksopb)
            sfp = 0.0
            if (ctoa.ge.0.0001) then
              sfp = ( fp - ( 
     &                alfa1 / ((1. + cstar1/ctoa)*(1. + cstar1/ctoa))
     &              + alfa2 / ((1. + cstar2/ctoa)*(1. + cstar2/ctoa))
     &                      ) / fcond ) * stoa / ctoa
            endif
            if (rh.lt.0.2) then
              sddm(ip,nspec+3) = 0.0
              if (ihetrxn(4).gt.0) sddm(ip,nspec+4) = 0.0
            elseif (rh.ge.0.4) then
              sddm(ip,nspec+3) = sfp * khon0
              if (ihetrxn(4).gt.0) sddm(ip,nspec+4) = sfp * khon20
            else
              sddm(ip,nspec+3) = sfp * khon0 * (rh - 0.2) / 0.2
              if (ihetrxn(4).gt.0) sddm(ip,nspec+4) = sfp * khon20
     &                                       * (rh - 0.2) / 0.2
            endif
          enddo
        endif
c
c======================== DDM End =======================
c
      endif
c
c-----Heterogeneous rxn of N2O5 / SO2 on aerosol surface area
c
      khet1 = 0.0
      khet2 = 0.0
      khet5 = 0.0

      k2fprime = beta * ( 1. - exp(-delta*molarPH2O) )      ! [1/s]
      cbarN2O5 = sqrt( 8.0*Rgc*tcell / (Pi*mwN2O5) )        ! [m/s]
      dc_N2O5 = (diff_N2O5/7.6e6) * (1013.25/pcell) * (tcell/296.)**1.75 ! [m2/s]

      cbarSO2 = sqrt( 8.0*Rgc*tcell / (Pi*mwSO2) )          ! [m/s]
      dc_SO2 = (diff_SO2/7.6e6) * (1013.25/pcell) * (tcell/296.)**1.75 ! [m2/s]
      if (rh.lt.0.5) then
        gamma_SO2 = g1_SO2
      else
        gamma_SO2 = g1_SO2 + (g2_SO2 - g1_SO2)*(rh - 0.5)/0.5
      endif
c
c-----CF PM chemistry
c
      if (aeropt.eq.'CF') then
        diam = sqrt(dcut(kph2o,1)*dcut(kph2o,2))            ! geometric mean dry diameter [um]
        !bk - limit hetero SO2 rxn only at dust particle surface
        saden = (6./diam) * ( con(kfprm)/roprt(kfprm)       ! dust particle surface area [m2/m3_air]
     &                      + con(kfcrs)/roprt(kfcrs) 
     &                      + con(kpfe) /roprt(kfcrs) 
     &                      + con(kpmn) /roprt(kfcrs) 
     &                      + con(kpca) /roprt(kfcrs) 
     &                      + con(kpmg) /roprt(kfcrs) 
     &                      + con(kpk)  /roprt(kfcrs) 
     &                      + con(kpal) /roprt(kfcrs) 
     &                      + con(kpsi) /roprt(kfcrs) 
     &                      + con(kpti) /roprt(kfcrs) )
        khet5 = saden / ( diam * 0.5e-6 / dc_SO2 +
     &                    4. / ( cbarSO2 * gamma_SO2 ) )    ! [1/s]
c
c-----Adjust dry diameter to wet diameter for N2O5 hydrolysis
c
        if (con(kph2o).le.bdnl(kph2o)) goto 800

        isempty = 1
        vt = 0.
        do l = ngas+1,nspec
          if (dcut(l,2).lt.(dcut(kph2o,2)+1.e-5)) then      ! fine
            if (l.ne.kph2o) then                            ! dry
              if (con(l).gt.bdnl(l)) isempty = 0
              vt = vt + con(l)/roprt(l)                     ! [ug/m3_air]/[g/m3]
            endif
          endif
        enddo
        if (isempty.eq.1) goto 800

        nden = vt / (diam*diam*diam*Pi/6.0)                 ! [ug/m3_air]/[um3-g/m3] = [#/m3_air] * 1.E-12
c
c======================== DDM Begin =======================
c
        if ( lddm ) then !bk - broken
        endif
c
c======================== DDM End =======================
c
        diam = diam * ( 1.0 +
     &         con(kph2o)/roprt(kph2o)/vt )**0.33333        ! wet diameter [um]
        saden = Pi*diam*diam*nden                           ! [m2/m3_air]
        r1 = k3k2b * (con(kph2o)/18.)/(con(kpno3)/62.)
        r2 = k4k2b * (con(kpcl)/35.)/(con(kpno3)/62.)
        rdenom = r1 + 1. + r2
        gamma_N2O5 = Afac * k2fprime * (1. - 1. / rdenom)
        y_clno2 = r2 / (r1 + r2)
        khet0 = saden / ( diam * 0.5e-6 / dc_N2O5 +
     &                    4. / ( cbarN2O5 * gamma_N2O5 ) )  ! [1/s]
        khet0 = MIN( khet0, rtau )                          ! check against upper limit [1/s]
        khet1 = (1. - y_clno2) * khet0
        khet2 = y_clno2 * khet0
c
c======================== DDM Begin =======================
c
        if ( lddm ) then !bk - broken
        endif
c
c======================== DDM End =======================
c
c-----CMU PM chemistry
c
      elseif (aeropt.eq.'CMU') then
        kwtr = (kph2o_1 - ngas)/nbin + 1
        if (nbin.eq.1) kwtr = kph2o_1 - ngas
        do isec = 1, nbin
          diam = sqrt(dcut(ngas+isec,1)*dcut(ngas+isec,2))
          !bk - limit hetero SO2 rxn only at dust particle surface
          saden = (6./diam) * con(kcrst_1+isec-1)/roprt(kcrst_1+isec-1)
          khet5 = khet5 + saden / ( diam * 0.5e-6 / dc_SO2 +
     &                              4. / ( cbarSO2 * gamma_SO2 ) )

          if (con(kph2o_1+isec-1).le.bdnl(kph2o_1+isec-1)) cycle
          isempty = 1
          vt = 0.
          do iaero = 1, naero
            if (iaero.ne.kwtr) then
              l = ngas + (iaero - 1)*nbin + isec
              if (con(l).gt.bdnl(l)) isempty = 0
              vt = vt + con(l)/roprt(l)
            endif
          enddo
          if (isempty.eq.1) cycle

          nden = vt / (diam*diam*diam*Pi/6.0)
          diam = diam * ( 1.0 +
     &           con(kph2o_1+isec-1)/roprt(kph2o_1+isec-1)/vt )**0.33333
          saden = Pi*diam*diam*nden

          r1 = k3k2b * (con(kph2o_1+isec-1)/18.)
     &               / (con(kpno3_1+isec-1)/62.)
          r2 = k4k2b * (con(kpcl_1+isec-1)/35.)
     &               / (con(kpno3_1+isec-1)/62.)
          rdenom = r1 + 1. + r2
          gamma_N2O5 = Afac * k2fprime * (1. - 1. / rdenom)
          y_clno2 = r2 / (r1 + r2)

          khet0 = saden / ( diam * 0.5e-6 / dc_N2O5 +
     &                      4. / ( cbarN2O5 * gamma_N2O5 ) )
          khet1 = khet1 + (1. - y_clno2) * khet0
          khet2 = khet2 + y_clno2 * khet0
        enddo
        rscl = MIN( 1., rtau / (khet1 + khet2) )
        khet1 = rscl * khet1
        khet2 = rscl * khet2
      else
        write(iout,'(//,a)') 'ERROR in KHETERO: invalid AEROPT'
        call camxerr()
      endif
c
c --- adjust rk(N2O5); rk(N2O5) will be reset by KTHERM every timestep
c
 800  continue
      rk(ihetrxn(1)) = rk(ihetrxn(1)) + khet1 * h2s / H2O   ! add gas-phase rate; [1/ppm-hr]
      if (ihetrxn(2).gt.0) rk(ihetrxn(2)) = khet2 * h2s / con(khcl)
      if (ihetrxn(5).gt.0) rk(ihetrxn(5)) = khet5 * h2s     ! [1/hr]

      return

      end
