      subroutine pigsampl(m1,m2,m3,ioff,joff,lrtrac,igrd,lgrd,
     &                    nsp,nav,nxs,nys,nrow,meshnst,inst1,jnst1,dt,
     &                    delx,dely,deltax,deltay,tempk,press,conc,cncsmp)
      use camxcom
      use chmstry
      use pigsty
      use tracer
      use rtracchm
      use node_mod
      implicit none
c
c----CAMx v7Beta6 190902
c
c     PIGSAMPL increments the time-averaged concentrations on PiG/RTRAC
c     sampling grids.
c
c     Copyright 1996 - 2018
c     Ramboll
c
c     Modifications: 
c        8/2/05       Generalized to sample regular model species or
c                     RTRAC species.
c        6/22/07      Temperature/pressure arrays interpolated to sampling
c                     grid for smoother results in complex terrain.
c       06/04/09      Added sigmax
c       10/29/15      Fixed MPI-related bugs, revised vertical sampling 
c                     criterion to be 10 m (rather than layer 1)
c
c     Input arguments:
c        m1                 number of columns in this slice
c        m2                 number of row in this slice
c        m3                 number of layers in this slice
c        ioff               X offset of slice in full domain
c        joff               Y offset of slice in full domain
c        lrtrac             Rtrac flag
c        igrd               nested grid index
c        lgrd               sampling grid index
c        nsp                number of species in input field
c        nav                number of species in output average field
c        nxs                number of sampling grid columns
c        nys                number of sampling grid rows
c        nrow               number of computational grid rows
c        meshnst            meshing factor for computational grid
c        inst1              starting I index of comp grid rel to master grid 
c        jnst1              starting J index of comp grid rel to master grid 
c        dt                 time step for computational grid concentration (s)
c        delx               master x-grid spacing in input units (km or deg)
c        dely               master y-grid spacing in input units (km or deg)
c        deltax             computational grid spacing in x (m)
c        deltay             computational grid spacing in y (m)
c        tempk              computational grid temperature field (K)
c        press              computational grid pressure field (mb)
c        conc               instantaneous gridded concentration (gas=umol/m3,
c                                                                other=ug/m3)
c        cncsmp             average sampling grid concentration (gas=ppm,
c                                                                other=ug/m3)
c
c     Output arguments:
c        cncsmp             average sampling grid concentration (gas=ppm,
c                                                                other=ug/m3)
c
c     Routines Called:
c        PIGCOORD
c        SMPL1PUF
c
c     Called by:
c        CAMx
c
      include "camx.prm"
      include "flags.inc"
c
      integer MXOUT
      parameter (MXOUT = MAX(MXSPEC,MXTRSP))

      integer m1
      integer m2
      integer m3
      integer ioff
      integer joff
      logical lrtrac
      integer igrd
      integer lgrd
      integer nsp
      integer nav
      integer nxs
      integer nys
      integer nrow
      integer meshnst
      integer inst1
      integer jnst1
      real    dt
      real    delx
      real    dely
      real    deltax(nrow)
      real    deltay
      real    tempk(m1,m2,m3)
      real    press(m1,m2,m3)
      real    conc(m1,m2,m3,nsp)
      real    cncsmp(nxs,nys,nav)
c
      real convfac(MXCOLSMP,MXROWSMP,MXOUT)
      real avcnc(MXCOLSMP,MXROWSMP,MXOUT)
      real pmass(MXOUT)
      real bndlow(MXOUT)
c
      integer l,mfac,jc,j,ic,i,n,nr,idum,jj,ng,ll,pgrd
      real dtfact,yc,xc,yf,xf,
     &     dcdx1,dcdx2,c1,c2,dcdy,xpig,ypig,fpigx,fpigy,bpigx,bpigy,
     &     xpmin,xpmax,ypmin,ypmax,xmax,ymax,xlen,sigxx,pdepth
      real deltaxs,deltays,tsmp,psmp
c
      common /compigsampl/ convfac, avcnc
c
c-----Entry point
c
      avcnc = 0.
      dtfact = dt/(dtout*60.)
      mfac = meshsmp(lgrd)/meshnst
      jj = (jsmp1(lgrd) - jnst1)*meshnst + 2
      deltaxs = deltax(jj)/mfac
      jj = jj-joff
      deltays = deltay/mfac
c
c-----Interpolate gridded concentration fields from host computational 
c     grid to sampling grid, and define conc units conversion factor 
c
      jc = (jsmp1(lgrd) - jnst1)*meshnst + 2
      jc = jc-joff
      yc = 0.5
      do 10 j = 1,nys
        ic = (ismp1(lgrd) - inst1)*meshnst + 2
        ic = ic-ioff
        xc = 0.5
        yf = (j - 0.5)/float(mfac)
        if (yc.le.yf) then
          jc = jc + 1
          yc = yc + 1.
        endif
        do 20 i = 1,nxs
          xf = (i - 0.5)/float(mfac)
          if (xc.le.xf) then
            ic = ic + 1
            xc = xc + 1.
          endif
c
c  --- calculate the conversion factor for only cells in 
c      this slice ----
c
          if( ic .GE. ia .AND. ic .LE. iz .AND. jc .GE. ja
     &                                    .AND. jc .LE. jz ) then
             dcdx1 = (tempk(ic,jc-1,1) - tempk(ic-1,jc-1,1))
             dcdx2 = (tempk(ic,jc,1) - tempk(ic-1,jc,1))
             c1 = tempk(ic,jc-1,1) - dcdx1*(xc - xf)
             c2 = tempk(ic,jc,1) - dcdx2*(xc - xf)
             dcdy = (c2 - c1)
             tsmp = c2 - dcdy*(yc - yf)
             dcdx1 = (press(ic,jc-1,1) - press(ic-1,jc-1,1))
             dcdx2 = (press(ic,jc,1) - press(ic-1,jc,1))
             c1 = press(ic,jc-1,1) - dcdx1*(xc - xf)
             c2 = press(ic,jc,1) - dcdx2*(xc - xf)
             dcdy = (c2 - c1)
             psmp = c2 - dcdy*(yc - yf)
             do l = 1,nav
               if (lrtrac) then
                 ng = nrtgas
                 ll = l
               else
                 ng = ngas
                 ll = lavmap(l)
               endif
               if (ll.gt.ng) then 
                 convfac(i,j,l) = 1.
               else
                 convfac(i,j,l) = densfac*(273./tsmp)*(psmp/1013.)
               endif
c
c  --- add background for cells within the computation domain
c      of this slice ----
c
               if( ic .GE. ia .AND. ic .LE. iz .AND. jc .GE. ja 
     &                                         .AND. jc .LE. jz ) then
                 if (lbckgrd) then
                   dcdx1 = (conc(ic,jc-1,1,ll) - conc(ic-1,jc-1,1,ll))
                   dcdx2 = (conc(ic,jc,1,ll) - conc(ic-1,jc,1,ll))
                   c1 = conc(ic,jc-1,1,ll) - dcdx1*(xc - xf)
                   c2 = conc(ic,jc,1,ll) - dcdx2*(xc - xf)
                   dcdy = (c2 - c1)
                   avcnc(i,j,l) = (c2 - dcdy*(yc - yf))/convfac(i,j,l)
                 else
                   avcnc(i,j,l) = 0.
                 endif
               endif
             enddo
          endif
 20     continue
 10   continue
c
c-----Start puff loop
c
      do 50 n = 1,npig
c
c-----Locate the pig in the grid; consider only puffs extending below 10 m
c
        xpig = (xpigf(n) + xpigb(n))/2.
        ypig = (ypigf(n) + ypigb(n))/2.
        call pigcoord(xpig,ypig,i,j,pgrd)
        if (igrd.ne.pgrd) goto 50
        if( puffbot(n) .GT. 10. ) goto 50
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 50
        if( j .LT. ja .OR. j .GT. jz ) goto 50
c
c-----Determine puff location relative to origin of sampling grid
c     and convert to meters
c
        fpigx = xpigf(n) - (ismp1(lgrd) - 1)*delx
        fpigy = ypigf(n) - (jsmp1(lgrd) - 1)*dely
        bpigx = xpigb(n) - (ismp1(lgrd) - 1)*delx
        bpigy = ypigb(n) - (jsmp1(lgrd) - 1)*dely

        fpigx = fpigx/delx*meshnst*deltax(j+joff)
        fpigy = fpigy/dely*meshnst*deltay
        bpigx = bpigx/delx*meshnst*deltax(j+joff)
        bpigy = bpigy/dely*meshnst*deltay
c
c-----Check if puff is not in or near the current sampling grid
c
        xpmax =  3.*sigx(n) + amax1(fpigx,bpigx)
        xpmin = -3.*sigx(n) + amin1(fpigx,bpigx)
        ypmax =  3.*sigy(n) + amax1(fpigy,bpigy)
        ypmin = -3.*sigy(n) + amin1(fpigy,bpigy)
        xmax = nxs*deltaxs
        ymax = nys*deltays
        if (xpmin .GT. xmax .OR. xpmax .LT. 0. .OR.
     &      ypmin .GT. ymax .OR. ypmax .LT. 0.) goto 50
c
c-----We should sample this one; load additional parameters for sampling
c
        xlen  = sqrt((fpigx - bpigx)**2 + (fpigy - bpigy)**2)
        sigxx = xlen/3. + sigx(n)
        pdepth = pufftop(n)
        do l = 1,nav
          pmass(l) = 0.
          do nr = 1,nreactr
            if (lrtrac) then
              pmass(l) = pmass(l) + puffrt(l,nr,n)
              bndlow(l) = rtlbnd(l)
            else
              ll = lavmap(l)
              pmass(l) = pmass(l) + puffmass(ll,nr,n)
              bndlow(l) = bdnl(ll)
            endif
          enddo
        enddo
c
        call smpl1puf(lbckgrd,MXCOLSMP,MXROWSMP,MXOUT,nxs,nys,
     &                nav,ismp1(lgrd),jsmp1(lgrd),mfac,meshnst,
     &                inst1,jnst1,ia,iz,ja,jz,ioff,joff,deltaxs,
     &                deltays,sigxx,sigy(n),fpigx,bpigx,fpigy,bpigy,
     &                pdepth,pmass,bndlow,convfac,avcnc)
  50  continue
c
c-----Increment running average 
c
      do l = 1,nav
        do 60 j = 1,nys
          do i = 1,nxs
            cncsmp(i,j,l) = cncsmp(i,j,l) + dtfact*avcnc(i,j,l)
          enddo
  60    continue
      enddo
c 
      return
      end
