      program mrgspc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   Copyright (C) 2008-2010  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  DESCRIPTION:
c
c     MRGSPC merges together multiple CAMx avrg files which have
c     different sets of species.
c
c  REVISON HISTORY:
c
c     created by bkoo (06/03/2008)
c     modified by bkoo (08/07/2008) - added file type check
c                                   - enhanced error reporting
c     modified by bkoo (01/16/2008) - OUTPUT now copies F01 header
c                                   - added scaling factors
c
      implicit none

      integer, parameter :: MXFILE = 99 ! maximum # of input files
      integer, parameter :: MXNSP = 999 ! maximum # of species
      integer, parameter :: jo = 10
      integer :: ji


      character(4):: name(10),note(60),lspec(10),name2(10),note2(60)
      integer :: nseg,nspec(MXFILE),ibdate,iedate,ibdate2,iedate2
      integer :: iutm,nx,ny,nz,iutm2,nx2,ny2,nz2
      integer :: nzlowr,nzuppr,ixseg,iyseg,nxseg,nyseg,iseg
      real :: btime,etime,btime2,etime2
      real :: xorg,yorg,delx,dely,xorg2,yorg2,delx2,dely2
      real :: refx,refy,htsur,htlow,htupp

      character(4), allocatable :: mspec(:,:)
      real, allocatable :: cin(:,:), cout(:,:,:,:)
      integer :: istat

      character(256) :: line
      character(10) :: tmpnam,spcnam(MXNSP)
      real :: sfac(MXFILE)
      integer :: idx(MXNSP,MXFILE)

      integer :: ifile,nfile,isp,nsp
      integer :: i,j,k,l,l2,n,len

c
c     open input file
c
      ifile = 0
 10   read(*,'(20x,a)') line
      line = ADJUSTL(line)
      if (TRIM(line).eq.'END') goto 100
      ifile = ifile + 1
      if (ifile.gt.MXFILE) then
        write(*,*) 'ERROR: exceeded maximum # of files - ' //
     &             'increase MXFILE to ',ifile
        stop
      endif
      sfac(ifile) = 1.0
      len = INDEX(line,',')
      if ( len.gt.0 ) then
        read(line(:len-1),*) sfac(ifile)
        line = ADJUSTL(line(len+1:))
      endif
      ji = jo + ifile
      write(*,'(a,i2.2,2a,2x,f10.4)')
     &           ' Opening input file (F',ifile,'): ',
     &                                     TRIM(line),sfac(ifile)
      open(unit=ji,file=line,status='old',form='UNFORMATTED')
c
c     read header
c
      read(ji) name,note,nseg,nspec(ifile),ibdate,btime,iedate,etime
      read(ji) refx,refy,iutm,xorg,yorg,delx,dely,nx,ny,nz,
     &         nzlowr,nzuppr,htsur,htlow,htupp
      read(ji) ixseg,iyseg,nxseg,nyseg

      write(tmpnam,'(10a1)') (name(n),n=1,10)
      if (tmpnam.ne.'AVERAGE   ' .and.
     &    tmpnam.ne.'EMISSIONS ' .and.
     &    tmpnam.ne.'AIRQUALITY' .and.
     &    tmpnam.ne.'INSTANT   ') then
        write(*,*) 'ERROR: unsupported file type - ',tmpnam
        stop
      else if (tmpnam.eq.'EMISSIONS ' .and. nz.ne.1) then
        write(*,*) 'WARNING: surface emissions file with invalid NZ -',
     &             nz
        write(*,*) '         NZ is set to 1'
        nz = 1
      endif

      allocate ( mspec(10,nspec(ifile)), stat = istat )
      if (istat.ne.0) stop'ERROR: memory allocation failed - mspec'
      read(ji) ((mspec(n,l),n=1,10),l=1,nspec(ifile))
c
c     check input consistency
c
      if (ifile.eq.1) then
        name2 = name
        note2 = note
        ibdate2 = ibdate
        btime2  = btime
        iedate2 = iedate
        etime2  = etime
        iutm2 = iutm
        xorg2 = xorg
        yorg2 = yorg
        delx2 = delx
        dely2 = dely
        nx2   = nx
        ny2   = ny
        nz2   = nz

        isp = nspec(ifile)
        if (isp.gt.MXNSP) then
          write(*,*) 'ERROR: exceeded maximum # of species - ' //
     &               'increase MXNSP to ',isp
          stop
        endif
        do l = 1, nspec(ifile)
          write(tmpnam,'(10a1)') (mspec(n,l),n=1,10)
          spcnam(l) = tmpnam
          idx(l,ifile) = l
        enddo
      else
        if (ibdate.ne.ibdate2 .or.
     &      iedate.ne.iedate2 .or.
     &      ABS(btime-btime2).gt.0.001 .or.
     &      ABS(etime-etime2).gt.0.001) then
          write(*,*) 'WARNING: date/time mismatch (start,end)'
          write(*,'(a,2(i10,f10.0))') '  F01 -',ibdate2,btime2,iedate2,etime2
          write(*,'(a,i2.2,a,2(i10,f10.0))') '  F',ifile,' -',
     &                                        ibdate,btime,iedate,etime
cbk          stop
        endif
        if (iutm.ne.iutm2) then
          write(*,*) 'WARNING: utm zone mismatch'
          write(*,'(a,i5)') '  F01 -',iutm2
          write(*,'(a,i2.2,a,i5)') '  F',ifile,' -',iutm
cbk          stop
        endif
        if (xorg.ne.xorg2 .or. yorg.ne.yorg2) then
          write(*,*) 'ERROR: grid origin mismatch (Xorg,Yorg)'
          write(*,'(a,2F10.4)') '  F01 -',xorg2,yorg2
          write(*,'(a,i2.2,a,2F10.4)') '  F',ifile,' -',xorg,yorg
          stop
        endif
        if (delx.ne.delx2 .or. dely.ne.dely2) then
          write(*,*) 'ERROR: grid cell size mismatch (delX,delY)'
          write(*,'(a,2F10.4)') '  F01 -',delx2,dely2
          write(*,'(a,i2.2,a,2F10.4)') '  F',ifile,' -',delx,dely
          stop
        endif
        if (nx.ne.nx2. or. ny.ne.ny2) then
          write(*,*) 'ERROR: grid dimension mismatch (Column,Row)'
          write(*,'(a,2i5)') '  F01 -',nx2,ny2
          write(*,'(a,i2.2,a,2i5)') '  F',ifile,' -',nx,ny
          stop
        endif
        if (nz.ne.nz2) then
          write(*,*) 'ERROR: number of layers mismatch'
          write(*,'(a,i5)') '  F01 -',nz2
          write(*,'(a,i2.2,a,i5)') '  F',ifile,' -',nz
          stop
        endif

        do l = 1, nspec(ifile)
          write(tmpnam,'(10a1)') (mspec(n,l),n=1,10)
          do l2 = 1, isp
            if (tmpnam.ne.spcnam(l2)) CYCLE
            idx(l,ifile) = l2
            goto 50
          enddo
          isp = isp + 1
          if (isp.gt.MXNSP) then
            write(*,*) 'ERROR: exceeded maximum # of species - ' //
     &                 'increase MXNSP to ',isp
            stop
          endif
          spcnam(isp) = tmpnam
          idx(l,ifile) = isp
 50       continue
        enddo
      endif

      deallocate(mspec)

      goto 10

 100  continue

      nfile = ifile
      nsp   = isp
      write(*,*) 'NFILE = ',nfile
      write(*,*) 'NSP   = ',nsp
c
c     merged species list
c
      allocate ( mspec(10,nsp), stat = istat )
      if (istat.ne.0) stop'ERROR: memory allocation failed - mspec'

      write(*,'(a,$)') ' ID  SPECIES   '
      do ifile = 1, nfile
        write(*,'(a3,i2.2,$)') '  F',ifile
      enddo
      write(*,*)
      do isp = 1, nsp
        read(spcnam(isp),'(10a1)') (mspec(n,isp),n=1,10)
        write(*,'(i3,2x,10a1,$)') isp,(mspec(n,isp),n=1,10)
        do ifile = 1, nfile
          do l = 1, nspec(ifile)
            if (idx(l,ifile).eq.isp) then
              write(*,'(i5,$)') l
              goto 110
            endif
          enddo
          write(*,'(a5,$)') '     '
 110      continue
        enddo
        write(*,*)
      enddo
c
c     open output file & write header
c
      read(*,'(20x,a)') line
      write(*,*)'Opening output file: ',TRIM(line)
      open(unit=jo,file=line,status='new',form='UNFORMATTED')

      write(jo) name2,note2,nseg,nsp,ibdate2,btime2,iedate2,etime2
      write(jo) refx,refy,iutm2,xorg2,yorg2,delx2,dely2,nx2,ny2,nz2,
     &          nzlowr,nzuppr,htsur,htlow,htupp
      write(jo) ixseg,iyseg,nxseg,nyseg

      write(jo) ((mspec(n,isp),n=1,10),isp=1,nsp)
c
c     time-variant portion
c
      allocate ( cin(nx,ny), stat = istat )
      if (istat.ne.0) stop'ERROR: memory allocation failed - cin'
      allocate ( cout(nx,ny,nz,nsp), stat = istat )
      if (istat.ne.0) stop'ERROR: memory allocation failed - cout'

 200  ji = jo + 1
      read(ji,end=900) ibdate,btime,iedate,etime
      write(*,'(2(i10,f10.0))') ibdate,btime,iedate,etime
      write(jo) ibdate,btime,iedate,etime
      do ifile = 2, nfile
        ji = jo + ifile
        read(ji) ibdate,btime,iedate,etime
      enddo

      cout = 0.0
      do ifile = 1, nfile
        ji = jo + ifile
        do l = 1, nspec(ifile)
          do k = 1, nz
            read(ji) iseg,(lspec(n),n=1,10),((cin(i,j),i=1,nx),j=1,ny)
            cout(:,:,k,idx(l,ifile)) = cout(:,:,k,idx(l,ifile)) + 
     &                                 cin * sfac(ifile)
          enddo
        enddo
      enddo

      do isp = 1, nsp
        do k = 1, nz
          write(jo) iseg,(mspec(n,isp),n=1,10),((cout(i,j,k,isp),
     &                                                 i=1,nx),j=1,ny)
        enddo
      enddo

      goto 200
 900  write(*,*) 'End Of File'
c
c     close files
c
      close(jo)
      do ifile = 1, nfile
        ji = jo + ifile
        close(ji)
      enddo

      end

