      PROGRAM PTSRKCELLCUT

      parameter (bart=1400000,big=1400000)

      CHARACTER*4 IFILE(10), NOTE(60), MSPEC(10,50), MNAME(10)
      CHARACTER*4  SPNAME(10,1500)
      DIMENSION EMISS(big)

      character*150 INFILE
      character*400 fname
c      character*100 fname
      character*10 cname
      character*20 fac(100)

      REAL*4     DUMX(big),DUMY(big),IDUMX(big),IDUMY(big)
      REAL       DUMXX(big),DUMYY(big),emscut(big)
      REAL       FLOW(big),PLUMHT(big)     
      REAL       ncut
      INTEGER    KCELL(big),ICELL(big),JCELL(big),stack(100)

      character*20 cstack(bart),cfacil(bart),ajunk3,ajunk4
      character*20 ajunk,ajunk2,scc(bart)
      character*4 zname(10),bmatch(10),bmiss(10),sic(bart)
      integer st(bart),cn(bart),ic(bart),jc(bart),state,county
      integer sta(100),cou(100),reg(100),zz,ii
      integer uamid,gmap(150,150),mask(big),nxc,nyc,iii
      real clat(bart),clon(bart)
      real xorg,yorg,lx
      real rlon0,rlat0,tlat1,tlat2,rlon,rlat
      data zname/'A','V','E','R','A','G','E',' ',' ',' '/
      data bmatch/'B','A','R','T',' ',' ',' ',' ',' ',' '/
      data bmiss/'M','I','S','S',' ',' ',' ',' ',' ',' '/
      data ione /1/
      data izero /0/
      data zero /0./

c      read(*,'(10x,*)') xorg,yorg,lx,nxc,nyc,rlon0,rlat0,tlat1,tlat2
c      read(*,'(10x,a)') fname
c      read(fname,*) xorg,yorg,lx,nxc,nyc,rlon0,rlat0,tlat1,tlat2
c      write(*,*)'QA Grid output specs:',xorg,yorg,lx,nxc,nyc
c      write(*,*)'Projection specs:',rlon0,rlat0,tlat1,tlat2


      read(*,'(10x,a)') infile
      open(7,file=infile,form='unformatted',status='old',err=999)
      write(*,*) 'Opened emissions file: ',infile

      read(*,'(10x,a)') infile
      open(9,file=infile,form='unformatted')
      write(*,*) 'Opened emissions adjusted file: ',infile

      read(*,'(10x,a)') infile
      open(20,file=infile,form='formatted',status='new')
      write(*,*) 'Opened ASCII Q/A file: ',infile

      read(*,'(10x,a)') fname
      read(fname,*) ncut
      write(*,*)'Emissions adjustment factor:',ncut
      
      read(*,'(10x,a)') fname
      read(fname,*) nstack
      write(*,*)'Number of source groups:',nstack

c      do i = 1, nstack
  
      read(*,'(10x,a)') fname
      write(*,*) fname
      read(fname,*) (stack(i),i=1,nstack)
      write(*,*)'source group to cut:',(stack(i),i=1,nstack)

c      enddo


 
 1001 format(i10,1x,i3,1x,i3,1x,a20,1x,a20,1x,f4.3,f4.3,a50)
 1002 format(i10,1x,i3,1x,i3,1x,a20,1x,a20,1x,a50)
 1050 format(i8,1x,i3,1x,i3,1x,a15,1x,a15,1x,4(f15.5,1x),2i4)


C   READ AND WRITE OUT FILE HEADER INFORMATION

      READ (7) IFILE,NOTE,NSEG,NSPECS,IDATE,BEGTIM,JDATE,ENDTIM
      READ (7) ORGX,ORGY,IZONE,UTMX,UTMY,DELTAX,DELTAY,NX,NY,
     $ NZ,NZLOWR,NZUPPR,HTSUR,HTLOW,HTUPP

      WRITE(9) IFILE,NOTE,NSEG,NSPECS,IDATE,BEGTIM,JDATE,ENDTIM
      WRITE(9) ORGX,ORGY,IZONE,UTMX,UTMY,DELTAX,DELTAY,NX,NY,
     $ NZ,NZLOWR,NZUPPR,HTSUR,HTLOW,HTUPP

      WRITE (*,1007) IDATE, BEGTIM, JDATE, ENDTIM 
      WRITE (*,*)'  '

C      IF (NSEG.EQ.0) GO TO 11    !no segments? 

C  SEGMENT ORIGIN AND SEGMENT SIZE  

      DO I=1,NSEG
      READ  (7) IX,IY,NXCLL,NYCLL
      WRITE (9) IX,IY,NXCLL,NYCLL
      WRITE (*,*) IX,IY,NXCLL,NYCLL 
      ENDDO

      IF (NSPECS.LE.0) STOP !no species? 

C  NAMES OF ALL SPECIES PRESENT  

      READ  (7) ((MSPEC(I,J),I=1,10),J=1,NSPECS)
      WRITE (9) ((MSPEC(I,J),I=1,10),J=1,NSPECS)
      WRITE (*,1013) ((MSPEC(I,J),I=1,10),J=1,NSPECS)

      IF (NSEG.LE.0) STOP !no segments, no emissions

C  TIME INVARIANT DATA
 
      DO I=1,NSEG
      READ (7) ISEGM,NPMAX
      WRITE(9) ISEGM,NPMAX
      write(*,*)'Total # of stacks',NPMAX
      IF (NPMAX.LE.0) STOP !no points--no emissions

C  FOR EACH POINT SOURCE LOCATION AND STACK PARAMETERS 

      READ (7) (DUMX(II),DUMY(II),IDUMX(II),IDUMY(II),DUMXX(II),
     &   DUMYY(II),II=1,NPMAX)
      WRITE(9) (DUMX(II),DUMY(II),IDUMX(II),IDUMY(II),DUMXX(II),
     &   DUMYY(II),II=1,NPMAX)

      ENDDO !end loop through segments

C    LOOP FOR 24 HOURS OF PT SOURCE DATA 

      DO IH = 1,24

          READ  (7) IBGDAT,BEGTIM,IENDAT,ENDTIM 
          WRITE (9) IBGDAT,BEGTIM,IENDAT,ENDTIM      
          WRITE(*,1007) IBGDAT,BEGTIM,IENDAT,ENDTIM,IH

C    GO THROUGH NUMBER OF SEGMENTS 

          IF (NSEG.LE.0) STOP
   
          DO J=1,NSEG

             READ (7) ISEGNM,NUMPTS
             WRITE(9) ISEGNM,NUMPTS

             IF (NUMPTS.LE.0) STOP

             READ(7) (ICELL(II),JCELL(II),KCELL(II),FLOW(II), 
     &          PLUMHT(II),II=1,NUMPTS) 


             do ii=1,NUMPTS
 
             if(nstack.eq.1) then 

               if(KCELL(ii).eq.-14) then  !approach 1
                  KCELL(ii) = -9
               elseif(KCELL(ii).eq.-17) then 
                  KCELL(ii) = -2
               elseif(KCELL(ii).eq.-22) then 
                  KCELL(ii) = -3
               elseif(KCELL(ii).eq.-35) then 
                  KCELL(ii) = -4
               elseif(KCELL(ii).eq.-38) then 
                  KCELL(ii) = -5
               elseif(KCELL(ii).eq.-43) then 
                  KCELL(ii) = -6
               elseif(KCELL(ii).eq.-46) then 
                  KCELL(ii) = -7
               elseif(KCELL(ii).eq.-48) then 
                  KCELL(ii) = -8
               else
                  KCELL(ii) = -1
               endif

              elseif(nstack.eq.2) then
               if(KCELL(ii).eq.-2) then  !approach 2
                  KCELL(ii) = -9
               elseif(KCELL(ii).eq.-4) then 
                  KCELL(ii) = -2
               elseif(KCELL(ii).eq.-10) then 
                  KCELL(ii) = -3
               elseif(KCELL(ii).eq.-13) then 
                  KCELL(ii) = -4
               elseif(KCELL(ii).eq.-18) then 
                  KCELL(ii) = -5
               elseif(KCELL(ii).eq.-24) then 
                  KCELL(ii) = -6
               elseif(KCELL(ii).eq.-25) then 
                  KCELL(ii) = -7
               elseif(KCELL(ii).eq.-36) then 
                  KCELL(ii) = -8
               else
                  KCELL(ii) = -1 
               endif
 
              elseif(nstack.eq.3) then
               if(KCELL(ii).eq.-7) then  !approach 3
                  KCELL(ii) = -9
               elseif(KCELL(ii).eq.-8) then 
                  KCELL(ii) = -2
               elseif(KCELL(ii).eq.-16) then 
                  KCELL(ii) = -3
               elseif(KCELL(ii).eq.-20) then 
                  KCELL(ii) = -4
               elseif(KCELL(ii).eq.-30) then 
                  KCELL(ii) = -5
               elseif(KCELL(ii).eq.-32) then 
                  KCELL(ii) = -6
               elseif(KCELL(ii).eq.-42) then 
                  KCELL(ii) = -7
               elseif(KCELL(ii).eq.-49) then 
                  KCELL(ii) = -8
               else
                  KCELL(ii) = -1
               endif
  
             elseif(nstack.eq.4) then
              if(KCELL(ii).eq.-15) then
                 KCELL(ii) = -2
              else
                 KCELL(ii) = -1
              endif
 
             endif  

             if(IH.eq.1) then
             WRITE(20,1058)ii,KCELL(II),emscut(ii)
             endif

             enddo

 1058        format(i8,1x,i4,1x,f5.3)
 1059        format(i8,1x,i4,1x,i4)


             WRITE(9)(ICELL(II),JCELL(II),KCELL(II),FLOW(II),
     &          PLUMHT(II),II=1,NUMPTS)

               
             IF (NSPECS.LE.0) STOP  

              DO K=1,NSPECS  
                READ (7) ISEGNM,(SPNAME(II,K),II=1,10),
     &           (EMISS(II),II=1,NUMPTS)  
                WRITE(9) ISEGNM,(SPNAME(II,K),II=1,10),
     &           (EMISS(II),II=1,NUMPTS)  

              ENDDO !end nspecs loop 
            ENDDO !end NSEG loop
          ENDDO !end loop over 24 hours

 
  999 STOP
 
 1003 FORMAT(10A1)
 1004 FORMAT(4I10)
 1005 FORMAT(6D20.7)
 1006 FORMAT(10A1,E18.9) 
 1007 FORMAT(2(I10,F10.2),I10)
 1008 FORMAT(3I12,2E18.9) 
 1009 FORMAT(I10,10A1)
 1010 FORMAT(9E14.9)  
 1013 FORMAT(1X,10A1) 
 1021 FORMAT(2I10)  
 2001 FORMAT(F10.1,1X,F10.1,1X,I3,F10.1,1X,F10.1,1X,2F20.10,5I4,3F7.0)
 2100 FORMAT(10A1,60A1,/,I2,1X,I2,1X,I6,F6.0,I6,F6.0)     
      END 

      subroutine lcpgeo(iway,phic,xlonc,truelat1,truelat2,xloc,yloc,
     &                  xlon,ylat)
c
c     LCPGEO performs Lambert Conformal to geodetic (lat/lon) translation
c
c     Code based on the TERRAIN preprocessor for MM5 v2.0,
c     developed by Yong-Run Guo and Sue Chen, National Center for
c     Atmospheric Research, and Pennsylvania State University
c     10/21/1993
c
c     Input arguments:
c        iway                Conversion type
c                            0 = geodetic to Lambert Conformal
c                            1 = Lambert Conformal to geodetic
c        phic                Central latitude (deg, neg for southern hem)
c        xlonc               Central longitude (deg, neg for western hem)
c        truelat1            First true latitute (deg, neg for southern hem)
c        truelat2            Second true latitute (deg, neg for southern hem)
c        xloc/yloc           Projection coordinates (km)
c        xlon/ylat           Longitude/Latitude (deg)
c
c     Output arguments:
c        xloc/yloc           Projection coordinates (km)
c        xlon/ylat           Longitude/Latitude (deg)
c
      data conv/57.29578/, a/6370./
c
c-----Entry Point
c
      if (phic.lt.0) then
        sign = -1.
      else
        sign = 1.
      endif
      pole = 90.
      if (abs(truelat1).gt.90.) then
        truelat1 = 60.
        truelat2 = 30.
        truelat1 = sign*truelat1
        truelat2 = sign*truelat2
      endif
      xn = alog10(cos(truelat1/conv)) - alog10(cos(truelat2/conv))
      xn = xn/(alog10(tan((45. - sign*truelat1/2.)/conv)) -
     &         alog10(tan((45. - sign*truelat2/2.)/conv)))           
      psi1 = 90. - sign*truelat1
      psi1 = psi1/conv
      if (phic.lt.0.) then
        psi1 = -psi1
        pole = -pole
      endif
      psi0 = (pole - phic)/conv
      xc = 0.
      yc = -a/xn*sin(psi1)*(tan(psi0/2.)/tan(psi1/2.))**xn
c
c-----Calculate lat/lon of the point (xloc,yloc)
c
      if (iway.eq.1) then
        xloc = xloc + xc
        yloc = yloc + yc
        if (yloc.eq.0.) then
          if (xloc.ge.0.) flp = 90./conv
          if (xloc.lt.0.) flp = -90./conv
        else
          if (phic.lt.0.) then
            flp = atan2(xloc,yloc)
          else
            flp = atan2(xloc,-yloc)
          endif
        endif
        flpp = (flp/xn)*conv + xlonc
        if (flpp.lt.-180.) flpp = flpp + 360.
        if (flpp.gt. 180.) flpp = flpp - 360. 
        xlon = flpp 
c
        r = sqrt(xloc*xloc + yloc*yloc)
        if (phic.lt.0.) r = -r
        cell = (r*xn)/(a*sin(psi1))
        rxn  = 1.0/xn
        cel1 = tan(psi1/2.)*cell**rxn
        cel2 = atan(cel1)
        psx  = 2.*cel2*conv
        ylat = pole - psx
c
c-----Calculate x/y from lat/lon
c
      else
        ylon = xlon - xlonc
        if (ylon.gt. 180.) ylon = ylon - 360.
        if (ylon.lt.-180.) ylon = ylon + 360.
        flp = xn*ylon/conv
        psx = (pole - ylat)/conv
        r = -a/xn*sin(psi1)*(tan(psx/2.)/tan(psi1/2.))**xn
        if (phic.lt.0.) then
          xloc = r*sin(flp)
          yloc = r*cos(flp)
        else
          xloc = -r*sin(flp)
          yloc =  r*cos(flp)
        endif
      endif
c
      xloc = xloc - xc
      yloc = yloc - yc
c
      return
      end

      subroutine utmgeo(iway,iutmzon,rx4,ry4,rlon4,rlat4)
c  
c-----UTMGEO performs UTM to geodetic (lat/lon) translation:
c
c     This is a Fortran version of the BASIC program "Transverse Mercator
c     Conversion", Copyright 1986, Norman J. Berls (Stefan Musarra, 2/94)
c     Based on algorithm taken from "Map Projections Used by the USGS"
c     by John P. Snyder, Geological Survey Bulletin 1532, USDI.
c
c     Input arguments:  
c        iway                Conversion type
c                            0 = geodetic to UTM 
c                            1 = UTM to geodetic
c        iutmzon             UTM zone
c        rx4                 UTM easting (km) 
c        ry4                 UTM northing (km) 
c        rlon4               Longitude (deg, negative for W)
c        rlat4               Latitude (deg)
c              
c     Output arguments:  
c        rx4                 UTM easting (km) 
c        ry4                 UTM northing (km) 
c        rlon4               Longitude (deg)
c        rlat4               Latitude (deg)
c              
      implicit real*8 (a-h,o-z)
c
      real   rx4,ry4,rlon4,rlat4
      real*8 north
c      
      parameter(pi=3.14159265358979)
      parameter(degrad=pi/180., raddeg=1./degrad)
      parameter(semimaj=6378206.4, semimin=6356583.8)
c     parameter(e2=1.0-(semimin/semimaj)**2.)
c     parameter(e4=e2*e2, e6=e2*e4, ep2=e2/(1.-e2))
      parameter(scfa=.9996)
      parameter(north=0., east=500000.)
c
      e2=1.0-(semimin/semimaj)**2.0
      e4=e2*e2
      e6=e2*e4
      ep2=e2/(1.-e2)
c
c-----Entry point
c
c-----Convert inputs from single to double precision
c
      if (iway.eq.1) then
        xx = 1000.*rx4
        yy = 1000.*ry4
      else
        dlat = rlat4
        dlon = rlon4
      endif
c
c-----Set Zone parameters
c
      zone = iutmzon
      cm = zone*6.0 - 183.0
      cmr = cm*degrad
c
c-----Lat/Lon to UTM conversion
c
      if (iway.eq.0) then
	rlat = degrad*dlat
	rlon = degrad*dlon

	delam = dlon - cm
	if (delam.lt.-180.) delam = delam + 360.
	if (delam.gt.180.) delam = delam - 360.
	delam = delam*degrad
	
        f1 = (1. - e2/4. - 3.*e4/64. - 5.*e6/256)*rlat 
        f2 = 3.*e2/8. + 3.*e4/32. + 45.*e6/1024. 
        f2 = f2*sin(2.*rlat) 
        f3 = 15.*e4/256.*45.*e6/1024. 
        f3 = f3*sin(4.*rlat) 
        f4 = 35.*e6/3072. 
        f4 = f4*sin(6.*rlat) 
        rm = semimaj*(f1 - f2 + f3 - f4) 
        if (dlat.eq.90. .or. dlat.eq.-90.) then 
          xx = 0. 
          yy = scfa*rm 
        else 
          rn = semimaj/sqrt(1. - e2*sin(rlat)**2) 
          t = tan(rlat)**2 
          c = ep2*cos(rlat)**2 
          a = cos(rlat)*delam 
           
          f1 = (1. - t + c)*a**3/6. 
          f2 = 5. - 18.*t + t**2 + 72.*c - 58.*ep2 
          f2 = f2*a**5/120. 
          xx = scfa*rn*(a + f1 + f2) 
          f1 = a**2/2. 
          f2 = 5. - t + 9.*c + 4.*c**2 
          f2 = f2*a**4/24. 
          f3 = 61. - 58.*t + t**2 + 600.*c - 330.*ep2 
          f3 = f3*a**6/720. 
          yy = scfa*(rm + rn*tan(rlat)*(f1 + f2 + f3)) 
        endif
	xx = xx + east
	yy = yy + north
c
c-----UTM to Lat/Lon conversion
c
      else
        xx = xx - east 
        yy = yy - north 
        e1 = sqrt(1. - e2) 
        e1 = (1. - e1)/(1. + e1) 
        rm = yy/scfa 
        u = 1. - e2/4. - 3.*e4/64. - 5.*e6/256. 
        u = rm/(semimaj*u) 
         
        f1 = 3.*e1/2. - 27.*e1**3./32. 
        f1 = f1*sin(2.*u) 
        f2 = 21.*e1**2/16. - 55.*e1**4/32. 
        f2 = f2*sin(4.*u) 
        f3 = 151.*e1**3./96. 
        f3 = f3*sin(6.*u) 
        rlat1 = u + f1 + f2 + f3 
        dlat1 = rlat1*raddeg 
        if (dlat1.ge.90. .or. dlat1.le.-90.) then 
          dlat1 = dmin1(dlat1,dble(90.) ) 
          dlat1 = dmax1(dlat1,dble(-90.) ) 
          dlon = cm 
        else 
          c1 = ep2*cos(rlat1)**2. 
          t1 = tan(rlat1)**2. 
          f1 = 1. - e2*sin(rlat1)**2. 
          rn1 = semimaj/sqrt(f1) 
          r1 = semimaj*(1. - e2)/sqrt(f1**3) 
          d = xx/(rn1*scfa) 
           
          f1 = rn1*tan(rlat1)/r1 
          f2 = d**2/2. 
          f3 = 5.*3.*t1 + 10.*c1 - 4.*c1**2 - 9.*ep2 
          f3 = f3*d**2*d**2/24. 
          f4 = 61. + 90.*t1 + 298.*c1 + 45.*t1**2. - 252.*ep2 - 3.*c1**2
          f4 = f4*(d**2)**3./720. 
          rlat = rlat1 - f1*(f2 - f3 + f4) 
          dlat = rlat*raddeg 
           
          f1 = 1. + 2.*t1 + c1 
          f1 = f1*d**2*d/6. 
          f2 = 5. - 2.*c1 + 28.*t1 - 3.*c1**2 + 8.*ep2 + 24.*t1**2. 
          f2 = f2*(d**2)**2*d/120. 
          rlon = cmr + (d - f1 + f2)/cos(rlat1) 
          dlon = rlon*raddeg 
          if (dlon.lt.-180.) dlon = dlon + 360. 
          if (dlon.gt.180.) dlon = dlon - 360. 
        endif 
      endif
c
c-----Convert precision of outputs
c
      if (iway.eq.1) then
        rlat4 = REAL(dlat)
        rlon4 = REAL(dlon)
      else
        rx4 = REAL(xx/1000.)
        ry4 = REAL(yy/1000.)
      endif
c
      return
      end
