c*** REPARTSA.F
c
      subroutine repartsa(numcol,numrow,numlay,igrid,
     &                                  icell,jcell,kcell,modcon,delcon)
      use grid
      use chmstry
      use tracer
c
c----CAMx v7Beta6 190902
c
c-----------------------------------------------------------------------
c   Description:
c     This routine repartitions the particulate and gaseous portions
c     nitric acid and ammonia tracers based on the ratio found in the
c     regular model.
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       modcon  R  array of model species concentrations
c       delcon  R  array of change in model species concentrations
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/05/07   --gwilson--  Added code to make an adjustment to
c                             tracer species if there is stray from 
c                             the regular - added to handle mass errors
c                             in regular model chemistry routines
c     08/25/16   --bkoo--     Updated for new SOAP
c     09/16/16   --bkoo--     Updated for in-cloud SOA formation
c     10/16/17   --bkoo--     Updated for SOA photolysis
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
      real    modcon(*)
      real    delcon(7,MXSPEC+1)
c
c-----------------------------------------------------------------------
c    Local parameters:
c-----------------------------------------------------------------------
c
c    FUZZ - fuzz value for comparing against lower bound
c
      real FUZZ
c
      parameter( FUZZ = 10.0 )
c
c-----------------------------------------------------------------------
c    Local variables:
c-----------------------------------------------------------------------
c
      integer*8 idxcel, idxhno3, idxpno3, idxnh3, idxpnh4
      integer*8 idxpo1, idxpo2, idxpo3, idxpo4
      integer*8 idxcg1, idxcg2, idxcg3, idxcg4
      integer*8 idxppa, idxppb, idxisp
      integer   i, ipno3, ipnh4, ipo1, ipo2, ipo3, ipo4
      integer   ippa, ippb
      real      wthno3, wtnh3, conhno3, conpno3, connh3, conpnh4
      real      concg1, concg2, concg3, concg4
      real      conpo1, conpo2, conpo3, conpo4, conppa, conppb
      real      wtcg1, wtcg2, wtcg3, wtcg4
      real      ratio, sumtrac, contrac, conadjust
      real      pratio1, pratio2
      real      sumisp, delppb
c
c-----------------------------------------------------------------------
c    Data statements:
c     Note: MWs for CG species should be consistent with those defined
c           in SOAP.INC
c-----------------------------------------------------------------------
c
      data wtnh3  /18.0/
      data wthno3 /62.0/
      data wtcg1 /150.0/
      data wtcg2 /150.0/
      data wtcg3 /180.0/
      data wtcg4 /180.0/
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 * numrow *(kcell-1))
c
      if( lnitrate ) then
c
c  --- calculate the ratio of gas to  particluate nitric acid ---
c
         conhno3 = modcon(khno3) * wthno3
         conpno3 = modcon(kpno3)
         ratio = 0.
         if( conhno3+conpno3 .GT. 0. )
     &                    ratio = conhno3 / (conpno3 + conhno3)
c
c  --- get the sum of all tracer species ---
c
         contrac = 0.0
         ipno3 = iptcls(idxipt(ITRPN3)) - 1
         do i=iptcls(idxipt(ITRHN3)),nptcls(idxipt(ITRHN3))
            idxhno3 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                  DBLE(numlay) * DBLE(i-1)
            ipno3 = ipno3 + 1
            idxpno3 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                              DBLE(numlay) * DBLE(ipno3-1)
            sumtrac = ptconc(idxhno3)*wthno3 + ptconc(idxpno3)
            if( sumtrac .GT. FUZZ*BNDLPT ) contrac = contrac + sumtrac
         enddo
         if( contrac .GT. 0. ) then
             conadjust = (conhno3+conpno3) / contrac
         else
             conadjust = 1.0
         endif
c
c  --- adjust the nitric acid tracers ----
c
         ipno3 = iptcls(idxipt(ITRPN3)) - 1
         do i=iptcls(idxipt(ITRHN3)),nptcls(idxipt(ITRHN3))
            idxhno3 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                  DBLE(numlay) * DBLE(i-1)
            ipno3 = ipno3 + 1
            idxpno3 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                              DBLE(numlay) * DBLE(ipno3-1)
c
c  --- adjust for mass imbalance, if needed ---
c
            sumtrac = ptconc(idxhno3)*wthno3 + ptconc(idxpno3)
            if( sumtrac .GT. FUZZ*BNDLPT ) sumtrac = conadjust * sumtrac
c
c  --- redistribute the gas and particulates ---
c
            ptconc(idxhno3) = (ratio * sumtrac ) / wthno3
            ptconc(idxhno3) = MAX(ptconc(idxhno3),BNDLPT)
            ptconc(idxpno3) = sumtrac - ptconc(idxhno3) * wthno3
            ptconc(idxpno3) = MAX(ptconc(idxpno3),BNDLPT)
         enddo
c
c  --- calculate the ratio of gas to particluate ammonia ---
c
         connh3 = modcon(knh3) * wtnh3
         conpnh4 = modcon(kpnh4)
         ratio = 0.
         if( connh3+conpnh4 .GT. 0. )
     &                    ratio = connh3 / (conpnh4 + connh3)
c
c  --- get the sum of all tracer species ---
c
         contrac = 0.0
         ipnh4 = iptcls(idxipt(ITRPN4)) - 1
         do i=iptcls(idxipt(ITRNH3)),nptcls(idxipt(ITRNH3))
            idxnh3 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                 DBLE(numlay) * DBLE(i-1)
            ipnh4 = ipnh4 + 1
            idxpnh4 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                             DBLE(numlay) * DBLE(ipnh4-1)
            sumtrac = ptconc(idxnh3)*wtnh3 + ptconc(idxpnh4)
            if( sumtrac .GT. FUZZ*BNDLPT ) contrac = contrac + sumtrac
         enddo
         if( contrac .GT. 0. ) then
             conadjust = (connh3+conpnh4) / contrac
         else
             conadjust = 1.0
         endif
c
c  --- adjust the ammonia tracers ----
c
         ipnh4 = iptcls(idxipt(ITRPN4)) - 1
         do i=iptcls(idxipt(ITRNH3)),nptcls(idxipt(ITRNH3))
            idxnh3 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                     DBLE(numlay) * DBLE(i-1)
            ipnh4 = ipnh4 + 1
            idxpnh4 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                 DBLE(numlay) * DBLE(ipnh4-1)	
c        
c  --- adjust for mass imbalance, if needed ---
c  
            sumtrac = ptconc(idxnh3)*wtnh3 + ptconc(idxpnh4)
            if( sumtrac .GT. FUZZ*BNDLPT ) sumtrac = conadjust * sumtrac
c
c  --- redistribute the gas and particulates ---
c
            ptconc(idxnh3) = ratio * sumtrac / wtnh3
            ptconc(idxnh3) = MAX(ptconc(idxnh3),BNDLPT)
            ptconc(idxpnh4) = sumtrac - ptconc(idxnh3) * wtnh3
            ptconc(idxpnh4) = MAX(ptconc(idxpnh4),BNDLPT)
         enddo
      endif
c
c  --- secondary aerosols ---
c
      if( lsoa ) then
c
c  --- apportion photolyzed SOPA & SOPB first
c
         conppa = modcon(ksopa) - delcon(5,ksopa) ! get equil con by reverting the changes
     &                          - delcon(6,ksopa) ! due to SOA polymerization & photolysis
         pratio1 = 0.
         if ( conppa .GT. 0. ) pratio1 = -delcon(6,ksopa) / conppa

         do i=iptcls(idxipt(ITRPPA)),nptcls(idxipt(ITRPPA))
           idxppa = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                   DBLE(numlay) * DBLE(i-1)
           ptconc(idxppa) = ptconc(idxppa) - pratio1*ptconc(idxppa)
         enddo

         conppb = modcon(ksopb) - delcon(5,ksopb) ! get equil con by reverting the changes
     &                          - delcon(6,ksopb) ! due to SOA polymerization, photolysis,
     &                          - delcon(7,ksopb) ! and in-cloud SOA formation
         pratio1 = 0.
         if ( conppb .GT. 0. ) pratio1 = -delcon(6,ksopb) / conppb

         do i=iptcls(idxipt(ITRPPB)),nptcls(idxipt(ITRPPB))
           idxppb = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                     DBLE(numlay) * DBLE(i-1)
           ptconc(idxppb) = ptconc(idxppb) - pratio1*ptconc(idxppb)
         enddo
c
c  --- calculate the ratio of CG1 to SOA1 ---
c
         concg1 = modcon(kcg1) * wtcg1
         conpo1 = modcon(ksoa1) - delcon(5,ksoa1) ! get equil con by reverting the changes
     &                          - delcon(6,ksoa1) ! due to SOA polymerization & photolysis
         ratio = 0.
         if( concg1+conpo1 .GT. 0. )
     &                    ratio = concg1 / (conpo1 + concg1)
         pratio1 = 0.
         if ( conpo1 .GT. 0. ) pratio1 = -delcon(6,ksoa1) / conpo1
         conpo1 = conpo1 + delcon(6,ksoa1)
         pratio2 = 0.
         if ( conpo1 .GT. 0. ) pratio2 = -delcon(5,ksoa1) / conpo1
c
c  --- adjust the CG1/PO1 ratio ---
c
         ipo1 = iptcls(idxipt(ITRPO1)) - 1
         ippa = iptcls(idxipt(ITRPPA)) - 1
         do i=iptcls(idxipt(ITRCG1)),nptcls(idxipt(ITRCG1))
            idxcg1 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                     DBLE(numlay) * DBLE(i-1)
            ipo1 = ipo1 + 1
            ippa = ippa + 1
            idxpo1 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                   DBLE(numlay) * DBLE(ipo1-1)
            idxppa = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                   DBLE(numlay) * DBLE(ippa-1)
            sumtrac = ptconc(idxcg1)*wtcg1 + ptconc(idxpo1)
            ptconc(idxcg1) = (ratio * sumtrac ) / wtcg1              ! by partitioning
            ptconc(idxpo1) = (1.0-ratio) * sumtrac                   ! by partitioning
            ptconc(idxpo1) = (1.0-pratio1) * ptconc(idxpo1)          ! by photolysis
            ptconc(idxppa) = ptconc(idxppa) + pratio2*ptconc(idxpo1) ! by polymerization
            ptconc(idxpo1) = (1.0-pratio2) * ptconc(idxpo1)          ! by polymerization
            ptconc(idxcg1) = MAX(ptconc(idxcg1),BNDLPT)
            ptconc(idxpo1) = MAX(ptconc(idxpo1),BNDLPT)
            ptconc(idxppa) = MAX(ptconc(idxppa),BNDLPT)
         enddo
c
c  --- calculate the ratio of CG2 to SOA2 ---
c
         concg2 = modcon(kcg2) * wtcg2
         conpo2 = modcon(ksoa2) - delcon(5,ksoa2) ! get equil con by reverting the changes
     &                          - delcon(6,ksoa2) ! due to SOA polymerization & photolysis
         ratio = 0.
         if( concg2+conpo2 .GT. 0. )
     &                    ratio = concg2 / (conpo2 + concg2)
         pratio1 = 0.
         if ( conpo2 .GT. 0. ) pratio1 = -delcon(6,ksoa2) / conpo2
         conpo2 = conpo2 + delcon(6,ksoa2)
         pratio2 = 0.
         if ( conpo2 .GT. 0. ) pratio2 = -delcon(5,ksoa2) / conpo2
c
c  --- adjust the CG2/PO2 ratio ---
c
         ipo2 = iptcls(idxipt(ITRPO2)) - 1
         ippa = iptcls(idxipt(ITRPPA)) - 1
         do i=iptcls(idxipt(ITRCG2)),nptcls(idxipt(ITRCG2))
            idxcg2 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                DBLE(numlay) * DBLE(i-1)
            ipo2 = ipo2 + 1
            ippa = ippa + 1
            idxpo2 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                  DBLE(numlay) * DBLE(ipo2-1)
            idxppa = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                  DBLE(numlay) * DBLE(ippa-1)
            sumtrac = ptconc(idxcg2)*wtcg2 + ptconc(idxpo2)
            ptconc(idxcg2) = (ratio * sumtrac ) / wtcg2              ! by partitioning
            ptconc(idxpo2) = (1.0-ratio) * sumtrac                   ! by partitioning
            ptconc(idxpo2) = (1.0-pratio1) * ptconc(idxpo2)          ! by photolysis
            ptconc(idxppa) = ptconc(idxppa) + pratio2*ptconc(idxpo2) ! by polymerization
            ptconc(idxpo2) = (1.0-pratio2) * ptconc(idxpo2)          ! by polymerization
            ptconc(idxcg2) = MAX(ptconc(idxcg2),BNDLPT)
            ptconc(idxpo2) = MAX(ptconc(idxpo2),BNDLPT)
            ptconc(idxppa) = MAX(ptconc(idxppa),BNDLPT)
         enddo
c
c  --- calculate the ratio of CG3 to SOA3 ---
c
         concg3 = modcon(kcg3) * wtcg3
         conpo3 = modcon(ksoa3) - delcon(5,ksoa3) ! get equil con by reverting the changes
     &                          - delcon(6,ksoa3) ! due to SOA polymerization & photolysis
         ratio = 0.
         if( concg3+conpo3 .GT. 0. )
     &                    ratio = concg3 / (conpo3 + concg3)
         pratio1 = 0.
         if ( conpo3 .GT. 0. ) pratio1 = -delcon(6,ksoa3) / conpo3
         conpo3 = conpo3 + delcon(6,ksoa3)
         pratio2 = 0.
         if ( conpo3 .GT. 0. ) pratio2 = -delcon(5,ksoa3) / conpo3
c
c  --- adjust the CG3/PO3 ratio ---
c
         ipo3 = iptcls(idxipt(ITRPO3)) - 1
         ippb = iptcls(idxipt(ITRPPB)) - 1
         do i=iptcls(idxipt(ITRCG3)),nptcls(idxipt(ITRCG3))
            idxcg3 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                 DBLE(numlay) * DBLE(i-1)
            ipo3 = ipo3 + 1
            ippb = ippb + 1
            idxpo3 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                DBLE(numlay) * DBLE(ipo3-1)
            idxppb = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                DBLE(numlay) * DBLE(ippb-1)
            sumtrac = ptconc(idxcg3)*wtcg3 + ptconc(idxpo3)
            ptconc(idxcg3) = (ratio * sumtrac ) / wtcg3              ! by partitioning
            ptconc(idxpo3) = (1.0-ratio) * sumtrac                   ! by partitioning
            ptconc(idxpo3) = (1.0-pratio1) * ptconc(idxpo3)          ! by photolysis
            ptconc(idxppb) = ptconc(idxppb) + pratio2*ptconc(idxpo3) ! by polymerization
            ptconc(idxpo3) = (1.0-pratio2) * ptconc(idxpo3)          ! by polymerization
            ptconc(idxcg3) = MAX(ptconc(idxcg3),BNDLPT)
            ptconc(idxpo3) = MAX(ptconc(idxpo3),BNDLPT)
            ptconc(idxppb) = MAX(ptconc(idxppb),BNDLPT)
         enddo
c
c  --- calculate the ratio of CG4 to SOA4 ---
c
         concg4 = modcon(kcg4) * wtcg4
         conpo4 = modcon(ksoa4) - delcon(5,ksoa4) ! get equil con by reverting the changes
     &                          - delcon(6,ksoa4) ! due to SOA polymerization & photolysis
         ratio = 0.
         if( concg4+conpo4 .GT. 0. )
     &                    ratio = concg4 / (conpo4 + concg4)
         pratio1 = 0.
         if ( conpo4 .GT. 0. ) pratio1 = -delcon(6,ksoa4) / conpo4
         conpo4 = conpo4 + delcon(6,ksoa4)
         pratio2 = 0.
         if ( conpo4 .GT. 0. ) pratio2 = -delcon(5,ksoa4) / conpo4
c
c  --- adjust the CG4/PO4 ratio ---
c
         ipo4 = iptcls(idxipt(ITRPO4)) - 1
         ippb = iptcls(idxipt(ITRPPB)) - 1
         do i=iptcls(idxipt(ITRCG4)),nptcls(idxipt(ITRCG4))
            idxcg4 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                DBLE(numlay) * DBLE(i-1)
            ipo4 = ipo4 + 1
            ippb = ippb + 1
            idxpo4 = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                               DBLE(numlay) * DBLE(ipo4-1)
            idxppb = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                               DBLE(numlay) * DBLE(ippb-1)
            sumtrac = ptconc(idxcg4)*wtcg4 + ptconc(idxpo4)
            ptconc(idxcg4) = (ratio * sumtrac ) / wtcg4              ! by partitioning
            ptconc(idxpo4) = (1.0-ratio) * sumtrac                   ! by partitioning
            ptconc(idxpo4) = (1.0-pratio1) * ptconc(idxpo4)          ! by photolysis
            ptconc(idxppb) = ptconc(idxppb) + pratio2*ptconc(idxpo4) ! by polymerization
            ptconc(idxpo4) = (1.0-pratio2) * ptconc(idxpo4)          ! by polymerization
            ptconc(idxcg4) = MAX(ptconc(idxcg4),BNDLPT)
            ptconc(idxpo4) = MAX(ptconc(idxpo4),BNDLPT)
            ptconc(idxppb) = MAX(ptconc(idxppb),BNDLPT)
         enddo
c
c  --- get the sum of the tracer species for ISP ---
c
         sumisp = 0.
         do i=iptcls(idxipt(ITRISP)),nptcls(idxipt(ITRISP))
            idxisp = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                DBLE(numlay) * DBLE(i-1)
            sumisp = sumisp + ptconc(idxisp)
         enddo
c
c  --- adjust PPB for in-cloud SOA formation based on distribution of ISP ---
c
         delppb = delcon(7,ksopb)
         ippb = iptcls(idxipt(ITRPPB)) - 1
         do i=iptcls(idxipt(ITRISP)),nptcls(idxipt(ITRISP))
            idxisp = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                    DBLE(numlay) * DBLE(i-1)
            ippb = ippb + 1
            idxppb = idxcel + DBLE(numcol) * DBLE(numrow) * 
     &                                 DBLE(numlay) * DBLE(ippb-1)
            ptconc(idxppb) = ptconc(idxppb) + delppb * ptconc(idxisp) / sumisp
            ptconc(idxppb) = MAX(ptconc(idxppb),BNDLPT)
         enddo
      endif
c
c-----------------------------------------------------------------------
c    Return point:
c-----------------------------------------------------------------------
c
      return
      end
