      program cmaq2uam
      implicit none
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   Copyright (C) 2003-2006  ENVIRON
c
c
c   This program is free software; you can redistribute it and/or
c   modify it under the terms of the GNU General Public License
c   as published by the Free Software Foundation; either version 2
c   of the License, or (at your option) any later version.
c
c   This program is distributed in the hope that it will be useful,
c   but WITHOUT ANY WARRANTY; without even the implied warranty of
c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c   GNU General Public License for more details.
c
c   To obtain a copy of the GNU General Public License
c   write to the Free Software Foundation, Inc.,
c   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
c
c
c   For comments and questions, send to bkoo@environcorp.com
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     CMAQ2UAM converts CMAQ input files (I/O-API) into corresponding
c     CAMx input files (UAM-IV format). It only converts file formats
c     and units if needed. No species mapping.
c
c     This code is a sub-program of CMAQ2CAMx v2.
c
c     INPUT:
c       One of the followings via environmental variable
c         CMAQICON      - logical name for CMAQ IC input file
c         CMAQBCON      - logical name for CMAQ BC input file
c         CMAQEMIS      - logical name for CMAQ Emissions input file
c
c       Following inputs via standard input
c         File Type     - input file type (ICON, BCON, or EMIS)
c         OutFile1      - first output file name (CAMx IC, BC, or AR)
c         OutFile2      - second output file name (CAMx TopC or PT)
c         Output Timezone - timezone for output (8 for PST, etc.)
c
c     OUTPUT:
c       OutFile1
c       OutFile2
c
c     HISTORY:
c       modified by bkoo (06/28/2006)     - modified to process 2-D
c                                           emissions files
c       created by bkoo (04/14/2006)
c
      include 'PARMS3.EXT'
      include 'IODECL3.EXT'
      include 'FDESC3.EXT'

      integer :: LOGUNIT
      integer :: JUNIT,SEC2TIME,TIME2SEC
      integer :: JDATE,JTIME,SIZE

      character(16), parameter :: ICNAME = 'CMAQICON'
      character(16), parameter :: BCNAME = 'CMAQBCON'
      character(16), parameter :: EMNAME = 'CMAQEMIS'
      character(16), parameter :: PGNAME = 'CMAQ2UAM'
      character(16) :: INFILE
      character(256) :: MESG

      real, allocatable :: BUFGRD(:,:,:,:), BUFBND(:,:,:)

      character(4), dimension(10) :: name
      character(4), dimension(60) :: note
      integer :: nspec,ibdate,iedate,nx,ny,nz,nstk
      integer, dimension(4) :: ncell
      real :: btime,etime,xorg,yorg,delx,dely,topc
      integer, parameter :: nseg = 1, iseg = 1, ixseg = 0, iyseg = 0
      integer, parameter :: iutm = 0 ! Always LCP
      integer, parameter :: nzlowr = 0, nzuppr = 0, izero = 0
      real, parameter :: refx = 0.0, refy = 0.0
      real, parameter :: htsur = 0.0, htlow = 0.0, htupp = 0.0
      character(10), parameter :: nmcamxic = 'AIRQUALITY'
      character(10), parameter :: nmcamxbc = 'BOUNDARY  '
      character(10), parameter :: nmcamxar = 'EMISSIONS '
      character(10), parameter :: nmcamxpt = 'PTSOURCE  '

      character(4), allocatable :: mspec(:,:)
      integer, allocatable :: icell(:,:)
      integer, allocatable :: lstk(:,:,:)
      integer, allocatable :: istk(:), jstk(:), kstk(:)
      real, allocatable :: xstk(:), ystk(:)
      real, allocatable :: conc(:,:,:,:)
      integer :: istat

      character(4) :: typ
      integer :: iout1,iout2
      integer :: itzone,nstep
      integer, parameter :: sdat0 = 1, edat0 = 99365
      real, parameter :: stim0 = 0.0, etim0 = 24.0
      real, parameter :: hstk = 1.0, dstk = 1.0, tstk = 273.0
      real, parameter :: vstk = 1.0, flow = 1.0, plmh = 1.0
      real, parameter :: emsfac = 3600. ! [*/s] * emsfac = [*/hr]

      integer :: i,j,k,l,n,ihr
      logical :: lpt = .TRUE.


11    format(20x,a)
12    format(20x,i)
c
c     Initialize I/O-API
c
      LOGUNIT = INIT3()
      write(*,*) ' '
c
c     Read file type
c
      write(*,*) 'Enter File Type to process (ICON, BCON, or EMIS):'
      read(*,11) typ
      write(*,*) typ
      if (typ.eq.'ICON') then
        INFILE = ICNAME
      elseif (typ.eq.'BCON') then
        INFILE = BCNAME
      elseif (typ.eq.'EMIS') then
        INFILE = EMNAME
      else
        write(*,*) 'ERROR: file type must be ICON, BCON, or EMIS'
        stop
      endif
c
c     Open CMAQ input file
c
      write(*,*) ' '
      if (.not.OPEN3(INFILE,FSREAD3,PGNAME)) then
        MESG = 'Cannot open ' // TRIM(INFILE)
        call M3EXIT(PGNAME,0,0,MESG,1)
      endif
      write(*,*) ' '
c
c     Get description of CMAQ input file
c
      if (.not.DESC3(INFILE)) then
        MESG = 'Cannot get description of ' // TRIM(INFILE)
        call M3EXIT(PGNAME,0,0,MESG,1)
      endif
c
c     Check header info
c
      if (GDTYP3D.ne.LAMGRD3) then
        MESG = 'Grid type of ' // TRIM(INFILE) // ' must be LAMGRD3'
        call M3EXIT(PGNAME,0,0,MESG,2)
      endif

      if (typ.eq.'ICON') then

        if (FTYPE3D.ne.GRDDED3) then
          MESG = 'Data type of ' // TRIM(INFILE) // ' must be GRDDED3'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        if (TSTEP3D.ne.0) then
          MESG = 'Time step of ' // TRIM(INFILE) // ' must be 0'
c          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        write(*,*) 'Enter name of CAMx IC file (OutFile1):'
        read(*,11) MESG
        write(*,*) TRIM(MESG)
        iout1 = JUNIT()
        open(iout1,file=MESG,status='NEW',form='UNFORMATTED')
        write(*,*) 'Enter name of CAMx TopC file (OutFile2):'
        read(*,11) MESG
        write(*,*) TRIM(MESG)
        iout2 = JUNIT()
        open(iout2,file=MESG,status='NEW',form='FORMATTED')
        do n = 1, 10
          name(n) = nmcamxic(n:n)
        enddo

      elseif (typ.eq.'BCON') then

        if (FTYPE3D.ne.BNDARY3) then
          MESG = 'Data type of ' // TRIM(INFILE) // ' must be BNDARY3'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        if (NTHIK3D.ne.1) then
          MESG = 'Boundary thickness must be 1'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        write(*,*) 'Enter name of CAMx BC file (OutFile1):'
        read(*,11) MESG
        write(*,*) TRIM(MESG)
        iout1 = JUNIT()
        open(iout1,file=MESG,status='NEW',form='UNFORMATTED')
        write(*,*) 'Enter blank (this line will be ignored)'
        read(*,11) MESG ! Second output file is ignored if BCON
        do n = 1, 10
          name(n) = nmcamxbc(n:n)
        enddo

      elseif (typ.eq.'EMIS') then

        if (FTYPE3D.ne.GRDDED3) then
          MESG = 'Data type of ' // TRIM(INFILE) // ' must be GRDDED3'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        write(*,*) 'Enter name of CAMx AR file (OutFile1):'
        read(*,11) MESG
        write(*,*) TRIM(MESG)
        iout1 = JUNIT()
        open(iout1,file=MESG,status='NEW',form='UNFORMATTED')
        if (NLAYS3D.eq.1) lpt = .FALSE.
        if (lpt) then
          write(*,*) 'Enter name of CAMx PT file (OutFile2):'
          read(*,11) MESG
          write(*,*) TRIM(MESG)
          iout2 = JUNIT()
          open(iout2,file=MESG,status='NEW',form='UNFORMATTED')
        else
          write(*,*) '2-D Emissions (CAMx PT file is not generated):'
          read(*,11) MESG
        endif
        do n = 1, 10
          name(n) = nmcamxar(n:n)
        enddo

      endif
c
c     Read output timezone
c
      write(*,*) 'Enter Output Timezone:'
      read(*,12) itzone
      if (TSTEP3D.eq.0) then
        write(*,*) 'Timezone is ignored for time-independent files'
        nstep = 1
        ibdate = sdat0
        iedate = edat0
        btime = stim0
        etime = etim0

c        if(etime.eq.0.and.btime.eq.0) etime = 1.

      else
        write(*,*) itzone
        if (MOD(STIME3D,10000).ne.0) then
          write(*,*)'ERROR: This program requires the CMAQ input to' //
     &              ' start at an exact hour'
          MESG = 'NO SUPPORT'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        nstep = MXREC3D - 1 ! CAMx will have one less records than CMAQ has
        JDATE = SDATE3D
        JTIME = STIME3D
        call NEXTIME(JDATE,JTIME,SEC2TIME(nstep*TIME2SEC(TSTEP3D)))
        ibdate = SDATE3D
        iedate = JDATE
        btime = real(STIME3D/10000-itzone) ! Always assume exact hr
        etime = real(JTIME/10000-itzone)   ! Always assume exact hr

c        if(etime.eq.0.and.btime.eq.0) etime = 1.

        call corrdate(ibdate,btime)
        call corrdate(iedate,etime)
      endif
      write(*,*) 'Output period (start date/time & end date/time):'
      write(*,'(2(i,f))') ibdate,btime,iedate,etime
c
c     Set CAMx species name
c
      allocate (mspec(10,NVARS3D), stat = istat)
      if (istat.ne.0) then
        MESG = 'Memory allocation failed: MSPEC'
        call M3EXIT(PGNAME,0,0,MESG,2)
      endif

      write(*,'(a4,1x,a)')' No.','CAMx Species Name'
      do l = 1, NVARS3D
        read(VNAME3D(l),'(10a1)') (mspec(n,l),n=1,10)
        write(*,'(i3,2x,10a1)') l,(mspec(n,l),n=1,10)
      enddo
c
c     Write CAMx header - set header info
c
      do n = 1, 60
        note(n) = FDESC3D(1)(n:n) ! First line of CMAQ file description
      enddo
      nspec = NVARS3D
      xorg = XORIG3D
      yorg = YORIG3D
      delx = XCELL3D
      dely = YCELL3D
      nx = NCOLS3D
      ny = NROWS3D
      nz = NLAYS3D
      if (typ.eq.'EMIS') nz = 1 ! For CAMx AR file
c
c     First Output (CAMx IC/BC/AR) - header
c
      write(iout1)name,note,nseg,nspec,ibdate,btime,iedate,etime
      write(iout1)refx,refy,iutm,xorg,yorg,delx,dely,nx,ny,nz,
     &            nzlowr,nzuppr,htsur,htlow,htupp
      write(iout1)ixseg,iyseg,nx,ny
      write(iout1)((mspec(n,l),n=1,10),l=1,nspec)
c
c     CAMx BC - additional header
c
      if (typ.eq.'BCON') then
        allocate (icell(max(nx,ny),4), stat = istat)
        if (istat.ne.0) then
          MESG = 'Memory allocation failed: ICELL'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif

        ncell(1) = ny ! West
        icell(1,1)  = 0 ! First row
        icell(ny,1) = 0 ! Last row
        ncell(2) = ny ! East
        icell(1,2)  = 0 ! First row
        icell(ny,2) = 0 ! Last row
        do i = 2, ny-1
          icell(i,1) = 2    ! West
          icell(i,2) = nx-1 ! East
        enddo
        ncell(3) = nx ! South
        icell(1,3)  = 0 ! First column
        icell(nx,3) = 0 ! Last column
        ncell(4) = nx ! North
        icell(1,4)  = 0 ! First column
        icell(nx,4) = 0 ! Last column
        do i = 2, nx-1
          icell(i,3) = 2    ! South
          icell(i,4) = ny-1 ! North
        enddo

        do j = 1, 4
          write(iout1)iseg,j,ncell(j),(icell(i,j),izero,izero,izero
     &                                           ,i=1,ncell(j))
        enddo

        deallocate (icell)
      endif
c
c     CAMx PT - header
c
      if (typ.eq.'EMIS' .and. lpt) then
        do n = 1, 10
          name(n) = nmcamxpt(n:n)
        enddo
        nz = NLAYS3D

        write(iout2)name,note,nseg,nspec,ibdate,btime,iedate,etime
        write(iout2)refx,refy,iutm,xorg,yorg,delx,dely,nx,ny,nz,
     &              nzlowr,nzuppr,htsur,htlow,htupp
        write(iout2)ixseg,iyseg,nx,ny
        write(iout2)((mspec(n,l),n=1,10),l=1,nspec)

        allocate (lstk(nx,ny,nz), stat = istat)
        if (istat.ne.0) then
          MESG = 'Memory allocation failed: LSTK'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        lstk = 0 ! Initialization
      endif
c
c     Read CMAQ input - hourly
c
      JDATE = SDATE3D
      JTIME = STIME3D
      
      if (typ.eq.'BCON') then
        SIZE = 2*NCOLS3D + 2*NROWS3D + 4
        allocate (BUFBND(SIZE,NLAYS3D,NVARS3D), stat = istat)
        if (istat.ne.0) then
          MESG = 'Memory allocation failed: BUFBND'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif

        if (.not.READ3(INFILE,ALLVAR3,ALLAYS3,JDATE,JTIME,BUFBND)) then
          MESG = 'Cannot read data from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,JDATE,JTIME,MESG,1)
        endif
      else
        allocate (BUFGRD(NCOLS3D,NROWS3D,NLAYS3D,NVARS3D), stat = istat)
        if (istat.ne.0) then
          MESG = 'Memory allocation failed: BUFGRD'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif

        if (.not.READ3(INFILE,ALLVAR3,ALLAYS3,JDATE,JTIME,BUFGRD)) then
          MESG = 'Cannot read data from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,JDATE,JTIME,MESG,1)
        endif
      endif
c
c     Time-independent file
c
      if (TSTEP3D.eq.0.or.typ.eq.'ICON') then

        write(iout1)ibdate,btime,iedate,etime

        if (typ.eq.'ICON') then
          do l = 1, nspec
            do k = 1, nz
              write(iout1)iseg,(mspec(n,l),n=1,10),
     &                               ((BUFGRD(i,j,k,l),i=1,nx),j=1,ny)
            enddo

            topc = 0.0
            do i = 1, nx
            do j = 1, ny
              topc = topc + BUFGRD(i,j,nz,l)
            enddo
            enddo
            topc = topc / real(nx*ny)
            write(iout2,'(10a1,e10.4)')(mspec(n,l),n=1,10),topc
          enddo
        elseif (typ.eq.'BCON') then
          do l = 1, nspec
            write(iout1)iseg,(mspec(n,l),n=1,10),1, ! West
     &                 ((BUFBND(2*nx+ny+4+i,k,l),k=1,nz),i=1,ncell(1))
            write(iout1)iseg,(mspec(n,l),n=1,10),2, ! East
     &                 ((BUFBND(nx+1+i,k,l),k=1,nz),i=1,ncell(2))
            write(iout1)iseg,(mspec(n,l),n=1,10),3, ! South
     &                 ((BUFBND(i,k,l),k=1,nz),i=1,ncell(3))
            write(iout1)iseg,(mspec(n,l),n=1,10),4, ! North
     &                 ((BUFBND(nx+ny+3+i,k,l),k=1,nz),i=1,ncell(4))
          enddo
        elseif (typ.eq.'EMIS') then
          do l = 1, nspec
            write(iout1)iseg,(mspec(n,l),n=1,10),
     &                        ((BUFGRD(i,j,1,l)*emsfac,i=1,nx),j=1,ny)
            if (lpt) call setlstk(nx,ny,nz,lstk,BUFGRD(1,1,1,l)) ! Set virtual stacks
          enddo

          if (lpt) then
            call cntstk(nx*ny*nz,lstk,nstk) ! Get nstk
            write(*,*) 'Number of virtual stacks: '
            write(*,*) nstk

            allocate (istk(nstk), stat = istat)
            if (istat.ne.0) then
              MESG = 'Memory allocation failed: ISTK'
              call M3EXIT(PGNAME,0,0,MESG,2)
            endif
            allocate (jstk(nstk), stat = istat)
            if (istat.ne.0) then
              MESG = 'Memory allocation failed: JSTK'
              call M3EXIT(PGNAME,0,0,MESG,2)
            endif
            allocate (kstk(nstk), stat = istat)
            if (istat.ne.0) then
              MESG = 'Memory allocation failed: KSTK'
              call M3EXIT(PGNAME,0,0,MESG,2)
            endif
            allocate (xstk(nstk), stat = istat)
            if (istat.ne.0) then
              MESG = 'Memory allocation failed: XSTK'
              call M3EXIT(PGNAME,0,0,MESG,2)
            endif
            allocate (ystk(nstk), stat = istat)
            if (istat.ne.0) then
              MESG = 'Memory allocation failed: YSTK'
              call M3EXIT(PGNAME,0,0,MESG,2)
            endif

            call setstkloc(nx,ny,nz,nstk,lstk,istk,jstk,kstk,xstk,ystk,
     &                     xorg,yorg,delx,dely) ! Set location of virtual stacks

            write(iout2)iseg,nstk
            write(iout2)(xstk(n),ystk(n),hstk,dstk,tstk,vstk,n=1,nstk)
            write(iout2)ibdate,btime,iedate,etime
            write(iout2)iseg,nstk
            write(iout2)(istk(n),jstk(n),kstk(n),flow,plmh,n=1,nstk)

            do l = 1, nspec
              write(iout2)iseg,(mspec(n,l),n=1,10),
     &             (BUFGRD(istk(n),jstk(n),kstk(n),l)*emsfac,n=1,nstk)
            enddo
          endif
        endif
c
c     Time-dependent file
c
      else

        if (typ.eq.'BCON') then
          allocate (conc(max(nx,ny),nz,4,nspec), stat = istat)
          if (istat.ne.0) then
            MESG = 'Memory allocation failed: CONC'
            call M3EXIT(PGNAME,0,0,MESG,2)
          endif
c     Initialize conc
          do l = 1, nspec
            do k = 1, nz
              do i = 1, ncell(1) ! West
                conc(i,k,1,l) = 0.5 * BUFBND(2*nx+ny+4+i,k,l)
              enddo
              do i = 1, ncell(2) ! East
                conc(i,k,2,l) = 0.5 * BUFBND(nx+1+i,k,l)
              enddo
              do i = 1, ncell(3) ! South
                conc(i,k,3,l) = 0.5 * BUFBND(i,k,l)
              enddo
              do i = 1, ncell(4) ! North
                conc(i,k,4,l) = 0.5 * BUFBND(nx+ny+3+i,k,l)
              enddo
            enddo
          enddo
        elseif (typ.eq.'EMIS') then
          allocate (conc(nx,ny,nz,nspec), stat = istat)
          if (istat.ne.0) then
            MESG = 'Memory allocation failed: CONC'
            call M3EXIT(PGNAME,0,0,MESG,2)
          endif
c     Initialize conc
          do l = 1, nspec
            do j = 1, ny
            do i = 1, nx
              conc(i,j,1,l) = 0.5 * BUFGRD(i,j,1,l)
            enddo
            enddo
            if (lpt) call setlstk(nx,ny,nz,lstk,BUFGRD(1,1,1,l)) ! Set virtual stacks
          enddo
        endif

        do ihr = 1, nstep
          ibdate = JDATE
          btime = real(JTIME/10000-itzone) ! Always assume exact hr
          call corrdate(ibdate,btime)

          call NEXTIME(JDATE,JTIME,TSTEP3D)

          iedate = JDATE
          etime = real(JTIME/10000-itzone) ! Always assume exact hr
          call corrdate(iedate,etime)

          write(iout1)ibdate,btime,iedate,etime
          write(*,'(2(i,f))')ibdate,btime,iedate,etime

          if (typ.eq.'BCON') then
            if (.not.
     &          READ3(INFILE,ALLVAR3,ALLAYS3,JDATE,JTIME,BUFBND)) then
              MESG = 'Cannot read data from ' // TRIM(INFILE)
              call M3EXIT(PGNAME,JDATE,JTIME,MESG,1)
            endif

            do l = 1, nspec
c     Get average between hrs
              do k = 1, nz
                do i = 1, ncell(1) ! West
                  conc(i,k,1,l) = conc(i,k,1,l)
     &                                 + 0.5 * BUFBND(2*nx+ny+4+i,k,l)
                enddo
                do i = 1, ncell(2) ! East
                  conc(i,k,2,l) = conc(i,k,2,l)
     &                                 + 0.5 * BUFBND(nx+1+i,k,l)
                enddo
                do i = 1, ncell(3) ! South
                  conc(i,k,3,l) = conc(i,k,3,l)
     &                                 + 0.5 * BUFBND(i,k,l)
                enddo
                do i = 1, ncell(4) ! North
                  conc(i,k,4,l) = conc(i,k,4,l)
     &                                 + 0.5 * BUFBND(nx+ny+3+i,k,l)
                enddo
              enddo
c     Write hourly outputs
              do j = 1, 4
                write(iout1)iseg,(mspec(n,l),n=1,10),j,
     &                           ((conc(i,k,j,l),k=1,nz),i=1,ncell(j))
              enddo
c     Re-initialize conc at each hr
              do k = 1, nz
                do i = 1, ncell(1) ! West
                  conc(i,k,1,l) = 0.5 * BUFBND(2*nx+ny+4+i,k,l)
                enddo
                do i = 1, ncell(2) ! East
                  conc(i,k,2,l) = 0.5 * BUFBND(nx+1+i,k,l)
                enddo
                do i = 1, ncell(3) ! South
                  conc(i,k,3,l) = 0.5 * BUFBND(i,k,l)
                enddo
                do i = 1, ncell(4) ! North
                  conc(i,k,4,l) = 0.5 * BUFBND(nx+ny+3+i,k,l)
                enddo
              enddo
            enddo
          elseif (typ.eq.'EMIS') then
            if (.not.
     &          READ3(INFILE,ALLVAR3,ALLAYS3,JDATE,JTIME,BUFGRD)) then
              MESG = 'Cannot read data from ' // TRIM(INFILE)
              call M3EXIT(PGNAME,JDATE,JTIME,MESG,1)
            endif

            do l = 1, nspec
c     Get average between hrs
              do j = 1, ny
              do i = 1, nx
                conc(i,j,1,l) = conc(i,j,1,l) + 0.5 * BUFGRD(i,j,1,l)
              enddo
              enddo
c     Write hourly outputs
              write(iout1)iseg,(mspec(n,l),n=1,10),
     &                          ((conc(i,j,1,l)*emsfac,i=1,nx),j=1,ny)
c     Re-initialize conc at each hr
              do j = 1, ny
              do i = 1, nx
                conc(i,j,1,l) = 0.5 * BUFGRD(i,j,1,l)
              enddo
              enddo
              if (lpt) call setlstk(nx,ny,nz,lstk,BUFGRD(1,1,1,l)) ! Set virtual stacks
            enddo
          endif
        enddo ! ihr
c
c     Now process for CAMx PT
c
        if (typ.eq.'EMIS' .and. lpt) then

          call cntstk(nx*ny*nz,lstk,nstk) ! Get nstk
          write(*,*) 'Number of virtual stacks: '
          write(*,*) nstk

          allocate (istk(nstk), stat = istat)
          if (istat.ne.0) then
            MESG = 'Memory allocation failed: ISTK'
            call M3EXIT(PGNAME,0,0,MESG,2)
          endif
          allocate (jstk(nstk), stat = istat)
          if (istat.ne.0) then
            MESG = 'Memory allocation failed: JSTK'
            call M3EXIT(PGNAME,0,0,MESG,2)
          endif
          allocate (kstk(nstk), stat = istat)
          if (istat.ne.0) then
            MESG = 'Memory allocation failed: KSTK'
            call M3EXIT(PGNAME,0,0,MESG,2)
          endif
          allocate (xstk(nstk), stat = istat)
          if (istat.ne.0) then
            MESG = 'Memory allocation failed: XSTK'
            call M3EXIT(PGNAME,0,0,MESG,2)
          endif
          allocate (ystk(nstk), stat = istat)
          if (istat.ne.0) then
            MESG = 'Memory allocation failed: YSTK'
            call M3EXIT(PGNAME,0,0,MESG,2)
          endif

          call setstkloc(nx,ny,nz,nstk,lstk,istk,jstk,kstk,xstk,ystk,
     &                   xorg,yorg,delx,dely) ! Set location of virtual stacks

          write(iout2)iseg,nstk
          write(iout2)(xstk(n),ystk(n),hstk,dstk,tstk,vstk,n=1,nstk)

c     Re-read CMAQ hourly emissions for CAMx PT
          JDATE = SDATE3D
          JTIME = STIME3D
          if (.not.
     &        READ3(INFILE,ALLVAR3,ALLAYS3,JDATE,JTIME,BUFGRD)) then
            MESG = 'Cannot read data from ' // TRIM(INFILE)
            call M3EXIT(PGNAME,JDATE,JTIME,MESG,1)
          endif
c     Initialize conc
          do l = 1, nspec
            do k = 2, nz
            do j = 1, ny
            do i = 1, nx
              conc(i,j,k,l) = 0.5 * BUFGRD(i,j,k,l)
            enddo
            enddo
            enddo
          enddo

          do ihr = 1, nstep
            ibdate = JDATE
            btime = real(JTIME/10000-itzone) ! Always assume exact hr
            call corrdate(ibdate,btime)

            call NEXTIME(JDATE,JTIME,TSTEP3D)

            iedate = JDATE
            etime = real(JTIME/10000-itzone) ! Always assume exact hr
            call corrdate(iedate,etime)

            write(iout2)ibdate,btime,iedate,etime
            write(*,'(2(i,f))')ibdate,btime,iedate,etime
            write(iout2)iseg,nstk
            write(iout2)(istk(n),jstk(n),kstk(n),flow,plmh,n=1,nstk)

            if (.not.
     &          READ3(INFILE,ALLVAR3,ALLAYS3,JDATE,JTIME,BUFGRD)) then
              MESG = 'Cannot read data from ' // TRIM(INFILE)
              call M3EXIT(PGNAME,JDATE,JTIME,MESG,1)
            endif

            do l = 1, nspec
c     Get average between hrs
              do k = 2, nz
              do j = 1, ny
              do i = 1, nx
                conc(i,j,k,l) = conc(i,j,k,l) + 0.5 * BUFGRD(i,j,k,l)
              enddo
              enddo
              enddo
c     Write hourly outputs
              write(iout2)iseg,(mspec(n,l),n=1,10),
     &               (conc(istk(n),jstk(n),kstk(n),l)*emsfac,n=1,nstk)
c     Re-initialize conc at each hr
              do k = 2, nz
              do j = 1, ny
              do i = 1, nx
                conc(i,j,k,l) = 0.5 * BUFGRD(i,j,k,l)
              enddo
              enddo
              enddo
            enddo
          enddo ! ihr - CAMx PT

        endif

      endif ! TSTEP3D.eq.0
c
c     Close output files
c
      close(iout1)
      if (typ.ne.'BCON' .and. lpt) close(iout2)

      end



      subroutine setstkloc(n1,n2,n3,n4,lflag,iloc,jloc,kloc,xloc,yloc,
     &                     x0,y0,dx,dy)
      implicit none
c
c     SETSTKLOC sets the location parameters of virtual stacks.
c
      integer :: n1, n2, n3, n4, i, j, k, n
      integer, dimension(n1,n2,n3) :: lflag
      integer, dimension(n4) :: iloc, jloc, kloc
      real, dimension(n4) :: xloc, yloc
      real :: x0, y0, dx, dy

      n = 0
      do k = 2, n3
      do j = 1, n2
      do i = 1, n1
        if (lflag(i,j,k).eq.1) then
          n = n + 1
          iloc(n) = i
          jloc(n) = j
          kloc(n) = k
          xloc(n) = (real(i)-0.5)*dx + x0
          yloc(n) = (real(j)-0.5)*dy + y0
        endif
      enddo
      enddo
      enddo

      if (n.ne.n4) then
        write(*,*) 'ERROR: incorrect stack counting'
        write(*,*) '       NSTK - ', n4
        write(*,*) '       N    - ', n
        stop
      endif

      return
      end


      subroutine cntstk(narray,iarray,nstack)
      implicit none
c
c     CNTSTK counts the number of virtual stacks.
c
      integer :: narray, nstack, i
      integer, dimension(narray) :: iarray

      nstack = 0
      do i = 1, narray
        if (iarray(i).eq.1) nstack = nstack + 1
      enddo

      return
      end


      subroutine setlstk(n1,n2,n3,lstack,rarray)
      implicit none
c
c     SETLSTK sets the flag for virtual stacks (LSTK).
c
      integer :: n1, n2, n3, i, j, k
      integer, dimension(n1,n2,n3) :: lstack
      real, dimension(n1,n2,n3) :: rarray

      do k = 2, n3
      do j = 1, n2
      do i = 1, n1
        if (rarray(i,j,k).gt.0.0) lstack(i,j,k) = 1
      enddo
      enddo
      enddo

      return
      end


      subroutine corrdate (idate,rtime)
      implicit none
c
c     CORRDATE checks the time and adjusts it if it is off the range
c     of 0-24. The julian date is also adjusted accordingly.
c
c     INPUT:
c       idate     - Julian date (YYYYJJJ)
c       rtime     - time (hr)
c
c     OUTPUT:
c       idate     - Julian date (YYJJJ)
c       rtime     - time (hr)
c
      integer :: idate
      real :: rtime

      if (rtime .lt. 0.) then
        rtime = rtime + 24.
        idate = idate - 1
      endif
      if (rtime .ge. 24.) then
        rtime = rtime - 24.
        idate = idate + 1
      endif
      if (MOD(idate,1000) .eq. 0) then
        if (MOD(INT(idate/1000)-1,4) .eq. 0) then
          idate = (INT(idate/1000)-1)*1000 + 366
        else
          idate = (INT(idate/1000)-1)*1000 + 365
        endif
      endif
      if (MOD(idate,1000) .gt. 365) then
        if (MOD(INT(idate/1000),4) .eq. 0) then
          if (MOD(idate,1000) .eq. 367)
     &    idate = (INT(idate/1000)+1)*1000 + 1
        else
          idate = (INT(idate/1000)+1)*1000 + 1
        endif
      endif

      idate=idate-idate/100000*100000 ! YYYYJJJ -> YYJJJ

      return
      end

