      subroutine piginit(numprocs,dt,npts,nrads,nspc,pttrace,height,windu,
     &                                         windv,tempk,press)
      use filunit
      use grid
      use chmstry
      use pigsty
      use ptemiss
      use rtracchm
      use tracer
      use node_mod
      implicit none
c
c----CAMx v7Beta6 190902
c
c     PIGINIT initializes PiG puffs from flagged sources
c
c     Copyright 1996 - 2018
c     Ramboll
c          
c     Modifications: 
c        07/05/02    Added code for IRON-PiG
c        06/24/03    Removed thetapig
c        09/04/03    Removed adiabatic factor; improved plume rise
c                    entrainment; revamped puff initiation for "chained" puffs
c        11/12/03    Added RTRAC species injection
c        08/25/05    Added new PiG-OSAT pointers to allow different input
c                    point source lists across model restarts
c        02/02/06    Removed GREASD-PiG specific conditional code
c        12/20/06    Added plume height overrided
c        06/04/09    Added sigmax and multiple puff releases according to
c                    PIGMXDT
c        04/08/10    Minor mods for updated plume rise algorithm
c        11/06/12    Added check to make sure PiGged source does not
c                    trigger plume distribution override
c        04/16/13    Added update to Lslice when puff array space
c                    is recycled.
c
c     Input arguments:
c        numprocs            number of MPI slices
c        dt                  time step size (s)
c        npts                number of point sources in the array
c        nrads               number of radical species
c        nspc                number of model species
c        pttrace             point source emissions (mol/s)
c        height              layer height (m)
c        windu               wind speed in x-direction (m/s)
c        windv               wind speed in y-direction (m/s)
c        tempk               air temprature (deg.K)
c        press               air pressure (mb)
c
c     Output arguments:
c        none
c
c     Subroutines Called:
c        PLUMERIS
c
c     Called By:
c        EMISTRNS
c
      include  "camx.prm"
      include  "flags.inc"
c
      integer numprocs
      real    dt
      integer npts
      integer nrads
      integer nspc
      real    pttrace(npts,nspc)
      real    height(*)
      real    windu(*)
      real    windv(*)
      real    tempk(*)
      real    press(*)
c
      real hght1d(0:MXLAYER),wind1d(MXLAYER),tempk1d(MXLAYER),
     &     dtdz1d(MXLAYER),press1d(MXLAYER)
      real frctr(MXRECTR)
      real gamma,p0,t1,t2,frsum,xtmp,ytmp,w2,dz,
     &     dtheta,dstktmp,zstk,dtpuff,dtime,grav,wp,tp,zrise,
     &     trise,sigma,pwidth,rip,fp,qp2,cq1,cq2,rkp,xlmax
      integer nr,l,m1,lsrc,n,i,j,k,igrd0,ip,ic,igrd,ig,
     &        ii,jj,kk,n3d,n3d_full,npuff,ipuff,m,m0,is,iempty
      integer iprc, igr
      real xerf
c
      data gamma /0.286/
      data p0 /1000./
      data grav /9.8/
      data cq1,cq2 /0.4,3.0/
c
c-----Entry point
c
c-----Divide initial mass among PiG reactors: (1) partition according to 
c     ERF function; (2) correct for truncation of xerf to conserve mass
c
      frsum = 0.
      do nr = 1,nreactr
        t1 = nr*sqrt(2.)/nreactr
        t2 = (nr-1)*sqrt(2.)/nreactr
        frctr(nr) = xerf(t1) - xerf(t2)
        frsum = frsum + frctr(nr)
      enddo
      do nr = 1,nreactr
        frctr(nr) = frctr(nr)/frsum
      enddo
c
c-----Loop over all elevated point sources in domain (master grid)
c
      m1 = 1
      do 50 lsrc = 1,nosrc(1)
        n = idsrc(lsrc,1)
        if (.not.lpiglet(n)) goto 50
        iempty = 0
        do l = 1,nspc
          if( lemmap(l) .GT. 0 .AND. pttrace(n,l) .GT. 0.) iempty = 1
        enddo
c
c======================== Source Apportion Begin =======================
c
        if( ltrace .AND. (tectyp .EQ. RTRAC .OR.
     &                                      tectyp .EQ. RTCMC) ) then
           do is = 1,nrtrac
             if (sapnts(n,is).gt.0.) iempty = 1
           enddo
        endif
c
c========================= Source Apportion End ========================
c
        if (iempty.eq.0) goto 50
c
c-----Found a PiG point source with non-zero emissions; load local 1-D
c     vectors from finest grid
c
        i = isrc(lsrc,1)
        j = jsrc(lsrc,1)
        igrd0 = 1
        do ip = 1,ngrid
          do ic = 1,nchdrn(ip)
            igrd = idchdrn(ic,ip)
            ig = mapgrd(igrd)
            if (i.ge.inst1(ig) .and. i.le.inst2(ig) .and.             
     &          j.ge.jnst1(ig) .and. j.le.jnst2(ig)) 
     &      igrd0 = igrd
          enddo
        enddo

        if (igrd0.eq.1) then
          ii = i
          jj = j
        else
          xtmp = xstk(n,1) - (inst1(igrd0) - 1)*delx
          ytmp = ystk(n,1) - (jnst1(igrd0) - 1)*dely
          ii = 2 + FLOOR(xtmp/delx*meshold(igrd0))
          jj = 2 + FLOOR(ytmp/dely*meshold(igrd0))
        endif
c
c--- adjust cell index and skip if PiG not in this slice ---
c
        ii = ii-mi0(igrd0)
        jj = jj-mj0(igrd0)
        if( ii .LT. 0 .OR. ii .GT. mmxp(igrd0) ) goto 50
        if( jj .LT. 0 .OR. jj .GT. mmyp(igrd0) ) goto 50
c
        hght1d(0) = 0.
        do k = 1,nlay(igrd0)
          n3d = ii + mmxp(igrd0)*(jj - 1) +
     &                         mmxp(igrd0)*mmyp(igrd0)*(k - 1)
          n3d_full = ii + ncol(igrd0)*(jj - 1) +
     &                         ncol(igrd0)*nrow(igrd0)*(k - 1)
          hght1d(k) = height(iptr3d(igrd0)-1+n3d_full)
          tempk1d(k) = tempk(iptr3d(igrd0)-1+n3d)
          press1d(k) = press(iptr3d(igrd0)-1+n3d)
          w2 = windu(iptr3d(igrd0)-1+n3d)*windu(iptr3d(igrd0)-1+n3d) +
     &         windv(iptr3d(igrd0)-1+n3d)*windv(iptr3d(igrd0)-1+n3d)
          wind1d(k) = amax1(sqrt(w2),0.1) 
        enddo
        do k = 1,nlay(igrd0)
          if (k.lt.nlay(igrd0)) then 
            dz = (hght1d(k+1) - hght1d(k-1))/2.  
            dtheta = (tempk1d(k+1)*(p0/press1d(k+1))**gamma -  
     &                tempk1d(k)*(p0/press1d(k))**gamma)  
            dtdz1d(k) = dtheta/dz  
          else
            dtdz1d(k) = dtdz1d(k-1)  
          endif  
        enddo
c
c-----Disable the plume distribution override for PiG
c
        if( flowrat(n) .LT. 0. ) then
            write(iout,'(//,a)') 'ERROR in PIGINIT:'
            write(iout,*) 'Found negative flow rate for a PiGged source.'
            write(iout,*) 'The plume distribution override is not ',
     &                    'allowed for PiGged sources.'
            write(iout,*) 'Either reset the flow rate ',
     &                                     'on all PiGged sources'
            write(iout,*) 'Or turn off PiG treatment.'
            call camxerr()
        endif
c
c-----Calculate plume rise
c
        dstktmp = abs(dstk(n))
        if (effph(n) .lt. 0.) then
          zstk = abs(effph(n))
        else
          call plumeris(nlay(igrd0),hght1d,tempk1d,dtdz1d,wind1d,
     &                  hstk(n),dstktmp,tstk(n),vstk(n),zstk) 
        endif
        do kk = 1,nlay(igrd0)
          if (hght1d(kk).gt.zstk) goto 15
        enddo
        kk = nlay(igrd0)
        if (zstk.gt.hght1d(kk)) zstk = hght1d(kk)
  15    continue
c
c-----Divide the emissions into several puffs if the length exceeds half
c     the finest grid cell size; divide into several puffs if the current
c     timestep exceeds a specific max timestep
c
        xlmax = 1.e6
        do ip = 1,ngrid
          xlmax = amin1(xlmax,deltay(ip)/2.)
        enddo
        npuff = int(wind1d(kk)*dt/xlmax - 0.1) + 1
        dtpuff = amin1(dt,pigmxdt)
        npuff = max(npuff,int(dt/dtpuff-0.01)+1)
        dtpuff = dt/float(npuff)
        dtime = dt
        do 40 ipuff = 1,npuff
          dtime = dtime - dtpuff
c
c-----Calculate initial cross-puff sigma accounting for turbulent 
c     entrainment during release
c
          wp = amax1(1.,vstk(n)/2.)
          tp = (tempk1d(kk) + tstk(n))/2.
          zrise = zstk - hstk(n)
          trise = zrise/wp
          pwidth = sqrt(2.)*dstktmp
          sigma =  pwidth
          rip = grav*pwidth*abs(tp - tempk1d(kk))/(tempk1d(kk)*wp*wp)
          fp = 1. + 4.*rip
          qp2 = fp*wp*wp*(cq1 + cq2*wind1d(kk)*wind1d(kk)/
     &          (wind1d(kk)*wind1d(kk) + wp*wp))
          rkp = 0.15*pwidth*sqrt(qp2)
          sigma = sqrt(sigma*sigma + 2.*rkp*trise)
c
c-----Looking for a spot for the pig; add to the first empty location
c
          do m = m1,npig
            if (ingrd(m).eq.0) then
              m0 = m
              m1 = m0 + 1
              goto 20
            endif
          enddo
          npig = npig + 1
          m0 = npig
          m1 = npig + 1
  20      continue
          if (m0.gt.MXPIG) then
            write(iout,'(//,a)') 'ERROR in PIGINIT:'
            write(iout,*) 'PiG number exceeds maximum of ',MXPIG
            write(iout,*) 'Increase MXPIG and recompile'
            call camxerr()
          endif
c
c-----Load puff variables
c
          lnewt(m0) = .true.
          lnewg(m0) = .true.
          idpig(m0) = n
          ingrd(m0) = igrd0
          xpigf(m0) = xstk(n,1)
          xpigb(m0) = xstk(n,1)
          ypigf(m0) = ystk(n,1)
          ypigb(m0) = ystk(n,1)
          zpig(m0)  = zstk
          sigz(m0)  = amin1(sigma,amax1(zrise,1.))
          sigx(m0)  = sigma
          sigy(m0)  = sigma
          agepigf(m0) = dtime + dtpuff
          agepigb(m0) = dtime
          htfms(m0) = 0.
          htfmb(m0) = 0.
          vtfms(m0) = 0.
          vtfmb(m0) = 0.
          if( lmpi ) then
              do iprc=1,numprocs
                 do igr=1,ngrid
                    Lslice(igr,iprc,m0) = 0
                  enddo
              enddo
          endif
c
c======================== Source Apportion Begin =======================
c
          if( ltrace .AND. (tectyp .NE. RTRAC .OR.
     &                                      tectyp .EQ. RTCMC) ) then
             ipufmap(m0) = ipigmap(n)
             ipufgrp(m0) = ipiggrp(n)
          endif
c
c========================= Source Apportion End ========================
c
c
c-----Calculate 3-D puff dimensions and fill with concentrations
c
          axisz(m0) = 3.*sigz(m0)
          axisy(m0) = 3.*sigy(m0)
          fmspig(m0) = 1.
          pufftop(m0) = zstk + 1.5*sigz(m0)
          puffbot(m0) = zstk - 1.5*sigz(m0)
          if (puffbot(m0).lt.0.) then
            axisz(m0) = axisz(m0) + puffbot(m0)
            puffbot(m0) = 0.
          endif
          if (pufftop(m0).gt.hght1d(nlay(igrd0))) then
            axisz(m0) = axisz(m0) - 
     &                  (pufftop(m0) - hght1d(nlay(igrd0)))
            pufftop(m0) = hght1d(nlay(igrd0))
          endif
c
c-----Initialize all PiG concentrations to zero, then fill emissions
c     and initialize radicals
c
          do nr = 1,nreactr
            do l = nrads+1,nspc
              puffmass(l,nr,m0) = 0.
              is = lemmap(l)
              if( is .GT. 0 )
     &           puffmass(l,nr,m0) = pttrace(n,l)*
     &                                      dtpuff*1.e6*frctr(nr)
            enddo
            do l = 1,nrads
              puffmass(l,nr,m0) = 1.e-9 ! puffmass(radical) in ppm
            enddo
c
c======================== Source Apportion Begin =======================
c
            if( ltrace .AND. (tectyp .EQ. RTRAC .OR.
     &                        tectyp .EQ. RTCMC) ) then
               do is = 1,nrtrac
                 puffrt(is,nr,m0) = sapnts(n,is)*
     &                              dtpuff*1.e6*frctr(nr)
               enddo
            endif
c
c========================= Source Apportion End ========================
c
          enddo

  40    continue
  50  continue
c
      return
      end
