      subroutine pigdrive(m1,m2,m3,ioff,joff,
     &                    igrd,iptr2d,ncol,nrow,nlay,nspcs,nspdep,
     &                    itzon,dt,dx,dy,delx,dely,meshfac,mapscl,
     &                    height,rkv,tempk,tsurf,press,water,windu,
     &                    windv,cldtrns,cwc,pwr,pws,pwg,cph,cellat,
     &                    cellon,topo,sfcz0,albedo,vdep,conc,pigdump,
     &                    pgmserr,fluxes,depfld,ipsa2d,ipsa3d,ipsadep,
     &                    ipa_cel,iproc_id)
      use camxcom
      use filunit
      use chmstry
      use o3colmap
      use pigsty
      use tracer
      use procan
      use rtracchm
      use rtcmcchm
      use node_mod
      implicit none
c
c----CAMx v7Beta6 190902
c
c     PIGDRIVE is the driver program for the PiG submodel.
c     It performs chemistry, puff growth, and mass dumping.
c
c     Copyright 1996 - 2018
c     Ramboll
c
c     Modifications:
c        02/02/06        Renamed PIGDRIVE from IRONDRIV
c        8/31/06         Added map scale factor
c        01/08/07        Added Stage 3 criteria for Mech6 (CB05)
c                        Added ETOH,MTBE,MBUT rxns to Stage 3 criteria for Mech5
c        06/22/07        Removed termination check for puffs extending to
c                        surface layer
c        07/11/07        Added code for RTCMC probing tool
c        02/21/08        Added coupling time for puff aerosol chemistry
c        03/10/08        Added avgcnc and acgrad to calls to IRONCHEM
c        03/31/08        Added average puff age when killed
c        05/14/08        Improved handling of puff PM increments
c        08/04/08        Added IPR code for PiG changes (lost starting at v4.40)
c        11/02/08        Removed z0 calculation to SRFRUF routine
c        04/15/09        Modifications to support OSAT/PSAT deposition output
c        06/17/09        Restructured puff sigmas to improve impacts of
c                        vertical directional shear on puff growth
c        05/28/10        Added mechanism 7 (CB6)
c        07/14/10        Added code for in-line TUV cloud adjustment
c        04/02/12        Removed RADM cloud adjustment option, cloud/aerosol
c                        adjustments now always done with in-line TUV; AHO
c                        file is now just ozone column
c        04/04/12        Introduced 4-dim interpolation of J values over
c                        zenith, height, albedo, terrain
c        04/12/12        Added T and P adjustments to photolysis rates
c        07/20/12        Added Mechanism 1 (CB06 with Iodine chemistry)
c        10/08/12        Added Mechanism 2 (CB6r1)
c        07/31/13        Restructured to improve GREASD PiG chemistry,
c                        wet dep, and dumping
c        10/14/13        Added pressure to KHETERO call
c        06/27/14        Added Mechanism 3 (CB6r2h)
c        08/11/14        Albedo now includes effects of snow
c        11/11/14        Added Mechanism 4 (CB6r3)
c        12/07/14        Modified for VBS
c        02/10/16        Updated for OSAT3
c
c     Input arguments:
c        m1                  number of columns in this slice
c        m2                  number of rows in this slice
c        m3                  number of layers in this slice
c        ioff                X offset of this slice in complete domain
c        ioff                Y offset of this slice in complete domain
c        igrd                grid index 
c        iptr2d              pointers into vectors for 2-D fields
c        ncol                number of columns
c        nrow                number of rows
c        nlay                number of layers
c        nspcs               number of species
c        nspdep              number of species in deposition array
c        itzon               time zone 
c        dt                  time step (s)
c        dx                  cell size in x-direction (m)
c        dy                  cell size in y-direction (m)
c        delx                cell size in x-direction (km or deg)
c        dely                cell size in y-direction (km or deg)
c        meshfac             grid meshing factor relative to master grid
c        mapscl              map scale factor
c        height              gridded layer height (m)
c        rkv                 gridded vertical diffusivity (m2/s)
c        tempk               gridded temperature (K)
c        tsurf               gridded surface temperature (K)
c        press               gridded pressure (mb)
c        water               gridded water vapor concentration (ppm)
c        windu               gridded x-component windspeed (m/s)
c        windv               gridded y-component windspeed (m/s)
c        cldtrns             gridded energy transmission coefficient (fraction)
c        cwc                 gridded cloud water (g/m3)
c        pwr                 gridded precip rain (g/m3)
c        pws                 gridded precip snow (g/m3)
c        pwg                 gridded precip graupel (g/m3)
c        cph                 gridded cloud water pH
c        cellat              gridded cell centroid latitude (deg)
c        cellon              gridded cell centroid longitude (deg)
c        topo                gridded cell terrain altitude (m MSL)
c        sfcz0               gridded cell surface roughness (m)
c        albedo              gridded cell surface albedo
c        vdep                gridded species deposition velocity (m/s)
c        conc                species concentration in the grid (umol/m3,ug/m3)
c        ipa_cel             gridded array to identify if cell is
c                            in a IPRM sub-domain
c        iproc_id            process ID for this slice
c
c     Output arguments:
c        conc                species concentration in the grid (umol/m3,ug/m3)
c        pigdump             PiG dumped mass (umol)
c        pgmserr             PiG dumping mass error (umol)
c        fluxes              boundary mass fluxes (umol, ug)
c        depfld              2-D array of wet deposited mass (mol/ha, g/ha)
c                            and surface liquid concentrations (mol/l, g/l)
c
c     Subroutines Called:
c        PIGCOORD
c        VIRTDUMP
c        KTHERM
c        GETZNTH
c        KPHOTO
c        KHETERO
c        NOXCHEM
c        IRONCHEM
c        CALDATE
c        PIGGROW 
c
c     Called by:
c        PIGEVOL
c
      include "camx.prm"
      include "flags.inc"
      include "deposit.inc"
      include 'chmdbg.inc'
      include 'camx_aero.inc'
c
c
      integer m1
      integer m2
      integer m3
      integer ioff
      integer joff
      integer igrd
      integer iptr2d
      integer ncol
      integer nrow
      integer nlay
      integer nspcs
      integer nspdep
      integer itzon
      real    dt
      real    dx(nrow)
      real    dy
      real    delx
      real    dely
      integer meshfac
      real    mapscl(m1,m2)
      real    height(m1,m2,m3)
      real    rkv(m1,m2,m3)
      real    tempk(m1,m2,m3)
      real    tsurf(m1,m2)
      real    press(m1,m2,m3)
      real    water(m1,m2,m3)
      real    windu(m1,m2,m3)
      real    windv(m1,m2,m3)
      real    cldtrns(m1,m2,m3)
      real    cwc(m1,m2,m3)
      real    pwr(m1,m2,m3)
      real    pws(m1,m2,m3)
      real    pwg(m1,m2,m3)
      real    cph(m1,m2,m3)
      real    cellat(m1,m2)
      real    cellon(m1,m2)
      real    topo(m1,m2)
      real    sfcz0(m1,m2)
      real    albedo(m1,m2)
      real    vdep(m1,m2,nspcs)
      real    conc(m1,m2,m3,nspcs)
      real*8  pigdump(*)
      real*8  pgmserr(*)
      real*8  fluxes(nspcs,11)
      real    depfld(m1,m2,nspdep)
      integer ipsa2d
      integer*8 ipsa3d
      integer ipsadep
      integer iproc_id
c
c======================== Process Analysis Begin ===============================
c
      integer ipa_cel(ncol,nrow,nlay)
c
c========================= Process Analysis End ================================
c
      logical ldarkp,ldump,lvdmp(MXPIG),lkill(MXPIG),lstage3
      logical lchemp
      real*8  flxdry(MXSPEC), flxwet(MXSPEC), puff_flux(MXPIG,MXSPEC)
      real    puff_dep(MXPIG,2*MXSPEC)
      real    conpig(MXSPEC+1,MXRECTR), conback(MXSPEC+1)
      real    connew(MXSPEC+1)
      real    dumpmass(MXSPEC), pconc(MXLAYER)
      real    volvec(MXLAYER), wtfac(MXLAYER), pcvfac(MXLAYER)
      real    pcorrec(MXLAYER,MXSPEC), othrpuf(MXLAYER,MXSPEC)
      real    pbdnl(MXLAYER)
      real    depth(MXLAYER), rkpig(MXLAYER), vd(MXSPEC)
      real    dpfl(2*MXSPEC)
      real    delconc(MXLAYER,MXSPEC+1),avgcnc(MXSPEC+1)
      real    sendum(1,1)
c
      integer n,nj,nz,i,j,k,kk,kkk,l,ll,is,ierr,
     &        kpb,kpt,ij,iozon,nr,idum,
     &        idxcel,noxcnt,
     &        hinox(MXRECTR),ierr1,ierr2,aero_flag
      integer ipa_idx
      integer kohso2,kohco,kohno2
      real areamx,pufscl,delt,dtpig,tpuff,wpuff,ppuff,sumwt,
     &     deplyr,alb,trn,
     &     hght,ctrns,zenith,convfac,volpuff,uu,vv,wind,press0,
     &     delms,rctdmp,volfac,tmass,tvol,resid,negconc,
     &     z0,bmass,hshear,vshear,dudx,dudy,dudz,dvdx,
     &     dvdy,dvdz,avgup,avgum,avgvp,avgvm,dz,zzp,zzm,cwpuff,
     &     phpuff,dtchem,totgas
      real xpig,ypig,xdist,ydist,xlen,axisx,pi,test,atm,o2,ch4,h2
      real pufang,shrang
      real pufno(MXRECTR),pufno2(MXRECTR)
c
      common /pigdrivdat/ puff_flux, puff_dep
c
c======================== Source Apportion Begin =======================
c
      integer*8 idx
      real    rtcon(MXTRSP),rtdump(MXTRSP)
      real    dpflrt(MXTRSP),vdrt(MXTRSP)
      real*8  flxdryrt(MXTRSP)
c
c========================= Source Apportion End ========================
c
      pi = 3.1415927
      atm = 1.e6
      O2  = 2.095e5
      CH4 = 1.75
      H2  = 0.50
c
c-----Entry point
c
      do i=1,MXSPEC+1
        avgcnc(i) = 0.
      enddo
      do i=1,MXPIG
        do j=1,MXSPEC
          puff_flux(i,j) = 0.
          puff_dep(i,j) = 0.
          puff_dep(i,MXSPEC+j) = 0.
        enddo
      enddo 
c
c-----Is PiG to do chemistry?
c
      lchemp = .false.
      if (lchem .and. idmech.ne.10) lchemp = .true.
c
c-----Check for negative incoming CONC values and initialize puff overlap
c     array VPCONC
c
      do is = 1,nspec
        do k = 1,nlay
          do j = 1,m2
            do i = 1,m1
              vpconc(i,j,k,is) = 0.
              if (conc(i,j,k,is) .lt. 0.) then
                write(iout,'(/,a)') 
     &             'Start of PIGDRIVE ...negative value in grid conc'
                write(iout,*) 'species, grid, i, j, k: ',is,igrd,i,j,k
                call flush(iout)
                call camxerr()
              endif
            enddo
          enddo
        enddo
      enddo
      do n = 1,npig
        lkill(n) = .false.
      enddo
c
c-----Perform dry deposition
c       
      if (ldry) then
        do 10 n = 1,npig
        if( lmpi ) Lslice(igrd,iproc_id,n) = 0
c
c-----Skip puffs that are not in the current grid
c
          if (ingrd(n).ne.igrd) goto 10
c
c-----Perform deposition once for new puffs over their age
c
          if (lnewg(n)) then
            delt = (agepigf(n) + agepigb(n))/2.
c
c-----Perform deposition for old puffs over current grid's timestep
c
          else
            delt = dt
          endif
c
          xpig = (xpigf(n) + xpigb(n))/2.
          ypig = (ypigf(n) + ypigb(n))/2.
          call pigcoord(xpig,ypig,i,j,idum)
c
c--- adjust cell index and skip if PiG not in this slice ---
c
          i = i-ioff
          j = j-joff
          if( i .LT. ia .OR. i .GT. iz ) goto 10
          if( j .LT. ja .OR. j .GT. jz ) goto 10
          if( lmpi ) Lslice(igrd,iproc_id,n) = 1
c
          if(puffbot(n) .GE.
     &            amax1(10.,height(i,j,1)/2.)) goto 10
          dz = height(i,j,1) - puffbot(n)
          do l = 1,nspec
            vd(l) = vdep(i,j,l)
            flxdry(l) = 0.
            dpfl(l) = 0.
          enddo
c
          call irondry(nrad,nspec,nreactr,dx(j+joff),dy,dz,mapscl(i,j),
     &                 delt,puffmass(1,1,n),axisz(n),vd,dpfl,flxdry)
c
          do 11 l = 1,nspec
            fluxes(l,11) = fluxes(l,11) + flxdry(l)
            do ll = 1,ndepspc
              if (l .eq. ldepmap(ll)) then
                depfld(i,j,ll) = depfld(i,j,ll) + dpfl(l)
                goto 11
              endif
            enddo
 11       continue
c
c======================== Source Apportion Begin =======================
c
          if (ltrace ) then
             if ( tectyp .EQ. RTRAC .OR. tectyp .EQ. RTCMC) then
                do l = 1,nrtrac
                  idx = DBLE(ipsa2d-1) + 
     &                DBLE(i) + DBLE(m1)*DBLE(j-1) + 
     &                            DBLE(m1)*DBLE(m2)*DBLE(l-1)
                  vdrt(l) = vdeprt(idx)
                enddo
                call irondry(0,nrtrac,nreactr,dx(j+joff),dy,dz,
     &                       mapscl(i,j),delt,puffrt(1,1,n),axisz(n),
     &                       vdrt,dpflrt,flxdryrt)
             else
                if( lptdepout ) then
                   call pigdrysa(m1,m2,m3,nrad,nspec,notimespc,i,j,n,
     &                           dpfl,ptconc(ipsa3d),ptdryfld(ipsadep))
                endif
             endif
          endif
c
c========================= Source Apportion End ========================
c
 10     continue
      endif
c
c-----Perform chemistry and wet deposition
c
      if (.not.lchemp .AND. .not.lwet) goto 101
c
c-----Aerosol timing parameters:
c     time_aero : time to call aerosol routine (HHMM)
c     date_aero : date to call aerosol routine (YYJJJ)
c
      aero_flag = 0
      if( lchemp .AND. naero.GT.0 .AND. aeropt.EQ.'CF' ) then
         aero_flag = 1
      endif
c
c-----Treat overlapping puffs via a virtual dump of "large" puffs
c 
      do 20 n = 1,npig
        lvdmp(n) = .false.
        if (ingrd(n).ne.igrd .or. .not.loverlap) goto 20
        xpig = (xpigf(n) + xpigb(n))/2.
        ypig = (ypigf(n) + ypigb(n))/2.
        call pigcoord(xpig,ypig,i,j,idum)
c
c--- adjust cell index and skip if PiG not in this slice ---
c
        i = i-ioff
        j = j-joff
        if( i .LT. ia .OR. i .GT. iz ) goto 20
        if( j .LT. ja .OR. j .GT. jz ) goto 20
        if( lmpi ) Lslice(igrd,iproc_id,n) = 1
c
        xdist = (xpigf(n) - xpigb(n))/delx*dx(j+joff)*meshfac
        ydist = (ypigf(n) - ypigb(n))/dely*dy*meshfac
        xlen  = sqrt(xdist**2 + ydist**2)
        axisx = xlen + 3.*sigx(n)
        pufscl = pi*(axisx/2.)*(1.5*sigy(n))
c
        if (DXYMAX .gt. 0.) then
          areamx = DXYMAX*DXYMAX
        elseif (DXYMAX .lt. 0.) then
          areamx = MIN(DXYMAX*DXYMAX,dx(j+joff)*dy)
        else
          areamx = dx(j+joff)*dy
        endif
        if (pufscl .ge. FLEAK*areamx) then
          lvdmp(n) = .true.
          call virtdump(m1,m2,m3,i0,j0,n,ncol,nrow,nlay,nspec,i,j,
     &                  dx(j+joff),dy,mapscl,height,tempk,press,pcorrec)
          do is = nrad+1,nspec
            do kk = 1,nlay
              vpconc(i,j,kk,is) = vpconc(i,j,kk,is) + pcorrec(kk,is)
            enddo
          enddo
        endif
 20   continue
c
c-----CHEMISTRY + WET DEP loop
c
c$omp parallel default(shared)
c$omp&  private(n,delt,dtpig,xpig,ypig,i,j,idum,k,xdist,ydist,xlen,
c$omp&          axisx,volpuff,kpb,kpt,kk,pcorrec,is,othrpuf,
c$omp&          wtfac,sumwt,deplyr,tpuff,ppuff,wpuff,cwpuff,
c$omp&          phpuff,conback,convfac,
c$omp&          connew,nr,conpig,l,dpfl,flxwet,
c$omp&          noxcnt,hinox,ij,iozon,trn,alb,hght,
c$omp&          ctrns,zenith,ldarkp,lstage3,test,
c$omp&          ierr1,ierr2,rtcon,avgcnc,nj,
c$omp&          dtchem,nz,ierr,totgas,pufno,pufno2,kohso2,kohco,kohno2)
c
c$omp do schedule(dynamic)
c
      do 100 n = 1,npig
c
c-----Skip puffs that are not in the current grid
c
        if (ingrd(n) .ne. igrd) goto 100
c
c-----Perform chemistry once for new puffs over their age
c
        if (lnewg(n)) then
          delt = (agepigf(n) + agepigb(n))/2.
c
c-----Perform chemistry for old puffs over current grid's timestep
c
        else
          delt = dt 
        endif 
        dtpig = min( (agepigf(n) + agepigb(n))/(2.*3600.),
     &                aero_dt(igrd)+(dt/3600.) )
c
c-----Locate the puff in the grid
c
        xpig = (xpigf(n) + xpigb(n))/2.
        ypig = (ypigf(n) + ypigb(n))/2.
        call pigcoord(xpig,ypig,i,j,idum)
c
c--- adjust cell index and skip if PiG not in this slice ---
c
        i = i-ioff
        j = j-joff
        if( i .LT. ia .OR. i .GT. iz ) goto 100
        if( j .LT. ja .OR. j .GT. jz ) goto 100
        if( lmpi ) Lslice(igrd,iproc_id,n) = 1
c
        do k = 1,nlay
          if( height(i,j,k ) .GT. zpig(n)) goto 25
        enddo
        k = nlay
c
  25    xdist = (xpigf(n) - xpigb(n))/delx*dx(j+joff)*meshfac
        ydist = (ypigf(n) - ypigb(n))/dely*dy*meshfac
        xlen  = sqrt(xdist**2 + ydist**2)
        axisx = xlen + 3.*sigx(n)
        volpuff = pi*axisx*axisy(n)*axisz(n)/6.
c
c-----Find layers containing top and bottom of puff
c
        kpb = 1
        kpt = nlay
        do kk = 1,nlay
          if( height(i,j,kk) .GE. pufftop(n) ) then
            kpt = kk
            goto 26
          endif
        enddo
  26    continue
        do kk = nlay,1,-1
          if( height(i,j,kk) .LE. puffbot(n) ) then
            kpb = kk + 1
            goto 27
          endif
        enddo
  27    continue
c
c-----Get overlap vector for this puff if it contributed to the 
c     virtual dump above (subtract out its contribution to VPCONC).
c     Otherwise, do not apply any puff overlap
c
        if (lvdmp(n)) then
          call virtdump(m1,m2,m3,i0,j0,n,ncol,nrow,nlay,nspec,i,j,
     &                  dx(j+joff),dy,mapscl,height,tempk,press,pcorrec)
          do is = nrad+1,nspec
            do kk = kpb,kpt
              othrpuf(kk,is) = vpconc(i,j,kk,is) - pcorrec(kk,is)
            enddo
          enddo
        else
          do is = nrad+1,nspec
            do kk = kpb,kpt
              othrpuf(kk,is) = 0.
            enddo
          enddo
        endif
c 
c-----Average temp, press, water, and background concs (=overlap+gridded)
c     over layers contained in this puff; define layer weighting as the 
c     depth of puff coverage weighted by density
c
        if (kpb .eq. kpt) then
          wtfac(kpb) = 1.
        else
          sumwt = 0.
          deplyr = pufftop(n) - height(i,j,kpt-1)
          wtfac(kpt) = deplyr*press(i,j,kpt)/tempk(i,j,kpt)
          sumwt = sumwt + wtfac(kpt)
c
          deplyr = height(i,j,kpb) - puffbot(n)
          wtfac(kpb) = deplyr*press(i,j,kpb)/tempk(i,j,kpb)
          sumwt = sumwt + wtfac(kpb)
c
          do kk = kpb+1,kpt-1
            deplyr = height(i,j,kk) - height(i,j,kk-1)
            wtfac(kk) = deplyr*press(i,j,kk)/tempk(i,j,kk)
            sumwt = sumwt + wtfac(kk)
          enddo
c
          do kk = kpb,kpt
            wtfac(kk) = wtfac(kk)/sumwt
          enddo
        endif
c
        tpuff = 0.
        ppuff = 0.
        wpuff = 0.
        cwpuff = 0.
        phpuff = 0.
        do is = nrad+1,nspec+1
          conback(is) = 0.
        enddo
        do kk = kpb,kpt
          tpuff = tpuff + wtfac(kk)*tempk(i,j,kk)
          ppuff = ppuff + wtfac(kk)*press(i,j,kk)
          wpuff = wpuff + wtfac(kk)*water(i,j,kk)
          cwpuff = cwpuff + wtfac(kk)*cwc(i,j,kk)
          phpuff = phpuff + wtfac(kk)*cph(i,j,kk)
          do is = nrad+1,nspec
            conback(is) = conback(is) + wtfac(kk)*
     &                    (othrpuf(kk,is) + conc(i,j,kk,is))
          enddo
        enddo
        convfac = densfac*(273./tpuff)*(ppuff/1013.)
        do is = nrad+1,nspec
          if (is.le.ngas) conback(is) = conback(is)/convfac
          conback(is) = amax1(bdnl(is),conback(is))
          connew(is) = conback(is)
        enddo
        connew(nspec+1) = 0.
c
c-----Perform wet deposition
c
        if (lwet .AND. (pwr(i,j,1).ge.cwmin .or.
     &                  pws(i,j,1).ge.cwmin .or.
     &                  pwg(i,j,1).ge.cwmin)) then
c
          do nr = 1,nreactr
            do is = nrad+1,nspec
              conpig(is,nr) = puffmass(is,nr,n)/(volpuff/float(nreactr))
              if (is.le.ngas) conpig(is,nr) = conpig(is,nr)/convfac
            enddo
          enddo
c
          call ironwet(nspec,n,dx(j+joff),dy,mapscl(i,j),cwpuff,
     &                 pwr(i,j,1),pws(i,j,1),pwg(i,j,1),phpuff,tpuff,
     &                 ppuff,volpuff,conpig,conback,axisz(n),delt,dpfl,
     &                 flxwet)
c
          do l = nrad+1,nspec
            puff_flux(n,l) = puff_flux(n,l) + flxwet(l)
            puff_dep(n,l) = puff_dep(n,l) + dpfl(l)
            puff_dep(n,MXSPEC+l) = puff_dep(n,MXSPEC+l)
     &                                             + dpfl(nspec+l)
          enddo
c
c======================== Source Apportion Begin =======================
c
c  --- if doing tracer depostion, call the routine
c      to aadd to tracer arrays ---
c
          if( lptdepout ) then
             call pigwetsa(m1,m2,m3,nrad,nspec,notimespc,i,j,kpb,kpt,n,
     &                           dpfl,ptconc(ipsa3d),ptwetfld(ipsadep))
          endif
c
c======================== Source Apportion Begin =======================
c
        endif
c
c-----Perform chemistry
c
        if (.not.lchemp) goto 100
c
c-----Add background concentrations (CONBACK) to puff increments in
c     each reactor to get total concentration (CONPIG).
c     Background radicals taken from nominal pig center (no averaging)
c     Puff radicals taken from last time step to use as initial guess
c
        do nr = 1,nreactr
          totgas = 0.
          do is = nrad+1,nspec
            conpig(is,nr) = puffmass(is,nr,n)/(volpuff/float(nreactr))
            if (is.le.ngas) then
              conpig(is,nr) = conpig(is,nr)/convfac
              totgas = totgas + abs(conpig(is,nr))
            endif
            conpig(is,nr) = max(bdnl(is),conback(is) + conpig(is,nr))
          enddo
          if (totgas.lt.0.001) goto 100
          pufno(nr)  = conpig(kno,nr)  - conback(kno)
          pufno2(nr) = conpig(kno2,nr) - conback(kno2)
        enddo
        do nr = 1,nreactr
          conpig(nspec+1,nr) = 0.
        enddo
        if (lchemp) then
          do nr = 1,nreactr
            do l = 1,nrad
              conback(l) = conc(i,j,k,l)
              connew(l) = conback(l)
              conpig(l,nr) = puffmass(l,nr,n) ! puffmass(radical) in ppm
            enddo
          enddo
        endif
c
c-----Determine if any puff reactors have high NOx values
c
        noxcnt = 0
        do nr = 1,nreactr
          hinox(nr) = 0
          if (conpig(kno,nr) + conpig(kno2,nr) .gt. 0.2) then
            noxcnt = noxcnt + 1
            hinox(nr) = 1
          endif
        enddo
c
c-----Determine thermal rate constants
c
        call ktherm(tpuff,ppuff)
c
c-----Load local values of ozone, albedo, terrain ht, zenith angle, altitude
c
        call getznth(cellat(i,j),cellon(i,j),timec(igrd),datec(igrd),
     &               itzon,zenith,ldarkp)
c
        alb   = albedo(i,j)
        ij = i + (j - 1)*m1
        iozon = icdozn(iptr2d - 1 + ij)
c
        trn = topo(i,j)/1000.
        hght = height(i,j,k)/2000.
        if (k.gt.1) hght = (height(i,j,k) + height(i,j,k-1))/2000.
        ctrns = 0.
        if (.not.ldarkp) ctrns = cldtrns(i,j,k)
c
c-----Determine photolysis rates through interpolation of look-up table
c
        call kphoto(iozon,alb,trn,hght,zenith,ctrns,ldarkp,tpuff,ppuff)
c
c-----Determine if any GREASD puff reactors are in Stage 3 chemistry
c
        if (ipigflg.EQ.GRESPIG) then
          lstage3 = .false.
          do nr = 1,nreactr
            if (idmech.eq.5) then
              kohso2 = 44
              kohco  = 29
              kohno2 = 25
            elseif (idmech.eq.6) then
              kohso2 = 63
              kohco  = 66
              kohno2 = 28
            elseif (idmech.eq.2 .OR. idmech.eq.3) then
              kohso2 = 52
              kohco  = 123
              kohno2 = 45
            elseif (idmech.eq.4) then
              kohso2 = 52
              kohco  = 120
              kohno2 = 45
            else
              write(iout,'(/,a)')
     &         'In PIGDRIVE, selected mechanism not implemented for PiG'
              write(iout,'(a,i3)') 'Mechanism: ',idmech
              call flush(iout)
              call camxerr()
            endif
            test = 0.
            if( pufno2(nr) .NE. 0. ) 
     &            test = (rk(kohso2)*conpig(kSO2,nr) +
     &              rk(kohco)*conpig(kCO,nr))  /
     &             (rk(kohno2)*pufno2(nr))
            if (.not.ldarkp .AND. test.gt.1. .AND.
     &          pufno2(nr).gt.pufno(nr)) lstage3 = .true.
            if (pufno2(nr).lt.conback(kNO2)) lstage3 = .true.
c
            if (conpig(kNO2,nr).le.0. .OR. conpig(kNO,nr).le.0. .OR.
     &                                              lstage3) then
              lkill(n) = .true.
              nkill(5) = nkill(5) + 1
              goto 100
            endif
c
          enddo
        endif
c
c-----For IRONCHEM, perform chemistry for CONBACK if total NOx
c     (puff+background) is less than 200 ppb.
c     All reactors use same background, save results in CONNEW
c
        if (ipigflg .EQ. IRONPIG) then
          igrdchm = -n
          ichm    = i + ioff
          jchm    = j + joff
          kchm    = 0
          ierr1   = 0
          ierr2   = 0
          tchm = tpuff
          wchm = wpuff
          ldchm = ldarkp
c
c-----Calculate the rate constant for heterogeneous hydrolysis of N2O5
c
          call khetero(tpuff,ppuff,wpuff,conback,.FALSE.,1,sendum)
c
          if (noxcnt .ne. nreactr) then
            call ironchem(delt,dtpig,aero_flag,wpuff,tpuff,ppuff,cwpuff,
     &                    phpuff,ldarkp,connew,conback,avgcnc,convfac,
     &                    ierr1,ierr2)
            if (ierr1 .ne. 0) then
              lkill(n) = .true.
              nkill(3) = nkill(3) + 1
c             write(idiag,*)'Dumping puff ',n,': bad background LSODE'
              goto 100
            endif
          endif
        endif
c
c-----Loop over reactors, perform chemistry for CONPIG.
c     Check for NOx concentrations > 200 ppb; if found, use the 
c     reduced chemistry module
c
        do nr = 1,nreactr
          igrdchm = -n
          ichm    = i + ioff
          jchm    = j + joff
          kchm    = -nr
          ierr1   = 0
          ierr2   = 0
          tchm = tpuff
          wchm = wpuff
          ldchm = ldarkp
c
c-----Calculate the rate constant for heterogeneous hydrolysis of N2O5
c
          call khetero(tpuff,ppuff,wpuff,conpig(1,nr),.FALSE.,1,sendum)
c
c-----For GREASD PiG, remove background except for oxidants (O3,H2O2,FORM,CO)
c
          if (ipigflg.EQ.GRESPIG) then
            do is = nrad+1,nspec
              conpig(is,nr) = conpig(is,nr) - conback(is)
            enddo
            conpig(ko3,nr)   = MAX(bdnl(ko3),
     &                             conpig(ko3,nr) + conback(ko3))
            conpig(kh2o2,nr) = MAX(bdnl(kh2o2),
     &                             conpig(kh2o2,nr) + conback(kh2o2))
            conpig(kco,nr)   = MAX(bdnl(kco),
     &                             conpig(kco,nr) + conback(kco))
            if (idmech.eq.5) then
              conpig(khcho,nr)   = MAX(bdnl(khcho),
     &                             conpig(khcho,nr) + conback(khcho))
            else
              conpig(kform,nr)   = MAX(bdnl(kform),
     &                             conpig(kform,nr) + conback(kform))
            endif
          endif
c
          if (hinox(nr) .eq. 1) then
            call noxchem(delt,bdnl(ko3),bdnl(kno),conpig(kno,nr),
     &                   conpig(kno2,nr),conpig(ko3,nr),rk(ipigrxn))
          else
            call ironchem(delt,dtpig,aero_flag,wpuff,tpuff,ppuff,cwpuff,
     &                    phpuff,ldarkp,conpig(1,nr),conback,avgcnc,
     &                    convfac,ierr1,ierr2)
            if (ierr1 .ne. 0) then
              lkill(n) = .true.
              nkill(4) = nkill(4) + 1
c             write(idiag,*)'Dumping puff ',n,': bad puff LSODE'
              goto 100
            endif
            if (ierr2 .ne. 0) then
              lkill(n) = .true.
              nkill(9) = nkill(9) + 1
c             write(idiag,*)'Dumping puff ',n,': bad puff PM'
              goto 100
            endif
          endif
c
c======================== Source Apportion Begin =======================
c
c-----Original RTRAC
c
          if( ltrace .AND. tectyp .EQ. RTRAC .AND.
     &        (nrtherm .GT. 0 .OR. nrtphot .GT. 0) ) then
             do is = 1,nrtgas
               rtcon(is) = puffrt(is,nr,n)/(volpuff/float(nreactr))
               rtcon(is) = rtcon(is)/convfac
             enddo
             if (hinox(nr) .eq. 1) then
               avgcnc(koh)  = conpig(koh,nr)
               avgcnc(kno3) = conpig(kno3,nr)
               avgcnc(ko3)  = conpig(ko3,nr)
             endif
             dtchem = delt/3600.
             call pigchmrt(ppuff,tpuff,avgcnc(koh),avgcnc(ko3),
     &                     avgcnc(kno3),dtchem,rtcon)
             do is = 1,nrtgas
               rtcon(is) = rtcon(is)*convfac
               puffrt(is,nr,n) = rtcon(is)*volpuff/float(nreactr)
             enddo
          endif
c
c-----RTCMC
c
          if( ltrace .AND. tectyp .EQ. RTCMC ) then
             if( ktype .EQ. 1 ) then
                call krtc(tpuff,ppuff)
             else
                call ksci(tpuff,ppuff)
             endif
             if( njschm .GT. 0 ) then
                 do nj = 1,njschm
                    srkrtc(ijschm(nj)) = 0.
                    do nz = nzschm,1,-1
                       if( zenith .LE. zenschm(nz) )
     &                    srkrtc(ijschm(nj)) = rjschm(nz,nj)
                    enddo
                    rkrtc(ijschm(nj))  = DBLE( srkrtc(ijschm(nj)) )
                 enddo
             endif
             do is = 1,nrtgas
               rtcon(is) = puffrt(is,nr,n)/(volpuff/float(nreactr))
               rtcon(is) = rtcon(is)/convfac
             enddo
             if (hinox(nr) .eq. 1) then
               do is = 1,nspec
                 avgcnc(is)  = conpig(is,nr)
               enddo
             endif
             dtchem = delt/3600.
             if( isolv .EQ. 1) then
                call dlsdrv(dtchem,wpuff,atm,O2,CH4,H2,avgcnc,
     &                      rtcon,ierr,.true.)
             elseif( isolv .EQ. 2) then
                call slsdrv(dtchem,wpuff,atm,O2,CH4,H2,avgcnc,
     &                      rtcon,ierr,.true.)
             elseif( isolv .EQ. 3) then
                call rbkdrv(dtchem,wpuff,atm,O2,CH4,H2,avgcnc,
     &                      rtcon,ierr,.true.)
             endif
             do is = 1,nrtgas
               rtcon(is) = rtcon(is)*convfac
               puffrt(is,nr,n) = rtcon(is)*volpuff/float(nreactr)
             enddo
          endif
c
c========================= Source Apportion End ========================
c
c
c-----Subtract off the new background to get new increment.
c     CONPIG and PUFFMASS can be negative again
c
          if (ipigflg.EQ.IRONPIG) then
            do is = nrad+1,nspec
              if (hinox(nr).eq.1) then
                conpig(is,nr) = conpig(is,nr) - conback(is)
              else
                conpig(is,nr) = conpig(is,nr) - connew(is)
              endif
            enddo
          else
            conpig(kO3,nr)   = conpig(kO3,nr)   - connew(kO3)
            conpig(kh2o2,nr) = conpig(kh2o2,nr) - connew(kh2o2)
            conpig(kco,nr)   = conpig(kco,nr)   - connew(kco)
            if (idmech.eq.5) then
              conpig(khcho,nr) = conpig(khcho,nr) - connew(khcho)
            else
              conpig(kform,nr) = conpig(kform,nr) - connew(kform)
            endif
          endif
          do is = nrad+1,nspec
            if (is.le.ngas) conpig(is,nr) = conpig(is,nr)*convfac
            puffmass(is,nr,n) = conpig(is,nr)*volpuff/float(nreactr)
          enddo
          do l = 1,nrad
            puffmass(l,nr,n) = conpig(l,nr) ! puffmass(radical) in ppm
          enddo
        enddo
 100  continue
c
c  --- end of parallelized loop ---
c
c$omp end parallel
c
c  --- put the fluxes and wet deposition into the totals array ---
c
      do 110 n = 1,npig
        if( ingrd(n) .EQ. igrd ) goto 110
c
c-----Locate the puff in the grid
c
        xpig = (xpigf(n) + xpigb(n))/2.
        ypig = (ypigf(n) + ypigb(n))/2.
        call pigcoord(xpig,ypig,i,j,idum)
c
c--- adjust cell index and skip if PiG not in this slice ---
c
        i = i-ioff
        j = j-joff
        if( i .LT. ia .OR. i .GT. iz ) goto 110
        if( j .LT. ja .OR. j .GT. jz ) goto 110
c
c---- accumulate in totals array ---
c
        do 111 l=1,nspec
          fluxes(l,11) = fluxes(l,11) + puff_flux(n,l)
          do ll=1,ndepspc
             if( l .EQ. ldepmap(ll) ) then
                depfld(i,j,ndepspc+ll) = depfld(i,j,ndepspc+ll) + 
     &                                                    puff_dep(n,l)
                depfld(i,j,2*ndepspc+ll) = depfld(i,j,2*ndepspc+ll) +
     &                                              puff_dep(n,MXSPEC+l)
     &                                    
                goto 111
             endif
          enddo
  111   continue
  110 continue
c
c-----End CHEMISTRY loop
c
c-----Start GROWTH and DUMPING loop
c
 101  continue
      do 200 n = 1,npig
c
c-----Skip puffs that are not in the current grid
c
        if (ingrd(n) .ne. igrd) goto 200
c
c-----Perform growth once for new puffs over their age
c
        if (lnewg(n)) then
          delt = (agepigf(n) + agepigb(n))/2.
c
c-----Perform growth for old puffs over current grid's timestep
c
        else
          delt = dt 
        endif 
c
c-----Locate the pig in the grid
c
        xpig = (xpigf(n) + xpigb(n))/2.
        ypig = (ypigf(n) + ypigb(n))/2.
        call pigcoord(xpig,ypig,i,j,idum)
c
c--- adjust cell index and skip if PiG not in this slice ---
c
        i = i-ioff
        j = j-joff
        if( i .LT. ia .OR. i .GT. iz ) goto 200
        if( j .LT. ja .OR. j .GT. jz ) goto 200
        if( lmpi ) Lslice(igrd,iproc_id,n) = 1

        idxcel =  i + m1*(j-1)
        do k = 1,nlay
          if( height(i,j,k) .gt. zpig(n) ) goto 15
        enddo
        k = nlay
  15    xdist = (xpigf(n) - xpigb(n))/delx*dx(j+joff)*meshfac
        ydist = (ypigf(n) - ypigb(n))/dely*dy*meshfac
        xlen  = sqrt(xdist**2 + ydist**2)
        pufang = atan(ydist/(xdist+1.e-10))
c
c-----Find layers containing top and bottom of puff
c
        kpb = 1
        kpt = nlay
        do kk = 1,nlay
          if( height(i,j,kk) .GE. pufftop(n) ) then
            kpt = kk
            goto 151
          endif
        enddo
  151   continue
        do kk = nlay,1,-1
          if( height(i,j,kk) .LE. puffbot(n) ) then
            kpb = kk + 1
            goto 152
          endif
        enddo
  152   continue
c
c-----Get the drivers on puff growth
c
c-----Resolved shear metrics: use deformation for net horizontal shear
c
        dudx = (windu(i,j,k) - windu(i-1,j,k))/dx(j+joff)
        dvdy = (windv(i,j,k) - windv(i,j-1,k))/dy
        avgup = (windu(i,j,k) + windu(i,j+1,k) + windu(i-1,j,k) + 
     &           windu(i-1,j+1,k))/4.
        avgum = (windu(i,j,k) + windu(i,j-1,k) + windu(i-1,j,k) + 
     &           windu(i-1,j-1,k))/4.
        dudy = (avgup - avgum)/dy
        avgvp = (windv(i,j,k) + windv(i+1,j,k) + windv(i,j-1,k) + 
     &           windv(i+1,j-1,k))/4.
        avgvm = (windv(i,j,k) + windv(i-1,j,k) + windv(i,j-1,k) + 
     &           windv(i-1,j-1,k))/4.
        dvdx = (avgvp - avgvm)/dx(j+joff)
        hshear = sqrt((dudy + dvdx)**2 + (dudx - dvdy)**2)
c
        dudz = 0.
        dvdz = 0.
        dz = height(i,j,k)
        if (k.gt.1) dz = height(i,j,k) - height(i,j,k-1)
        if( axisz(n) .gt. dz .and. kpt-kpb .ge. 1 ) then
          zzp = (height(i,j,kpt) + height(i,j,kpt-1))/2.
          zzm = height(i,j,kpb)/2.
          if (kpb.gt.1) zzm = (height(i,j,kpb) + height(i,j,kpb-1))/2.
          avgup = (windu(i,j,kpt) + windu(i-1,j,kpt))/2.
          avgum = (windu(i,j,kpb) + windu(i-1,j,kpb))/2.
          dudz = (avgup - avgum)/(zzp - zzm)
          avgvp = (windv(i,j,kpt) + windv(i,j-1,kpt))/2.
          avgvm = (windv(i,j,kpb) + windv(i,j-1,kpb))/2.
          dvdz = (avgvp - avgvm)/(zzp - zzm)
        endif
        vshear = sqrt(dudz**2 + dvdz**2)
        shrang = atan(dvdz/(dudz+1.e-10))
c
c-----Surface layer met
c
        uu = 0.5*(windu(i,j,1) + windu(i-1,j,1))
        vv = 0.5*(windv(i,j,1) + windv(i,j-1,1))
        wind = sqrt(uu*uu + vv*vv)
        do kk = 1,nlay
          depth(kk) = height(i,j,kk)
          if (kk.gt.1) depth(kk) = height(i,j,kk) - height(i,j,kk-1)
          rkpig(kk) = rkv(i,j,kk)
        enddo
        press0 = press(i,j,1) -  depth(1)*(press(i,j,2) - 
     &                                   press(i,j,1))/height(i,j,2)
        z0 = sfcz0(i,j)
c
c-----Calculate new puff size and fraction of mass that has leaked (DELMS)
c
        if (DXYMAX .gt. 0.) then
          areamx = DXYMAX*DXYMAX
        elseif (DXYMAX .lt. 0.) then
          areamx = MIN(DXYMAX*DXYMAX,dx(j+joff)*dy)
        else
          areamx = dx(j+joff)*dy
        endif
        ldump = .false.
        call piggrow(n,nlay,k,delt,xlen,areamx,
     &               height(i,j,nlay),wind,
     &               tempk(i,j,1),tsurf(i,j),press(i,j,1),press0,rkpig,
     &               depth,z0,hshear,vshear,pufang,shrang,delms,ldump)
c
c-----Skip dumping if it's a new puff
c
        if (lnewg(n) .and. .not.lkill(n)) goto 40
c
c-----Find layers containing top and bottom of puff
c
        kpb = 1       
        kpt = nlay
        do kk = 1,nlay
          if (height(i,j,kk) .GE. pufftop(n)) then
            kpt = kk     
            goto 16
          endif       
        enddo
  16    continue
        do kk = nlay,1,-1
          if (height(i,j,kk) .LE. puffbot(n)) then
            kpb = kk + 1
            goto 17
          endif
        enddo
  17    continue
c
c-----Revise vertical averaging to account for growth; define layer weighting 
c     as the depth of puff coverage weighted by density
c
        if (kpb .eq. kpt) then
          wtfac(kpb) = 1.
        else
          sumwt = 0.
          deplyr = pufftop(n) - height(i,j,kpt-1)
          wtfac(kpt) = deplyr*press(i,j,kpt)/tempk(i,j,kpt)
          sumwt = sumwt + wtfac(kpt)

          deplyr = height(i,j,kpb) - puffbot(n)
          wtfac(kpb)= deplyr*press(i,j,kpb)/tempk(i,j,kpb)
          sumwt = sumwt + wtfac(kpb)
c
          do kk = kpb+1,kpt-1
            deplyr = height(i,j,kk) - height(i,j,kk-1)
            wtfac(kk) = deplyr*press(i,j,kk)/tempk(i,j,kk)
            sumwt = sumwt + wtfac(kk)
          enddo

          do kk = kpb,kpt
            wtfac(kk) = wtfac(kk)/sumwt
          enddo
        endif
c
        do kk = 1,nlay
          deplyr = height(i,j,1)
          if (kk.gt.1) deplyr = height(i,j,kk) - height(i,j,kk-1)
          volvec(kk) = deplyr*dx(j+joff)*dy/mapscl(i,j)**2
          pcvfac(kk) = densfac*(273./tempk(i,j,kk))*
     &                         (press(i,j,kk)/1013.)
        enddo
c
c-----First check for a full dump
c     Slaughter the puff if its horizontal size exceeds AREAMX
c     Slaughter the puff if remaining mass fraction < 0.1
c     Slaughter the puff according to chemistry flags
c
        axisx = xlen + 3.*sigx(n)
        pufscl = pi*(axisx/2.)*(1.5*sigy(n))
c
        if (lkill(n) .OR. pufscl.GT.areamx .OR. fmspig(n).lt.0.1) then
          if (.not.lkill(n)) then
            if (pufscl.gt.areamx) then
              nkill(1) = nkill(1) + 1
            elseif (fmspig(n).lt.0.1) then
              nkill(2) = nkill(2) + 1
            endif
          endif
          ingrd(n) = 0
          ldump = .true.
c
          do is = nrad+1,nspec
            dumpmass(is) = 0.
            do nr = 1,nreactr
              dumpmass(is) = dumpmass(is) + puffmass(is,nr,n)
              puffmass(is,nr,n) = 0.
            enddo
          enddo
          pigage(igrd) = pigage(igrd) + agepigf(n)
          nage(igrd) = nage(igrd) + 1
c
c======================== Source Apportion Begin =======================
c
          if( ltrace .AND. (tectyp .EQ. RTRAC .OR.
     &                      tectyp .EQ. RTCMC) ) then
             do is = 1,nrtrac
               rtdump(is) = 0.
               do nr = 1,nreactr
                 rtdump(is) = rtdump(is) + puffrt(is,nr,n)
                 puffrt(is,nr,n) = 0.
               enddo
             enddo
          endif
c
c========================= Source Apportion End ========================
c
          goto 60
        endif
c
c-----If no full dump, check for leakage
c     (LDUMP and DELMS were set in PIGGROW)
c
        if (ldump) then
          do is = nrad+1,nspec
            dumpmass(is) = 0.
            do nr = 1,nreactr
              rctdmp = puffmass(is,nr,n)*delms
              dumpmass(is) = dumpmass(is) + rctdmp
              puffmass(is,nr,n) = puffmass(is,nr,n) - rctdmp
            enddo
          enddo
c
c======================== Source Apportion Begin =======================
c
          if( ltrace .AND. (tectyp .EQ. RTRAC .OR.
     &                      tectyp .EQ. RTCMC) ) then
             do is = 1,nrtrac
               rtdump(is) = 0.
               do nr = 1,nreactr
                 rctdmp = puffrt(is,nr,n)*delms
                 rtdump(is) = rtdump(is) + rctdmp
                 puffrt(is,nr,n) = puffrt(is,nr,n) - rctdmp
               enddo
             enddo
          endif
c
c========================= Source Apportion End ========================
c
        endif
c
c-----Update total dumped mass 
c
 60     if (ldump) then
          do is = nrad+1,nspec
            pigdump(is) = pigdump(is) + DBLE(dumpmass(is))
          enddo
          delconc = 0.0 ! must be initialized
c
c-----Dump DUMPMASS into grid; allocate mass in proportion to layer densities
c     among layers contained within this puff
c
          do 300 is = nrad+1,nspec
            if (dumpmass(is).eq.0.) goto 300
            do kk = kpb,kpt
              volfac = wtfac(kk)/volvec(kk)
              pconc(kk) = conc(i,j,kk,is) + dumpmass(is)*volfac
              if (is.le.ngas) then
                pbdnl(kk) = bdnl(is)*pcvfac(kk)
              else
                pbdnl(kk) = bdnl(is)
              endif
c
c-----Woops! DUMPMASS leads to negative grid concentration
c
              if (pconc(kk) .lt. pbdnl(kk)) then
                tmass = 0.
                bmass = 0.
                tvol = 0.
                do kkk = kpb,kpt
                  tmass = tmass + conc(i,j,kkk,is)*volvec(kkk)
                  tvol  = tvol + volvec(kkk)
                  if (is.le.ngas) then
                    pbdnl(kkk) = bdnl(is)*pcvfac(kkk)
                  else
                    pbdnl(kkk) = bdnl(is)
                  endif
                  bmass = bmass + pbdnl(kkk)*volvec(kkk)
                enddo
                resid = tmass + dumpmass(is)
c
c-----Try spreading any positive residual mass evenly over layers contained 
c     within puff depth (if there is enough grid mass)
c
                if (resid .ge. bmass) then
                  do kkk = kpb,kpt
                    pconc(kkk) = resid/tvol
                  enddo
                  goto 301
c
c-----Not enough grid mass for a postive residual. Set grid conc to lower
c     bound and add negative mass increment to tracking array
c
                else
                  negconc = pconc(kk) - pbdnl(kk)
                  pgmserr(is) = pgmserr(is) + DBLE(negconc*volvec(kk))
                  pconc(kk) = pbdnl(kk)
                endif
              endif
            enddo
 301        do kk = kpb,kpt
              delconc(kk,is) = pconc(kk) - conc(i,j,kk,is)
              conc(i,j,kk,is) = pconc(kk)
c
c======================== Process Analysis Begin ====================================
c
              if( lipr ) then
                if( i .GE. ia .AND. i .LE. iz .AND. j .GE. ja
     &                                         .AND. j .LE. jz ) then
                     if( ipa_cel(i+i0,j+j0,kk) .GT. 0 ) then
                       ipa_idx = ipa_cel(i+i0,j+j0,kk)
                       cipr(IPR_PIGEMIS, ipa_idx, is) =
     &                   cipr(IPR_PIGEMIS, ipa_idx, is) + delconc(kk,is)
                     endif
                endif
              endif
c
c========================= Process Analysis End ================================
c
            enddo
 300      continue
c
c======================== Source Apportion Begin =======================
c
          if( ltrace ) then
c
c-----Update RTRAC tracer concentration in the grid ---
c
            if( tectyp .EQ. RTRAC .OR. tectyp .EQ. RTCMC ) then
              do is = 1,nrtrac
                do kk = kpb,kpt
                  volfac = wtfac(kk)/volvec(kk)
                  idx = DBLE(ipsa3d-1) + DBLE(i) + 
     &               DBLE(m1)*DBLE(j-1) + DBLE(m1)*DBLE(m2)*DBLE(kk-1) + 
     &                            DBLE(m1)*DBLE(m2)*DBLE(nlay)*DBLE(is-1)
                  ptconc(idx) = ptconc(idx) + rtdump(is)*volfac
                enddo
              enddo
c
c-----Call routine to update SA tracer concentration in the grid ---
c
            else
              call pigdumpsa(m1,m2,m3,nrad,nspec,ntotsp,
     &                       i,j,kpb,kpt,n,delconc,ptconc(ipsa3d) )
            endif
          endif
c
c========================= Source Apportion End ========================
c
        endif
c
c-----Reset ldump and set LNEWG to false
c
 40     continue
        ldump = .false.
        lnewg(n) = .false.
c
c-----End GROWTH and DUMPING loop
c
 200  continue
c
      return
      end
