      program ncf2uam
      implicit none
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Copyright (C) 2003-2016  Ramboll Environ
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 go to the Free Software Foundation at http://www.fsf.org.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     NCF2UAM converts CMAQ IO-API files to CAMx file format (UAM).
c     It does not alter any variables or species. If a CMAQ
c     variable name is longer than 10 characters, it truncates the 
c     name except for CMAQ's SGTOT species: SGTOT species will lose
c     'TOT' in their names.
c     The output can have either all or one of the input variables.
c
c     INPUT:
c       Via environmental variable:
c         CMAQOUT       - logical name for CMAQ output file
c
c       Via standard input:
c         OutFile       - CAMx-formatted output file name
c         OutVariable   - Output variable name
c                         'ALL' to output all input variables
c         OutTimeZone   - Output time zone (8 for PST, etc.)
c
c     OUTPUT:
c       OutFile
c
c     HISTORY:
c       created by bkoo (09/02/2004)
c       revised by bkoo (05/26/2006)
c
      include 'PARMS3.EXT'
      include 'IODECL3.EXT'
      include 'FDESC3.EXT'

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

      character(16), parameter :: INFILE = 'CMAQOUT'
      character(16), parameter :: PGNAME = 'NCF2UAM'
      character(256) :: MESG

      character(4), dimension(10) :: name
      character(4), dimension(60) :: note
      integer :: nspec,ibdate,iedate,nx,ny,nz
      real :: btime,etime,xorg,yorg,delx,dely
      integer, parameter :: nseg = 1, iseg = 1, ixseg = 0, iyseg = 0
      integer, parameter :: iutm = 0 ! Always LCP
      integer, parameter :: nzlowr = 0, nzuppr = 0
      real, parameter :: refx = 0.0, refy = 0.0
      real, parameter :: htsur = 0.0, htlow = 0.0, htupp = 0.0
      character(10), parameter :: namecamx = 'AVERAGE'

      character(4), allocatable :: mspec(:,:)
      real, allocatable :: BUFF(:,:,:,:)
      integer :: istat

      character(10) :: outvar
      integer :: iout,idx,itzone,nstep

      integer :: i,j,k,l,lq,n,ihr


 11    format(20x,a)
 12    format(20x,i5)
c
c     Initialize I/O-API
c
      LOGUNIT = INIT3()
c
c     Open CMAQ output file
c
      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 output 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

      write(*,*) 'Enter name of output file (OutFile):'
      read(*,11) MESG
      write(*,*) TRIM(MESG)
      iout = JUNIT()
      open(iout,file=MESG,status='NEW',form='UNFORMATTED')

      do n = 1, 10
        name(n) = namecamx(n:n)
      enddo
      do n = 1, 60
        note(n) = FDESC3D(1)(n:n) ! First line of CMAQ file description
      enddo
c
c     Read output variable name
c
      write(*,*) 'Enter Output Variable (OutVariable):'
      read(*,11) outvar
      idx = 0
      if (outvar.eq.'ALL') then
        nspec = NVARS3D
      else
        do lq = 1, NVARS3D
          if (outvar.ne.VNAME3D(lq)) cycle
          idx = lq
          exit
        enddo
        if (idx.eq.0) then
          MESG = TRIM(outvar) // ' is not found in ' // TRIM(INFILE)
          call M3EXIT(PGNAME,0,0,MESG,2)
        endif
        nspec = 1
      endif
c
c     Set output variable name
c
      allocate (mspec(10,nspec), stat = istat)
      if (istat.ne.0) then
        MESG = 'Memory allocation failed: MSPEC'
        call M3EXIT(PGNAME,0,0,MESG,2)
      endif

      do l = 1, nspec
        lq = l
        if (idx.gt.0) lq = idx
        if (VNAME3D(lq)(1:6).eq.'SGTOT_') then
          read(VNAME3D(lq),'(2a1)') (mspec(n,l),n=1,2)
          read(VNAME3D(lq)(6:),'(8a1)') (mspec(n,l),n=3,10)
        else
          read(VNAME3D(lq),'(10a1)') (mspec(n,l),n=1,10)
        endif
        write(*,'(i3,2x,10a1)') l,(mspec(n,l),n=1,10)
      enddo
c
c     Read output timezone
c
      write(*,*) 'Enter Output Time Zone (OutTimeZone):'
      read(*,12) itzone
      write(*,*) itzone

      nstep = MXREC3D ! CAMx will have the same number of records as CMAQ
      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)

      write(*,*) 'Output period (start date/time & end date/time):'
      write(*,'(2(i10,f10.0))') ibdate,btime,iedate,etime
c
c     Write CAMx header - set header info
c
      xorg = XORIG3D
      yorg = YORIG3D
      delx = XCELL3D
      dely = YCELL3D
      nx = NCOLS3D
      ny = NROWS3D
      nz = NLAYS3D

      write(iout)name,note,nseg,nspec,ibdate,btime,iedate,etime
      write(iout)refx,refy,iutm,xorg,yorg,delx,dely,nx,ny,nz,
     &           nzlowr,nzuppr,htsur,htlow,htupp
      write(iout)ixseg,iyseg,nx,ny
      write(iout)((mspec(n,l),n=1,10),l=1,nspec)
c
c     Read CMAQ input - hourly
c
      allocate (BUFF(NCOLS3D,NROWS3D,NLAYS3D,NVARS3D), stat = istat)
      if (istat.ne.0) then
        MESG = 'Memory allocation failed: BUFF'
        call M3EXIT(PGNAME,0,0,MESG,2)
      endif

      JDATE = SDATE3D
      JTIME = STIME3D

      do ihr = 1, nstep

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

        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(iout)ibdate,btime,iedate,etime
        write(*,'(a,2(i10,f10.0))')'Writing ',ibdate,btime,iedate,etime

        do l = 1, nspec
          lq = l
          if (idx.gt.0) lq = idx
          do k = 1, nz
            write(iout)iseg,(mspec(n,l),n=1,10),
     &                               ((BUFF(i,j,k,lq),i=1,nx),j=1,ny)
          enddo
        enddo

      enddo
c
c     Close output files
c
      close(iout)

      MESG = 'Successful completion of ' // PGNAME
      call M3EXIT(PGNAME,0,0,MESG,0)
	
      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

