program cmaq2asc
!----------------------------------------------------------------------
!   CMAQ2ASC 1.0 Convert 3D CMAQ conc to gridded ascii file
!----------------------------------------------------------------------

c *** Variable declarations
      implicit none

      include 'PARMS3.EXT'
      include 'IODECL3.EXT'
      include 'FDESC3.EXT'

      INTEGER :: NCOLS, NROWS, NLAYS, NVARS, JDATE, JTIME, RUNLEN
      INTEGER :: TSTEP, NSTEPS, I,V, VMAX, UMAX, DMAX,dayflag,fillflag
      CHARACTER*16    VNAME( MXVARS3 ) !  list of vble names, from user
      CHARACTER*16    UNITS( MXVARS3 ) !  list of vble units
      CHARACTER*80    VDESC( MXVARS3 ) !  list of vble descs
      CHARACTER*256   MESG    !  buffer for m3exit(), etc
      CHARACTER*16    INAME   !  logical name of the input file
      LOGICAL         EFLAG !flag: error has happened

      REAL    :: rlon0,rlat0,tlat1,tlat2,xloc,yloc,rlat,rlon,xpos,ypos
      INTEGER :: xorg,yorg,nx,ny,j,deltax,colrow,bigyear,bigflag
      integer    c, r, k, layer, hour, ihour, idate, n, nfiles, tday
      real cfac, lat, long,junk
      integer saveday(365),z,NRECS,tempcell,useryear,coordflag,numberdays

      REAL, ALLOCATABLE :: TOTS( :,:,: ),small(:,:),clat(:,:),clon(:,:)
      REAL, ALLOCATABLE :: GRID( :,:,:),ALTGRID(:,:),ozone(:,:,:)
      REAL, ALLOCATABLE :: nh4(:,:,:),no3(:,:,:),so4(:,:,:),ec(:,:,:)
      REAL, ALLOCATABLE :: pm25(:,:,:),cm(:,:,:),crustal(:,:,:),oc(:,:,:)

      CHARACTER(LEN=16) :: progname
      CHARACTER(LEN=16) :: met, outfile
      CHARACTER(LEN=180) :: ename(366),fname,land

      INTEGER :: istatus
      INTEGER :: TRIMLEN
      real       envreal

      EXTERNAL TRIMLEN, envreal

c *** Initialize variables
      progname = 'CMAQ2ASC'
      cfac = (1./1000.) * 0.001102293 * 3600. !convert from g to tons and sec to hour

      read(*,'(a)') fname
      read(fname,*) fillflag
      write(*,*) 'Fill in boundary cells with nearby cells (1=YES,0=NO):',fillflag

      read(*,'(a)') fname
      read(fname,*) bigflag
      write(*,*) 'Ouput large (complete) SMAT-CE input file (1=YES,0=NO):',bigflag

      read(*,'(a)') fname
      read(fname,*) useryear
      write(*,*) 'Year for SMAT-CE input file:',useryear

      read(*,'(a)') fname
      read(fname,*) rlon0,rlat0,tlat1,tlat2,nx,ny,xorg,yorg,deltax
      write(*,'(a,t20,3f10.0)') 'Projection:',rlon0,rlat0,tlat1,tlat2
      write(*,*) 'X and Y cells:',nx,ny
      write(*,*) 'x and y orig:',xorg,yorg,deltax

      read(*,'(a)') fname
      open(11,file=fname,recl=1000)
      write(*,*)'Openend output text file:',fname

      read(*,'(a)') fname
      open(21,file=fname,recl=1000)
      write(*,*)'Openend small output SMAT-CE file:',fname

      read(*,'(a)') fname
      open(12,file=fname,form='formatted')
      write(*,*)'Openend file with cells for small output SMAT-CE file:',fname

      read(*,'(a)') fname
      if(fname.eq.'NONE'.or.fname.eq.'FALSE'.or.fname.eq.'') then
       coordflag = 0
      else
       coordflag = 1
      open(14,file=fname,form='formatted')
      write(*,*)'Opened optional cell coordinates file:',fname
      endif

      read(*,'(a)') outfile
      write(*,*)'IOAPI output file:', outfile

      read(*,'(a)') fname
      read(fname,*) numberdays

      write(*,*)'Number of input days:', numberdays

      read(*,'(a)') fname
      read(fname,*) nfiles

      do n = 1, nfiles
      read(*,'(a)') ename(n)
      enddo

      do n = 1,nfiles
      write(*,*)n,ename(n)
      enddo

c----- header info for text output file

      if(bigflag.eq.1) then
      write(11,1998)
      write(11,1999)
      endif

      write(21,1998)
      write(21,1999)

 1998 format('Day')
 1999 format('_ID, _TYPE, LAT, LONG, DATE, CRUSTAL, NH4, SO4, EC, NO3, OC, PM25, CM')

c---- allocate and clear arrays

      allocate ( GRID(nx, ny, 1) )
      allocate ( ALTGRID(nx, ny) )
      allocate ( ozone(nx,ny,numberdays))
      allocate ( no3(nx,ny,numberdays))
      allocate ( so4(nx,ny,numberdays))
      allocate ( ec(nx,ny,numberdays))
      allocate ( crustal(nx,ny,numberdays))
      allocate ( cm(nx,ny,numberdays))
      allocate ( oc(nx,ny,numberdays))
      allocate ( nh4(nx,ny,numberdays))
      allocate ( pm25(nx,ny,numberdays))
      allocate ( TOTS(nx, ny ,500) )
      allocate ( small(nx,ny))
      allocate ( clat(nx,ny))
      allocate ( clon(nx,ny))


           DO R = 1, ny
            DO C = 1, nx
                small(C,R) = 0.
                do V = 1, 365
                ozone(C,R,V) = 0.
                enddo
             DO V = 1,500
                TOTS( C,R,V ) = 0.0
             ENDDO
            ENDDO
           ENDDO

c---- read in list of cells for small file


 100  read(12,'(i6)',end=200) tempcell
      i = int(tempcell/1000)
      j = tempcell - i*1000
c      write(*,*) tempcell,i,j
      if(i.gt.0.and.j.gt.0.and.i.le.nx.and.j.le.ny) then
      small(i,j) = 1.
      endif
      goto 100
 200  continue

c--- read in optional list of cell coordinates

      if(coordflag.ne.0) then

      read(14,'(a)')
      read(14,'(a)')
      read(14,'(a)') !skip 3 header lines

 101  read(14,'(a)',end=201) fname
      read(fname,*) j,i,lat,long

c 101  read(14,'(i3,i3,f10.6,f10.6)',end=201) j,i,lat,long
c      write(*,*) i,j,lat,long,nx,ny
      clat(i,j) = lat
      clon(i,j) = -(abs(long))
      goto 101
 201  continue

      endif

c---- open ioapi files

      dayflag = 0

      do N = 1, nfiles
      land = ename(N)

      if ( .not. open3( land, FSREAD3, progname ) ) THEN
         MESG = 'Could not open file "' //
     &   land( 1: TRIMLEN(land))
     &   // '" for input'
         CALL M3EXIT( progname, 0, 0, MESG, 2 )
      end if

      IF ( .NOT. DESC3(land))THEN
         MESG = 'Could not get description info for file "' //
     &            land( 1: TRIMLEN( land) ) //'"'
         CALL M3EXIT( progname, 0, 0, MESG, 2 )
      ENDIF

        NCOLS = NCOLS3D
        NROWS = NROWS3D
        NLAYS = NLAYS3D
        TSTEP = TSTEP3D
        NRECS = MXREC3D

        write(*,*) NCOLS,NROWS,NLAYS,TSTEP

C.......   Get max string-lengths for use in variables-listing:

        VMAX = TRIMLEN( VNAME3D( 1 ) )
        UMAX = TRIMLEN( UNITS3D( 1 ) )
        DMAX = TRIMLEN( VDESC3D( 1 ) )
        DO  I = 1, NVARS3D
            VMAX = MAX( VMAX , TRIMLEN( VNAME3D( I ) ) )
            UMAX = MAX( UMAX , TRIMLEN( UNITS3D( I ) ) )
            DMAX = MAX( DMAX , TRIMLEN( VDESC3D( I ) ) )
        END DO

        WRITE( *,92000 )
     &  ' ', 'The list of variables in this file is:', ' ',
     &  ( VNAME3D( I )( 1:VMAX ) // ' (' //
     &    UNITS3D( I )( 1:UMAX ) // '): ' //
     &    VDESC3D( I )( 1:DMAX ), I = 1, NVARS3D )

        write(*,*) NVARS3D

        IF ( NVARS3D .EQ. 1 ) THEN

            NVARS = 1
            VNAME( NVARS ) = VNAME3D( 1 )
            UNITS( NVARS ) = UNITS3D( 1 )
            VDESC( NVARS ) = VDESC3D( 1 )

            IF ( VTYPE3D( 1 ) .NE. M3REAL ) THEN
                MESG = 'Variable "' //
     &                  VNAME3D( 1 )( 1: TRIMLEN( VNAME3D( 1 ) ) )//
     &                 '" not of type REAL; ' //
     &                 'VERTOT processes REAL only'
                CALL M3EXIT( progname, 0, 0, MESG, 2 )
            END IF

        ELSE    !  else nvars3d > 1:

                    NVARS = NVARS3D
                    DO  V = 1, NVARS3D
                        VNAME( V ) = VNAME3D( V )
                        UNITS( V ) = UNITS3D( V )
                        VDESC( V ) = VDESC3D( V )
                        EFLAG = EFLAG .OR. ( VTYPE3D( V ) .NE. M3REAL )
                    END DO
        ENDIF !end loop over variables on file

c---- read in variables and do averaging

       do z = 1, NRECS

       JDATE = SDATE3D
       JTIME = STIME3D

       tday = JDATE - int(JDATE/1000)*1000
c       saveday(tday) = SDATE3D
       dayflag=dayflag + 1
       saveday(dayflag) = useryear*1000 + tday

       write(*,*) SDATE3D,STIME3D,tday,dayflag,saveday(dayflag)

        DO  V = 1, NVARS

            IF ( .NOT. READ3( land, VNAME( V ), ALLAYS3,
     &                        JDATE, JTIME, GRID ) ) THEN
                MESG = 'Read failure:  file ' // land //
     &                 ' variable ' // VNAME( V )
                CALL M3EXIT( 'VERTOT:VERSTEP', JDATE, JTIME,
     &                       MESG, 2 )
            END IF      !  if read3() failed

           DO R = 1, NROWS
            DO C = 1, NCOLS
                TOTS( C,R,V ) = TOTS(C,R,V) + GRID( C,R,1 )

                if(VNAME(V).eq.'PM25_NO3') then
                no3( C,R, dayflag ) =  GRID( C,R,1 )
                elseif(VNAME(V).eq.'PM25_NH4') then
                nh4( C,R, dayflag ) =  GRID( C,R,1 )
                elseif(VNAME(V).eq.'PM25_SO4') then
                so4( C,R, dayflag ) =  GRID( C,R,1 )
                elseif(VNAME(V).eq.'PM25_EC') then
                ec( C,R, dayflag) =  GRID( C,R,1 )
                elseif(VNAME(V).eq.'PM25_OM') then
                oc( C,R, dayflag) =  GRID( C,R,1 )
                elseif(VNAME(V).eq.'Fine_CRUSTAL') then
                crustal( C,R, dayflag ) =  GRID( C,R,1 )
                elseif(VNAME(V).eq.'PM25_bulk') then
                pm25( C,R, dayflag ) =  GRID( C,R,1 )
                elseif(VNAME(V).eq.'Coarse_CRUSTAL') then
                cm( C,R, dayflag ) =  GRID( C,R,1 )
                endif

            ENDDO
           ENDDO

        ENDDO !end loop over variables on file

        SDATE3D = SDATE3D + 1
       enddo


       write(*,*) 'got here'

      if ( .not. close3 ( land ) ) THEN
         MESG = 'Could not close file'
         CALL M3EXIT( progname, 0, 0, MESG, 2 )
      end if

        ENDDO !end loop over number of input files

c----- for old CAMx applications fill in missing data around the domain border

       if (fillflag.gt.0) then

       do tday = 1,numberdays

        do j = 1, NROWS

         crustal(1,j,tday) = crustal(2,j,tday)
         crustal(NCOLS,j,tday) = crustal((NCOLS-1),j,tday)
         so4(1,j,tday) = so4(2,j,tday)
         so4(NCOLS,j,tday) = so4((NCOLS-1),j,tday)
         nh4(1,j,tday) = nh4(2,j,tday)
         nh4(NCOLS,j,tday) = nh4((NCOLS-1),j,tday)
         ec(1,j,tday) = ec(2,j,tday)
         ec(NCOLS,j,tday) = ec((NCOLS-1),j,tday)
         no3(1,j,tday) = no3(2,j,tday)
         no3(NCOLS,j,tday) = no3((NCOLS-1),j,tday)
         oc(1,j,tday) = oc(2,j,tday)
         oc(NCOLS,j,tday) = oc((NCOLS-1),j,tday)
         cm(1,j,tday) = cm(2,j,tday)
         cm(NCOLS,j,tday) = cm((NCOLS-1),j,tday)
         pm25(1,j,tday) = pm25(2,j,tday)
         pm25(NCOLS,j,tday) = pm25((NCOLS-1),j,tday)

        enddo

        do i = 1, NCOLS

         crustal(i,1,tday) = crustal(i,2,tday)
         crustal(i,NROWS,tday) = crustal(i,(NROWS-1), tday)
         so4(i,1,tday) = so4(i,2,tday)
         so4(i,NROWS,tday) = so4(i,(NROWS-1), tday)
         nh4(i,1,tday) = nh4(i,2,tday)
         nh4(i,NROWS,tday) = nh4(i,(NROWS-1), tday)
         ec(i,1,tday) = ec(i,2,tday)
         ec(i,NROWS,tday) = ec(i,(NROWS-1), tday)
         no3(i,1,tday) = no3(i,2,tday)
         no3(i,NROWS,tday) = no3(i,(NROWS-1), tday)
         oc(i,1,tday) = oc(i,2,tday)
         oc(i,NROWS,tday) = oc(i,(NROWS-1), tday)
         cm(i,1,tday) = cm(i,2,tday)
         cm(i,NROWS,tday) = cm(i,(NROWS-1), tday)
         pm25(i,1,tday) = pm25(i,2,tday)
         pm25(i,NROWS,tday) = pm25(i,(NROWS-1), tday)

       enddo

       enddo

       endif


c------ output SMAT-CE input file(s)

       do tday = 1,numberdays
       do i = 1, NCOLS
        do j = 1, NROWS
         colrow = i*1000 + j

         if(coordflag.eq.0) then !estimate cell coordinates
          xpos = (XORG + (i - 0.5)*DELTAX)
          ypos = (YORG + (j - 0.5)*DELTAX)
          call lcpgeo(1,rlat0,rlon0,tlat1,tlat2,xpos,ypos,
     &                  long,lat)
         else !use optional input file with cell coordinates
          lat=clat(i,j)
          long=clon(i,j)
         endif

        bigyear = saveday(tday)
        call caldate(bigyear)

        if(bigflag.eq.1) then
        write(11,1117) colrow,lat,long,bigyear,crustal(i,j,tday),nh4(i,j,tday),so4(i,j,tday),ec(i,j,tday),
     &                no3(i,j,tday),oc(i,j,tday),pm25(i,j,tday),cm(i,j,tday)
        endif

        if(small(i,j).gt.0) then
        write(21,1117) colrow,lat,long,bigyear,crustal(i,j,tday),nh4(i,j,tday),so4(i,j,tday),ec(i,j,tday),
     &                no3(i,j,tday),oc(i,j,tday),pm25(i,j,tday),cm(i,j,tday)
        endif

        enddo
       enddo
       enddo
 1009  format(i3,',',i3,',D24HourMean',',QuarterlyMean,,"',3(f8.5,','),f8.5,'"')
 1007  format(i3,',',i3,',D24HourMean',',QuarterlyMean,Mean,"',f8.5,'"')
 1117  format(i6,',"",',f10.6,',',f11.6,',',i8,8(',',f12.6))


c------ create IOAPI output file

          DO V =1,NVARS
           DO R = 1, NROWS
            DO C = 1, NCOLS
                TOTS( C,R,V ) = TOTS(C,R,V) / nfiles
            ENDDO
          ENDDO
          ENDDO

c        IF ( TSTEP .EQ. 0 ) THEN
            JDATE  = 0
            JTIME  = 0
            NSTEPS = 1
c        ELSE

            SDATE3D = JDATE
            STIME3D = JTIME
            NLAYS3D = 1
            NVARS3D = NVARS
            DO I = 1, NVARS
                VNAME3D( I ) = VNAME( I )
                UNITS3D( I ) = UNITS( I )
                VDESC3D( I ) = VDESC( I )
                VTYPE3D( I ) = M3REAL
            END DO

      if ( .not. open3( outfile, FSUNKN3, progname ) ) THEN
          MESG = 'Could not open file "' //
     &     outfile( 1: TRIMLEN(outfile))
     &     // '" for output'
           CALL M3EXIT( progname, 0, 0, MESG, 2 )
      end if

          do V = 1,NVARS
           write(*,*)VNAME3D(V)

      IF ( .NOT. WRITE3( OUTFILE,VNAME(V),JDATE, JTIME, TOTS(1,1,V) )) THEN
                 MESG = 'Could not write totals to ' // VNAME(V)
                  CALL M3EXIT( 'VERTOT:VERSTEP', JDATE, JTIME,
     &                          MESG, 2 )
                END IF      !  if write3() failed

          enddo !end loop over variables


C******************  FORMAT  STATEMENTS   ******************************

C...........   Error and warning message formats..... 91xxx

91000   FORMAT ( //5X , '*** ERROR ABORT in program VERTOT ***',
     &            /5X , A ,
     &           // )        !  generic error message format

C...........   Informational (LOG) message formats... 92xxx

92000   FORMAT ( 5X , A )

92999   FORMAT ( //5X , A, // )

C...........   Formatted file I/O formats............ 93xxx

93000   FORMAT ( A16 )

C...........   Internal buffering formats............ 94xxx

C...........   Miscellaneous formats................. 95xxx

95000   FORMAT ( /5X , A , $ )          !  generic prompt format.


      stop

      END


c-----Start subroutines
      subroutine lcpgeo(iway,phic,xlonc,truelat1,truelat2,xloc,yloc,
     &                  xlon,ylat)
c
c     LCPGEO performs Lambert Conformal to geodetic (lat/lon) translation
c
c     Code based on the TERRAIN preprocessor for MM5 v2.0,
c     developed by Yong-Run Guo and Sue Chen, National Center for
c     Atmospheric Research, and Pennsylvania State University
c     10/21/1993
c
c     Input arguments:
c        iway                Conversion type
c                            0 = geodetic to Lambert Conformal
c                            1 = Lambert Conformal to geodetic
c        phic                Central latitude (deg, neg for southern hem)
c        xlonc               Central longitude (deg, neg for western hem)
c        truelat1            First true latitute (deg, neg for southern hem)
c        truelat2            Second true latitute (deg, neg for southern hem)
c        xloc/yloc           Projection coordinates (km)
c        xlon/ylat           Longitude/Latitude (deg)
c
c     Output arguments:
c        xloc/yloc           Projection coordinates (km)
c        xlon/ylat           Longitude/Latitude (deg)
c
      data conv/57.29578/, a/6370./
c
c-----Entry Point
c
      if (phic.lt.0) then
        sign = -1.
      else
        sign = 1.
      endif
      pole = 90.
      if (abs(truelat1).gt.90.) then
        truelat1 = 60.
        truelat2 = 30.
        truelat1 = sign*truelat1
        truelat2 = sign*truelat2
      endif
      xn = alog10(cos(truelat1/conv)) - alog10(cos(truelat2/conv))
      xn = xn/(alog10(tan((45. - sign*truelat1/2.)/conv)) -
     &         alog10(tan((45. - sign*truelat2/2.)/conv)))
      psi1 = 90. - sign*truelat1
      psi1 = psi1/conv
      if (phic.lt.0.) then
        psi1 = -psi1
        pole = -pole
      endif
      psi0 = (pole - phic)/conv
      xc = 0.
      yc = -a/xn*sin(psi1)*(tan(psi0/2.)/tan(psi1/2.))**xn
c
c-----Calculate lat/lon of the point (xloc,yloc)
c
      if (iway.eq.1) then
        xloc = xloc + xc
        yloc = yloc + yc
        if (yloc.eq.0.) then
          if (xloc.ge.0.) flp = 90./conv
          if (xloc.lt.0.) flp = -90./conv
        else
          if (phic.lt.0.) then
            flp = atan2(xloc,yloc)
          else
            flp = atan2(xloc,-yloc)
          endif
        endif
        flpp = (flp/xn)*conv + xlonc
        if (flpp.lt.-180.) flpp = flpp + 360.
        if (flpp.gt. 180.) flpp = flpp - 360.
        xlon = flpp
c
        r = sqrt(xloc*xloc + yloc*yloc)
        if (phic.lt.0.) r = -r
        cell = (r*xn)/(a*sin(psi1))
        rxn  = 1.0/xn
        cel1 = tan(psi1/2.)*cell**rxn
        cel2 = atan(cel1)
        psx  = 2.*cel2*conv
        ylat = pole - psx
c
c-----Calculate x/y from lat/lon
c
      else
        ylon = xlon - xlonc
        if (ylon.gt. 180.) ylon = ylon - 360.
        if (ylon.lt.-180.) ylon = ylon + 360.
        flp = xn*ylon/conv
        psx = (pole - ylat)/conv
        r = -a/xn*sin(psi1)*(tan(psx/2.)/tan(psi1/2.))**xn
        if (phic.lt.0.) then
          xloc = r*sin(flp)
          yloc = r*cos(flp)
        else
          xloc = -r*sin(flp)
          yloc =  r*cos(flp)
        endif
      endif
c
      xloc = xloc - xc
      yloc = yloc - yc
c
      return
      end

      subroutine caldate(idate)

      integer idate
      dimension nday(12)
      data nday/31,28,31,30,31,30,31,31,30,31,30,31/

c
c-----Entry point
c
c-----If it is already in calender date, return
c

      icent = int(idate/100000)
      iyear = int((idate - icent*100000)/1000)
      jday = idate - icent*100000 - iyear*1000

c     nday(2) = 28
c uncomment line below if base year is a leap year
c Do not use mod(year,4)=0. Future year will have 366 days if base year does.
      nday(2) = 29
      mday = 0
      do imonth = 1,12
        mday = mday + nday(imonth)
        if (mday.ge.jday) go to 20
      enddo
 20   iday = jday - (mday - nday(imonth))
      idate = icent*1000000 + iyear*10000 + imonth*100 + iday
c
c      write(*,'(4i6)') icent,iyear,imonth,iday

      return
      end