      program cmaq2uam
      implicit none
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   Copyright (C) 2006-2011  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 an environment 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       Optional input (environment variable)
c         STKGROUP      - logical name for CMAQ stack group file
c                         (required only if CMAQEMIS is CMAQ in-line
c                          point source emission file)
c
c       Following inputs via standard input
c         File Type     - input file type (ICON, BCON, EMIS or IPNT)
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       Flag for BCON interpolation (environment variable)
c         BCON_INTRP    - T or Y to interpolate BC; otherwise, F or N
c
c     OUTPUT:
c       OutFile1
c       OutFile2
c
c     HISTORY:
c       created by bkoo (04/14/2006)
c       modified by bkoo (06/28/2006)     - modified to process 2-D
c                                           emissions files
c       modified by bkoo (02/21/2008)     - added BC timestep intrp
c       modified by bkoo (06/30/2011)     - added option to process
c                                           in-line CMAQ point emiss
c
      include 'PARMS3.EXT'
      include 'IODECL3.EXT'
      include 'FDESC3.EXT'

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

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

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

      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(:), hstk(:), dstk(:)
      real, allocatable :: tstk(:), vstk(:), flow(:)
      real, allocatable :: conc(:,:,:,:)
      integer :: istat

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

      integer :: i,j,k,l,n,ihr,isub
      logical :: l2nd = .TRUE. ! OutFile2 is generated if TRUE
      logical :: lintrp = .FALSE. ! Hourly BC is interpolated if 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, EMIS or IPNT):'
      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
      elseif (typ.eq.'IPNT') then
        INFILE = SGNAME
      else
        write(*,*) 'ERROR: file type must be ICON, BCON, EMIS or IPNT'
        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'
          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')
        l2nd = .FALSE.
        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

        lintrp = ENVYN( 'BCON_INTRP', 'flag for BCON interpolation',
     &                  .false., istat )
        if (istat.gt.0) then
          MESG = 'Bad value for BCON_INTRP'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        if (lintrp) then
          nsubstep = TSTEP3D / 10000
        else
          nsubstep = 1
        endif

      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')
cKB        if (NLAYS3D.eq.1) l2nd = .FALSE.
        if (l2nd) 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

      elseif (typ.eq.'IPNT') then

        if (FTYPE3D.ne.GRDDED3) then
          MESG = 'Data type of ' // TRIM(INFILE) // ' must be GRDDED3'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        if (NCOLS3D.ne.1 .or. NLAYS3D.ne.1) then
          MESG = 'NCOLS & NLAYS of ' // TRIM(INFILE) // ' must be 1'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        if (MXREC3D.ne.1) then
          MESG = 'MXREC3D of ' // TRIM(INFILE) // ' must be 1'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif

        xorg = XORIG3D
        yorg = YORIG3D
        delx = XCELL3D
        dely = YCELL3D
        nstk = NROWS3D
        nz = 1
        nx = 0
        ny = 0
        do i = 1, MXDESC3
          j = INDEX(FDESC3D(i),'/NCOLS3D/')
          if (j.gt.0) read(FDESC3D(i)(10:),'(i)') nx
          j = INDEX(FDESC3D(i),'/NROWS3D/')
          if (j.gt.0) read(FDESC3D(i)(10:),'(i)') ny
        enddo
        if (nx.lt.1 .or. ny.lt.1) then
          MESG = 'Cannot find NCOLS & NROWS in the file description'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        write(*,*) ' '
        write(*,*) 'NX, NY     = ', nx, ny
        write(*,*) 'NSTK       = ', nstk
        write(*,*) 'XORG, YORG = ', xorg, yorg
        write(*,*) 'DELX, DELY = ', delx, dely

        allocate (xstk(nstk),ystk(nstk),hstk(nstk),dstk(nstk),
     &            tstk(nstk),vstk(nstk),flow(nstk), stat = istat)
        if (istat.ne.0) then
          MESG = 'Memory allocation failed: STK real parameters'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        allocate (istk(nstk),jstk(nstk),kstk(nstk), stat = istat)
        if (istat.ne.0) then
          MESG = 'Memory allocation failed: STK int parameters'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif

        if (.not.READ3(INFILE,'STKDM',1,SDATE3D,STIME3D,dstk)) then ! [m]
          MESG = 'Cannot read STKDM from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,SDATE3D,STIME3D,MESG,1)
        endif
        if (.not.READ3(INFILE,'STKHT',1,SDATE3D,STIME3D,hstk)) then ! [m]
          MESG = 'Cannot read STKHT from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,SDATE3D,STIME3D,MESG,1)
        endif
        if (.not.READ3(INFILE,'STKTK',1,SDATE3D,STIME3D,tstk)) then ! [K]
          MESG = 'Cannot read STKTK from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,SDATE3D,STIME3D,MESG,1)
        endif
        if (.not.READ3(INFILE,'STKVE',1,SDATE3D,STIME3D,vstk)) then ! [m/s]
          MESG = 'Cannot read STKVE from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,SDATE3D,STIME3D,MESG,1)
        endif
        vstk = vstk * 3600. ! [m/s] -> [m/hr]
        if (.not.READ3(INFILE,'STKFLW',1,SDATE3D,STIME3D,flow)) then ! [m^3/s]
          MESG = 'Cannot read STKFLW from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,SDATE3D,STIME3D,MESG,1)
        endif
        flow = flow * 3600. ! [m^3/s] -> [m^3/hr]
        if (.not.READ3(INFILE,'XLOCA',1,SDATE3D,STIME3D,xstk)) then ! [m]
          MESG = 'Cannot read XLOCA from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,SDATE3D,STIME3D,MESG,1)
        endif
        if (.not.READ3(INFILE,'YLOCA',1,SDATE3D,STIME3D,ystk)) then ! [m]
          MESG = 'Cannot read YLOCA from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,SDATE3D,STIME3D,MESG,1)
        endif
        if (.not.READ3(INFILE,'COL',1,SDATE3D,STIME3D,istk)) then
          MESG = 'Cannot read COL from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,SDATE3D,STIME3D,MESG,1)
        endif
        if (.not.READ3(INFILE,'ROW',1,SDATE3D,STIME3D,jstk)) then
          MESG = 'Cannot read ROW from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,SDATE3D,STIME3D,MESG,1)
        endif
        if (.not.READ3(INFILE,'ISTACK',1,SDATE3D,STIME3D,kstk)) then ! KB add source appt group ID
          MESG = 'Cannot read ROW from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,SDATE3D,STIME3D,MESG,1)
        endif
c        kstk = 0 ! No info from CMAQ stack group file
c
c     Open CMAQ in-line pnt emiss file
c
        INFILE = EMNAME

        if (.not.OPEN3(INFILE,FSREAD3,PGNAME)) then
          MESG = 'Cannot open ' // TRIM(INFILE)
          call M3EXIT(PGNAME,0,0,MESG,1)
        endif
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 (FTYPE3D.ne.GRDDED3) then
          MESG = 'Data type of ' // TRIM(INFILE) // ' must be GRDDED3'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        if (NCOLS3D.ne.1 .or. NLAYS3D.ne.1) then
          MESG = 'NCOLS & NLAYS of ' // TRIM(INFILE) // ' must be 1'
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        if (NROWS3D.ne.nstk) then
          MESG = 'Inconsistent NSTK in ' // TRIM(INFILE)
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        if (XORIG3D.ne.xorg .or. YORIG3D.ne.yorg) then
          MESG = 'Inconsistent XORG/YORG in ' // TRIM(INFILE)
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        if (XCELL3D.ne.delx .or. YCELL3D.ne.dely) then
          MESG = 'Inconsistent DELX/DELY in ' // TRIM(INFILE)
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif

        write(*,*) 'In-line PNT emiss file ',
     &             '(CAMx AR file is not generated):'
        read(*,11) MESG
        write(*,*) 'Enter name of CAMx PT file (OutFile2):'
        read(*,11) MESG
        write(*,*) TRIM(MESG)
        iout2 = JUNIT()
        open(iout2,file=MESG,status='NEW',form='UNFORMATTED')

      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
      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
        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
      nz = NLAYS3D
      if (typ.eq.'EMIS') nz = 1 ! For CAMx AR file
      if (typ.ne.'IPNT') then
        nx = NCOLS3D
        ny = NROWS3D
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)
      endif
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. l2nd) .or. typ.eq.'IPNT' ) 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)

        if (typ.ne.'IPNT') then
          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
      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),
     &            BUF2(SIZE,NLAYS3D,NVARS3D), stat = istat)
        if (istat.ne.0) then
          MESG = 'Memory allocation failed: BUFBND/BUF2'
          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) then

        if (typ.ne.'IPNT') 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 (l2nd) call setlstk(nx,ny,nz,lstk,BUFGRD(1,1,1,l)) ! Set virtual stacks
          enddo

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

            allocate (xstk(nstk),ystk(nstk),hstk(nstk),dstk(nstk),
     &                tstk(nstk),vstk(nstk),flow(nstk), stat = istat)
            if (istat.ne.0) then
              MESG = 'Memory allocation failed: STK real parameters'
              call M3EXIT(PGNAME,0,0,MESG,2)
            endif
            allocate (istk(nstk),jstk(nstk),kstk(nstk), stat = istat)
            if (istat.ne.0) then
              MESG = 'Memory allocation failed: STK int parameters'
              call M3EXIT(PGNAME,0,0,MESG,2)
            endif
            hstk = 1.0
            dstk = 1.0
            tstk = 273.0
            vstk = 1.0
            flow = 1.0

            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(n),dstk(n),
     &                                   tstk(n),vstk(n),n=1,nstk)
            write(iout2)ibdate,btime,iedate,etime
            write(iout2)iseg,nstk
            write(iout2)(istk(n),jstk(n),kstk(n),flow(n),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
        elseif (typ.eq.'IPNT') then
          write(iout2)iseg,nstk
          write(iout2)(xstk(n),ystk(n),hstk(n),dstk(n),
     &                                 tstk(n),vstk(n),n=1,nstk)
          write(iout2)ibdate,btime,iedate,etime
          write(iout2)iseg,nstk
          write(iout2)(istk(n),jstk(n),-(kstk(n)),flow(n),plmh,n=1,nstk)!set KSTK for source appt

          do l = 1, nspec
            write(iout2)iseg,(mspec(n,l),n=1,10),
     &                              (BUFGRD(1,n,1,l)*emsfac,n=1,nstk)
          enddo
        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 (l2nd) call setlstk(nx,ny,nz,lstk,BUFGRD(1,1,1,l)) ! Set virtual stacks
          enddo
        elseif (typ.eq.'IPNT') then
          allocate (conc(1,nstk,1,nspec), stat = istat)
          if (istat.ne.0) then
            MESG = 'Memory allocation failed: CONC'
            call M3EXIT(PGNAME,0,0,MESG,2)
          endif

          goto 100
        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)

          if (.not.lintrp) then
            write(iout1)ibdate,btime,iedate,etime
            write(*,'(2(i,f))')ibdate,btime,iedate,etime
          endif

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

            do isub = 1, nsubstep
            if (lintrp) then
              iedate = ibdate
              etime = btime + 1.0
              call corrdate(iedate,etime)
              write(iout1)ibdate,btime,iedate,etime
              write(*,'(2(i,f))')ibdate,btime,iedate,etime
              BUFBND = BUFBND * ( real(nsubstep-isub)
     &                               / real(nsubstep-isub+1) )
     &                 + BUF2 * ( 1. / real(nsubstep-isub+1) )
            else
              BUFBND = BUF2
            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
            if (lintrp) then
              ibdate = iedate
              btime = etime
            endif
            enddo ! isub
          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 (l2nd) 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
100     continue
        if ( (typ.eq.'EMIS' .and. l2nd) .or. typ.eq.'IPNT' ) then

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

            allocate (xstk(nstk),ystk(nstk),hstk(nstk),dstk(nstk),
     &                tstk(nstk),vstk(nstk),flow(nstk), stat = istat)
            if (istat.ne.0) then
              MESG = 'Memory allocation failed: STK real parameters'
              call M3EXIT(PGNAME,0,0,MESG,2)
            endif
            allocate (istk(nstk),jstk(nstk),kstk(nstk), stat = istat)
            if (istat.ne.0) then
              MESG = 'Memory allocation failed: STK int parameters'
              call M3EXIT(PGNAME,0,0,MESG,2)
            endif
            hstk = 1.0
            dstk = 1.0
            tstk = 273.0
            vstk = 1.0
            flow = 1.0

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

          write(iout2)iseg,nstk
          write(iout2)(xstk(n),ystk(n),hstk(n),dstk(n),
     &                                 tstk(n),vstk(n),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
          if (typ.eq.'IPNT') then
            do l = 1, nspec
              do n = 1, nstk
                conc(1,n,1,l) = 0.5 * BUFGRD(1,n,1,l)
              enddo
            enddo
          else
            do l = 1, nspec
              do k = 1, nz !KB
              do j = 1, ny
              do i = 1, nx
                conc(i,j,k,l) = 0.5 * BUFGRD(i,j,k,l)
              enddo
              enddo
              enddo
            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(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(n),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

            if (typ.eq.'IPNT') then
              do l = 1, nspec
                do n = 1, nstk
                  conc(1,n,1,l) = conc(1,n,1,l) + 0.5 * BUFGRD(1,n,1,l)
                enddo
                write(iout2)iseg,(mspec(n,l),n=1,10),
     &                                  (conc(1,n,1,l)*emsfac,n=1,nstk)
                do n = 1, nstk
                  conc(1,n,1,l) = 0.5 * BUFGRD(1,n,1,l)
                enddo
              enddo
            else
              do l = 1, nspec
c     Get average between hrs
                do k = 1, nz !KB
                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 = 1, nz !KB
                do j = 1, ny
                do i = 1, nx
                  conc(i,j,k,l) = 0.5 * BUFGRD(i,j,k,l)
                enddo
                enddo
                enddo
              enddo
            endif
          enddo ! ihr - CAMx PT

        endif

      endif ! TSTEP3D.eq.0
c
c     Close output files
c
      if (typ.ne.'IPNT') close(iout1)
      if (l2nd) 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 = 1, n3 !KBedit
      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

cKB      do k = 2, n3
      do k = 1, 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

