c*** PSATSA
c
      subroutine psatsa(numcol,numrow,numlay,igrid,icell,jcell,kcell,
     &                                         nspc,delcon,nprec,yfac)
      use grid
      use chmstry
      use tracer
      implicit none
c
c----CAMx v7Beta6 190902
c
c-----------------------------------------------------------------------
c   Description:
c     This routine makes the "chemistry" adjustments to the tracer 
c     species for PSAT.  The adjustments are based on the differences in 
c     concentrations of the regular model species before and after
c     the regular model chemistry.  This is essentially an adjustment
c     for production or decay.
c
c     Copyright 1996 - 2018
c     Ramboll 
c
c   Argument descriptions:
c     Inputs:
c       numcol  I number of columns in this slice
c       numrow  I number of rows in this slice
c       numlay  I number of layers in this slice
c       igrid     I  grid number
c       icell     I  the X grid location of current cell
c       jcell     I  the X grid location of current cell
c       kcell     I  the vertical grid location of current layer
c       nspc      I  number of model species
c       delcon    R  array of change in concentrations total
c       nprec     I  number of SOA precursors
c       yfac      R  array of weighting factors for SOA
c
c-----------------------------------------------------------------------
c   LOG:
c-----------------------------------------------------------------------
c
c     09/28/03   --gwilson--  Original development
c     12/29/06   --bkoo--     Revised for the updated SOA scheme
c     05/04/07   --gwilson--  Changed call to cyctpnsa so that conversion
c                             factor is set for lower bound values -
c                             the conversion is passed here from chemdriv
c     03/18/14   --bkoo--     Revised for benzene SOA
c     08/25/16   --bkoo--     Updated for new SOAP
c     10/16/17   --bkoo--     Fixed DELARO that was missing IVOA changes
c                             Updated for SOA photolysis
c     01/12/18   --bkoo--     Removed BNZA/TOLA/XYLA/ISP/TRP
c     01/09/19   --cemery--   Added DMS
c     06/18/19   --bkoo--     Updated delARO/delTRP/delSQT for SAPRC07T
c
c-----------------------------------------------------------------------
c   Include files:
c-----------------------------------------------------------------------
c
      include 'camx.prm'
c
c-----------------------------------------------------------------------
c    Argument declarations:
c-----------------------------------------------------------------------
c
      integer   numcol
      integer   numrow
      integer   numlay
      integer   igrid
      integer   icell
      integer   jcell
      integer   kcell
      integer   nspc
      real      delcon(7,nspc+1)
      integer   nprec
      real      yfac(3,nprec)
c
c-----------------------------------------------------------------------
c    Local variables:
c-----------------------------------------------------------------------
c
      integer   icg1, icg2, isopa, iaro, iisp, itrp, isqt
      integer   iso2, idms, ihg0, i
      integer*8 idxcel, idx, jdx, kdx, ldx, idx0, idx2
      real      delARO, delISP, delTRP, delSQT, delso2, delps4, deldms, delhg2
      real      delCG1, delCG2, delNV1
      real      sumaro, sumisp, sumtrp, sumsqt, sumso2, sumps4, sumdms
      real      sumhg0, sumhg2

      real, parameter :: wtsopb = 220.0 ! must be consistent with MWSOPB in SOAP.INC
c
c-----------------------------------------------------------------------
c    Entry point:
c-----------------------------------------------------------------------
c
c   --- calculate the index of the cell in the grid ---
c
      idxcel =  DBLE(ipsa3d(igrid)-1) + DBLE(icell) + 
     &            DBLE(numcol) * DBLE(jcell-1) + 
     &                 DBLE(numcol) * DBLE(numrow) * DBLE(kcell-1)
c
c----------------------------------------------------------------------
c  ---  SOA tracers ----
c----------------------------------------------------------------------
c
      if( lsoa ) then
         delARO = delcon(1,kbenz)                   ! [umol/m3] for gases
     &          + delcon(1,ktol)  + delcon(1,ktolu) + delcon(1,karo1)
     &          + delcon(1,kxyl)  + delcon(1,koxyl) + delcon(1,kmxyl)
     &          + delcon(1,kpxyl) + delcon(1,kb124) + delcon(1,karo2)
     &          + delcon(1,kivoa)
         delISP = delcon(1,kisop)
         delTRP = delcon(1,kterp) + delcon(1,kapin)
         delSQT = delcon(1,ksqt)  + delcon(1,ksesq)
         delCG1 = delcon(1,kcg1)
         delCG2 = delcon(1,kcg2)
         delNV1 = delcon(3,ksopa) - delcon(5,ksopa) ! adjust non-vol SOA delta [ug/m3] to exclude
     &                            - delcon(6,ksopa) ! changes due to SOA polymerization/photolysis
c
c  --- get the sum of the tracer species for ARO/CG1 yield ---
c
         icg1 = iptcls(idxipt(ITRCG1)) - 1
         sumaro = 0.
         do i=iptcls(idxipt(ITRARO)),nptcls(idxipt(ITRARO))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                              DBLE(numlay) * DBLE(i-1)
            icg1 = icg1 + 1
            sumaro = sumaro + ptconc(idx) *
     &                      ( yhrates(icg1)*yfac(1,1)
     &                      + ylrates(icg1)*yfac(2,1) )
         enddo
c
c  --- make adjustment to CG1 tracers based on change in CG1 and
c      yield rates for ARO ---
c
         iaro = iptcls(idxipt(ITRARO)) - 1
         do i=iptcls(idxipt(ITRCG1)),nptcls(idxipt(ITRCG1))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                 DBLE(numlay) * DBLE(i-1)
            iaro = iaro + 1
            jdx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                               DBLE(numlay) * DBLE(iaro-1)
            ptconc(idx) = ptconc(idx) + delCG1 * ( ptconc(jdx) *
     &                                  ( yhrates(i)*yfac(1,1)
     &                                  + ylrates(i)*yfac(2,1) ) )
     &                                                    / sumaro
            ptconc(idx) = MAX(BNDLPT,ptconc(idx))
         enddo
c
c  --- get the sum of the tracer species for ARO/CG2 yield ---
c
         icg2 = iptcls(idxipt(ITRCG2)) - 1
         sumaro = 0.
         do i=iptcls(idxipt(ITRARO)),nptcls(idxipt(ITRARO))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                    DBLE(numlay) * DBLE(i-1)
            icg2 = icg2 + 1
            sumaro = sumaro + ptconc(idx) *
     &                      ( yhrates(icg2)*yfac(1,1)
     &                      + ylrates(icg2)*yfac(2,1) )
         enddo
c
c  --- make adjustment to CG2 tracers based on change in CG2 and
c      yield rates for ARO ---
c
         iaro = iptcls(idxipt(ITRARO)) - 1
         do i=iptcls(idxipt(ITRCG2)),nptcls(idxipt(ITRCG2))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                DBLE(numlay) * DBLE(i-1)
            iaro = iaro + 1
            jdx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                             DBLE(numlay) * DBLE(iaro-1)
            ptconc(idx) = ptconc(idx) + delCG2 * ( ptconc(jdx) *
     &                                  ( yhrates(i)*yfac(1,1)
     &                                  + ylrates(i)*yfac(2,1) ) )
     &                                                    / sumaro
            ptconc(idx) = MAX(BNDLPT,ptconc(idx))
         enddo
c
c  --- get the sum of the tracer species for ARO/SOPA yield ---
c      SPECIAL CASE: non-volatile CG; ARO -> SOPA directly (skip CG)
c
         isopa = iptcls(idxipt(ITRPPA)) - 1
         sumaro = 0.
         do i=iptcls(idxipt(ITRARO)),nptcls(idxipt(ITRARO))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                DBLE(numlay) * DBLE(i-1)
            isopa = isopa + 1
            sumaro = sumaro + ptconc(idx) *
     &                      ( yhrates(isopa)*yfac(1,1)
     &                      + ylrates(isopa)*yfac(2,1) )
         enddo
c
c  --- make adjustment to SOPA tracers based on change in SOPA and
c      yield rates for ARO ---
c      SPECIAL CASE: non-volatile CG; ARO -> SOPA directly (skip CG)
c
         iaro = iptcls(idxipt(ITRARO)) - 1
         do i=iptcls(idxipt(ITRPPA)),nptcls(idxipt(ITRPPA))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                DBLE(numlay) * DBLE(i-1)
            iaro = iaro + 1
            jdx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                             DBLE(numlay) * DBLE(iaro-1)
            ptconc(idx) = ptconc(idx) + delNV1 * ( ptconc(jdx) *
     &                                  ( yhrates(i)*yfac(1,1)
     &                                  + ylrates(i)*yfac(2,1) ) )
     &                                                    / sumaro
            ptconc(idx) = MAX(BNDLPT,ptconc(idx))
         enddo
c
c  --- get the sum of the tracer species for ARO ---
c
         sumaro = 0.
         do i=iptcls(idxipt(ITRARO)),nptcls(idxipt(ITRARO))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                  DBLE(numlay) * DBLE(i-1)
            sumaro = sumaro + ptconc(idx)*wtkoh(i)
         enddo
c
c  --- make adjustment to ARO based on change in aromatic precursors and kOH ---
c
         do i=iptcls(idxipt(ITRARO)),nptcls(idxipt(ITRARO))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                  DBLE(numlay) * DBLE(i-1)
            ptconc(idx) = ptconc(idx) + delARO * 
     &                                  ptconc(idx)*wtkoh(i) / sumaro
            ptconc(idx) = MAX(BNDLPT,ptconc(idx))
         enddo
c
c  --- get the sum of the tracer species for ISP/TRP/SQT ---
c
         sumisp = 0.
         do i=iptcls(idxipt(ITRISP)),nptcls(idxipt(ITRISP))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                DBLE(numlay) * DBLE(i-1)
            sumisp = sumisp + ptconc(idx)
         enddo

         sumtrp = 0.
         do i=iptcls(idxipt(ITRTRP)),nptcls(idxipt(ITRTRP))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                              DBLE(numlay) * DBLE(i-1)
            sumtrp = sumtrp + ptconc(idx)
         enddo

         sumsqt = 0.
         do i=iptcls(idxipt(ITRSQT)),nptcls(idxipt(ITRSQT))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                              DBLE(numlay) * DBLE(i-1)
            sumsqt = sumsqt + ptconc(idx)
         enddo
c
c  --- make adjustment to CG3 tracers based on change in ISP/TRP/SQT and
c      distribution of ISP/TRP/SQT ----
c
         iisp = iptcls(idxipt(ITRISP)) - 1
         itrp = iptcls(idxipt(ITRTRP)) - 1
         isqt = iptcls(idxipt(ITRSQT)) - 1
         do i=iptcls(idxipt(ITRCG3)),nptcls(idxipt(ITRCG3))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                DBLE(numlay) * DBLE(i-1)
            iisp = iisp + 1
            jdx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                             DBLE(numlay) * DBLE(iisp-1)
            itrp = itrp + 1
            kdx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                             DBLE(numlay) * DBLE(itrp-1)
            isqt = isqt + 1
            ldx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                             DBLE(numlay) * DBLE(isqt-1)
            ptconc(idx) = ptconc(idx) - delISP * yfac(1,5) *
     &                                           ptconc(jdx) / sumisp
     &                                - delTRP * yfac(1,6) *
     &                                           ptconc(kdx) / sumtrp
     &                                - delSQT * yfac(1,7) *
     &                                           ptconc(ldx) / sumsqt
            ptconc(idx) = MAX(BNDLPT,ptconc(idx))
         enddo
c
c  --- make adjustment to CG4 tracers based on change in ISP/TRP/SQT and
c      distribution of ISP/TRP/SQT ----
c
         iisp = iptcls(idxipt(ITRISP)) - 1
         itrp = iptcls(idxipt(ITRTRP)) - 1
         isqt = iptcls(idxipt(ITRSQT)) - 1
         do i=iptcls(idxipt(ITRCG4)),nptcls(idxipt(ITRCG4))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                               DBLE(numlay) * DBLE(i-1)
            iisp = iisp + 1
            jdx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                            DBLE(numlay) * DBLE(iisp-1)
            itrp = itrp + 1
            kdx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                            DBLE(numlay) * DBLE(itrp-1)
            isqt = isqt + 1
            ldx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                            DBLE(numlay) * DBLE(isqt-1)
            ptconc(idx) = ptconc(idx) - delISP * yfac(2,5) *
     &                                           ptconc(jdx) / sumisp
     &                                - delTRP * yfac(2,6) *
     &                                           ptconc(kdx) / sumtrp
     &                                - delSQT * yfac(2,7) *
     &                                           ptconc(ldx) / sumsqt
            ptconc(idx) = MAX(BNDLPT,ptconc(idx))
         enddo
c
c  --- make adjustment to SOPB tracers based on change in ISP/TRP/SQT and
c      distribution of ISP/TRP/SQT ----
c      SPECIAL CASE: non-volatile CG; ISP/TRP/SQT -> SOPB directly (skip CG)
c
         iisp = iptcls(idxipt(ITRISP)) - 1
         itrp = iptcls(idxipt(ITRTRP)) - 1
         isqt = iptcls(idxipt(ITRSQT)) - 1
         do i=iptcls(idxipt(ITRPPB)),nptcls(idxipt(ITRPPB))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                  DBLE(numlay) * DBLE(i-1)
            iisp = iisp + 1
            jdx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                               DBLE(numlay) * DBLE(iisp-1)
            itrp = itrp + 1
            kdx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                               DBLE(numlay) * DBLE(itrp-1)
            isqt = isqt + 1
            ldx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                               DBLE(numlay) * DBLE(isqt-1)
            ptconc(idx) = ptconc(idx) - wtsopb * 
     &                                ( delISP * yfac(3,5) *
     &                                           ptconc(jdx) / sumisp
     &                                + delTRP * yfac(3,6) *
     &                                           ptconc(kdx) / sumtrp
     &                                + delSQT * yfac(3,7) *
     &                                           ptconc(ldx) / sumsqt )
            ptconc(idx) = MAX(BNDLPT,ptconc(idx))
         enddo
c
c  --- make adjustment to ISP based on change in ISP ---
c
         do i=iptcls(idxipt(ITRISP)),nptcls(idxipt(ITRISP))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                 DBLE(numlay) * DBLE(i-1)
            ptconc(idx) = ptconc(idx) + delISP * ptconc(idx) / sumisp
            ptconc(idx) = MAX(BNDLPT,ptconc(idx))
         enddo
c
c  --- make adjustment to TRP based on change in TRP ---
c
         do i=iptcls(idxipt(ITRTRP)),nptcls(idxipt(ITRTRP))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                     DBLE(numlay) * DBLE(i-1)
            ptconc(idx) = ptconc(idx) + delTRP * ptconc(idx) / sumtrp
            ptconc(idx) = MAX(BNDLPT,ptconc(idx))
         enddo
c
c  --- make adjustment to SQT based on change in SQT ---
c
         do i=iptcls(idxipt(ITRSQT)),nptcls(idxipt(ITRSQT))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                    DBLE(numlay) * DBLE(i-1)
            ptconc(idx) = ptconc(idx) + delSQT * ptconc(idx) / sumsqt
            ptconc(idx) = MAX(BNDLPT,ptconc(idx))
         enddo
      endif
c
c----------------------------------------------------------------------
c  ---  Sulfate tracers ---
c----------------------------------------------------------------------
c
      if( lsulfate ) then
         delso2 = delcon(3,kso2)
         delps4 = delcon(3,kpso4)
c
c  --- get the sum of the tracer species for SO2 ---
c
         sumso2 = 0.
         do i=iptcls(idxipt(ITRSO2)),nptcls(idxipt(ITRSO2))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                   DBLE(numlay) * DBLE(i-1)
            sumso2 = sumso2 + ptconc(idx)
         enddo
c
c  --- get the sum of the tracer species for DMS ---
c
         if( ldmschm ) then
            idms = iptcls(idxipt(ITRDMS)) - 1
            deldms = delcon(1,kdms)
            sumdms = 0.
            do i=iptcls(idxipt(ITRDMS)),nptcls(idxipt(ITRDMS))
               idx = idxcel + numcol * numrow * numlay * (i-1)
               sumdms = sumdms + ptconc(idx)
            enddo
         endif
c
c  --- if PSO4 change is positive, make adjustment to PS4 
c      tracers based on distribtion of SO2 ---
c
         if( delps4 .GT. 0. ) then
            iso2 = iptcls(idxipt(ITRSO2)) - 1
            do i=iptcls(idxipt(ITRPS4)),nptcls(idxipt(ITRPS4))
               idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                  DBLE(numlay) * DBLE(i-1)
               iso2 = iso2 + 1
               jdx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                               DBLE(numlay) * DBLE(iso2-1)
               ptconc(idx) = ptconc(idx) + delps4 * ptconc(jdx) / sumso2
               ptconc(idx) = MAX(BNDLPT,ptconc(idx))
            enddo
c
c  --- if PSO4 change is negative, make adjustment to PS4 
c      tracers based on distribtion of PS4 ---
c
         else
c
c  --- get the sum of the tracer species for PS4 ---
c
            sumps4 = 0.
            do i=iptcls(idxipt(ITRPS4)),nptcls(idxipt(ITRPS4))
               idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                   DBLE(numlay) * DBLE(i-1)
               sumps4 = sumps4 + ptconc(idx)
            enddo
            do i=iptcls(idxipt(ITRPS4)),nptcls(idxipt(ITRPS4))
               idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                    DBLE(numlay) * DBLE(i-1)
               ptconc(idx) = ptconc(idx) + delps4 * ptconc(idx) / sumps4
               ptconc(idx) = MAX(BNDLPT,ptconc(idx))
            enddo
         endif
c
c  --- make adjustment to SO2 tracers based on change in SO2 (and DMS) ---
c
         do i=iptcls(idxipt(ITRSO2)),nptcls(idxipt(ITRSO2))
            idx = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                  DBLE(numlay) * DBLE(i-1)
            ptconc(idx) = ptconc(idx) + delso2 * ptconc(idx) / sumso2
            if( ldmschm ) then
               idms = idms + 1
               jdx = idxcel + numcol * numrow * numlay * (idms-1)
               ptconc(idx) = ptconc(idx) - deldms * ptconc(jdx) / sumdms
            endif
            ptconc(idx) = MAX(BNDLPT,ptconc(idx))
         enddo
c
c  --- make adjustment to DMS tracers based on change in DMS ---
c
         if( ldmschm ) then
           do i=iptcls(idxipt(ITRDMS)),nptcls(idxipt(ITRDMS))
              idx = idxcel + numcol * numrow * numlay * (i-1)
              ptconc(idx) = ptconc(idx) + deldms * ptconc(idx) / sumdms
              ptconc(idx) = MAX(BNDLPT,ptconc(idx))
           enddo
         endif
      endif
c
c----------------------------------------------------------------------
c  ---  Mercury tracers ----
c----------------------------------------------------------------------
c
      if( lmercury ) then
         delhg2 = delcon(3,khg2)

c  --- get the sum of the tracer species ---
c
         sumhg0 = 0.
         do i=iptcls(idxipt(ITRHG0)),nptcls(idxipt(ITRHG0))
            idx0 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                   DBLE(numlay) * DBLE(i-1)
            sumhg0 = sumhg0 + ptconc(idx0)
         enddo
         sumhg2 = 0.
         do i=iptcls(idxipt(ITRHG2)),nptcls(idxipt(ITRHG2))
            idx2 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                   DBLE(numlay) * DBLE(i-1)
            sumhg2 = sumhg2 + ptconc(idx2)
         enddo
c
c  --- HG0 being oxidized to HG2 ----
c
         if( delhg2 .GT. 0 ) then
            ihg0 = iptcls(idxipt(ITRHG0)) - 1
            do i=iptcls(idxipt(ITRHG2)),nptcls(idxipt(ITRHG2))
               idx2 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                   DBLE(numlay) * DBLE(i-1)
               ihg0 = ihg0 + 1
               idx0 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                DBLE(numlay) * DBLE(ihg0-1)
               ptconc(idx2) = ptconc(idx2) + 
     &                                   delhg2 * ptconc(idx0) / sumhg0
               ptconc(idx2) = MAX(BNDLPT,ptconc(idx2))
               ptconc(idx0) = ptconc(idx0) - 
     &                                   delhg2 * ptconc(idx2) / sumhg2
               ptconc(idx0) = MAX(BNDLPT,ptconc(idx0))
            enddo
c
c  --- HG2 beign reduced to HG0 ----
c
         else
            ihg0 = iptcls(idxipt(ITRHG0)) - 1
            do i=iptcls(idxipt(ITRHG2)),nptcls(idxipt(ITRHG2))
               idx2 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                      DBLE(numlay) * DBLE(i-1)
               ihg0 = ihg0 + 1
               idx0 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                   DBLE(numlay) * DBLE(ihg0-1)
               ptconc(idx2) = ptconc(idx2) + 
     &                                  delhg2 * ptconc(idx2) / sumhg2
               ptconc(idx2) = MAX(BNDLPT,ptconc(idx2))
               ptconc(idx0) = ptconc(idx0) - 
     &                                  delhg2 * ptconc(idx0) / sumhg0
               ptconc(idx0) = MAX(BNDLPT,ptconc(idx0))
            enddo
         endif
      endif
c
c-----------------------------------------------------------------------
c    Return point:
c-----------------------------------------------------------------------
c
      return
      end
