c     program to total concentrations across facilities by HAP
      module main1
      implicit none
      
      character hap*20,runtyp*10,st*2,faclist*250,logfil*250,ptcmaq*250,
     +recdir*250,hapdir*250,xdir*250,gridst(4)*4
     
      integer*8 ncat,ngrids,nblocks,nzblocks,nmons,logunit,gridres,
     +nruns,runno,zblks(4,2)  !these are the first occurences of non-pop blocks
      
      real*8, allocatable, dimension(:,:) :: gridconc  ! gridded AERMOD conc, #of receptors,AERMOD # of categories +1 (total)
      
      real*8, allocatable, dimension(:,:) :: monconc  ! monitor AERMOD conc, #of receptors,AERMOD # of categories +1 (total)
      
      real*8, allocatable, dimension(:,:) :: blokconc  ! block AERMOD conc, #of receptors,AERMOD # of categories +1 (total)
      
      real*8, allocatable, dimension(:,:) :: zblokconc  ! non-pop block AERMOD conc, #of receptors,AERMOD # of categories +1 (total); only used in non-CONUS

 
      logical runpt,runap,runon,runnon,runrwc,runcmv,runnphi,runnplow,
     +runog,runag,akhiprvi,noncontr,lgo,lgrids(4)
      
c     all gridded receptors
      type gridrecs
          character recid*15
          integer col,row
          integer*8 irec
      end type gridrecs
      
c     all block receptors
      type blkrecs
          character recid*15
          integer col,row
          integer*8 irec
      end type blkrecs

c     non-CONUS non-populated receptors
      type zblkrecs
          character recid*15
          integer*8 col,row
          integer*8 irec
      end type zblkrecs
c     all monitors
      type monrecs
          character recid*15
          integer col,row
          integer*8 irec
      end type monrecs

      type(gridrecs), dimension(:),allocatable:: gridrec
      type(blkrecs), dimension(:),allocatable:: blockrec
      type(zblkrecs), dimension(:),allocatable:: zblockrec
      type(monrecs), dimension(:),allocatable:: monrec
      data gridst /'12US','09AK','03HI','03PR'/
      parameter(logunit=12)
      end module main1
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      program totalconc
      use main1
      implicit none
      
      
c     get HAP, # of categories, runtype, and state (if non-CONUS gridded)
      call readargs
      
c     get ancillary directories/files needed for processing
      call getdir

c     get appropriate receptor arrays
      call rec_matrix
      
c     read AERMOD files and calculate total concentrations
      call read_post
      
c     write final output
      call writeout
      
      end
      
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc    
c     subroutine to read command line arguments
      subroutine readargs
      use main1
      implicit none
c     i:      loop counter
c     ilen:   length of argument
      integer i,ilen,istat,n1,n2,ndash
      
      character acat*2,runtyp1*12,logfile*250,agrid*2,noncon*1,arun*6,
     +arunno*6

      runpt=.false.
      runap=.false.
      runcmv=.false.
      runon=.false.
      runnon=.false.
      runag=.false.
      runog=.false.
      runrwc=.false.
      runnphi=.false.
      runnplow=.false.
      noncontr=.false.
      lgrids=.false.
      
      call getarg(1,hap)
      call getarg(2,acat)
      read(acat,*)ncat
      call getarg(3,runtyp1)
      call getarg(4,arunno)
      call getarg(5,arun)
      call getarg(6,st,istat)
      read(arun,*)nruns
      read(arunno,*)runno
      
      
      ilen=len_trim(hap)
      do i=1,ilen
          call upcase(hap(i:i))
      enddo
      
      ndash=0 
      ilen=len_trim(runtyp)
      do i=1,ilen
          call upcase(runtyp1(i:i))
          if (runtyp1(i:i) .eq. '_') ndash=ndash+1
      enddo
      
      if (ndash .eq. 0) then
          runtyp=trim(adjustl(runtyp1))
      elseif (ndash .eq. 1) then
          n1=index(runtyp1,'_')
          runtyp=trim(adjustl(runtyp1(1:n1-1)))
c          read(runtyp1(n1+1:ilen),*)gridres
          read(runtyp1(n1+1:ilen),*)agrid
          call checkgrid(agrid)
          if (gridres .ne. 4 .and. gridres .ne. 12 .and. gridres .ne. 3 
     +    .and. gridres .ne. 9) then
              write(*,5)agrid
              stop
          endif
      elseif (ndash .eq. 2) then
          n1=index(runtyp1,'_')
          n2=index(runtyp1(n1+1:ilen),'_')+n1
          runtyp=trim(adjustl(runtyp1(1:n1-1)))
c          read(runtyp1(n1+1:n2-1),*)gridres
          read(runtyp1(n1+1:n2-1),*)agrid
          call checkgrid(agrid)
          if (gridres .ne. 4 .and. gridres .ne. 12 .and. gridres .ne. 3
     +    .and. gridres .ne. 9) then
              write(*,5)agrid
              stop
          endif
          read(runtyp1(n2+1:ilen),'(a)')noncon
          ilen=len_trim(noncon)
          do i=1,ilen
              call upcase(noncon(i:i))
          enddo
          if (noncon .eq. 'Y') then
              noncontr=.true.
              akhiprvi=.true.
          elseif (noncon .eq. 'N') then
              noncontr=.false.
          else
              write(*,5)noncon
              stop
          endif
      else
          write(*,5)runtyp1
          stop
      endif
      
      if (runtyp .eq. 'PT') then
          runpt=.true.
      elseif (runtyp .eq. 'AP') then
          runap=.true.
      elseif (runtyp .eq. 'PORT' .or. runtyp .eq. 'UNDER' .or. 
     +runtyp .eq. 'CMV') then
          runcmv=.true.
      elseif (runtyp .eq. 'HDON' .or. runtyp .eq. 'HDOFF' .or. 
     +    runtyp .eq. 'LDON' .or. runtyp .eq. 'LDOFF' .or. 
     +    runtyp .eq. 'ONOFF' .or. runtyp .eq. 'HOTEL') then 
          runon=.true.
      elseif (runtyp .eq. 'NONRD') then
          runnon=.true.
      elseif (runtyp .eq. 'NPHI') then
          runnphi=.true.
      elseif (runtyp .eq. 'NPLO') then
          runnplow=.true.
      elseif (runtyp .eq. 'RWC') then
          runrwc=.true.
      elseif (runtyp .eq. 'OILGAS') then
          runog=.true.
      elseif (runtyp .eq. 'AG') then
          runag=.true.
      else
          write(*,5)runtyp
          stop
      endif
      
      if (istat .eq. -1) then
          write(logfil,'(4(a))')trim(adjustl(hap)),'_',
     +    trim(adjustl(runtyp)),'_totalconc.log'
          st=''
      else
          write(logfil,'(5(a))')trim(adjustl(hap)),'_',
     +    trim(adjustl(runtyp)),trim(adjustl(st)),'_totalconc.log'
      endif
      
      open(unit=logunit,file=logfil,status='replace')
      write(logunit,*)'Starting program version 18010'
      call datetime
 
5     format(a,' is not valid, stopping program')
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine getdir
      use main1
      implicit none
      
c     eof:        end of file indicator
c     n:          index variable
c     i:          loop counter
      integer*8 eof,n,nout,inpunit

c     line:       line of text read from input_directories.dat
c     tag:        source parameter file tag, point or area
c                 or airport runway or non-runway
c                 or port or underway for CMV
c     fname:      source parameter filename written to paramfil array
c     acol:       character string of CMAQ column
c     arow:       character string of CMAQ row
c     ares:       character string of grid resolution in km
      character line*300,infil*250,ffil*250,ares*2

c     lloc:       denotes if locations folder keyword found or locations file exists
c     lparam:     2D array denotes if parameter file(s) exist
c                 1D = point, 2D=area
c     lcmaq:      denotes if facility grid-cell xref folder keyword found
c     lmet:       denotes if met directory keyword found
c     lrec:       denotes if receptor directory keyword found
c     lout:       denotes if output directory keyword found
c     lparam1:    denotes if parameter folder keyword found 
c     ltemp1:     denotes if temporal folder keyword found 
c     lout:       denotes if output folder keyword found
c     lexist:     denotes that input_directories.dat is in directory
      logical lcmaq,lrec,lexist,lxwalk,lhap,lffile
      
c     initialize variables
      eof=0
      lcmaq=.false.
      lrec=.false.
      lxwalk=.false.
      lhap=.false.
c      lxfile=.false.
c     xfile1='NA'
      ffil='NA'
c      nx=0
c      nxfile=0
      lffile=.false.
c      paramfil='NA'
c      temfil='NA'
c      locfil='NA'
      
c     create facility ID for gridded sources based on column and row
!      if (runon .or. runnon .or. runnphi .or. runnplow .or. runrwc) then
!          write(acol,'(i3)')col
!          write(arow,'(i3)')row
!          write(facid,'(4a)')'G',trim(adjustl(acol)),'R',
!     +    trim(adjustl(arow))
!      endif
      
c     check to see if input_directories.dat is in directory
c     version 17037, now file has rungroup in name
      if (runpt .or. runcmv .or. runap) then
          write(infil,'(3(a))')'input_',trim(adjustl(runtyp)),
     +    '_directories.dat'
      else
          write(ares,'(i2)')gridres
          if (.not. noncontr) then
              write(infil,'(5(a))')'input_',trim(adjustl(runtyp)),
     +        '_',trim(adjustl(ares)),'_directories.dat'
          else
              write(infil,'(7(a))')'input_',trim(adjustl(runtyp)),
     +        '_',trim(adjustl(ares)),'_directories_',trim(adjustl(st)),
     +        '.dat'
          endif
      endif
      
      inpunit=11
      inquire(file=infil,exist=lexist)
      if (lexist) then
          open(unit=inpunit,file=infil,status='old')
      else
          write(logunit,*)trim(adjustl(infil)),' does not exist'
          write(logunit,50)
          stop
      endif
      
      
      read(11,5,iostat=eof)line
 10   if (eof .eq. 0) then
          if (index(line,'PTCMAQ ') .gt. 0) then
              n=index(line,'PTCMAQ ')
              ptcmaq=line(n+7:300)
              lcmaq=.true.
          endif
          if (index(line,'RECDIR ') .gt. 0) then
              n=index(line,'RECDIR ')
              recdir=line(n+7:300)
              lrec=.true.
          endif
          if (index(line,'XWALK ') .gt. 0) then  !emissions directory
              n=index(line,'XWALK ')
              xdir=line(n+6:300)
              lxwalk=.true.
          endif
          if (index(line,'HAPDIR ') .gt. 0) then
              n=index(line,'HAPDIR ')
              hapdir=line(n+7:300)
              lhap=.true.
          endif
c          if (index(line,'XFILE ') .gt. 0) then
c              n=index(line,'XFILE ')
c              nx=nx+1
c              xfile1(nx)=line(n+6:300)
c              lxfile=.true.
c          endif
          if (index(line,'FACFIL ') .gt. 0) then
              n=index(line,'FACFIL ')
              ffil=line(n+6:300)
              lffile=.true.
          endif
          read(11,5,iostat=eof)line
          goto 10
      endif
      
      write(faclist,'(2(a))')trim(adjustl(ptcmaq)),trim(adjustl(ffil))
c      nxfile=nx
c      allocate(xfile(nxfile))
c      do n=1,nxfile
c          write(xfile(n),'(2(a))')trim(adjustl(xdir)), 
c     +    trim(adjustl(xfile1(n)))
c      enddo
       
c      noutdir=nout
c      allocate(outdir(noutdir))
c      do n=1,noutdir
c          outdir(n)=outfil(n)
c      enddo
      
c     if any keyword not found, stop program
      if (.not. lcmaq) write(logunit,30)
      if (.not. lrec) write(logunit,40)
      if (.not. lxwalk) write(logunit,46)
      if (.not. lhap) write(logunit,47) 
      if(.not. lcmaq .or. .not. lrec .or. .not. lxwalk .or. .not. lhap)
     + then
          write(logunit,50)
          stop
      endif
      
 5    format(a300)
 30   format('PTCMAQ keyword not found')
 40   format('RECDIR keyword not found')
 46   format('XWALK keyword not found')
 47   format('HAPDIR keyword not found')
 50   format('Aborting run..')
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     create national matrix of receptors
      subroutine rec_matrix
      use main1
      implicit none
      
      integer*8 eof,irec,iunit,c,r,gr,br,mr,cr,grec,br1,ii,igrid
      real*8 x,y,barea
      character rfile*300,rid*15,fcode*5
      
      write(logunit,*)'Getting gridded receptors'
      zblks=0 
      iunit=60
c     first read in gridded receptors file
      
c      recdir='/work/ROMO/2014_NATA/receptor_files/'
      if (trim(adjustl(st)) .eq. 'AK') then
          igrid=2
      elseif (trim(adjustl(st)) .eq. 'HI') then
          igrid=3
      elseif (trim(adjustl(st)) .eq. 'PR') then
          igrid=4
      else
          igrid=1
      endif
      
      if (igrid .eq. 1) then
          write(rfile,'(2(a))')trim(adjustl(recdir)),
     +    'gridded_receptors.csv'
      else
          write(rfile,'(3(a))')trim(adjustl(recdir)),
     +    trim(adjustl(gridst(igrid))),'_gridded_receptors.csv'
      endif
      open(unit=iunit,file=rfile,status='old')
      
      eof=0
      gr=0
      read(iunit,*,iostat=eof)rid
      read(iunit,*,iostat=eof)c
 5    if (eof .eq. 0) then
          gr=gr+1
          read(iunit,*,iostat=eof)c
          goto 5
      endif
      rewind(iunit)
      
      ngrids=gr
      allocate(gridconc(ngrids,ncat))
      allocate(gridrec(ngrids))
      gridrec%col=0
      gridrec%row=0
      gridrec%irec=0
      gridrec%recid='NA'
      gridconc=0.0
      
c     now fill in array
      eof=0
      gr=0
      read(iunit,*,iostat=eof)rid
      read(iunit,*,iostat=eof)c,r,irec,x,y
 10   if (eof .eq. 0) then
          cr=c*1000+r
          gr=gr+1
c         version 20002, receptor ID now includes gridst
          write(rid,'(a4,i6.6,i4.4)')
     +    trim(adjustl(gridst(igrid))),cr,irec
          gridrec(gr)%recid=rid
          gridrec(gr)%col=c
          gridrec(gr)%row=r
          gridrec(gr)%irec=gr
          read(iunit,*,iostat=eof)c,r,irec,x,y
          goto 10
      endif      
      close(iunit)
      
c     now blocks
      write(logunit,*)'Getting block receptors'
  
      write(rfile,'(2(a))')trim(adjustl(recdir)),'block_receptors.csv'
      open(unit=iunit,file=rfile,status='old')
      eof=0
      br=0
      read(iunit,*,iostat=eof)rid
      read(iunit,*,iostat=eof)rid
 15   if (eof .eq. 0) then
          br=br+1
          read(iunit,*,iostat=eof)rid
          goto 15
      endif
      rewind(iunit)
      nblocks=br
      br=0
      eof=0

      allocate(blokconc(nblocks,ncat))
      allocate(blockrec(nblocks))

      
      blokconc=0.0
      blockrec%col=0
      blockrec%row=0
      blockrec%irec=0
      blockrec%recid='NA'
      
      read(iunit,*,iostat=eof)rid
      read(iunit,*,iostat=eof)rid,x,y,c,r,barea
 20   if (eof .eq. 0) then
c          cr=c*1000+r
          br=br+1
          blockrec(br)%recid=rid
          blockrec(br)%col=c
          blockrec(br)%row=r
          blockrec(br)%irec=br
          read(iunit,*,iostat=eof)rid,x,y,c,r,barea
          goto 20
      endif
      close(iunit)      

c     now non-populated blocks, only want those that are in non-CONUS
      write(logunit,*)'Getting non-CONUS zero pop block receptors'
  
      write(rfile,'(2(a))')trim(adjustl(recdir)),
     +'zero_pop_block_receptors.csv'
      open(unit=iunit,file=rfile,status='old')
      eof=0
      br=0
c     ak1=0
c     hi1=0
c     pr1=0
c     vi1=0
      br1=0
      read(iunit,*,iostat=eof)rid
      read(iunit,*,iostat=eof)rid
 25   if (eof .eq. 0) then
          br1=br1+1
          if (rid(1:2) .eq. '02' .or. rid(1:2) .eq. '15' .or. 
     +    rid(1:2) .eq. '72' .or. rid(1:2) .eq. '78')then
             if (rid(1:2) .eq. '02') then
                ii=1
             elseif (rid(1:2) .eq. '15') then
                ii=2
             elseif (rid(1:2) .eq. '72') then
                ii=3
             else
                ii=4
             endif
             br=br+1
c            if (rid(1:2) .eq. '02' .and. ak1 .eq. 0) ak1=br 
             if (zblks(ii,1) .eq. 0) zblks(ii,1)=br1
             zblks(ii,2)=zblks(ii,2)+1
c            if (rid(1:2) .eq. '15' .and. hi1 .eq. 0) hi1=br
c            if (rid(1:2) .eq. '72' .and. pr1 .eq. 0) pr1=br
c            if (rid(1:2) .eq. '78' .and. vi1 .eq. 0) vi1=br
          endif
          read(iunit,*,iostat=eof)rid
          goto 25
      endif
      rewind(iunit)
      nzblocks=br
      br=0
      eof=0

      allocate(zblokconc(nzblocks,ncat))
      allocate(zblockrec(nzblocks))

      
      zblokconc=0.0
      zblockrec%col=0
      zblockrec%row=0
      zblockrec%irec=0
c     zblockrec%irec1=0
      zblockrec%recid='NA'
      br1=0
      read(iunit,*,iostat=eof)rid
      read(iunit,*,iostat=eof)rid,x,y,c,r,barea
 30   if (eof .eq. 0) then
          br1=br1+1
          if (rid(1:2) .eq. '02' .or. rid(1:2) .eq. '15' .or. 
     +    rid(1:2) .eq. '72' .or. rid(1:2) .eq. '78') then
              cr=c*1000+r
              br=br+1
              zblockrec(br)%recid=rid
              zblockrec(br)%col=c
              zblockrec(br)%row=r
c             zblockrec(br)%irec1=br
              zblockrec(br)%irec=br1 !master record
          endif
          read(iunit,*,iostat=eof)rid,x,y,c,r,barea
          goto 30
      endif
      close(iunit)  
      
c     write(*,*)zblks(1,1),zblks(1,2)
c     write(*,*)zblks(2,1),zblks(2,2)
c     write(*,*)zblks(3,1),zblks(3,2)
c     write(*,*)zblks(4,1),zblks(4,2)


c     now monitors
      write(logunit,*)'Getting monitor receptors'
      
      write(rfile,'(2(a))')trim(adjustl(recdir)),'monitor_receptors.csv'
      open(unit=iunit,file=rfile,status='old')
      eof=0
      mr=0
      read(iunit,*,iostat=eof)rid
      read(iunit,*,iostat=eof)rid,fcode
 35   if (eof .eq. 0) then
         mr=mr+1
          read(iunit,*,iostat=eof)rid
          goto 35
      endif
      rewind(iunit)
      nmons=mr
      mr=0
      eof=0
      allocate(monconc(nmons,ncat))
      allocate(monrec(nmons))
      
      monconc=0.0
      monrec%col=0
      monrec%row=0
      monrec%irec=0
      monrec%recid='NA'
      
      read(iunit,*,iostat=eof)rid
      read(iunit,*,iostat=eof)rid,fcode,x,y,c,r
 40   if (eof .eq. 0) then
          mr=mr+1
          monrec(mr)%col=c
          monrec(mr)%row=r
          monrec(mr)%irec=mr
          monrec(mr)%recid=rid
          read(iunit,*,iostat=eof)rid,fcode,x,y,c,r
          goto 40
      endif
      close(iunit)      
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     subroutine to loop through file list and calculate totals
      subroutine read_post
      use main1
      implicit none
      integer eof
      integer*8 isrc,nsrc,nsrc1,j,jj
      character fname*300,cas*20,facid*15,zipline*310,ares*2
      logical lexist,lfound,lstop,lexist1
      
      integer*8, allocatable, dimension(:,:) :: igrps
      
      write(logunit,*)'Reading HAP-facility list and concentrations'
      call datetime
      if (.not. runpt .and. .not. runap .and. .not. runcmv)
     +write(ares,'(i2.2)')gridres
c     get a list of files for the HAP
c      write(listfil,'()')trim(adjustl(hap)),
c     +'_list.txt'
c      write(str,'(a,1x,4(a))')'ls',trim(adjustl(hapdir)),
c     +trim(adjustl(hap)),'_PT_*_concs.txt >',trim(adjustl(listfil))
c      call system(copyline)
      inquire(file=faclist,exist=lexist)
      if (.not. lexist) then
          write(logunit,*)'no facility file, stop'
          stop
      endif
      
      eof=0
      open(unit=14,file=faclist,status='old')
      
c     first get # of facilities for the requested HAP
      lfound=.false.
      lstop=.false.
      isrc=0
      do while(eof .eq. 0 .and. .not. lstop)
          read(14,*,iostat=eof)cas,facid
          if (trim(adjustl(cas)) .eq. trim(adjustl(hap))) then
              lfound=.true.
              isrc=isrc+1
          else
              if (lfound) lstop=.true.
          endif
c         read(14,*,iostat=eof)cas,facid
      enddo
      
      nsrc1=isrc
      nsrc=nsrc1/nruns
      
      if (mod(real(nsrc1),real(nruns)) .ne. 0.0) nsrc=nsrc+1
      
      allocate(igrps(nruns,2))
      do j=1,nruns
          if (j .eq. 1) then
              igrps(j,1)=1
          else
              igrps(j,1)=igrps(j-1,2)+1
          endif
          igrps(j,2)=j*nsrc
      enddo
      igrps(nruns,2)=nsrc1
      jj=igrps(runno,2)-igrps(runno,1)+1

c     now read POST files
      rewind(14)
      lfound=.false.
      lstop=.false.
      isrc=0
      eof=0
      do while(eof .eq. 0 .and. .not. lstop)
          read(14,*,iostat=eof)cas,facid
          if (trim(adjustl(cas)) .eq. trim(adjustl(hap))) then
              lfound=.true.
              isrc=isrc+1
              if (isrc .ge. igrps(runno,1) .and. isrc .le. 
     +            igrps(runno,2)) then
c                 first check to see if a zipped version of the file exists
                  if (runpt .or. runap .or. runcmv) then
                    write(fname,'(7(a))')trim(adjustl(hapdir)),
     +              trim(adjustl(hap)),'_',trim(adjustl(runtyp)),'_',
     +              trim(adjustl(facid)),'_1_concs.txt.gz'
                  else
                    write(fname,'(9(a))')trim(adjustl(hapdir)),
     +              trim(adjustl(hap)),'_',trim(adjustl(runtyp)),'_',
     +              trim(adjustl(facid)),'_',trim(adjustl(ares)),
     +              '_1_concs.txt.gz'
                  endif
                  inquire(file=fname,exist=lexist1)
                  if (lexist1) then  
c                     unzip the file
                      write(zipline,'(a,1x,a)')'gunzip -f',
     +                trim(adjustl(fname))
                      call system(zipline)
c                     now rewrite filename without .gz extension
                      if (runpt .or. runap .or. runcmv) then
                        write(fname,'(7(a))')trim(adjustl(hapdir)),
     +                  trim(adjustl(hap)),'_',trim(adjustl(runtyp)),'_'
     +                  ,trim(adjustl(facid)),'_1_concs.txt'
                      else
                        write(fname,'(9(a))')trim(adjustl(hapdir)),
     +                  trim(adjustl(hap)),'_',trim(adjustl(runtyp)),'_'
     +                  ,trim(adjustl(facid)),'_',trim(adjustl(ares)),
     +                  '_1_concs.txt'
                      endif
                      call read_file(fname)
                  else
c                     check for unzipped file
                      if (runpt .or. runap .or. runcmv) then
                        write(fname,'(7(a))')trim(adjustl(hapdir)),
     +                  trim(adjustl(hap)),'_',trim(adjustl(runtyp)),'_' 
     +                  ,trim(adjustl(facid)),'_1_concs.txt'
                      else
                        write(fname,'(9(a))')trim(adjustl(hapdir)),
     +                  trim(adjustl(hap)),'_',trim(adjustl(runtyp)),'_'
     +                  ,trim(adjustl(facid)),'_',trim(adjustl(ares)),
     +                  '_1_concs.txt'
                      endif
                      inquire(file=fname,exist=lexist1)
                      if (lexist1) then
                          call read_file(fname)
                      else
                          write(logunit,5)trim(adjustl(facid)),
     +                    trim(adjustl(cas))
                      endif
                  endif
              endif
          else
              if (lfound) lstop=.true.
          endif
c         read(14,*,iostat=eof)cas,facid
      enddo
      
      if (.not. lfound) then
         lgo=.false. 
       else
         lgo=.true.
       endif
 5    format('File does not exist for ',a,1x,a,' skipping')
      close(14)
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine read_file(fname)
      use main1
      implicit none
      integer*8 eof,itype,ic,c,r,findrec,irec,irec1,ii,jj
      character rid*15,fname*300,line*1000,zipline*310
      logical lexist,lfound

      real*8, allocatable, dimension(:) :: conc

      allocate(conc(ncat))

      write(logunit,*)'Reading ',trim(adjustl(fname))
      call datetime
      eof=0
      open(unit=13,file=fname,status='old')
      read(13,'(a1000)')line
      read(13,*,iostat=eof)rid,irec,itype,c,r,(conc(ic),ic=1,ncat)
 10   if (eof .eq. 0) then
          if (itype .eq. 1) then
              do ic=1,ncat
                  gridconc(irec,ic)=gridconc(irec,ic)+conc(ic)
              enddo
          endif
          if (itype .eq. 2) then
              do ic=1,ncat
                  blokconc(irec,ic)=blokconc(irec,ic)+conc(ic)
              enddo
          endif
          if (itype .eq. 4) then
              if (rid(1:2) .eq. '02') then  
                 ii=1
                 jj=0
              elseif (rid(1:2) .eq. '15') then
                 ii=2
                 jj=zblks(ii-1,2)
              elseif (rid(1:2) .eq. '72') then
                 ii=3 
                 jj=zblks(ii-1,2)+zblks(ii-2,2)
              else
                 ii=4
                 jj=zblks(ii-1,2)+zblks(ii-2,2)+zblks(ii-3,2)
              endif
c             lfound=.false.
c             do while (irec1 .le. nzblocks .and. .not. lfound)
c                if (zblockrec(irec1)%irec .eq. irec) then
c                    lfound=.true.
c                else
c                    irec1=irec1+1
c                endif
c             enddo  
              irec1=irec-(zblks(ii,1)-1)+jj
              do ic=1,ncat
c                 zblokconc(irec,ic)=zblokconc(irec,ic)+conc(ic)
                  zblokconc(irec1,ic)=zblokconc(irec1,ic)+conc(ic)
              enddo
          endif
          if (itype .eq. 3) then
              do ic=1,ncat
                  monconc(irec,ic)=monconc(irec,ic)+conc(ic)
              enddo
          endif
          read(13,*,iostat=eof)rid,irec,itype,c,r,
     +    (conc(ic),ic=1,ncat)
          goto 10
      endif
      close(13)
      
c     gzip file
      write(zipline,'(a,1x,a)')'gzip -f',trim(adjustl(fname))
      call system(zipline)
      deallocate(conc)
      call datetime
      write(logunit,*)'Done with ',trim(adjustl(fname))
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     subroutine to write output files
      subroutine writeout
      use main1
      implicit none
      integer*8 r,c,rt
      character form1*45,acom*1,acom1,fname2*250,zipline*300
      allocatable :: acom(:)

      if (.not. lgo) then
         write(logunit,*)'No files found..'
         goto 95
      endif
      write(logunit,*)'Writing output files..'

      if (nruns .eq. 1) then
          write(fname2,'(6(a))')trim(adjustl(hapdir)),'total/',
     +    trim(adjustl(hap)),'_',trim(adjustl(runtyp)),
     +    '_total_concs_1.txt'
      else
          write(fname2,'(6(a),i3.3,a)')trim(adjustl(hapdir)),'total/',
     +    trim(adjustl(hap)),'_',trim(adjustl(runtyp)),
     +    '_total_concs_1_',runno,'.txt'
      endif
      open(unit=14,file=fname2,status='replace')

c      write(form1,'(a,i3,a)')'(a15,1x,i1,',ncats,'(1x,e13.6))'
c      write(form1,'(a,i3,a)')'(a,a1,i1,a1,',ncats,'(e12.6,a1))'
c      write(form1,'(a,i3,a)')'(a,a1,i1,a1,2(i3.3,a1),',ncats,
c     +'(e12.6,a1))'
      write(form1,'(a,i3,a)')'(a,a1,i7,a1,i1,a1,2(i3.3,a1),',ncat,
     +'(e12.6,a1))'


      allocate(acom(ncat))
      acom1=','
      acom=','
      acom(ncat)=' '  !set last one to blank
      do r=1,ngrids !gridded
          rt=1
          write(14,form1)trim(adjustl(gridrec(r)%recid)),acom1,
     +    gridrec(r)%irec,acom1,rt,acom1,gridrec(r)%col,acom1,
     +    gridrec(r)%row,acom1,(gridconc(r,c),acom(c),c=1,ncat)
      enddo
      do r=1,nblocks !blocks
          rt=2
          write(14,form1)trim(adjustl(blockrec(r)%recid)),
     +    acom1,blockrec(r)%irec,acom1,rt,acom1,blockrec(r)%col,
     +    acom1,blockrec(r)%row,acom1,(blokconc(r,c),
     +    acom(c),c=1,ncat)
      enddo
      do r=1,nzblocks !blocks
          rt=4
          write(14,form1)trim(adjustl(zblockrec(r)%recid)),
     +    acom1,zblockrec(r)%irec,acom1,rt,acom1,zblockrec(r)%col,
     +    acom1,zblockrec(r)%row,acom1,(zblokconc(r,c),
     +    acom(c),c=1,ncat)
      enddo
      do r=1,nmons !monitors
          rt=3
          write(14,form1)trim(adjustl(monrec(r)%recid)),acom1,
     +    monrec(r)%irec,acom1,rt,acom1,monrec(r)%col,acom1,
     +    monrec(r)%row,acom1,(monconc(r,c),acom(c),c=1,ncat)
      enddo

      close(14)
c     zip the file
      write(zipline,'(a,1x,a)')'gzip -f',trim(adjustl(fname2))
      call system(zipline)  !zip file
      write(logunit,*)'Output files written'
      call datetime
 95   if (allocated(gridconc)) deallocate(gridconc)
      if (allocated(gridrec)) deallocate(gridrec)
      if (allocated(blokconc)) deallocate(blokconc)
      if (allocated(blockrec)) deallocate(blockrec)
      if (allocated(zblokconc)) deallocate(zblokconc)
      if (allocated(zblockrec)) deallocate(zblockrec)
      if (allocated(monconc)) deallocate(monconc)
      if (allocated(monrec)) deallocate(monrec)
      
     
      close(logunit)
     
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine upcase(flag)
c     convert lowercase variable to upper case variable

c flag:       variable to checked
c lower:      string of lowercase letters
c upper:      string of uppercase letters
c i:          index of flag in lower
      integer*8 i
      character flag*1
      character lower*26,upper*26
      i=0
      lower='abcdefghijklmnopqrstuvwxyz'
      upper='ABCDEFGHIJKLMNOPQRSTUVWXYZ'
      
      i=index(lower,flag)

      if (i .gt. 0) flag=upper(i:i)
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     subroutine to make sure grid resolution from input argument is valid
      subroutine checkgrid(agrid)
      use main1
      implicit none
      integer i,ii,ilen
      character agrid*2
      logical lbad
      
      gridres=0
      
     
      ilen=len_trim(agrid)
      i=1
      lbad=.false.
      do while (i .le. ilen .and. .not. lbad)
          ii=ichar(agrid(i:i))
          if (ii .lt. 48 .or. ii .gt. 57) lbad=.true.
          i=i+1
      enddo
      if (.not. lbad) read(agrid,*)gridres
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine datetime
c subroutine to write starting and ending times to screen and log file
c subroutine calls the intrinsic FORTRAN subroutine date_and_time
c inputs are iflag which tells the subroutine if the program starts
c(iflag=1)
c or ends (iflag not equal to 1)

c variables
c idattim:  array of date/time returned by date_and_time with 8 elements
c           elements are:
c           (1) 4-digit year
c           (2) month of year (1 to 12)
c           (3) day of month
c           (4) time difference with respect to GMT in minutes
c           (5) hour of day (0 to 23) in local time
c           (6) minute of time
c           (7) second of time
c           (8) milliseconds of time
c hr:       integer hour of day
c cdate:    character variable of date
c ctime:    character variable of time
c czone:    character variable of time zone from GMT
c ampm:     character variable with value of 'AM' or 'PM'     
    
      use main1
      implicit none
      integer*8 idattim(8),hr
      character cdate*8,ctime*10,czone*5,ampm*2,amonths(12)*10

      data amonths /'January','February','March','April','May','June',
     &  'July','August','September','October','November','December'/
      
      call date_and_time(cdate,ctime,czone,idattim)
      
      if (idattim(5) .le. 11) then
        hr=idattim(5)
        ampm='AM'
      else
        hr=12+(idattim(5)-12)
        ampm='PM'
      endif
      write(logunit,10)trim(adjustl(amonths(idattim(2)))),idattim(3),
     + idattim(1),hr,idattim(6),idattim(7),ampm

 10   format(1x,a,1x,i2.2,', ',i4,2x,i2,':',
     +  i2.2,':',i2.2,1x,a2)
      return
      end