      program pthtq
      implicit none
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   Copyright (C) 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     PTHTQ adds "effective plume height" to a CAMx elevated point source
c     emissions file (CAMx PT) that was created by CMAQ2UAM processor.
c     The effective plume height is assigned to the center heihgt of the
c     cell where the virtual stack is located. The cell center height is
c     read from CMAQ METCRO3D file (ZH). The grid definition and time
c     interval of the METCRO3D file is assumed to be consistent with the
c     CAMx PT file.
c
c     This code is a sub-program of CMAQ2CAMx v2.
c
c     INPUT: via standard input
c       Input CAMx PT   - Input CAMx PT file name
c       Output CAMx PT  - Output CAMx PT file name
c
c     INPUT: via environment variable
c       CMAQ_METCRO3D   - CMAQ METCRO3D file name
c
c     OUTPUT:
c       CAMx PT file with "effective plume height"
c
c     HISTORY:
c       created by bkoo (02/10/2011)      - modified from PTHT
c
      include 'PARMS3.EXT'
      include 'IODECL3.EXT'
      include 'FDESC3.EXT'

      integer :: LOGUNIT
      integer :: JDATE,JTIME

      character(16), parameter :: INFILE = 'CMAQ_METCRO3D'
      character(16), parameter :: PGNAME = 'PTHTQ'


      integer, parameter :: jin1 = 10, jout = 20

      character(256) :: line, MESG
      character(4), dimension(10) :: name,lspec
      character(4), dimension(60) :: note
      integer :: nseg,nstk,nspec,iutm,nx,ny,nz,ibdate,iedate
      integer :: iseg,ixseg,iyseg,nxseg,nyseg,nzlowr,nzuppr
      real :: xorg,yorg,delx,dely,btime,etime
      real :: refx,refy,htsur,htlow,htupp

      character(4), allocatable :: mspec(:,:)
      integer, allocatable :: ijkstk(:,:)
      real, allocatable :: parstk(:,:)
      real, allocatable :: ptems(:)
      real, allocatable :: plmht(:)
      real, allocatable :: ZH(:,:,:)
      integer :: istat

      integer :: i,j,k,l,n
      logical :: lfirst
c
c     Initialize I/O-API
c
      LOGUNIT = INIT3()
      write(*,*) ' '
c
c     Open CMAQ METCRO3D file
c
      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
      write(*,*) ' '
c
c     Check header info
c
      if (FTYPE3D.ne.GRDDED3) then
        MESG = 'Data type of ' // TRIM(INFILE) // ' must be GRDDED3'
        call M3EXIT(PGNAME,0,0,MESG,2)
      endif
      write(*,*) ' '


11    format(20x,a)
c
c     Read input CAMx PT file name and open it
c
      write(*,*) 'Enter name of input CAMx PT file:'
      read(*,11) line
      write(*,*) TRIM(line)
      open(jin1,file=line,status='OLD',form='UNFORMATTED')
c
c     Read output CAMx PT file name and open it
c
      write(*,*) 'Enter name of output CAMx PT file:'
      read(*,11) line
      write(*,*) TRIM(line)
      open(jout,file=line,status='NEW',form='UNFORMATTED')
c
c     Read CAMx PT header portion
c
      read(jin1) name,note,nseg,nspec,ibdate,btime,iedate,etime
      read(jin1) refx,refy,iutm,xorg,yorg,delx,dely,nx,ny,nz,
     &           nzlowr,nzuppr,htsur,htlow,htupp
      read(jin1) ixseg,iyseg,nxseg,nyseg
c
c     Check grid definition
c
      if ( nx.NE.NCOLS3D .OR. ny.NE.NROWS3D ) then
        write(*,*) 'ERROR: grid dimension mismatch'
        write(*,*) ' Emis Cols & Rows = ',nx,ny
        write(*,*) '  Met Cols & Rows = ',NCOLS3D,NROWS3D
        stop
      endif
      if ( nz.GT.NLAYS3D ) then
        write(*,*) 'ERROR: insufficient met layers'
        write(*,*) ' Emis Layers = ',nz
        write(*,*) '  Met Layers = ',NLAYS3D
        stop
      endif
      if ( delx.NE.XCELL3D .OR. dely.NE.YCELL3D ) then
        write(*,*) 'ERROR: grid cell size mismatch'
        write(*,*) ' Emis delX & delY = ',delx,dely
        write(*,*) '  Met delX & delY = ',XCELL3D,YCELL3D
        stop
      endif
      if ( xorg.NE.XORIG3D .OR. yorg.NE.YORIG3D ) then
        write(*,*) 'ERROR: grid origin mismatch'
        write(*,*) ' Emis Xorg & Yorg = ',xorg,yorg
        write(*,*) '  Met Xorg & Yorg = ',XORIG3D,YORIG3D
        stop
      endif


      allocate (mspec(10,nspec), stat = istat)
      if (istat.ne.0) stop'ERROR: memory allocation failed - MSPEC'

      read(jin1) ((mspec(n,l),n=1,10),l=1,nspec)
      read(jin1) iseg,nstk

      allocate (parstk(6,nstk), stat = istat)
      if (istat.ne.0) stop'ERROR: memory allocation failed - PARSTK'
      allocate (ijkstk(3,nstk), stat = istat)
      if (istat.ne.0) stop'ERROR: memory allocation failed - IJKSTK'
      allocate (ptems(nstk), stat = istat)
      if (istat.ne.0) stop'ERROR: memory allocation failed - PTEMS'
      allocate (ZH(NCOLS3D,NROWS3D,NLAYS3D), stat = istat)
      if (istat.ne.0) stop'ERROR: memory allocation failed - ZH'
      allocate (plmht(nstk), stat = istat)
      if (istat.ne.0) stop'ERROR: memory allocation failed - PLMHT'


      read(jin1) ((parstk(k,n),k=1,6),n=1,nstk)
c
c     Write CAMx PT header portion
c
      write(jout)name,note,nseg,nspec,ibdate,btime,iedate,etime
      write(jout)refx,refy,iutm,xorg,yorg,delx,dely,nx,ny,nz,
     &           nzlowr,nzuppr,htsur,htlow,htupp
      write(jout)ixseg,iyseg,nxseg,nyseg
      write(jout)((mspec(n,l),n=1,10),l=1,nspec)
      write(jout)iseg,nstk
      write(jout)((parstk(k,n),k=1,6),n=1,nstk)
c
c     Time-variant portion
c
      JDATE = SDATE3D
      JTIME = STIME3D

      lfirst = .true.

100   read(jin1,end=900) ibdate,btime,iedate,etime
      write(jout)        ibdate,btime,iedate,etime
      write(*,'(2(i,f))')ibdate,btime,iedate,etime

      read(jin1) iseg,nstk
      write(jout)iseg,nstk

      read(jin1) ((ijkstk(j,n),j=1,3),(parstk(k,n),k=1,2),n=1,nstk)

      if (lfirst) then
        if (.NOT.READ3(INFILE,'ZH',ALLAYS3,JDATE,JTIME,ZH)) then
          MESG = 'Cannot read data from ' // TRIM(INFILE)
          call M3EXIT(PGNAME,JDATE,JTIME,MESG,1)
        endif
        do n = 1, nstk
          if (ijkstk(3,n).eq.1) then
c            write(*,*) 'ERROR: elevated source (#',n,')'
c            write(*,*) '       is put in the ground level'
c            stop
          plmht(n) = 1.0
          else
c          endif
          plmht(n) = ZH(ijkstk(1,n),ijkstk(2,n),ijkstk(3,n))
          endif 
        enddo
        lfirst = .false.
      endif ! lfirst

      call NEXTIME(JDATE,JTIME,TSTEP3D)

      if (.NOT.READ3(INFILE,'ZH',ALLAYS3,JDATE,JTIME,ZH)) then
        MESG = 'Cannot read data from ' // TRIM(INFILE)
        call M3EXIT(PGNAME,JDATE,JTIME,MESG,1)
      endif
      do n = 1, nstk
          if(ijkstk(3,n).eq.1) then
          parstk(2,n) = -1.0
          else
        parstk(2,n) = plmht(n)
        plmht(n) = ZH(ijkstk(1,n),ijkstk(2,n),ijkstk(3,n))
        parstk(2,n) = -0.5 * ( parstk(2,n) + plmht(n) ) ! "Negative" time avg
          endif
      enddo

      write(jout)((ijkstk(j,n),j=1,3),(parstk(k,n),k=1,2),n=1,nstk)

      do l = 1, nspec
        read(jin1) iseg,(lspec(j),j=1,10),(ptems(n),n=1,nstk)
        write(jout)iseg,(lspec(j),j=1,10),(ptems(n),n=1,nstk)
      enddo

      goto 100

900   write(*,*) 'End of CAMx PT file'
c
c     Close files
c
      close(jin1)
      close(jout)

      end

