      module main1
c     character variables
c     runtyp:     AERMOD runtype see subroutine readargs for values
c     facid:      facility ID
c     state2:     2-character state abbreviation
c     srcnam:     source name used for AERMOD.INP file title
c     ptloc:      directory containing AERMOD include location files
c     ptparam:    directory containing AERMOD include source parameter files
c     ptcmaq:     
c     pttemp:     directory containing AERMOD include temporal factor files
c                 or hourly emission files
c     metdir:     directory containing AERMOD meteorology files
c     recdir:     directory containing receptor files.  these are not
c                 in AERMOD format
c     outdir:     directory for AERMOD output
c     locsrc:     array of source IDs from location files
c     parsc:      array of source IDs from parameter files
c     temsrc:     array of source IDs from temporal files
c     locfil:     facility specific AERMOD source locations file
c     paramfil:   facility specific AERMOD source parameter file(s)
c                 variable is 2D.  1 dim is point sources filename
c                 2nd dim is area sources filename
c                 if no point or no area sources, then paramfil
c                 is set to NA
c     temsrc:     facility specific AERMOD temporal factor file or
c                 hourly emissions file
c     singlhap:   specific hap to process    
c     gridst:     grids
      character runtyp*7,facid*15,state2*2,ptloc*250,ptparam*250,
     +pttemp*250,ptcmaq*250,metdir*250,recdir*250,outdir*250,
     +hapdir*250,haps*15,xdir*250,rtype*15,cats*100,singlhap*15,
     +haplist*15,hap1*15,fac1*15,st2*2,singlfac*15,singlst*2,xfile*250,
     +states(54)*2,states1(54)*2,versn*5,logfil*50,ncst*2,gridst(4)*4!,grecid*15,
c     +monid*15,blockid*15,mfips*5

c      allocatable :: locsrc(:)
c      allocatable :: parsrc(:)
c      allocatable :: temsrc(:)
      allocatable :: haps(:)
      allocatable :: haplist(:)
      allocatable :: facid(:)
      allocatable :: state2(:)
      allocatable :: rtype(:)
      allocatable :: cats(:)
      allocatable :: hap1(:)
      allocatable :: fac1(:)
      allocatable :: st2(:)
      allocatable :: outdir(:)
      allocatable :: xfile(:)
c      allocatable :: grecid(:)
c      allocatable :: monid(:)
c     allocatable :: blockid(:)
c    allocatable :: mfips(:)
      
c     lparam:     2-D array denoting parameter files found
c                 1 dim is point parameter file
c                 2 dim is area parameter file for point sources
c     ltemp:      3-D array denoting temporal files found
c                 1 dim is emissions factor file
c                 2 dim is hourly emissions file
c                 3 dim is both have been found
c     lcnty:      logical variable denoting county for grid cell has been
c                 found in RWC CMAQ to county cross reference
c     runpt:      point source being processed
c     runap:      airport source being processed
c     runon:      onroad sources being processed
c     runnon:     nonroad sources being processed
c     runrwc:     RWC sources being processed
c     runcmv:     CMV (port or underway) sources being processed
c     runnphi:    nonpoint 10 m sources being processed
c     runnplow:   low-level nonpoint sources being processed
c     akhiprvi:   variable denoting that AK, HI, PR, or VI being processed
c     lgrids:     denotes which grid files to read      
      logical runpt,runap,runon,runnon,runrwc,runcmv,runnphi,runnplow,
     +runog,runag,akhiprvi,lsinghap,writefac,noncontr,conus(2),lsingfac,
     +lsingst,writetot,lgrids(5),readzero,writheid  !write file for HEID purposes
      
      logical, allocatable, dimension(:) :: lhaps  !when processing a set number of HAPs, to see if HAP is in rungroup
c     d2r:        degrees to radians conversion factor
c     blkdist:    distance to search from source for block receptors
c     grdddist:   distance to search from source for gridded receptors
c     profbase:       profbase written to AERMOD.INP
      real*8 d2r,blkdist,grddist,profbase
      
      real*8, allocatable, dimension(:,:,:) :: gridconc  !# of categories, # of receptors
      real*8, allocatable, dimension(:,:,:) :: monconc  !# of categories, # of receptors
      real*8, allocatable, dimension(:,:,:) :: blokconc  !# of categories, # of receptors
      real*8, allocatable, dimension(:,:,:) :: zblokconc  !# of categories, # of receptors for nonpopulated blocks in non-CONUS version 17236
      
      real*8, allocatable, dimension(:,:,:) :: modgc !modeled gridded concentrations
      real*8, allocatable, dimension(:,:,:) :: modbc !modeled block concentrations
      real*8, allocatable, dimension(:,:,:) :: modzbc !modeled nonpopulated non-CONUS block concentrations version 17236
      real*8, allocatable, dimension(:,:,:) :: modmc !modeled monitor concentrations
      real*8, allocatable, dimension(:,:,:) :: intbc !modeled block concentrations
      real*8, allocatable, dimension(:,:,:) :: intmc !modeled monitor concentrations
      
      real*8, allocatable, dimension(:,:) :: polyx !cmvx
      real*8, allocatable, dimension(:,:) :: polyy !cmvy
      
      real*8, allocatable, dimension(:,:) :: memis !seasonal emissions by group
      
c      real*8, allocatable, dimension(:) :: gridx  !cmaq grid coordinates
c      real*8, allocatable, dimension(:) :: gridy  !cmaq grid coordinates
c      real*8, allocatable, dimension(:) :: monx  !cmaq grid coordinates
c      real*8, allocatable, dimension(:) :: mony  !cmaq grid coordinates
c      real*8, allocatable, dimension(:) :: blockx  !cmaq grid coordinates
c      real*8, allocatable, dimension(:) :: blocky  !cmaq grid coordinates
c     col:        CMAQ column
c     row:        CMAQ row
c     inpunit:    input_directories.txt file unit, lists directories (11)
c     locunit:    locations file unit (12)
c     parunit:    parameter file unit (13)
c     temunit:    temporal file unit (14)
c     recunit:    receptor file unit (16)
c     aerunit:    AERMOD.INP file unit (27)
c     metunit:    file of AK and HI met stations (32)
c     ngrps:      # of source groups, usually equal to number of sources
c     utmzone:    UTM zone of source
c     hrunit:     hourly file unit (28)
c     gridres:    grid resolution in km
c     cntyunit:   county file unit (29)
c     fips:       County FIPS code
c     wrfunit:    WRF elevations file unit (31)
c     days:       array of number of days/month
c     ncats:      number of categories output by source group
      integer*8 col,row,inpunit,recunit,days(12),modyear,
     +ngrps,utmzone,gridres,fips,nlines,nhaps,nsrc,ngrids,nhours(12),
     +nblocks,nmons,nmodg,nmodb,nmodm,nmodz,nzblocks,nintm,
     +nintb,pltunit,ncats,nmonths,ngrps1,nzblocks1,seashrs(4,24),
     +hapunit,intunit,ngrids1,nblocks1,nmons1,noutdir,nxfile,
     +nruns,runno,nhap1,nhap2,ncolrow,istate,errunit,logunit
    
      integer*8, allocatable, dimension(:) :: srccr !array of source column/rows used for subsetting receptor arrays
      integer*8, allocatable, dimension(:) :: srccr1 !array of source column/rows used for subsetting receptor arrays,includes duplicates
      integer*8, allocatable, dimension(:) :: stgrid !array of grid identifiers for srccr
      integer*8, allocatable, dimension(:) :: stgrid1 !array of grid identifiers for srccr1
      integer*8, allocatable, dimension(:) :: nsrcs !number of sources at facility
      
c      integer, allocatable, dimension(:) :: blockcr 
c      integer, allocatable, dimension(:) :: moncr
      
c     AERMOD point source inputs
      type facinfo1
          character srcid*12,srctyp*8
          real*8 utmx,utmy,xdim,ydim,angle
      end type facinfo1
c     AERMOD point source emission info  
      type facinfo
c          character cat*100,facid*12,cas*15,state*2,srcgrp*8
          character cat*100,facid*15,cas*15,state*2,srcgrp*12
          integer*8 factype,procid,unitid,relid  !process id, unit id, release id
          real*8 emis(1)
      end type facinfo

c     AERMOD CMV source inputs
      type cmvinfo1
          character srcid*12,srctyp*8
          real*8 utmx1,utmy1
          integer*8 nvert
      end type cmvinfo1
c     AERMOD CMV source emission info  
      type cmvinfo
          character cat*100,facid*15,cas*15,state*2,srcgrp*12
          real*8 emis(1)
      end type cmvinfo
      
c     AERMOD tract source inputs
      type tractinfo1
          character srcid*12,srctyp*8
          real*8 utmx1,utmy1
          integer*8 nvert
      end type tractinfo1
      
c     AERMOD tract source emission info  
      type tractinfo
          character cat*100,facid*15,cas*15,state*2,srcgrp*12
          real*8 emis(13)
      end type tractinfo
      
c     AERMOD airport inputs 
      type apinfo1 
          character srcid*12,srctyp*8
          real*8 area,utmx1,utmx2,utmy1,utmy2,frac,xdim,ydim,angle
      end type apinfo1
      
c     airport facility emissions info      
      type apinfo 
c          character facid*12,cas*15,state*2,cat*100,srcgrp*8
          character facid*15,cas*15,state*2,cat*100,srcgrp*12
          real*8 emis(1)
      end type apinfo
   
      type gridinfo1
          character srcid*12
          real*8 utmx1,utmx2,utmx3,utmx4,utmy1,utmy2,utmy3,utmy4
          integer*8 nvert
      end type gridinfo1
      
c     onroad, nonroad
      type gridinfo
c          character facid*12,cas*15,fips*5,cat*100,srcgrp*8
          character facid*15,cas*15,fips*5,cat*100,srcgrp*12
          real*8 emis(13) !monthly emissions and annual. annual may not be used
      end type gridinfo
      
c     all gridded receptors
      type gridrecs
          character recid*15,rtype*15
          real*8 x,y  !x,y are cmaqx,cmaqy
          integer*8 colrow,iflag,col,row,irec,irec1
      end type gridrecs

      type gridrecs1
          character recid*15
          integer*8 irec
      end type gridrecs1
c     modeled gridded receptors
      type modgrids
          character recid*15
          real*8 ux,uy,x,y,dist !ux,uy are utm, x,y are cmaqx,cmaqy
          integer*8 colrow,iflag,mastrec   !mastrec is the index of the receptor in national gridded network
      end type modgrids
      
c     all block receptors
      type blkrecs
          character recid*15,rtype*15
          real*8 x,y !x,y are cmaqx,cmaqy
          integer*8 colrow,iflag,col,row,irec,irec1
      end type blkrecs
      type blkrecs1
          character recid*15
          integer*8 irec
      end type blkrecs1
c     modeled block receptors
      type modblks
          character recid*15
          real*8 ux,uy,x,y,dist !ux,uy are utm, x,y are cmaqx,cmaqy
          integer*8 colrow,iflag,mastrec   !mastrec is the index of the receptor in national block network
      end type modblks
c     all nonpopulated non-CONUS block receptors
      type zblkrecs
          character recid*15,rtype*15
          real*8 x,y !x,y are cmaqx,cmaqy
          integer*8 colrow,iflag,col,row,irec,irec1
      end type zblkrecs
      type zblkrecs1
          character recid*15
          integer*8 irec
      end type zblkrecs1
c     modeled nonpopulated non-CONUSblock receptors
      type zmodblks
          character recid*15
          real*8 ux,uy,x,y,dist !ux,uy are utm, x,y are cmaqx,cmaqy
          integer*8 colrow,iflag,mastrec   !mastrec is the index of the receptor in national block network
      end type zmodblks
c     interpolated blocks
      type intblks
          character recid*15
          real*8 ux,uy,x,y,dist !ux,uy are utm, x,y are cmaqx,cmaqy
          integer*8 colrow,mastrec   !mastrec is the index of the receptor in national block network
      end type intblks
      
c     all monitors
      type monrecs
          character recid*15,fip*5,rtype*5
          real*8 x,y !x,y are cmaqx,cmaqy
          integer*8 colrow,iflag,col,row,irec,irec1
      end type monrecs
      type monrecs1
          character recid*15
          integer*8 irec
      end type monrecs1
c     modeled monitor receptors
      type modmons
          character recid*15
          real*8 ux,uy,x,y,dist !ux,uy are utm, x,y are cmaqx,cmaqy
          integer*8 colrow,iflag,mastrec   !mastrec is the index of the receptor in national monitor network
      end type modmons

c     interpolated monitors
      type intmons
          character recid*15
          real*8 ux,uy,x,y,conc,dist !ux,uy are utm, x,y are cmaqx,cmaqy
          integer*8 colrow,mastrec   !mastrec is the index of the receptor in national monitor network
      end type intmons
      
      type(facinfo1), dimension(:),allocatable:: ptsrcs
      type(facinfo), dimension(:),allocatable:: facility
      type(cmvinfo1), dimension(:),allocatable:: cmvsrcs
      type(cmvinfo), dimension(:),allocatable:: cmv
      type(tractinfo1), dimension(:),allocatable:: tractsrcs
      type(tractinfo), dimension(:),allocatable:: tracts
      type(apinfo), dimension(:),allocatable:: airport
      type(gridinfo), dimension(:),allocatable:: gridded
c     type(mobinfo), dimension(:),allocatable:: mobile
      type(apinfo1), dimension(:),allocatable:: apsrcs
      type(gridinfo1), dimension(:),allocatable:: gridsrcs
      type(gridrecs), dimension(:),allocatable:: gridrec
      type(gridrecs1), dimension(:),allocatable:: gridrec1
      type(modgrids), dimension(:),allocatable:: modgrid
      type(blkrecs), dimension(:),allocatable:: blockrec
      type(blkrecs1), dimension(:),allocatable:: blockrec1
      type(modblks), dimension(:),allocatable:: modblock
      type(intblks), dimension(:),allocatable:: intblock
      type(zblkrecs), dimension(:),allocatable:: zblockrec
      type(zblkrecs1), dimension(:),allocatable:: zblockrec1
      type(zmodblks), dimension(:),allocatable:: zmodblock
      type(monrecs), dimension(:),allocatable:: monrec
      type(monrecs1), dimension(:),allocatable:: monrec1
      type(modmons), dimension(:),allocatable:: modmon
      type(intmons), dimension(:),allocatable:: intmon
c     parameters
      parameter(inpunit=11,recunit=16,xwalk=38,pltunit=39,
     +hapunit=41,intunit=26,errunit=27,d2r=3.141592654/180.0,
     +grddist=50000.0,blkdist=10000.0,logunit=37)
c     assign number of days per month.  February set to 28
c                 J  F  M  A  M  J  J  A S  O  N  D
      data days /31,28,31,30,31,30,31,31,30,31,30,31/  
      data versn /'20016'/
      data states /'AL','AK','AZ','AR','CA','CO','CT','DE','DC','FL',
     +'GA','HI','ID','IL','IN','IA','KS','KY','LA','ME','MD','MA','MI',
     +'MN','MS','MO','MT','NE','NV','NH','NJ','NM','NY','NC','ND','OH',
     +'OK','OR','PA','RI','SC','SD','TN','TX','UT','VT','VA','WA','WV',
     +'WI','WY','PR','VI','TB'/
      data states1 /'01','02','04','05','06','08','09','10','11','12',
     +'13','15','16','17','18','19','20','21','22','23','24','25','26',
     +'27','28','29','30','31','32','33','34','35','36','37','38','39',
     +'40','41','42','44','45','46','47','48','59','50','51','53','54',
     +'55','56','72','78','88'/
      data gridst /'12US','09AK','03HI','03PR'/
      end module main1
c##########################################################################
      
      program postaermod
      use main1
      implicit none
c     read command line arguments to get run type
      call readargs      
c     get directories
      call getdir
c     get list of HAPs to process
      call gethaps
c     read crosswalk
      lgrids=.false.
      if (runap) then
          call read_air
      elseif (runpt) then
          call read_point
      elseif (runcmv) then
          call read_cmv
c      elseif (noncontr:q:q) then
c          call read_tract
      else
          call read_grid
      endif
      
c     determine if source array has mix of CONUS and non-CONUS sources and get column/rows of sources
      call src_cr
c     get grid cells that are covered by
c     sources being processed to help
c     cut down number of receptors in master array
c     call haps_sources  !now called in read_ subroutines
c     create matrix of all blocks, gridded and monitors
      call rec_matrix
      
c     read in modeled sources and receptors
      call read_model
      
c     write output
      if (writetot) call writeout
      close(logunit)
      end

c##########################################################################
c     subroutine to read command line arguments
      subroutine readargs
      use main1
      implicit none
      
c     i:          loop counter
c     ilen:       length of argument
c     agrid:      character string of grid resolution (km)
c     awrite:     character string to denonte write facility specific output (yes or no)
      
      integer*8 i,ilen,istat,ndash,n1,n2,get_state
      character agrid*2,awrite*3,aruns*10,arun*10,noncon*1,ahap1*3,
     +ahap2*3,runtyp1*12,awrite1*3,awrite2*3

c      write(*,'(2(a))')'Starting program version: ',versn
c      call datetime
     
c     initialize logical flags for run types
      runpt=.false.
      runap=.false.
      runcmv=.false.
      runon=.false.
      runnon=.false.
      runrwc=.false.
      runnphi=.false.
      runnplow=.false.
      runog=.false.
      runag=.false.
c     initialize variable denoting if AK, HI, PR, or VI
c     being processed
      akhiprvi=.false.
      noncontr=.false.
      
c     get the requested run type 
c      call getarg(1,runtyp)
      call getarg(1,runtyp1)
c     these are subject to change
c     convert to upper case
c      ilen=len_trim(runtyp)
      ilen=len_trim(runtyp1)
      ndash=0
      do i=1,ilen
c          call upcase(runtyp(i:i))
          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) 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
c         allow first letter of state to be equivalent of  yes
          if (noncon .eq. 'Y' .or. noncon .eq. 'A' .or. noncon .eq. 'H'
     +    .or .noncon .eq. 'P' .or. noncon .eq. 'V') then  
              noncontr=.true.
              if (noncon .eq. 'A') then
                  ncst='AK'
              elseif (noncon .eq. 'H') then
                  ncst='HI'
              elseif (noncon .eq. 'P') then
                  ncst='PR'
              elseif (noncon .eq. 'VI') then
                  ncst='VI'
               else
                   ncst='NA'
                endif
          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.
          nmonths=1
      elseif (runtyp .eq. 'AP') then
          runap=.true.
          nmonths=1
      elseif (runtyp .eq. 'PORT' .or. runtyp .eq. 'UNDER' .or. runtyp 
     +.eq. 'CMV') then
          runcmv=.true.
          nmonths=1
      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.
          nmonths=12
      elseif (runtyp .eq. 'NONRD') then
          runnon=.true.
          nmonths=12
c          if (runnon .and. noncontr) nmonths=1 !reset nonroad to 1 season for tract sources
      elseif (runtyp .eq. 'NPHI') then
          runnphi=.true.
          nmonths=1
      elseif (runtyp .eq. 'NPLO') then
          runnplow=.true.
          nmonths=1
      elseif (runtyp .eq. 'RWC') then
          runrwc=.true.
          nmonths=1
      elseif (runtyp .eq. 'OILGAS') then
          runog=.true.
          nmonths=1
      elseif (runtyp .eq. 'AG') then
          runag=.true.
          nmonths=1
      else
          write(*,5)runtyp
          stop
      endif
      
c     single or all haps
      call getarg(2,singlhap,istat)
      ilen=len_trim(singlhap)
      do i=1,ilen
          call upcase(singlhap(i:i))
      enddo
      if (singlhap .eq. 'ALL') then
          lsinghap=.false.
      else
          lsinghap=.true.
          if (trim(adjustl(singlhap)) .eq. 'BENZN') singlhap='BENZENE' !make consistent across run groups
      endif
      
c     HAP #s
      call getarg(3,ahap1,istat)
      if (lsinghap) then
          nhap1=-1
          nhap2=-1
      else
          read(ahap1,*)nhap1
          call getarg(4,ahap2)
          read(ahap2,*)nhap2
      endif

c     single facility or state
      call getarg(5,singlfac,istat)
      if (istat .eq. -1) then
          lsingfac=.false.
          lsingst=.false.
      else
          ilen=len_trim(singlfac)
          do i=1,ilen
              call upcase(singlfac(i:i))
          enddo
          if (ilen .eq. 2) then !state
              lsingst=.true.
              singlst=trim(adjustl(singlfac))
              istate=get_state(singlst)  !get index of state
              lsingfac=.false.
          else
              lsingst=.false.
              if (singlfac .eq. 'ALL') then
                  lsingfac=.false.
              else
                  lsingfac=.true.
              endif
          endif
      endif
c     write facility level output
      call getarg(6,awrite,istat)
      if (istat .eq. -1) then
          writefac=.false.
      else
          ilen=len_trim(awrite)
          do i=1,ilen
              call upcase(awrite(i:i))
          enddo
          if (awrite .eq. 'YES') then
              writefac=.true.
          else
              writefac=.false.
          endif
      endif
      
c     # of runs
      call getarg(7,aruns,istat)
      if (istat .eq. -1 .or. lsingfac) then
          nruns=1
      else
          read(aruns,*)nruns
      endif
      
c     run number
      call getarg(8,arun,istat)
      if (istat .eq. -1 .or. lsingfac) then
          runno=1
      else
          read(arun,*)runno
      endif
          
c     write total output across all processed facilities
c     may not want to do when just one facility
      call getarg(9,awrite1,istat)
      if (istat .eq. -1) then
          writetot=.false.
      else
          ilen=len_trim(awrite1)
          do i=1,ilen
              call upcase(awrite1(i:i))
          enddo
          if (awrite1 .eq. 'YES') then
              writetot=.true.
          else
              writetot=.false.
          endif
      endif
      
c     write facility level data in format for HEID purpsoes
      call getarg(10,awrite2,istat)
      if (istat .eq. -1) then
          writheid=.false.
      else
          ilen=len_trim(awrite2)
          do i=1,ilen
              call upcase(awrite2(i:i))
          enddo
          if (awrite2 .eq. 'YES') then
              writheid=.true.
          else
              writheid=.false.
          endif
      endif
      
      if (runno .gt. nruns) then
          write(*,10)runno,nruns
          stop
      endif
      if ((nhap1 .gt. 0 .and. nhap1 .gt. nhap2) .or. nhap1 .eq. 0) then
          write(*,*)'Number of HAP range invalid'
          stop
      endif
      if (lsingst .and. istate .eq. 0) then
          write(*,15)trim(adjustl(singlst))
          stop
      endif
      
      if (.not. lsingst) istate=1  !set to 1 for later logic
      
c     write filename for logfile based on arguments
      if (lsingfac) then
          if (lsinghap) then
              write(logfil,'(4(a))')trim(adjustl(singlfac)),'_',
     +        trim(adjustl(singlhap)),'.log'
          else
              write(logfil,'(2(a),i3.3,a,i3.3,a)')
     +        trim(adjustl(singlfac)),'_',nhap1,'_',nhap2,'.log'
          endif
      elseif (lsingst) then
          if (lsinghap) then
              write(logfil,'(4(a))')trim(adjustl(singlst)),'_',
     +        trim(adjustl(singlhap)),'.log'
          else
              write(logfil,'(2(a),i3.3,a,i3.3,a)')
     +        trim(adjustl(singlst)),'_',nhap1,'_',nhap2,'.log'
          endif
      else
          if (lsinghap) then
              write(logfil,'(4(a))')trim(adjustl(runtyp)),'_',
     +        trim(adjustl(singlhap)),'.log'
          else
              write(logfil,'(2(a),i3.3,a,i3.3,a)')
     +        trim(adjustl(runtyp)),'_',nhap1,'_',nhap2,'.log'
          endif
      endif
      open(unit=logunit,file=logfil,status='replace')
      write(logunit,'(2(a))')'Starting program version: ',versn
      call datetime
      
 5    format(a,'is not valid, stopping program')
 10   format('Run number: ',i10/'is greater than number of runs: ',i10/
     +'stopping program')
 15   format('State identifier ',a,' is not valid, stopping program')
      return
      end      
c##########################################################################
c     subroutine to get pathnames from input_directories.dat
c     edit on sol once moved over
      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,nx

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,ares*2,infil*250,outfil(2)*250,xfile1(4)*250

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 lloc,lparam2(2),lcmaq,lmet,lrec,lparam1,
     +ltemp1,lout,lexist,lxwalk,lhap,lxfile
      
c     initialize variables
      eof=0
      lloc=.false.
c      lparam=.false.
c      ltemp=.false.
      lcmaq=.false.
      lmet=.false.
      lrec=.false.
      lparam1=.false.
      ltemp1=.false.
      lout=.false.
c      lcnty=.false.
      lxwalk=.false.
      lhap=.false.
      outfil='NA'
      nout=0
      lxfile=.false.
      xfile1='NA'
      nx=0
      nxfile=0
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_',ncst,'.dat'
          endif
      endif
      
      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,'PTLOC ') .gt. 0) then
              n=index(line,'PTLOC ')
              ptloc=line(n+6:300)
              lloc=.true.
          endif
          if (index(line,'PTPARAM ') .gt. 0) then
              n=index(line,'PTPARAM ')
              ptparam=line(n+8:300)
              lparam1=.true.
          endif
          if (index(line,'PTTEMP ') .gt. 0) then
              n=index(line,'PTTEMP ')
              pttemp=line(n+7:300)
              ltemp1=.true.
          endif
          if (index(line,'PTCMAQ ') .gt. 0) then
              n=index(line,'PTCMAQ ')
              ptcmaq=line(n+7:300)
              lcmaq=.true.
          endif
          if (index(line,'METDIR ') .gt. 0) then
              n=index(line,'METDIR ')
              metdir=line(n+7:300)
              lmet=.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
          if (index(line,'OUTPUT ') .gt. 0) then
              n=index(line,'OUTPUT ')
              nout=nout+1
              outfil(nout)=line(n+7:300)
              lout=.true.
          endif
          if (index(line,'XFILE ') .gt. 0) then
              n=index(line,'XFILE ')
              nx=nx+1
              xfile1(nx)=line(n+6:300)
              lxfile=.true.
          endif
          read(11,5,iostat=eof)line
          goto 10
      endif
      
      nxfile=nx
      allocate(xfile(nxfile))
      do n=1,nxfile
          write(xfile(n),'(2(a))')trim(adjustl(xdir)), 
     +    trim(adjustl(xfile1(n)))
      enddo
       
      noutdir=nout
      allocate(outdir(noutdir))
      do n=1,noutdir
          outdir(n)=outfil(n)
      enddo
      
c     if any keyword not found, stop program
      if (.not. lloc) write(logunit,15)
      if (.not. lparam1) write(logunit,20)
      if (.not. ltemp1) write(logunit,25)
      if (.not. lcmaq) write(logunit,30)
      if (.not. lmet) write(logunit,35)
      if (.not. lrec) write(logunit,40)
      if (.not. lout) write(logunit,45)
      if (.not. lxwalk) write(logunit,46)
      if (.not. lhap) write(logunit,47)
      if (.not. lloc .or. .not. lparam1 .or. .not. ltemp1 .or. 
     +.not. lcmaq .or. .not. lmet .or. .not. lrec .or. .not. lxwalk 
     +.or. .not. lout .or. .not. lhap) then
          write(logunit,50)
          stop
      endif

      

      
 5    format(a300)
 15   format('PTLOC keyword not found')
 20   format('PTPARAM keyword not found')
 25   format('PTEMP keyword not found')
 30   format('PTCMAQ keyword not found')
 35   format('METDIR keyword not found')
 40   format('RECDIR keyword not found')
 45   format('OUTDIR keyword not found')
 46   format('XWALK keyword not found')
 47   format('HAPDIR keyword not found')
 50   format('Aborting run..')
 55   format('location file not found for')
 60   format('parameter file(s) not found for')
 65   format('temporal file(s) not found for')   
 70   format('County file not found for')
      return
      end
c##########################################################################
c     subroutine to get haps
      subroutine gethaps
      use main1
      implicit none
      integer eof,i,j
      character ahap*15,hapfile*250
      logical lexist
      
      if (nhap1 .gt. 0) then
          write(hapfile,'(2(a))')trim(adjustl(ptcmaq)),
     +    'hap_priority_list.csv'
          inquire(file=hapfile,exist=lexist)
          if (lexist) then
              nhaps=nhap2-nhap1+1
              allocate(haps(nhaps))
              allocate(lhaps(nhaps))
              open(unit=hapunit,file=hapfile,status='old')
              eof=0
              i=0
              j=0
              read(hapunit,*,iostat=eof)ahap
              read(hapunit,*,iostat=eof)ahap
   5          if (eof .eq. 0) then
                  if (ahap .eq. 'BENZN') ahap='BENZENE'
                  i=i+1
                  if (i .ge. nhap1 .and. i .le. nhap2) then
                    j=j+1
                    haps(j)=trim(adjustl(ahap))
                    lhaps(j)=.false.
                  endif 
                  read(hapunit,*,iostat=eof)ahap
                  goto 5
              endif
              close(hapunit)
          else
              write(logunit,*)'HAP list does not exist'
              write(logunit,*)'use single hap or all haps'
          endif
      else !single hap or all HAPs
          if (lsinghap) then
              nhaps=1
              allocate(haps(nhaps))
              allocate(lhaps(nhaps))
              haps(1)=trim(adjustl(singlhap))
              lhaps(1)=.false.
          else
              write(hapfile,'(2(a))')trim(adjustl(ptcmaq)),
     +        'hap_priority_list.csv'
              inquire(file=hapfile,exist=lexist)
              if (lexist) then
                  open(unit=hapunit,file=hapfile,status='old')
                  eof=0
                  i=0
                  read(hapunit,*,iostat=eof)ahap
                  read(hapunit,*,iostat=eof)ahap
   10            if (eof .eq. 0) then
                      if (ahap .eq. 'BENZN') ahap='BENZENE'
                      i=i+1
                      read(hapunit,*,iostat=eof)ahap
                      goto 10
                  endif
                  rewind(hapunit)
                  nhaps=i
                  allocate(haps(nhaps))
                  allocate(lhaps(nhaps))
                  eof=0
                  i=0
                  read(hapunit,*,iostat=eof)ahap
                  read(hapunit,*,iostat=eof)ahap
   15             if (eof .eq. 0) then
                      if (ahap .eq. 'BENZN') ahap='BENZENE'
                      i=i+1
                      haps(i)=trim(adjustl(ahap))
                      lhaps(i)=.false.
                      read(hapunit,*,iostat=eof)ahap
                      goto 15
                  endif
                  close(hapunit)
              else
                  write(logunit,*)'HAP list does not exist'
              endif
          endif
      endif
      return
      end   
c##########################################################################
c     subroutine to read airport pollutant/category/AERMOD source cross-walk
c     for airports
      subroutine read_air
      use main1
      implicit none
      
      integer*8 eof,ilines,i,ilen,ihap,findhap,nlines1,findfac,is,ityp
      character facid2*15,relid*8,scc1*8,cas1*15,src*12,
     +junk*1,st2a*2,junk2*20
      real*8 ann
      logical lexist,lgo

      write(logunit,*)'Reading emissions file'

     
     
      
c     check for file existence
      inquire(file=xfile(1),exist=lexist)
      if (.not. lexist) then
          write(logunit,5)trim(adjustl(xfile(1)))
          stop
      endif

c     read in cross walk and assign variables to appropriate arrays
      lgo=.false.
      open(unit=xwalk,file=xfile,status='old')
      ilines=0
      read(xwalk,*)junk
      ncats=1
      read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann
 10   if (eof .eq. 0) then
          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          if (ihap .gt. 0) then
              if ((lsingst .and. (trim(adjustl(st2a)) .eq.  
     +        trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) 
     +        .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. 
     +        trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. 
     +        (.not. lsingst .and. .not. lsingfac)) then
                      ilines=ilines+1
                      lhaps(ihap)=.true.
                      lgo=.true.
              endif
          endif
          read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann
          goto 10
      endif
      
      if (.not. lgo) then
          close(xwalk)
          write(logunit,*)'Requested not HAPs found'
          stop
      endif
      
      rewind(xwalk)
c      nlines=ilines
      nlines1=ilines
  
c      allocate(hap1(nlines1))
      allocate(st2(nlines1))
      allocate(fac1(nlines1))
      fac1='NA'
c      allocate(airport(nlines))
c      airport%srcgrp='ALL'
c      airport%cat='NR-Point-Airports'
      eof=0
      ilines=0
      read(xwalk,*)junk
      read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann
 15   if (eof .eq. 0) then
          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          if (ihap .gt. 0) then
              if ((lsingst .and. (trim(adjustl(st2a)) .eq.  
     +        trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) 
     +        .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. 
     +        trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. 
     +        (.not. lsingst .and. .not. lsingfac)) then
                  ilines=ilines+1
c                 hap1(ilines)=cas1
                  st2(ilines)=st2a
                  fac1(ilines)=facid2
              endif
          endif
          read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann
          goto 15
      endif
 

c     now get unique list of haps/sources. sources based on run#
      call haps_sources(nlines1)
      
      rewind(xwalk)
      eof=0
      ilines=0
      read(xwalk,*)junk
      read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann
 20   if (eof .eq. 0) then
          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          
          if (ihap .gt. 0)then
              is=findfac(facid2)
              if (is .gt. 0) ilines=ilines+1
          endif
          read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann
          goto 20
      endif           
      
      nlines=ilines

      allocate(airport(nlines))
      airport%srcgrp='ALL'
      airport%cat='NR-Point-Airports'
c     now assign to the array based on HAPs and sources
      rewind(xwalk)
      ilines=0
      eof=0
      read(xwalk,*)junk
      read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann
 25   if (eof .eq. 0) then
          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          if (ihap .gt. 0 )then
              is=findfac(facid2)
              if (is .gt. 0) then
                  ilines=ilines+1
                  airport(ilines)%state=st2a
                  if (st2a .eq. '02') then
                      lgrids(2)=.true.
                  elseif (st2a .eq. '15') then
                      lgrids(3)=.true.
                  elseif (st2a .eq. '72' .or. st2a .eq. '78') then
                      lgrids(4)=.true.
                  else
                      lgrids(1)=.true.
                  endif
                  
                  airport(ilines)%facid=facid2
                  airport(ilines)%cas=cas1
                  airport(ilines)%emis=ann
              endif
          endif
          read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann
          goto 25
      endif      
c      write(*,*)nlines
      close(xwalk)
c      do ilines=1,nlines
c         write(107,'(3(a,1x),e13.6)')
c    +    airport(ilines)%state,airport(ilines)%facid,
c    +    airport(ilines)%cas,airport(ilines)%emis
c     enddo
 5    format(a,/' cross-walk file does not exist...stopping program')
      
      return
      end
c##########################################################################
c     subroutine to read pollutant/category/AERMOD source cross-walk
c     for point
      subroutine read_point
      use main1
      implicit none
      
      integer*8 eof,ilines,ftype,i,ilen,ihap,findhap,nlines1,findfac,is,
     + prid,uid,rid,n1,n2,nn,nn1
      character facid2*15,relid*8,scc1*8,cas1*15,src*12,
     +junk*1,st2a*2,junk2*20,ares*2,eline*200
      real*8 ann
      logical lexist,lgo,lfix

      write(logunit,*)'Reading emissions file'
      
c     check for file existence
      inquire(file=xfile(1),exist=lexist)
      if (.not. lexist) then
          write(logunit,5)trim(adjustl(xfile(1)))
          stop
      endif

c     write the cross walk to a temporary cross walk and put in quotes for the 
c     facility name
c     close the existing cross walk and open the new one using the same file unit
c     open(unit=xwalk,file=xfile,status='old')
c     open(unit=105,file='temp_emis.txt',status='replace')
c     eof=0
c     do while (eof .eq. 0) 
c         read(xwalk,'(a200)',iostat=eof)eline
c         if (eof .eq. 0) then
c             nn=1
c             n1=0
c             n2=0
c             nn1=0
c             lfix=.false.
c             if (index(eline,'"') .eq. 0) then
c               find 2nd and 3rd commas
c               do while (nn .le. 200 .and. .not. lfix)
c                   if (eline(nn:nn) .eq. ',') then
c                       nn1=nn1+1
c                       if (nn1 .eq. 2)n1=nn
c                       if (nn1 .eq. 3)n2=nn
c                   endif
c                   if (n1 .gt. 0 .and. n2 .gt. 0)lfix=.true.
c                   nn=nn+1
c               enddo
c               if (eline(n1+1:n1+1) .ne. '"')eline(n1+1:n1+1)='"'
c               if (eline(n2-1:n2-1) .ne. '"') eline(n2-1:n2-1)='"'
c             endif
c             write(105,'(a200)')eline
c         endif
c     enddo
      
c     close(xwalk)
c     close(99)

      lgo=.false.
c     read in cross walk and assign variables to appropriate arrays
       open(unit=xwalk,file=xfile,status='old')
c     open(unit=xwalk,file='temp_emis.txt',status='old')
      ilines=0
      read(xwalk,*)junk
c     set initial value of ftype to 999 before reading in
c     because some records are missing the value
      ftype=999
      read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,rid,
     +cas1,ann
 10   if (eof .eq. 0) then

          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          if (ihap .gt. 0) then
               if ((lsingst .and. (trim(adjustl(st2a)) .eq.  
     +        trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) 
     +        .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. 
     +        trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. 
     +        (.not. lsingst .and. .not. lsingfac)) then
                  ilines=ilines+1
                  lhaps(ihap)=.true.
                  lgo=.true.
              endif
          endif
          ftype=999
          read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,
     +    rid,cas1,ann
          goto 10
      endif
      
      if (.not. lgo) then
          write(logunit,*)'Requested HAPs not found'
          stop
      endif
      
      rewind(xwalk)
c      nlines=ilines
      nlines1=ilines
c      allocate(hap1(nlines1))
      allocate(st2(nlines1))
      allocate(fac1(nlines1))
c      allocate(facility(nlines))
      
      eof=0
      ilines=0
      read(xwalk,*)junk
      ftype=999
      read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,rid,
     +cas1,ann
 15   if (eof .eq. 0) then
          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          if (ihap .gt. 0) then
              if ((lsingst .and. (trim(adjustl(st2a)) .eq.  
     +        trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) 
     +        .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. 
     +        trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. 
     +        (.not. lsingst .and. .not. lsingfac)) then
                  ilines=ilines+1
c                 hap1(ilines)=cas1
                  st2(ilines)=st2a
                  fac1(ilines)=facid2
              endif
          endif
          ftype=999
          read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,
     +    rid,cas1,ann
          goto 15
      endif
     
 
c     now get unique list of haps/sources. sources based on run#
      call haps_sources(nlines1)
      rewind(xwalk)
      eof=0
      ilines=0
      ftype=999
      read(xwalk,*)junk
      read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,rid,
     +cas1,ann
 20   if (eof .eq. 0) then
          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          if (ihap .gt. 0) then
              is=findfac(facid2)
              if (is .gt. 0) ilines=ilines+1
          endif
          ftype=999
          read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,
     +    rid,cas1,ann
          goto 20
      endif
      rewind(xwalk)
      eof=0
      nlines=ilines
      allocate(facility(nlines))
      ilines=0
      read(xwalk,*)junk
      ftype=999
      read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,rid,
     +cas1,ann
 25   if (eof .eq. 0) then
          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          if (ihap .gt. 0) then
              is=findfac(facid2)
              if (is .gt. 0) then
                  ilines=ilines+1
                  facility(ilines)%state=st2a
                  if (st2a .eq. '02') then
                      lgrids(2)=.true.
                  elseif (st2a .eq. '15') then
                      lgrids(3)=.true.
                  elseif (st2a .eq. '72' .or. st2a .eq. '78') then
                      lgrids(4)=.true.
                  else
                      lgrids(1)=.true.
                  endif
                  facility(ilines)%facid=facid2
                  facility(ilines)%cas=cas1
                  facility(ilines)%emis=ann
                  facility(ilines)%factype=ftype
                  facility(ilines)%srcgrp=src
                  facility(ilines)%procid=prid
                  facility(ilines)%unitid=uid
                  facility(ilines)%relid=rid
                  if (ftype .eq. 151) then
                      facility(ilines)%cat='NR-Point_Railyards'
                  else
                      facility(ilines)%cat='Point'
                  endif 
              endif
          endif
          ftype=999
          read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,
     +    rid,cas1,ann
          goto 25
      endif
      close(xwalk)
      do ilines=1,nlines
c         write(107,'(4(a,1x),i3,1x,e13.6)')
c    +    facility(ilines)%state,facility(ilines)%facid,
c    +    facility(ilines)%srcgrp,facility(ilines)%cas,
c    +    facility(ilines)%factype,facility(ilines)%emis
      enddo
    
 5    format(a,/' cross-walk file does not exist...stopping program')
      
      return
      end      
c##########################################################################
c     subroutine to read pollutant/category/AERMOD source cross-walk
c     for CMV sources
      subroutine read_cmv
      use main1
      implicit none
      
      integer*8 eof,ilines,i,ilen,ihap,findhap,is,findfac,nlines1,nn
      character facid2*15,st2a*2,rungrp*3,src*12,cat1*100,cas1*15,
     +tag*8,junk*1
      real*8 ann
      logical lexist,lgo
      
      write(logunit,*)'Reading emissions file'

c     if (trim(adjustl(runtyp)) .eq. 'PORT') then
c         tag='ports'
c     else
c         tag='underway'
c     endif
      
c     write(xfile,'(3(a))')trim(adjustl(xdir)),trim(adjustl(tag)),
c    +'_emis.csv'
      
c     check for file existence
      inquire(file=xfile(1),exist=lexist)
      if (.not. lexist) then
          write(logunit,5)trim(adjustl(xfile(1)))
          stop
      endif
      lgo=.false.
c     read in cross walk and assign variables to appropriate arrays
      open(unit=xwalk,file=xfile,status='old')
      ilines=0
      read(xwalk,*)junk
      read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann
 10   if (eof .eq. 0) then
          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          if (ihap .gt. 0) then
              if ((lsingst .and. (trim(adjustl(st2a)) .eq.  
     +        trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) 
     +        .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. 
     +        trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. 
     +        (.not. lsingst .and. .not. lsingfac)) then
                  ilines=ilines+1
                  lhaps(ihap)=.true.
                  lgo=.true.
              endif
          endif
          read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann
          goto 10
      endif
      
      if (.not. lgo) then
          write(logunit,*)'Requested HAPs not found'
          stop
      endif
      
      rewind(xwalk)
      nlines1=ilines
c      nlines=ilines
c      allocate(cmv(nlines))
c      allocate(hap1(nlines1))
      allocate(st2(nlines1))
      allocate(fac1(nlines1))
      
      eof=0
      ilines=0
      read(xwalk,*)junk
      read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann
 15   if (eof .eq. 0) then
          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          if (ihap .gt. 0) then
              if ((lsingst .and. (trim(adjustl(st2a)) .eq.  
     +        trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) 
     +        .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. 
     +        trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. 
     +        (.not. lsingst .and. .not. lsingfac)) then
                  ilines=ilines+1
c                 hap1(ilines)=cas1
                  st2(ilines)=st2a
                  fac1(ilines)=facid2
              endif
          endif
          read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann
          goto 15
      endif

c     now get unique list of haps/sources. sources based on run#
      call haps_sources(nlines1)
      rewind(xwalk)
      eof=0
      ilines=0
      
      read(xwalk,*)junk
      read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann
 20   if (eof .eq. 0) then
          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          if (ihap .gt. 0) then
              is=findfac(facid2)
              if (is .gt. 0) ilines=ilines+1
          endif
          read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann
          goto 20
      endif
      
      rewind(xwalk)
      eof=0
      nlines=ilines
      allocate(cmv(nlines))
      ilines=0
      read(xwalk,*)junk
      read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann
 25   if (eof .eq. 0) then
          ilen=len_trim(cas1)
          do i=1,ilen
              call upcase(cas1(i:i))
          enddo
          if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
          ihap=findhap(cas1)
          if (ihap .gt. 0) then
              is=findfac(facid2)
              if (is .gt. 0) then
                  ilines=ilines+1
                  cmv(ilines)%state=st2a
                  if (st2a .eq. '02') then
                      lgrids(2)=.true.
                  elseif (st2a .eq. '15') then
                      lgrids(3)=.true.
                  elseif (st2a .eq. '72' .or. st2a .eq. '78') then
                      lgrids(4)=.true.
                  else
                      lgrids(1)=.true.
                  endif
                  cmv(ilines)%facid=facid2  
                  cmv(ilines)%cas=cas1
                  cmv(ilines)%emis=ann
                  nn=index(src,'c')
                  write(cmv(ilines)%srcgrp,'(2(a1))')'C',src(nn+1:nn+1)
c                 cmv(ilines)%srcgrp=src
                  if (runtyp .eq. 'PORT' .or. runtyp .eq. 'UNDER') then
                      write(cmv(ilines)%cat,'(2(a))')'NR-CMV_',
     +                trim(adjustl(cat1))
                  else
                      cmv(ilines)%cat=cat1
                  endif
              endif
          endif
          read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann
          goto 25
      endif
      close(xwalk)
      do ilines=1,nlines
          write(107,'(4(a,1x),1x,e13.6)')
     +    cmv(ilines)%state,cmv(ilines)%facid,
     +    cmv(ilines)%srcgrp,cmv(ilines)%cas,
     +    cmv(ilines)%emis
      enddo
      
 5    format(a,/' cross-walk file does not exist...stopping program')
      
      return
      end
c##########################################################################
c     subroutine to read pollutant/category/AERMOD source cross-walk
c     for tract sources
      subroutine read_tract
      use main1
      implicit none
      
      integer*8 eof,ilines,i,ilen,is,ihap,findhap,findfac,nlines1,is1,m
      character st2a*2,rungrp*3,src*12,cat1*100,cas1*15,
     +junk*10,stfips1(4)*2,ares*2,fips6*6,poly*4,
     +tract1*11,facid2*16
      real*8 ann,memis1(12)!,win,spr,sum,fall
      logical lexist(4),lgo
      
      
      write(logunit,*)'Reading emissions file'
      
      stfips1(1)='02'
      stfips1(2)='15'
      stfips1(3)='72'
      stfips1(4)='78'
      write(ares,'(i2)')gridres
      lgo=.false.
      
      do is=1,4        
c         check for file existence
          inquire(file=xfile(is),exist=lexist(is))
      enddo
c     check to see if any files exist
      if (.not. lexist(1) .and. .not. lexist(2) .and. .not. lexist(3) 
     +.and. .not. lexist(4)) then
          write(logunit,*)' no x-walk files exist...stopping program'
          stop
      endif
      
      ilines=0
      do is=1,4
          if (lexist(is)) then
c             read in cross walk and assign variables to appropriate arrays
              open(unit=xwalk,file=xfile(is),status='old')
              read(xwalk,*)junk
              read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly,
     +        facid2,src,cat1,cas1,ann
 10           if (eof .eq. 0) then
                  ilen=len_trim(cas1)
                  do i=1,ilen
                      call upcase(cas1(i:i))
                  enddo
                  if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
                  ihap=findhap(cas1)
                  if (ihap .gt. 0) then
                      if ((lsingst .and. (trim(adjustl(st2a)) .eq.  
     +                trim(adjustl(states(istate))) .or. 
     +                trim(adjustl(st2a)) .eq. 
     +                trim(adjustl(states1(istate))))) .or. (lsingfac 
     +                .and. trim(adjustl(src)) .eq. 
     +                trim(adjustl(singlfac))) .or. (.not. lsingst .and. 
     +                .not. lsingfac)) then
                          ilines=ilines+1
                          lhaps(ihap)=.true.
                          lgo=.true.
                      endif
                  endif
                  read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly,
     +            facid2,src,cat1,cas1,ann
                  goto 10
              endif
              close(xwalk)
          endif
      enddo
      
      if (.not. lgo) then
          write(logunit,*)'Requested HAPs not found'
          stop
      endif  
          
      nlines1=ilines
c      allocate(hap1(nlines1))
      allocate(st2(nlines1))
      allocate(fac1(nlines1))
c      nlines=ilines
c      allocate(tracts(nlines))
      ilines=0
      do is=1,4
          if (lexist(is)) then
              open(unit=xwalk,file=xfile(is),status='old')
              eof=0
              read(xwalk,*)junk
              read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly,
     +        facid2,src,cat1,cas1,ann
 15           if (eof .eq. 0) then
                  ilen=len_trim(cas1)
                  do i=1,ilen
                      call upcase(cas1(i:i))
                  enddo
                  if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
                  ihap=findhap(cas1)
                  if (ihap .gt. 0) then
                      if ((lsingst .and. (trim(adjustl(st2a)) .eq.  
     +                trim(adjustl(states(istate))) .or. 
     +                trim(adjustl(st2a)).eq. 
     +                trim(adjustl(states1(istate))))) .or. (lsingfac 
     +                .and. trim(adjustl(src)) .eq. 
     +                trim(adjustl(singlfac))).or. (.not. lsingst .and. 
     +                .not. lsingfac)) then
                          ilines=ilines+1
c                         hap1(ilines)=cas1
                          st2(ilines)=st2a
                          fac1(ilines)=src
                      endif
                  endif
                  read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly,
     +            facid2,src,cat1,cas1,ann
                  goto 15
              endif
 
              close(xwalk)
          endif
      enddo
c     now get unique list of haps/sources. sources based on run#
      call haps_sources(nlines1)
      
      ilines=0
      do is=1,4
          if (lexist(is)) then
              open(unit=xwalk,file=xfile(is),status='old')
              eof=0
              read(xwalk,*)junk
              read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly,
     +        facid2,src,cat1,cas1,ann
              
 20           if (eof .eq. 0) then
                  ilen=len_trim(cas1)
                  do i=1,ilen
                      call upcase(cas1(i:i))
                  enddo
                  if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
                  ihap=findhap(cas1)
                  if (ihap .gt. 0) then
                      is1=findfac(src)
                      if (is1 .gt. 0) ilines=ilines+1
                  endif
                  read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly,
     +            facid2,src,cat1,cas1,ann
                  goto 20
              endif
 
              close(xwalk)
          endif
      enddo
      

      nlines=ilines
      ilines=0
      allocate(tracts(nlines))
      do is=1,4
          if (lexist(is)) then
              open(unit=xwalk,file=xfile(is),status='old')
              eof=0
              read(xwalk,*)junk
              if (runon) then
                  read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly,
     +            facid2,src,cat1,cas1,ann,(memis1(m),m=1,12) !,win,spr,sum,fall
              else
                  read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly,
     +            facid2,src,cat1,cas1,ann
              endif
              
 25           if (eof .eq. 0) then
                  ilen=len_trim(cas1)
                  do i=1,ilen
                      call upcase(cas1(i:i))
                  enddo
                  if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
                  ihap=findhap(cas1)
                  if (ihap .gt. 0) then
                      is1=findfac(src)
                      if (is1 .gt. 0) then
                          ilines=ilines+1
                          tracts(ilines)%state=st2a
                          tracts(ilines)%facid=src  !for tract, the AERMOD source ID is the facility ID
                          tracts(ilines)%cas=cas1
                          if (runon .or. runnon) then
                              do m=1,12
                                  tracts(ilines)%emis(m)=memis1(m)
                              enddo
                          endif
                          tracts(ilines)%emis(13)=ann
                          tracts(ilines)%srcgrp=src
                          tracts(ilines)%cat=cat1
                      endif
                  endif
                  if (runon) then
                      read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,
     +                poly,facid2,src,cat1,cas1,ann,(memis1(m),m=1,12) !win,spr,sum,fall
                  else
                      read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,
     +                poly,facid2,src,cat1,cas1,ann
                  endif
                  goto 25
              endif
 
              close(xwalk)
          endif
      enddo
      
      do ilines=1,nlines
          write(107,'(4(a,1x),1x,e13.6)')
     +    tracts(ilines)%state,tracts(ilines)%facid,
     +    tracts(ilines)%srcgrp,tracts(ilines)%cas,
     +    tracts(ilines)%emis(5)
      enddo
          
          
c      do ilines=1,nlines
c          write(105,*)tracts(ilines)%cas
c      enddo
      
 5    format(a,/' cross-walk file does not exist...stopping program')
      
      return
      end

c##########################################################################
c     subroutine to readpollutant/category/AERMOD source cross-walk
c     for gridded sources
      subroutine read_grid
      use main1
      implicit none
      
      integer*8 eof,ilines,i,ilen,ihap,findhap,isrc,findfac,ngsrc,nrows,
     +ncols,r,c,nsrc1,j,jj,k,singc,singr,m
      integer*8, allocatable, dimension(:,:) :: igrps
      character facid2a*15,relid*8,scc1*8,cas1*15,src*12,
     +junk*1,st2a*2,junk2*20,ares*2,afips*5,cat1*100,facid2*15,fid*15
      real*8 ann,memis1(12) !win,spr,sum,fall
      logical lexist,lf,lgo

      allocatable :: facid2(:)
      
      write(logunit,*)'Reading emissions file'

      write(ares,'(i2)')gridres
      
c     check for file existence
      inquire(file=xfile(1),exist=lexist)
      if (.not. lexist) then
          write(logunit,5)trim(adjustl(xfile(1)))
          stop
      endif
      
c     if single state
      if (lsingst) then
          open(unit=xwalk,file=xfile(1),status='old')
          ilines=0
          read(xwalk,*)junk
          read(xwalk,*,iostat=eof)junk2,afips,facid2a   
 20       if (eof .eq. 0) then
              st2a=afips(1:2)  !may change for v2, may be 5 characters instead of 6
              if (trim(adjustl(st2a)) .eq. trim(adjustl(states(istate))) 
     +        .or. trim(adjustl(st2a)) .eq. 
     +        trim(adjustl(states1(istate)))) isrc=1
              if (isrc .gt. 0) then
                  ilen=len_trim(cas1)
                  do i=1,ilen
                      call upcase(cas1(i:i))
                  enddo
                  if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
                  ihap=findhap(cas1)
                  if (ihap .gt. 0) then
                      read(facid2a,'(2(a1,i3))')junk,c,junk,r
                      call gridsrc(r,c,fid,lf)
                      if (lf) then
                          ilines=ilines+1
                          lhaps(ihap)=.true.
                          lgo=.true.
                      endif
                  endif
              endif
              read(xwalk,*,iostat=eof)junk2,afips,facid2a
              goto 20
          endif
          if (.not. lgo) then
              write(logunit,40)trim(adjustl(states(istate)))
              stop
          endif
          nsrc1=ilines
          allocate(facid2(nsrc1))
          rewind(xwalk)
          read(xwalk,*)junk
          read(xwalk,*,iostat=eof)junk2,afips,facid2a
 25       if (eof .eq. 0) then
              st2a=afips(1:2)
              if (trim(adjustl(st2a)) .eq. trim(adjustl(states(istate))) 
     +        .or. trim(adjustl(st2a)) .eq. 
     +        trim(adjustl(states1(istate)))) isrc=1
              if (isrc .gt. 0) then
                  ilen=len_trim(cas1)
                  do i=1,ilen
                      call upcase(cas1(i:i))
                  enddo
                  if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
                  ihap=findhap(cas1)
                  if (ihap .gt. 0) then
                      read(facid2a,'(2(a1,i3))')junk,c,junk,r
                      call gridsrc(r,c,fid,lf)
                      if (lf) then
                          ilines=ilines+1
                          facid2(i)=fid
                      endif
                  endif
              endif
              read(xwalk,*,iostat=eof)junk2,afips,facid2a
              goto 25
          endif 
          close(xwalk)
      endif        
      
c     get column/row of the single grid cell from command line
      if (lsingfac) then
          read(singlfac,'(2(a1,i3))')junk,singc,junk,singr
          call gridsrc(singr,singc,fid,lf)
          if (lf) then
              nsrc1=1
              allocate(facid2(nsrc1))
              facid2(1)=singlfac
          else
              write(logunit,45)trim(adjustl(singlfac))
              stop
          endif
      endif
      
c     not single state or grid cell
      if (.not. lsingfac .or. lsingst) then
          ngsrc=0
          nrows=246
          ncols=396
          do r=1,nrows
              do c=1,ncols
                  call gridsrc(r,c,fid,lf)
                  if (lf)ngsrc=ngsrc+1
              enddo
          enddo
          nsrc1=ngsrc
          ngsrc=0
          allocate(facid2(nsrc1))
          do r=1,nrows
              do c=1,ncols
                  call gridsrc(r,c,fid,lf)
                  if (lf) then
                      ngsrc=ngsrc+1
                      facid2(ngsrc)=fid
                  endif
              enddo
          enddo
      endif
      
c     common to scenarios, single state, single facility, or multiple facilities
      nsrc=nsrc1/nruns
      if (mod(real(nsrc1),real(nruns)) .ne. 0.0) nsrc=nsrc+1
      allocate(igrps(nruns,2))

c      write(*,*)'a ',nsrc1,nruns,runno,nsrc
c     pause   !delete later 
      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
      allocate(facid(nsrc))
      facid='NA'
      
      k=0
      do i=1,nsrc1
          if (i .ge. igrps(runno,1) .and. i .le. 
     +    igrps(runno,2)) then
              k=k+1
              facid(k)=facid2(i)
          endif
      enddo
      
      write(logunit,'(a,1x,i7)')'Total sources:',nsrc1
      write(logunit,'(a,1x,i7)')'Sources to process:',nsrc
      
c     read in cross walk and assign variables to appropriate arrays
      open(unit=xwalk,file=xfile,status='old')
      ilines=0
      read(xwalk,*)junk
      read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1,cas1,ann    
 10   if (eof .eq. 0) then
          isrc=findfac(facid2a)
          if (isrc .gt. 0) then
              ilen=len_trim(cas1)
              do i=1,ilen
                  call upcase(cas1(i:i))
              enddo
              if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
              ihap=findhap(cas1)
              if (ihap .gt. 0) then
                  ilines=ilines+1
                  lhaps(ihap)=.true.
                  lgo=.true.
              endif
          endif
          read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1,cas1,ann
          goto 10
      endif
      
      if (.not. lgo) then
          write(logunit,*)'Requested HAPs not found'
          close(xwalk)
          stop
      endif
      
      rewind(xwalk)
      nlines=ilines  
  
      call haps_sources(nlines)
      allocate(gridded(nlines))
      eof=0
      ilines=0
      read(xwalk,*)junk
      if (runnplow .or. runnphi .or. runrwc .or. runog .or. runag) then
          read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1,cas1,ann
      else
          read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1,cas1,
     +    ann,(memis1(m),m=1,12) !win,spr,sum,fall
      endif
 15   if (eof .eq. 0) then
          isrc=findfac(facid2a)
          if (isrc .gt. 0) then
              ilen=len_trim(cas1)
              do i=1,ilen
                  call upcase(cas1(i:i))
              enddo
              if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE'  !consistent across run groups
              call fixcat(cat1)
              ihap=findhap(cas1)
              if (ihap .gt. 0) then
                  call fixcat(cat1)
                  ilines=ilines+1
                  gridded(ilines)%fips=afips
                  if (afips(1:2) .eq. '02') then
                      lgrids(2)=.true.
                  elseif (afips(1:2) .eq. '15') then
                      lgrids(3)=.true.
                  elseif (afips(1:2) .eq. '72' .or. st2a .eq. '78') then
                      lgrids(4)=.true.
                  else
                      lgrids(1)=.true.
                  endif
                  gridded(ilines)%facid=facid2a
                  gridded(ilines)%cas=cas1
                  gridded(ilines)%emis(13)=ann
                  gridded(ilines)%srcgrp=src
                  gridded(ilines)%cat=cat1
                  if (runon .or. runnon) then
                      do m=1,12
                          gridded(ilines)%emis(m)=memis1(m)
                      enddo
                  endif
              endif
          endif
          if (runnplow .or. runnphi .or. runrwc .or. runog 
     +   .or. runag) then
              read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1,
     +        cas1,ann
          else
              read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1,
     +        cas1,ann,(memis1(m),m=1,12) !win,spr,sum,fall
          endif
          goto 15
      endif
 
      
      close(xwalk)
      deallocate(facid2)
      deallocate(igrps)      
 5    format(a,/' cross-walk file does not exist...stopping program')
 40   format(a,/'AERMOD output not found or no sources for state ',a)
 45   format(a,/'AERMOD output not found for ',a,
     +' or not listed in cross-walk')
      return
      end
c##########################################################################
c     subroutine to replace blanks and commas with underscores in
c     category names
      subroutine fixcat(cat1)
      use main1
      implicit none
      integer ilen,i
      character cat1*100
      
      ilen=len_trim(cat1)
      do i=1,ilen
          if (cat1(i:i) .eq. ' ' .or. cat1(i:i) .eq. ',') cat1(i:i)='_'
      enddo
      return
      end
      
      
c##########################################################################
c     subroutine to get unique list of HAPs,categories, and sources to be used when writing
c     HAP specific output files
      subroutine haps_sources(nlines1)
      use main1
      implicit none
      
      integer*8 i,ihap,n,ihap1,nsrc1,j,k,jj,ihap2,nlines1
      integer*8, allocatable, dimension(:,:) :: igrps
      character h*15,
     st*2,rt*15,facid2*15,st2a*2
      allocatable :: h(:)
      allocatable :: facid2(:)
      allocatable :: st2a(:)
      logical ldup,lreset
      
      write(logunit,*)'Getting HAPs and sources'
   
c NR-C1C2_ports
c NR-C3_ports
c NR-C1C2C3_underway
c NP-AgriculturalLivestock
   
c     determine if any requsted HAPS need to be removed from list for rungroup
c      if (nhap1 .gt. 0) then
      lreset=.false.
      ihap1=nhaps
      ihap2=nhaps
      do ihap=1,nhaps
          if (.not. lhaps(ihap)) then
              ihap1=ihap1-1
              lreset=.true.
              write(logunit,'(2(a))')'Removing ',
     +        trim(adjustl(haps(ihap)))
          endif
      enddo
      if (lreset) then
          allocate(h(nhaps))
          h=haps
          deallocate(haps)
          nhaps=ihap1
          allocate(haps(nhaps))
          ihap1=0
          do ihap=1,ihap2
              if (lhaps(ihap)) then
                  ihap1=ihap1+1
                  haps(ihap1)=h(ihap)
              endif
          enddo
          deallocate(h)
      endif
c      endif
      
c     temporary array of haps, sources, duplicates included
      ihap1=nlines1


      write(logunit,'(a,1x,i3)')'Number of HAPs to process:',nhaps
      
c     do ihap=1,nhaps
c         write(102,*)haps(ihap)
c     enddo
c      write(*,*)'haps done'
c     eliminate duplicate facid and state if applicable
c      write(*,*)'start sources'
      if (runcmv .or. runap .or. runpt) then ! .or. noncontr) then 
          ihap=1  
          do i=2,ihap1
              ldup=.false.
c             n=1
              if (runcmv .or. runap .or. runpt) then
                  n=i-1
              else
                  n=1
              endif
              do while (n .lt. i .and. .not. ldup)
                  if (fac1(i) .eq. fac1(n)) ldup=.true.
c              write(*,*)fac1(i),fac1(n),i,n,ldup,ihap
c              pause
                  n=n+1
              enddo 
              if (.not. ldup) ihap=ihap+1
          enddo
          nsrc=ihap
          nsrc1=ihap
c      write(*,*)'eliminated duplicates'
c      allocate(facid(nsrc))
c      if (runpt .or. runap .or. runcmv) allocate(state2(nsrc))
          allocate(facid2(nsrc1))
          if (allocated(st2)) allocate(st2a(nsrc1))
      
c      facid(1)=fac1(1)
c     if (runpt .or. runap .or. runcmv) state2(1)=st2(1)
          facid2(1)=fac1(1)
          if (allocated(st2a)) st2a(1)=st2(1)
          ihap=1
          do i=2,ihap1
              ldup=.false.
c          n=1
              if (runcmv .or. runap .or. runpt) then
                  n=i-1
              else
                  n=1
              endif
          
              do while (n .lt. i .and. .not. ldup)
                  if (fac1(i) .eq. fac1(n)) ldup=.true.
                  n=n+1
              enddo 
              if (.not. ldup) then
                  ihap=ihap+1
c              facid(ihap)=fac1(i)
c              if (runpt .or. runap .or. runcmv) state2(ihap)=st2(i)
                  facid2(ihap)=fac1(i)
                  if (allocated(st2a)) st2a(ihap)=st2(i)
              endif
          enddo     
      
c     break up 
          nsrc=nsrc1/nruns
          if (nsrc1 .lt. nruns) then
              write(logunit,5)nruns,nsrc1
              stop
          endif
      
          if (mod(real(nsrc1),real(nruns)) .ne. 0.0) nsrc=nsrc+1
      
          allocate(igrps(nruns,2))

c      write(*,*)'a ',nsrc1,nruns,runno,nsrc
      
          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
c      do j=1,nruns
c          write(*,*)'b ',j,igrps(j,1),igrps(j,2)
c      enddo
      
          jj=igrps(runno,2)-igrps(runno,1)+1
c      write(*,*)jj,igrps(runno,2),igrps(runno,2)
      
          allocate(facid(nsrc))
          facid='NA'
          if (allocated(st2a)) then
              allocate(state2(nsrc))
              state2='NA'
          endif
      
          k=0
          do i=1,nsrc1
              if (i .ge. igrps(runno,1) .and. i .le. 
     +        igrps(runno,2)) then
                  k=k+1
                  facid(k)=facid2(i)
                  if (allocated(st2a))state2(k)=st2a(i)
              endif
          enddo
      
          write(logunit,'(a,1x,i7)')'Total sources:',nsrc1
c         write(*,'(a,1x,i7)')'Sources to process:',nsrc
      
      endif
c     write(*,'(a,1x,i7)')'Total sources:',nsrc1
      write(logunit,'(a,1x,i7)')'Sources to process:',nsrc

c     do ihap=1,nsrc
c         write(104,*)facid(ihap)
c     enddo
      
c     write(*,*)'sources done'
c     duplicates for categories if applicable
c      version 17030, now assign categories, don't rely
c     on using xwalk because some HAPs don't have all
c     categories and want consistency between HAPs
      if (runap) then
          allocate(cats(1))
          cats='NR-Point-Airports'
          ncats=1
      elseif (runpt) then
          allocate(cats(2))
          cats(1)='Point'
          cats(2)='NR-Point_Railyards'
          ncats=2
      elseif (runcmv) then
          if (runtyp .eq. 'PORT') then
              allocate(cats(1))
              cats(1)='NR-CMV_ports'
              ncats=1
          elseif (runtyp .eq. 'UNDER') then
              allocate(cats(1))
              cats(1)='NR-CMV_underway'
              ncats=1
          else
              allocate(cats(3))
              cats(1)='NR-C1C2_ports'
              cats(2)='NR-C3_ports'
              cats(3)='NR-C1C2C3_underway'
              ncats=3
          endif
      elseif (runnphi) then
          allocate(cats(3))
          cats(1)='NP-FuelCombustion_not_RWC'
          cats(2)='NP-Industrial'
          cats(3)='NP-WasteDisposal'
          ncats=3
      elseif (runnplow) then
          allocate(cats(5))
          cats(1)='NP-CommercialCooking'
c          cats(2)='NP-FuelCombustion_not_RWC'
          cats(2)='NP-MiscellaneousNonindustrial'
          cats(3)='NP-SolventsCoatings'
          cats(4)='NP-StorageTransfer_BulkTerminals_GasStage1'
          cats(5)='NR-Locomotives'
          ncats=5
      elseif (runnon) then
          allocate(cats(7))
          cats(1)='NR-Agriculture'
          cats(2)='NR-CommercialEquipment'
          cats(3)='NR-CommercialLawnGarden'
          cats(4)='NR-Construction'
          cats(5)='NR-Recreational-inc-PleasureCraft'
          cats(6)='NR-ResidentialLawnGarden'
          cats(7)='NR-AllOther'
          ncats=7
      elseif (runog) then
          allocate(cats(1))
          cats(1)='NP-OilGas'
          ncats=1
      elseif (runon) then
          if (runtyp .eq. 'HDOFF') then
              allocate(cats(2))
              cats(1)='OR-HeavyDuty-OffNetwork-Diesel'
              cats(2)='OR-HeavyDuty-OffNetwork-Gas'
              ncats=2
          elseif (runtyp .eq. 'LDOFF') then
              allocate(cats(2))
              cats(1)='OR-LightDuty-OffNetwork-Diesel'
              cats(2)='OR-LightDuty-OffNetwork-Gas'
              ncats=2
          elseif (runtyp .eq. 'HDON') then
              allocate(cats(2))
              cats(1)='OR-HeavyDuty-OnNetwork-Diesel'
              cats(2)='OR-HeavyDuty-OnNetwork-Gas'
              ncats=2
          elseif (runtyp .eq. 'LDON') then
              allocate(cats(3))
              cats(1)='OR-LightDuty-OnNetwork-Diesel'
              cats(2)='OR-LightDuty-OnNetwork-Gas'
              cats(3)='OR-Refueling'
              ncats=3
          else  !hotelling
              allocate(cats(1))
              cats(1)='OR-HeavyDuty-Hoteling'
              ncats=1
          endif
      elseif (runag) then
          allocate(cats(1))
          cats(1)='NP-AgriculturalLivestock'
          ncats=1
      else !rwc
          allocate(cats(1))
          cats(1)='NP-ResidentialWoodCombustionRWC'
          ncats=1
      endif
      
c          ihap=1  
c          do i=2,ihap1
c              ldup=.false.
c              n=1
c              do while (n .lt. i .and. .not. ldup)
c                  if (cats1(i) .eq. cats1(n)) ldup=.true.
c                  n=n+1
c              enddo 
c              if (.not. ldup) ihap=ihap+1
c          enddo
c          ncats=ihap

c          allocate(cats(ncats))
          
      
c          cats(1)=cats1(1)
      
c          ihap=1
c          do i=2,ihap1
c              ldup=.false.
c              n=1
c              do while (n .lt. i .and. .not. ldup)
c                  if (cats1(i) .eq. cats1(n)) ldup=.true.
c                  n=n+1
c              enddo 
c              if (.not. ldup) then
c                  ihap=ihap+1
c                  cats(ihap)=cats1(i)
c              endif
c          enddo              
c      endif
      
c      do i=1,ncats
c          write(103,*)cats(i),ncats
c      enddo
c      write(*,*)'cats done'
      
c      deallocate(hap1)
      if (allocated(h)) deallocate(h)
      if (allocated(fac1)) deallocate(fac1)
      if (allocated(facid2)) deallocate(facid2)
      if (allocated(igrps)) deallocate(igrps)
      if (allocated(st2)) then
          deallocate(st2)
          deallocate(st2a)
      endif
      
 5    format('# of runs: ',i10,/'exceeds # of sources: ',i10/
     +'stopping program')
      return
      end   
c##########################################################################
c     subroutine to get gridded sources ids instead of checking emissions files
      subroutine gridsrc(r,c,fid,lf)
      use main1
      implicit none
      integer*8 r,c
      logical lf
      character fid*15,ares*2,inpfil*300
      write(fid,'(a1,i3.3,a1,i3.3)')'G',c,'R',r
      
      write(ares,'(i2.2)')gridres
      write(inpfil,'(8(a))')trim(adjustl(outdir(1))),  
     +'aermod_',trim(adjustl(runtyp)),'_',
     +trim(adjustl(ares)),'_',trim(adjustl(fid)),'.out'
      inquire(file=inpfil,exist=lf)
      return
      end
c##########################################################################
      subroutine src_cr
      use main1
      implicit none
      
      integer*8 cr,ns,i,c,r,ig
      character fid*15,junk*1
      
      
c     conus(1)=source is CONUS source
c     conus(2)=source is non-CONUS source
      conus=.false.
      readzero=.false.
c     determine if any sources are conus or nonconus, this is used in subroutine rec_matrix as well
      if (runcmv .or. runpt. or. runap) then
          ns=1
c         see if any sources are non-CONUS sources
          do while (ns .le. nsrc .and. .not. conus(2))
              if (state2(ns) .eq. 'AK' .or. state2(ns) .eq. 'HI' 
     +        .or. state2(ns) .eq. 'PR' .or. state2(ns) .eq. 'VI' 
     +        .or. state2(ns) .eq. '02' .or. state2(ns) .eq. '15' 
     +        .or. state2(ns) .eq. '72' .or. state2(ns) .eq. '78')
     +        conus(2)=.true.
              ns=ns+1
          enddo
        
c         see if any sources are CONUS sources
          ns=1
          do while (ns .le. nsrc .and. .not. conus(1))
              if (state2(ns) .ne. 'AK' .and. state2(ns) .ne. 'HI'
     +        .and. state2(ns) .ne. 'PR' .and. state2(ns) .ne. 'VI'
     +        .and. state2(ns) .ne. '02' .and. state2(ns) .ne. '15'
     +        .and. state2(ns) .ne. '72' .and. state2(ns) .ne. '78')
     +        conus(1)=.true.
              ns=ns+1
          enddo
      
      else if (noncontr) then
          conus(2)=.true.  !conus(1) is false
      else
          conus(1)=.true. !conus(2) is false
      endif
      
c     if conus(1) = true then get column/row of source
c     version 20002 now get column/row for all sources, not just CONUS

      if (runap .or. runcmv .or. runpt) then
          allocate(srccr1(nsrc))
          allocate(stgrid1(nsrc))
          srccr1=0
      else
          allocate(srccr(nsrc))
          allocate(stgrid(nsrc))
          srccr=0
          ncolrow=nsrc
      endif
      do ns=1,nsrc
          if (runap .or. runcmv .or. runpt) then
              call read_co(ns,cr)
              srccr1(ns)=cr
              if (srccr1(ns) .le. 0)readzero=.true.
              if (trim(adjustl(state2(ns))) .eq. '02') then
                  ig=2
              elseif (trim(adjustl(state2(ns))) .eq. '15') then
                  ig=3
              elseif (trim(adjustl(state2(ns))) .eq. '72' .or. 
     +        trim(adjustl(state2(ns))) .eq. '78') then
                  ig=4
              else
                  ig=1
              endif
              stgrid1(ns)=ig
          else
              fid=facid(ns)
              read(fid,'(a1,i3,a1,i3)')junk,c,junk,r
              if (c .le. 0 .or. r .le. 0) then
                  readzero=.true.
                  srccr(ns)=-1*(abs(c)*1000+abs(r))
              else
                  srccr(ns)=c*1000+r
              endif
!             20002, find grid for gridded source
              do ig=1,4
                  if (lgrids(ig)) stgrid(ns)=ig
              enddo
              
          endif
      enddo
      if (runap .or. runpt .or. runcmv) call sort_co

            
      return
      end

c##########################################################################  
c     subroutine to read CO pathway of AERMOD output file to get column and row
      subroutine read_co(ns,cr)
      use main1
      implicit none
      integer*8 cr,nc,nr,i,eof,nu,ns
      character inpfil*300,line*250,rt(2)*10
      logical lexist1,lgo,lexist2
c     borrow from subroutine src_recs for non-gridded sources
c     for point check to see if non-egu or egu or airport is
c     runway or non-runway
    
      if (runap) then 
          i=1
          lexist1=.false.
          do while (i .le. 2 .and. .not. lexist1) !check for runway or non-runway
              write(inpfil,'(5(a))')trim(adjustl(outdir(i))), 
     +        trim(adjustl(state2(ns))),'/aermod_',
     +        trim(adjustl(facid(ns))),'.out'
              inquire(file=inpfil,exist=lexist1)
              i=i+1
          enddo
          if (.not. lexist1) then
              write(logunit,5)trim(adjustl(inpfil))
              lgo=.false.
              goto 95
c              stop
          endif
      elseif (runpt) then
          i=1
          lexist1=.false.
          do while (i .le. 1 .and. .not. lexist1) !check for EGU or non-EGU
              write(inpfil,'(5(a))')trim(adjustl(outdir(i))), 
     +        trim(adjustl(state2(ns))),'/aermod_',
     +        trim(adjustl(facid(ns))),'.out'
              inquire(file=inpfil,exist=lexist1)
              i=i+1
          enddo
          if (.not. lexist1) then  !need to see if source run was split
              i=1
              lexist2=.false.
              do while (i .le. 2 .and. .not. lexist2) !check for EGU or non-EGU
                  write(inpfil,'(5(a))')trim(adjustl(outdir(i))), 
     +            trim(adjustl(state2(ns))),'/aermod_',
     +            trim(adjustl(facid(ns))),'_01.out'
                  inquire(file=inpfil,exist=lexist2)
                  i=i+1
              enddo
              if (.not. lexist2) then
                  write(logunit,5)trim(adjustl(inpfil))
                  lgo=.false.
                  goto 95
              endif
c              stop
          endif
      else 
          write(inpfil,'(5(a))')trim(adjustl(outdir(1))),
     +    trim(adjustl(state2(ns))),'/aermod_',
     +    trim(adjustl(facid(ns))),'.out'
          inquire(file=inpfil,exist=lexist1)
          if (.not. lexist1) then
              write(logunit,5)trim(adjustl(inpfil))
              lgo=.false.
              goto 95
          endif
      endif
      
      eof=0
      open(unit=inpunit,file=inpfil,status='old')
      read(inpunit,10,iostat=eof)line
      if (eof .eq. 0) then
          nc=index(line,'COL')
          nr=index(line,'ROW')
          nu=index(line,'UTM ZONE')
          if (index(line(nc+4:nr-1),'*') .gt. 0) then
              col=0
          else
              read(line(nc+4:nr-1),*)col
          endif
          if (index(line(nr+4:nu-1),'*') .gt. 0) then
              row=0
          else
              read(line(nr+4:nu-1),*)row
          endif
          cr=col*1000+row
       
      endif
      close(inpunit)
 5    format('READ CO AERMOD file ',a,' not found'/'skip')  
 10   format(//a200)
 95   return
      end

      
c##########################################################################
c     subroutine to sort source column/row
      subroutine sort_co
      use main1
      implicit none
      integer*8 i,n,ic
      logical ldup
      
c     20002, now check state as well since using grids for all areas
      ic=1
      do i=2,nsrc
          ldup=.false.
          n=1
          do while (n .lt. i .and. .not. ldup)
              if (srccr1(i) .eq. srccr1(n) .and. stgrid1(i) .eq. 
     +        stgrid1(n)) ldup=.true. !duplicate if col/row and state match
              n=n+1
          enddo
          if (.not. ldup) ic=ic+1
      enddo
      
      ncolrow=ic
      ic=1
      allocate(srccr(ncolrow))
      allocate(stgrid(ncolrow))
      srccr=0
      stgrid=0
      srccr(1)=srccr1(1)
      stgrid(1)=stgrid1(1)
      do i=2,nsrc
          ldup=.false.
          n=1
          do while (n .lt. i .and. .not. ldup)
              if (srccr1(i) .eq. srccr1(n) .and. stgrid1(i) .eq. 
     +        stgrid1(n)) ldup=.true.
              n=n+1
          enddo
          if (.not. ldup) then
              ic=ic+1
              srccr(ic)=srccr1(i)
              stgrid(ic)=stgrid1(i)
          endif
      enddo
      return
      end
      
c##########################################################################
c     create national matrix of receptors
      subroutine rec_matrix
      use main1
      implicit none
      
      integer*8 eof,irec,iunit,c,r,gr,br,mr,cr,br1,mr1,gr1,i
      real*8 x,y,barea
      character rfile*300,rid*15,fcode*5
      logical luse
      
c     version 17186
c     allocate gridded network if in CONUS area
c     allocate only non-CONUS blocks/monitors when only non-CONUS sources
c     allocate only CONUS blocks/monitors when only CONUS areas
c     both if a mix of CONUS/non-CONUS
      
      write(logunit,*)'getting receptor arrays'
      iunit=60
      
c     version 20002
c     now have gridded receptors for all areas
      write(logunit,*)'Getting gridded receptors'
      call datetime
      
c     first read in gridded receptors file(s)
c     loop through the potential receptor files
      do i=1,4
          if (lgrids(i)) then
              iunit=100*(60+i)
              write(rfile,'(3(a))')trim(adjustl(recdir)),
     +        trim(adjustl(gridst(i))),'_gridded_receptors.csv'
              open(unit=iunit,file=rfile,status='old')
      
              eof=0
              gr=0
              read(iunit,*,iostat=eof)rid
              read(iunit,*,iostat=eof)c,r
 5            if (eof .eq. 0) then
                  call checkrec(c,r,i,luse)
                  if (luse) gr=gr+1
                  read(iunit,*,iostat=eof)c,r
                  goto 5
              endif
              rewind(iunit)
          endif
      enddo
      ngrids=gr
      allocate(gridconc(nhaps,ngrids,ncats))
      allocate(gridrec(ngrids))
      gridrec%x=0.0
      gridrec%y=0.0
      gridrec%colrow=0
      gridrec%col=0
      gridrec%row=0
      gridrec%irec=0
      gridrec%irec1=0
      gridrec%iflag=0
      gridrec%recid='NA'
      gridrec%rtype='NA'
      gridconc=0.0

      
c     now fill in array
      do i=1,4
          if (lgrids(i)) then
              iunit=100*(60+i)
              eof=0
              gr=0
              gr1=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
                  gr1=gr1+1
                  call checkrec(c,r,i,luse)
                  if (luse) then
                      gr=gr+1
c                 version 20002, receptor ID now includes gridst
                      write(rid,'(a4,i6.6,i4.4)')
     +                trim(adjustl(gridst(i))),cr,irec
c              write(rid,'(i6,i4.4)')cr,irec
                      gridrec(gr)%recid=rid
                      gridrec(gr)%x=x
                      gridrec(gr)%y=y
                      gridrec(gr)%rtype='GRIDDED'
                      gridrec(gr)%colrow=cr
                      gridrec(gr)%col=c
                      gridrec(gr)%row=r
                      gridrec(gr)%irec=gr1 !master
                      gridrec(gr)%irec1=gr 
                  endif
                  read(iunit,*,iostat=eof)c,r,irec,x,y
                  goto 10
              endif  
              close(iunit)
          endif
      enddo
      
      

c     now blocks
c     version 20002, not all blocks in non-CONUS modeled, some will be interpolated
      iunit=60
      write(logunit,*)'Getting block receptors'
      call datetime
      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,x,y,c,r,barea
 15   if (eof .eq. 0) then
          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 (conus(2)) then
                  if (rid(1:2) .eq. '02') then
                      i=2
                  elseif (rid(1:2) .eq. '15') then
                      i=3
                  else
                      i=4
                  endif
                  
                  if (c .eq. 0 .and. r .eq. 0) then !keep blocks not in a grid
                      luse=.true.
                  else
                      call checkrec(c,r,i,luse)
                  endif
                  if (luse) br=br+1
              endif
          else
              if (conus(1)) then
                  i=1
                  call checkrec(c,r,i,luse)
                  if (luse) br=br+1
              endif
          endif
          read(iunit,*,iostat=eof)rid,x,y,c,r,barea
          goto 15
      endif
      rewind(iunit)
      nblocks=br
  
      br=0
      br1=0
      eof=0
!      allocate(blockid(nblocks))
      if (nblocks .gt. 0) then
          allocate(blokconc(nhaps,nblocks,ncats))
          allocate(blockrec(nblocks))
      
!      allocate(blockx(nblocks))
!      allocate(blocky(nblocks))
!      allocate(blockcr(nblocks))
          blokconc=0.0
          blockrec%x=0.0
          blockrec%y=0.0
          blockrec%colrow=0
          blockrec%col=0
          blockrec%row=0
          blockrec%irec=0
          blockrec%irec1=0
          blockrec%iflag=0
          blockrec%recid='NA'
          blockrec%rtype='NA'
      
          read(iunit,*,iostat=eof)rid
          read(iunit,*,iostat=eof)rid,x,y,c,r
 20       if (eof .eq. 0) then
              cr=c*1000+r
c             br=br+1
              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 (conus(2)) then
                      if (rid(1:2) .eq. '02') then
                          i=2
                      elseif (rid(1:2) .eq. '15') then
                          i=3
                      else
                          i=4
                      endif
                      if (c .eq. 0 .and. r .eq. 0) then
                          luse=.true.
                      else
                          call checkrec(c,r,i,luse)
                      endif
                      if (luse) br=br+1
                  else
                      luse=.false.
                  endif
              else
                  if (conus(1)) then
                      i=1
                      call checkrec(c,r,i,luse)
                      if (luse) br=br+1
                  else
                      luse=.false.
                  endif
              endif
              if (luse) then
                  blockrec(br)%recid=rid
                  blockrec(br)%rtype='BLOCK'
                  blockrec(br)%x=x
                  blockrec(br)%y=y
                  blockrec(br)%colrow=cr
                  blockrec(br)%col=c
                  blockrec(br)%row=r
                  blockrec(br)%irec=br1 !master record #
                  blockrec(br)%irec1=br 
              endif
!             blockid(br)=rid
!             blockx(br)=x
!             blocky(br)=y
!             blockcr(br)=cr
c             write(108,*)blockrec(br)%recid,blockrec(br)%x,blockrec(br)%y,
c     +       blockrec(br)%colrow
              read(iunit,*,iostat=eof)rid,x,y,c,r
              goto 20
          endif
      endif
      close(iunit)
  
c     now nonpopulated blocks in non-conus areas
c     version 20002, now only if source is not in CMAQ domain
c     only in AK
      if (readzero) then
          write(logunit,*)'Getting nonpoulated block receptors'
          call datetime
          write(rfile,'(2(a))')trim(adjustl(recdir)),
     +    'zero_pop_block_receptors.csv'
          open(unit=iunit,file=rfile,status='old')
          eof=0
          br=0
          read(iunit,*,iostat=eof)rid
          read(iunit,*,iostat=eof)rid,x,y,c,r,barea
 35       if (eof .eq. 0) then
              if (rid(1:2) .eq. '02' .and. c .eq. 0 .and. r .eq. 0 )
     +        br=br+1
              read(iunit,*,iostat=eof)rid,x,y,c,r,barea
              goto 35
          endif
          rewind(iunit)
          nzblocks=br
  
          br=0
          br1=0
          eof=0
!      allocate(blockid(nblocks))
          if (nzblocks .gt. 0) then
              allocate(zblokconc(nhaps,nzblocks,ncats))
              allocate(zblockrec(nzblocks))
      
!      allocate(blockx(nblocks))
!      allocate(blocky(nblocks))
!      allocate(blockcr(nblocks))
              zblokconc=0.0
              zblockrec%x=0.0
              zblockrec%y=0.0
              zblockrec%colrow=0
              zblockrec%col=0
              zblockrec%row=0
              zblockrec%irec=0
              zblockrec%irec1=0
              zblockrec%iflag=0
              zblockrec%recid='NA'
              zblockrec%rtype='NA'
      
              read(iunit,*,iostat=eof)rid
              read(iunit,*,iostat=eof)rid,x,y,c,r
 40           if (eof .eq. 0) then
                  cr=c*1000+r
c                 br=br+1
                  br1=br1+1
                  if (rid(1:2) .eq. '02' .and. c .eq. 0 .and. 
     +            r .eq. 0) then
                      br=br+1
                      luse=.true.
                  else
                      luse=.false.
                  endif
                  if (luse) then
                      zblockrec(br)%recid=rid
                      zblockrec(br)%rtype='BLOCK'
                      zblockrec(br)%x=x
                      zblockrec(br)%y=y
                      zblockrec(br)%colrow=cr
                      zblockrec(br)%col=c
                      zblockrec(br)%row=r
                      zblockrec(br)%irec=br1 !master record #
                      zblockrec(br)%irec1=br 
                  endif
!                 blockid(br)=rid
!                 blockx(br)=x
!                 blocky(br)=y
!                 blockcr(br)=cr
c                 write(108,*)blockrec(br)%recid,blockrec(br)%x,blockrec(br)%y,
c     +           blockrec(br)%colrow
                  read(iunit,*,iostat=eof)rid,x,y,c,r
                  goto 40
              endif
          endif
          close(iunit)
      endif
      
c     now monitors
      write(logunit,*)'Getting monitor receptors'
      call datetime
      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,x,y,c,r
 25   if (eof .eq. 0) then
          if (fcode(1:2) .eq. '02' .or. fcode(1:2) .eq. '15' .or. 
     +    fcode(1:2) .eq. '72' .or. fcode(1:2) .eq. '78') then
              if (conus(2)) then
                  if (fcode(1:2) .eq. '02') then
                      i=2
                  elseif (fcode(1:2) .eq. '15') then
                      i=3
                  else
                      i=4
                  endif
                  call checkrec(c,r,i,luse)
                  if (luse) mr=mr+1
              endif
          else
              if (conus(1)) then
                  i=1
                  call checkrec(c,r,i,luse)
                  if (luse) mr=mr+1
              endif
          endif
          read(iunit,*,iostat=eof)rid,fcode,x,y,c,r
          goto 25
      endif
      rewind(iunit)
      nmons=mr
      mr=0
      mr1=0
      eof=0
!      allocate(monid(nmons))
!     allocate(mfips(nmons))
      if (nmons .gt. 0) then
          allocate(monconc(nhaps,nmons,ncats))
          allocate(monrec(nmons))
!         allocate(monx(nmons))
!         allocate(mony(nmons))
!         allocate(moncr(nmons))
      
          monconc=0.0
          monrec%x=0.0
          monrec%y=0.0
          monrec%colrow=0
          monrec%col=0
          monrec%row=0
          monrec%irec=0
          monrec%irec1=0
          monrec%iflag=0
          monrec%recid='NA'
          monrec%rtype='NA'
          monrec%fip='NA'
      
          read(iunit,*,iostat=eof)rid
          read(iunit,*,iostat=eof)rid,fcode,x,y,c,r
 30       if (eof .eq. 0) then
              cr=c*1000+r
              mr1=mr1+1
              if (fcode(1:2) .eq. '02' .or. fcode(1:2) .eq. '15' .or. 
     +        fcode(1:2) .eq. '72' .or. fcode(1:2) .eq. '78') then
                 if (conus(2)) then
                     if (fcode(1:2) .eq. '02') then
                         i=2
                     elseif (fcode(1:2) .eq. '15') then
                         i=3
                     else
                         i=4
                     endif
                     call checkrec(c,r,i,luse)
                     if (luse) mr=mr+1
                     
                 else
                     luse=.false.
                 endif
              else
                  if (conus(1)) then
                      i=1
                      call checkrec(c,r,i,luse)
                      if (luse) mr=mr+1
                  else
                      luse=.false.
                  endif
              endif
              if (luse) then
                  monrec(mr)%x=x
                  monrec(mr)%y=y
                  monrec(mr)%colrow=cr
                  monrec(mr)%col=c
                  monrec(mr)%row=r
                  monrec(mr)%irec=mr1
                  monrec(mr)%irec1=mr
                  monrec(mr)%recid=rid
                  monrec(mr)%fip=fcode
                  monrec(mr)%rtype='MONITOR'
              endif
!             monid(mr)=rid
!             mfips(mr)=fcode
!             monx(mr)=x
!             mony(mr)=y
!             moncr(mr)=cr
c             write(109,*)monrec(mr)%recid,monrec(mr)%x,monrec(mr)%y,
c     +       monrec(mr)%fip,monrec(mr)%colrow
              read(iunit,*,iostat=eof)rid,fcode,x,y,c,r
              goto 30
          endif
      endif

      write(logunit,*)'number of receptors ',ngrids,nblocks,nmons
c     write(*,'(a,3(1x,i8.0))')'Number of receptors',
c    + ngrids,nblocks,nmons
      
      close(iunit)      
      return
      end
c##########################################################################
c     subroutine to check if a receptor is within a specified block
c     of grid cells of source cells
      subroutine checkrec(c,r,ig,luse)
      use main1
      implicit none
      integer*8 c,r,modcol,modrow,i,cr,col1,row1,col2,row2,scol,srow,ig
      logical luse
      
      if (.not. runcmv) then
c          modcol=5
c          modrow=5
          if (ig .eq. 1) then !US
              modcol=8
              modrow=8
          elseif (ig .eq. 2) then !AK
              modcol=7
              modrow=7
          else !HI OR PR/VI
              modcol=20
              modrow=20
          endif
      else
          modcol=40
          modrow=40
      endif
      
      
      
      cr=c*1000+r
      i=1
      luse=.false.
c     version 20002, make sure looking at correct grid
      do while (i .le. ncolrow .and. .not. luse)
          scol=int(srccr(i)/1000)
          srow=srccr(i)-scol*1000
          col1=scol-modcol
          col2=scol+modcol
          row1=srow-modrow
          row2=srow+modrow
          if (c .ge. col1 .and. c .le. col2 .and. r .ge. row1 .and. 
     +    r .le. row2 .and. ig .eq. stgrid(i)) luse=.true.
          i=i+1
      enddo
      return
      end
c##########################################################################
c     subroutine to loop through sources and read model output
      subroutine read_model
      use main1
      implicit none
      integer*8 ns
      logical lgo
      do ns=1,nsrc
          if (trim(adjustl(facid(ns))) .ne. 'NA') then
              if (runpt .or. runap .or. runcmv) then
                  write(logunit,'(2(a),1x,a)')'Start ',
     +            trim(adjustl(state2(ns))),trim(adjustl(facid(ns)))
              else
                  write(logunit,'(2(a))')'Start ',
     +            trim(adjustl(facid(ns)))    
              endif
              call datetime
              if (runcmv .or. runpt. or. runap) then ! .or. noncontr) then
                  if (state2(ns) .eq. 'AK' 
     +            .or. state2(ns) .eq. 'HI' .or. state2(ns) .eq. 'PR' 
     +            .or. state2(ns) .eq. 'VI' .or. state2(ns) .eq. '02' 
     +            .or. state2(ns) .eq. '15' .or. state2(ns) .eq. '72' 
     +            .or. state2(ns) .eq. '78') then
                      akhiprvi=.true.
                  else
                      akhiprvi=.false.
                  endif
              else
                  if (noncontr) then
                      akhiprvi=.true.
                  else
                      akhiprvi=.false.
                  endif
              endif
              if (writefac) call cleanfiles(ns)  !check and delete HAP facility files before processing
              call src_recs(ns,lgo)   !get source locations and modeled receptors
              if (lgo) then
                  call int_recs(ns)   !get blocks and monitors to interpolate 
              
c                 if (.not. akhiprvi) then
                  call interp(ns) !interpolate blocks and monitors
c                 endif
                  call hapconc(ns)        !calculate HAP specific concentrations and output results
                  call cleanup        !deallocate arrays specific to the source
              else
                  write(logunit,*)'No processing for ',
     +            trim(adjustl(facid(ns)))
              endif
          
              call datetime
              if (runpt .or. runcmv .or. runpt) then
                  write(logunit,'(2(a),1x,a)')'End ',
     +            trim(adjustl(state2(ns))),trim(adjustl(facid(ns)))
              else
                  write(logunit,'(2(a))')'End ',trim(adjustl(facid(ns)))
              endif
          
          endif
      enddo            

      return
      end
c##########################################################################
c     subroutine to clean HAP-facility output files before new processing
      subroutine cleanfiles(ns)
      use main1
      implicit none
      integer*8 ns,h
      character ofile*250
      logical lexist 
      
      write(logunit,*)'Deleting old HAP-facility output files'
      do h=1,nhaps
          write(ofile,'(5(a))')trim(adjustl(hapdir)),
     +    trim(adjustl(haps(h))),'_',
     +    trim(adjustl(facid(ns))),'_concs.txt'
          inquire(file=ofile,exist=lexist)
          if (lexist) then
              open(unit=99,file=ofile,status='old')
              close(99,status='delete')
          endif
      enddo
      
      return
      end
      
c##########################################################################
c     subroutine to read in aermod output and  plot file
      subroutine src_recs(ns,lgo)
      use main1
      implicit none
      
      integer*8 eof,igrid,iblock,isrc,irec,nc,nr,nu,imon,n,ns,i,
     +imonth,nv,izblock,findsrc,m1,iyr,idir,ist1,get_state,imonth1
      real*8 x,y,conc,z1,z2,z3,emis(12)
c     character block*15,pltfil*300,inpfil*300,line*250,ave*6,srcgrp*8,
      character pltfil*300,inpfil*300,line*250,ave*6,
     +srcgrp*12,rt(2)*10,ares*2,errfil*300,srcgrp1*12
      logical lexist1,lrecstart,lrecend,lexist2,lexist3,lgo,lmonth,
     +lsplit,lexist,lexist1a,lexist2a
      write(logunit,*)'Reading source and modeled receptors'
      
      lsplit=.false.
c     for point check to see if non-egu or egu or airport is
c     runway or non-runway
      if (runap) then 
          i=1
          lexist1=.false.
          lexist2=.false.
          do while (i .le. 2 .and. .not. lexist1)
              write(inpfil,'(5(a))')trim(adjustl(outdir(i))), 
     +        trim(adjustl(state2(ns))),
     +        '/aermod_',trim(adjustl(facid(ns))),'.out'
              inquire(file=inpfil,exist=lexist1)
              write(pltfil,'(5(a))')trim(adjustl(outdir(i))),
     +        trim(adjustl(state2(ns))),
     +        '/',trim(adjustl(facid(ns))),'.plt'
              inquire(file=pltfil,exist=lexist2)
              i=i+1
          enddo
      elseif (runpt) then
          i=1
          lexist1=.false.
          lexist2=.false.
          do while (i .le. 1 .and. .not. lexist1)
              write(inpfil,'(5(a))')trim(adjustl(outdir(i))), 
     +        trim(adjustl(state2(ns))),
     +        '/aermod_',trim(adjustl(facid(ns))),'.out'
              inquire(file=inpfil,exist=lexist1)
              write(pltfil,'(5(a))')trim(adjustl(outdir(i))),
     +        trim(adjustl(state2(ns))),
     +        '/',trim(adjustl(facid(ns))),'.plt'
              inquire(file=pltfil,exist=lexist2)
              i=i+1
          enddo    
          if (.not. lexist1 .and. .not. lexist2) then
              i=1
              lexist1a=.false.
              lexist2a=.false.
              do while (i .le. 2 .and. .not. lexist1a)
                  write(inpfil,'(5(a))')trim(adjustl(outdir(i))), 
     +            trim(adjustl(state2(ns))),
     +            '/aermod_',trim(adjustl(facid(ns))),'_01.out'
                  inquire(file=inpfil,exist=lexist1a)
                  write(pltfil,'(5(a))')trim(adjustl(outdir(i))),
     +            trim(adjustl(state2(ns))),
     +            '/',trim(adjustl(facid(ns))),'_01.plt'
                  inquire(file=pltfil,exist=lexist2a)
                  i=i+1
              enddo
              if (lexist1a .or. lexist2a) lsplit=.true.
              idir=i-1
          endif    
      elseif (runcmv) then 
          i=1
          write(inpfil,'(6(a))')trim(adjustl(outdir(i))),
     +    trim(adjustl(state2(ns))),'/aermod_',
     +    trim(adjustl(facid(ns))),'.out'
          write(pltfil,'(6(a))')trim(adjustl(outdir(i))),
     +    trim(adjustl(state2(ns))),'/',
     +    trim(adjustl(facid(ns))),'.plt'
c     elseif (noncontr) then
c         ist1=get_state(ncst)
c         write(ares,'(i2.2)')gridres
c         write(inpfil,'(9(a))')trim(adjustl(outdir(i))),  
c    +    trim(adjustl(states1(ist1))),'/aermod_',trim(adjustl(runtyp)),
c    +    '_',trim(adjustl(ares)),'_',trim(adjustl(facid(ns))),'.out'
c         if (runon .or. runnon) then
c             write(pltfil,'(9(a))')trim(adjustl(outdir(i))),
c    +        trim(adjustl(states1(ist1))),'/',trim(adjustl(runtyp)),'_'
c    +        ,trim(adjustl(ares)),'_',trim(adjustl(facid(ns))),
c    +        '_month.pst'
c             write(errfil,'(5(a),i2.2,3(a))')trim(adjustl(outdir(i))),
c    +        trim(adjustl(states1(ist1))),'/',trim(adjustl(runtyp)),'_'
c    +        ,gridres,'_',trim(adjustl(facid(ns))),'_errors.out'
c         else
c             write(pltfil,'(9(a))')trim(adjustl(outdir(i))),
c    +        trim(adjustl(states1(ist1))),'/',trim(adjustl(runtyp)),'_'
c    +        ,trim(adjustl(ares)),'_',trim(adjustl(facid(ns))),'.plt'
c         endif
      else   !gridded, check for seasonal files for mobile, otherwise, annual files
          i=1
          write(ares,'(i2.2)')gridres
          write(inpfil,'(8(a))')trim(adjustl(outdir(i))),  
     +    'aermod_',trim(adjustl(runtyp)),'_',
     +    trim(adjustl(ares)),'_',trim(adjustl(facid(ns))),'.out'
          if (runon .or. runnon) then
              write(pltfil,'(7(a))')trim(adjustl(outdir(i))),
     +        trim(adjustl(runtyp)),'_',trim(adjustl(ares)),'_',
     +        trim(adjustl(facid(ns))),'_month.pst'
              write(errfil,'(3(a),i2.2,3(a))')trim(adjustl(outdir(i))),
     +        trim(adjustl(runtyp)),'_',gridres,'_',
     +        trim(adjustl(facid(ns))),'_errors.out'
          else
              write(pltfil,'(7(a))')trim(adjustl(outdir(i))),
     +        trim(adjustl(runtyp)),'_',trim(adjustl(ares)),'_',
     +        trim(adjustl(facid(ns))),'.plt'
          endif
      endif
      inquire(file=inpfil,exist=lexist1)
      inquire(file=pltfil,exist=lexist2)
      if (runon .or. runnon) then
          inquire(file=errfil,exist=lexist3)
          if (.not. lexist1 .or. .not. lexist2 .or. .not. lexist3) then
              if (.not. lexist1) write(logunit,5)trim(adjustl(inpfil))
              if (.not. lexist2) write(logunit,5)trim(adjustl(pltfil))
              if (.not. lexist3) write(logunit,5)trim(adjustl(errfil))
              lgo=.false.
              goto 95
          endif 
      else
          if (.not. lexist1 .or. .not. lexist2) then
              if (.not. lexist1) write(logunit,5)trim(adjustl(inpfil))
              if (.not. lexist2) write(logunit,5)trim(adjustl(pltfil))
              lgo=.false.
              goto 95
          endif
      endif
c seasonal files: runtyp_res_facid_season.plt res is 04 or 12
c     read aermod.out file to get # of sources, # of gridded receptors, # of block receptors, and # of monitors
c     that were explicitly modeled (this number for later writing of HAP specific concentrations)

      igrid=0
      iblock=0
      izblock=0
      eof=0
      isrc=0
      imon=0
      lrecstart=.false.
      lrecend=.false.
      lgo=.true.
      modyear=0
      open(unit=inpunit,file=inpfil,status='old')
      read(inpunit,10,iostat=eof)line
 15   if (eof .eq. 0) then
          if (index(line,'CO TITLEONE') .gt. 0) 
     +    read(line(13:16),'(i4)')modyear
          if (index(line,'CO TITLETWO') .gt. 0) then
              nc=index(line,'COL')
              nr=index(line,'ROW')
              nu=index(line,'UTM ZONE')
              if (index(line(nc+4:nr-1),'*') .gt. 0) then
                  col=0
              else
                  read(line(nc+4:nr-1),*)col
              endif
              if (index(line(nr+4:nu-1),'*') .gt. 0) then
                  row=0
              else
                  read(line(nr+4:nu-1),*)row
              endif
c             read(line(nc+4:nr-1),*)col
c             read(line(nr+4:nu-1),*)row
              read(line(nu+8:nu+18),*)utmzone
           
          endif
          if (index(line,'SO LOCATION') .gt. 0 .or. 
     +    index(line,'** LOCATION') .gt. 0)isrc=isrc+1
          if (index(line,'** NUMBER OF BLOCKS') .gt. 0) then
              n=index(line,'** NUMBER OF BLOCKS')
              read(line(n+20:250),*)iblock
          endif
          if (index(line,'** NUMBER OF ZERO BLOCKS') .gt. 0) then
              n=index(line,'** NUMBER OF ZERO BLOCKS')
              read(line(n+25:250),*)izblock
          endif
          if (index(line,'** NUMBER OF GRIDDED RECEPTORS') .gt. 0) then
              n=index(line,'** NUMBER OF GRIDDED RECEPTORS')
              read(line(n+31:250),*)igrid
          endif
          if (index(line,'** NUMBER OF MONITORS') .gt. 0) then
              n=index(line,'** NUMBER OF MONITORS')
              read(line(n+22:250),*)imon
          endif
          if (index(line,'AERMOD Finishes UN-successfully') .gt. 0)
     +    lgo=.false.
          read(inpunit,10,iostat=eof)line
          goto 15
      endif
      
c     read errors file to get # of valid hours per month
      if (runon .or. runnon) call readerr(errfil)
      
      ngrps=isrc
      if (runap) then
          ngrps1=1
      else
          ngrps1=ngrps
      endif
      
c      if (runpt) then
c          ngrps=isrc
c      else
c          ngrps=1
c      endif
      
      nmodg=igrid
      nmodb=iblock
      nmodm=imon
      nmodz=izblock
c       write(*,*)nmodb,nmodg,nmodm,ngrps
c      pause
      if (lgo) then
          if (runap) then
              allocate(apsrcs(ngrps))
              apsrcs%utmx1=-999.0
              apsrcs%utmy1=-999.0
              apsrcs%utmx2=-999.0
              apsrcs%utmy2=-999.0
              apsrcs%frac=-999.0
              apsrcs%xdim=-999.0
              apsrcs%ydim=-999.0
              apsrcs%area=-999.0
              apsrcs%angle=-999.0
              apsrcs%srctyp='NA'
              apsrcs%srcid='NA'
          elseif (runpt) then
              allocate(ptsrcs(ngrps))
              ptsrcs%srcid='NA'
              ptsrcs%srctyp='NA'
              ptsrcs%utmx=-999.0
              ptsrcs%utmy=-999.0
              ptsrcs%xdim=-999.0
              ptsrcs%ydim=-999.0
              ptsrcs%angle=-999.0
          elseif (runcmv) then
              allocate(cmvsrcs(ngrps))
              cmvsrcs%srcid='NA'
              cmvsrcs%srctyp='NA'
              cmvsrcs%utmx1=-999.0
              cmvsrcs%utmy1=-999.0
              cmvsrcs%nvert=-9
c          elseif (noncontr) then
c              allocate(tractsrcs(ngrps))
c              tractsrcs%srcid='NA'
c              tractsrcs%srctyp='NA'
c              tractsrcs%utmx1=-999.0
c              tractsrcs%utmy1=-999.0
c              tractsrcs%nvert=-9
          else
              allocate(gridsrcs(ngrps))
              gridsrcs%utmx1=-999.0
              gridsrcs%utmy1=-999.0
              gridsrcs%utmx2=-999.0
              gridsrcs%utmy2=-999.0
              gridsrcs%utmx3=-999.0
              gridsrcs%utmy3=-999.0
              gridsrcs%utmx4=-999.0
              gridsrcs%utmy4=-999.0
              gridsrcs%nvert=-9
              gridsrcs%srcid='NA'
          endif
      else
          goto 95
      endif
      
      allocate(memis(ngrps,13))
      memis=0.0
      emis=0.0
      if (nmonths .eq. 1) then
          do isrc=1,ngrps
              memis(isrc,13)=10000.
          enddo
      endif
      
      
      if (nmodg .gt. 0) then
          allocate(modgc(ngrps1,nmodg,nmonths))
          allocate(modgrid(nmodg))
          modgc=0.0
          modgrid%ux=0.0
          modgrid%uy=0.0
          modgrid%x=0.0
          modgrid%y=0.0
          modgrid%colrow=0
          modgrid%iflag=0
          modgrid%recid='NA'
          modgrid%mastrec=0
          modgrid%dist=0.0
      endif
      
      if (nmodb .gt. 0) then
          allocate(modbc(ngrps1,nmodb,nmonths))
          allocate(modblock(nmodb))
          modbc=0.0
          modblock%ux=0.0
          modblock%uy=0.0
          modblock%x=0.0
          modblock%y=0.0
          modblock%recid='NA'
          modblock%iflag=0
          modblock%mastrec=0
          modblock%dist=-9.0
      endif
      if (nmodz .gt. 0) then
          allocate(modzbc(ngrps1,nmodz,nmonths))
          allocate(zmodblock(nmodz))
          modzbc=0.0
          zmodblock%ux=0.0
          zmodblock%uy=0.0
          zmodblock%x=0.0
          zmodblock%y=0.0
          zmodblock%recid='NA'
          zmodblock%iflag=0
          zmodblock%mastrec=0
          zmodblock%dist=-9.0
      endif
      if (nmodm .gt. 0) then
          allocate(modmc(ngrps1,nmodm,nmonths))
          allocate(modmon(nmodm))
          modmc=0.0
          modmon%ux=0.0
          modmon%uy=0.0
          modmon%x=0.0
          modmon%y=0.0
          modmon%iflag=0.0
          modmon%recid='NA'
          modmon%mastrec=0
          modmon%dist=-9.0
      endif
      
      call subsetrecs(ns)
      rewind(inpunit)
c     get gridded receptors Lambert coordinates,column/row
c     get block ID and column/row
      eof=0
      igrid=0
      iblock=0
      izblock=0
      imon=0
      irec=0
      lrecstart=.false.
      lrecend=.false.
      isrc=0
      lmonth=.false.
      read(inpunit,10,iostat=eof)line
 20   if (eof .eq. 0) then
          if (index(line,'SO LOCATION') .gt. 0 .or. 
     +    index(line,'** LOCATION') .gt. 0) then
              isrc=isrc+1
              call source_info(line,isrc)
          endif
          
          if (index(line,'SO SRCPARAM') .gt. 0 .or. 
     +    index(line,'** SRCPARAM')) then
              call source_info(line,isrc)
              if (runcmv) then ! .or. noncontr) then
                  nv=cmvsrcs(isrc)%nvert
c                  if (runcmv) then
c                      nv=cmvsrcs(isrc)%nvert
c                  else
c                      nv=tractsrcs(isrc)%nvert
c                  endif
                  if (isrc .eq. 1) then
                      allocate(polyx(ngrps,nv))
                      allocate(polyy(ngrps,nv))
                      polyx=-999.0
                      polyy=-999.0
                  endif
              endif
              
          endif
          if (index(line,'SO AREAVERT') .gt. 0) then
              call source_info(line,isrc)
          endif
c         seasonal emissions version 17243
          if (nmonths .eq. 12) then
              if (index(line,'** MONTHEMIS') .gt. 0) then
                  if (.not. lmonth) then
                      lmonth=.true.
                      m1=1
                  else
                      lmonth=.false.
                      m1=7
                  endif
                  
                  read(line(14:250),*)srcgrp,(emis(i),i=m1,m1+5)
                  isrc=findsrc(srcgrp)
                  if (isrc .gt. 0) then
                      do i=1,13
                          if (i .le. 12) then
                              memis(isrc,i)=emis(i)
                          else
                              memis(isrc,13)=10000.
                          endif
                      enddo
                  endif
              endif
          endif
      
          if (index(line,'RE STARTING') .gt. 0) lrecstart=.true.
          if (index(line,'RE FINISHED') .gt. 0) lrecend=.true.
          if (lrecstart .and. .not. lrecend) then
              if (index(line,'**') .gt. 0) then
                  irec=irec+1
                  if (index(line,'**COLROW=') .gt. 0) then
                      igrid=igrid+1  !gridded receptor
                      call recinfo(ns,line,igrid,irec)
                  endif
                  if (index(line,'**BLOCK=') .gt. 0) then
                      iblock=iblock+1 !block receptor
                      call recinfo(ns,line,iblock,irec)   
                  endif
                  if (index(line,'**ZEROBLOCK=') .gt. 0) then
                      izblock=izblock+1 !block receptor
                      call recinfo(ns,line,izblock,irec)   
                  endif
                  if (index(line,'**MONITOR=') .gt. 0) then
                      imon=imon+1 !block receptor
                      call recinfo(ns,line,imon,irec)   
                  endif
              endif
          endif
          read(inpunit,10,iostat=eof)line
          goto 20
      endif
      close(inpunit)
      
   
c      do i=1,tractsrcs(isrc)%nvert !JAT
c          write(101,*)polyx(i),polyy(i)
c      enddo
      
c     now read the plt file
c     order of receptors by source group is block,gridded,monitor
      eof=0
      irec=0
      write(logunit,*)'Reading model concentrations'
      if (lsplit) then
          i=1
          lexist=.true.
          do while (lexist)
              write(pltfil,'(5(a),i2.2,a)')trim(adjustl(outdir(idir))),
     +        trim(adjustl(state2(ns))),'/',
     +        trim(adjustl(facid(ns))),'_',i,'.plt'
              inquire(file=pltfil,exist=lexist)
              if (lexist) then
                  write(logunit,*)'reading ',trim(adjustl(pltfil))
                  open(unit=pltunit,file=pltfil,status='old')
                  srcgrp1='A'
                  read(pltunit,10,iostat=eof)line
 40               if (eof .eq. 0) then
                      if (line(1:1) .eq. '*') then
                          irec=0
                          if (nmonths .eq. 12) then
                             imonth1=0
                          else
                             imonth1=1
                          endif
                      else
                          if (nmonths .eq. 12) then
                             read(line,30)x,y,conc,z1,z2,z3,ave,srcgrp,
     +                       iyr,imonth
                          else
                             read(line,25)x,y,conc,z1,z2,z3,ave,srcgrp
                             imonth=1
                          endif
                          if (trim(adjustl(srcgrp)) .ne.
     +                    trim(adjustl(srcgrp1))) then
                             irec=0
                             srcgrp1=srcgrp
                          endif
                          if (imonth .ne. imonth1) then
                             irec=0
                             imonth1=imonth
                          endif
                          irec=irec+1
                          call assignmod(x,y,conc,irec,imonth,srcgrp)
                      endif
                      read(pltunit,10,iostat=eof)line
                      goto 40
                  endif 
                  close(pltunit)
              endif
              i=i+1
          enddo
      else
          open(unit=pltunit,file=pltfil,status='old')
          read(pltunit,10,iostat=eof)line
 35       if (eof .eq. 0) then
              if (line(1:1) .eq. '*') then
                  irec=0
                  if (nmonths .eq. 12) then
                     imonth1=0
                  else
                     imonth1=1
                   endif
              else
                  if (nmonths .eq. 12) then
                      read(line,30)x,y,conc,z1,z2,z3,ave,srcgrp,iyr,
     +                imonth
                  else
                      read(line,25)x,y,conc,z1,z2,z3,ave,srcgrp
                      imonth=1
                  endif
                  if (trim(adjustl(srcgrp)) .ne. 
     +            trim(adjustl(srcgrp1))) then                  
                      irec=0
                      srcgrp1=srcgrp
                  endif 
                  if (imonth .ne. imonth1) then
                      irec=0
                      imonth1=imonth
                  endif
                  irec=irec+1
                  call assignmod(x,y,conc,irec,imonth,srcgrp)
              endif
              read(pltunit,10,iostat=eof)line
              goto 35
          endif 

          close(pltunit)
      endif
     
5     format('SRC_REC AERMOD file ',a,' not found'/'skip')
 10   format(a200)
 25   format(2(1X,F13.5),1X,E13.6,3(1X,F8.2),2X,A6,2X,A8)
c30   format(2(1X,F13.5),1X,E13.6,3(1X,F7.2),2X,A8,2X,2(i2))
 30   format(2(1X,F13.5),1X,E13.6,3(1X,F8.2),2X,A6,2X,A8,2X,2(i2.2))
 95   return
      end
c##########################################################################
      subroutine readerr(efile)
      use main1
      implicit none
      character efile*300,line*100
      integer*8 eof,m,nh
      logical lleap
      
c     determine if year is leap year
      call leapyr(lleap)
      
c     calculate max # of hours per month
      do m=1,12
          if (m .eq. 2 .and. lleap) then
              nhours(m)=(days(m)+1)*24
          else
              nhours(m)=days(m)*24
          endif
      enddo
      
      nh=0
c     read the errors file and subtract 1 from the hourly total for the appropriate month
c     when a calm or missing hour message is encountered
      open(unit=errunit,file=efile,status='old')
      read(errunit,5,iostat=eof)line
 10   if (eof .eq. 0) then
          if (index(line,'MX I440') .gt. 0 .or. 
     +    index(line,'MX 1460') .gt. 0) then
              read(line(88:89),'(i2)')m
              nhours(m)=nhours(m)-1
          endif
          read(errunit,5,iostat=eof)line
          goto 10
      endif
      close(errunit)
      write(logunit,15)
      do m=1,12
        nh=nh+nhours(m)
        write(logunit,20)m,nhours(m)
      enddo
      write(logunit,25)nh
 5    format(a100)
 15   format(/'Number of valid hours per month'/)
 20   format(i2.2,1x,i3)
 25   format(/'Total valid hours:',1x,i4/)
      return
      end
c##########################################################################
c      subroutine leapyr(iyear,leap)
      subroutine leapyr(lleap)
      use main1
      implicit none
      
c iyear:    integer 4-digit year of processed year
c ileap4:   integer result of MOD(iyear,4)
c ileap100: integer result of MOD(iyear,100)
c ileap400: integer result of MOD(iyear,400)
      integer ileap4,ileap100,ileap400
      logical lleap
      
      ileap4=mod(modyear,4)   
      ileap100=mod(modyear,100)
      ileap400=mod(modyear,400)

c if year is not divisible by 4 or year is divisible by 100 but not divisible by 400 then year 
c is not a leap year
c 2001 is not a leap year because it is not divisible by 4
c 2000 is is a leap year because it is divisible by 400      
      if (ileap4 .ne. 0 .or. 
     +  (ileap100 .eq. 0 .and. ileap400 .ne. 0)) then 
        lleap=.false.  
      else
        lleap=.true.  
      endif
      return
      end
c##########################################################################
      subroutine source_info(line,is)
      use main1
      implicit none
      integer*8 is,nn
      real*8 e,ht,j3,x1,y1,x2,y2
      character line*250,src*12,stype*8,j1*2,j2*8
      
      if (index(line,'SO LOCATION') .gt. 0 .or. 
     +index(line,'** LOCATION') .gt. 0) then
          if (runap) then
              read(line,*)j1,j2,src,stype
              if (stype .eq. 'LINE') then
                  read(line,*)j1,j2,apsrcs(is)%srcid,apsrcs(is)%srctyp,
     +            apsrcs(is)%utmx1,apsrcs(is)%utmy1,apsrcs(is)%utmx2,
     +            apsrcs(is)%utmy2
              else
                  read(line,*)j1,j2,apsrcs(is)%srcid,apsrcs(is)%srctyp,
     +            apsrcs(is)%utmx1,apsrcs(is)%utmy1
              endif        
          elseif (runpt) then
              read(line,*)j1,j2,ptsrcs(is)%srcid,ptsrcs(is)%srctyp,
     +        ptsrcs(is)%utmx,ptsrcs(is)%utmy
          elseif (runcmv) then
c             read(line,*)j1,j2,cmvsrcs(is)%srcid,cmvsrcs(is)%srctyp,
              read(line,*)j1,j2,src,cmvsrcs(is)%srctyp,
     +        cmvsrcs(is)%utmx1,cmvsrcs(is)%utmy1
              nn=index(src,'c')
              write(cmvsrcs(is)%srcid,'(a1,a1)')'C',src(nn+1:nn+1)
c          elseif (noncontr) then
c              read(line,*)j1,j2,tractsrcs(is)%srcid,
c     +        tractsrcs(is)%srctyp,tractsrcs(is)%utmx1,
c     +        tractsrcs(is)%utmy1
          else
              read(line,*)j1,j2,gridsrcs(is)%srcid,stype,
     +        gridsrcs(is)%utmx1,gridsrcs(is)%utmy1
          endif
      endif
      if (index(line,'SO SRCPARAM') .gt. 0 .or. 
     +index(line,'** SRCPARAM') .gt. 0) then
          if (runap) then
                  read(line,*)j1,j2,src,e,ht,apsrcs(is)%xdim,j3
                  if (apsrcs(is)%srctyp .eq. 'AREA') apsrcs(is)%ydim=j3
          elseif (runpt) then
              if(ptsrcs(is)%srctyp .eq. 'AREA')
     +        read(line,*)j1,j2,src,e,ht,ptsrcs(is)%xdim,ptsrcs(is)%ydim
          elseif (runcmv) then
              read(line,*)j1,j2,src,e,ht,cmvsrcs(is)%nvert
c          elseif (noncontr) then
c              read(line,*)j1,j2,src,e,ht,tractsrcs(is)%nvert
          else
              read(line,*)j1,j2,src,e,ht,gridsrcs(is)%nvert
          endif
      endif
      if (index(line,'SO AREAVERT') .gt. 0) then
c          if ((runon .or. runnon .or. runrwc .or. runnplow 
c     +    .or. runnphi .or. runog) .and. .not. noncontr) then
          if (.not. runcmv) then
              read(line,*)j1,j2,src,x1,y1,x2,y2
              if (gridsrcs(is)%utmx2 .eq. -999.0) then
c                  gridsrcs(is)%utmx1=x1
c                  gridsrcs(is)%utmy1=y1
                  gridsrcs(is)%utmx2=x2
                  gridsrcs(is)%utmy2=y2
              else
                  gridsrcs(is)%utmx3=x1
                  gridsrcs(is)%utmy3=y1
                  gridsrcs(is)%utmx4=x2
                  gridsrcs(is)%utmy4=y2
              endif
          else
              call get_vertices(line,is)
          endif
          
      endif
c      write(107,*)ptsrcs(is)%srcid,ptsrcs(is)%srctyp,ptsrcs(is)%utmx,
c     + ptsrcs(is)%utmy
      return
      end
c##########################################################################
c     subroutine to fill in CMV vertices array
      subroutine get_vertices(line,is)
      use main1
      implicit none
      integer*8 iv,i,ilen,j,j1,j2,v,nv,is,nvmaster
      character line*250,line1*250
      logical lmiss
      
      real*8, allocatable, dimension(:) :: x
      real*8, allocatable, dimension(:) :: y
c     check to see where the vertices array has been filled in

      
      v=1
      lmiss=.false.
      if (runcmv) then
          nvmaster=cmvsrcs(is)%nvert
      else
          nvmaster=tractsrcs(is)%nvert
      endif
      
      do while (v .le. nvmaster .and. .not. lmiss)
          if (polyx(is,v) .eq. -999.0) then
              lmiss=.true.
          else
              v=v+1
          endif
      enddo
      
c     subset line to part of line past source ID
      line1=line(25:250)
      ilen=len_trim(line1)
c     loop through line1 to see how many pairs of numbers there will be on the line
c      most is 3 pairs
      iv=0
      do i=1,ilen
          j=ichar(line1(i:i))
          if (j .eq. 32 .and. i .gt. 1) then  !ignore if first character is blank
              j1=ichar(line1(i-1:i-1)) 
              j2=ichar(line1(i+1:i+1)) 
              if (j1 .ge. 48 .and. j1 .le. 57 .and. j2 .ge. 48 .and. 
     +        j2 .le. 57)iv=iv+1
          endif
      enddo
      
      allocate(x(iv))
      allocate(y(iv))
      
      nv=iv
c     now fill in array
      read(line1,*)(x(iv),y(iv),iv=1,nv)
      do iv=1,nv
          polyx(is,v+iv-1)=x(iv)
          polyy(is,v+iv-1)=y(iv)
      enddo
      
      deallocate(x)
      deallocate(y)
      return
      end
      
c##########################################################################
c     subroutine to fill in receptor arrays
      subroutine recinfo(is,line,i,j)
      use main1
      implicit none
      
      integer*8 i,j,ncr,nx,ny,nblk,nmon,cr,nf,nr,iflag,irec,
     +findrec,recflag,nzblk,is,igrid
      real*8 x,y,ux,uy,dist,distance,zerod,maxd
      character line*200,recid*15
      logical lwrite
      
      zerod=0.0
      nx=index(line,'CMAQX=')
      ny=index(line,'CMAQY=')
      ncr=index(line,'COLROW=')
      nblk=index(line,'**BLOCK=')
      nzblk=index(line,'**ZEROBLOCK=')
      nmon=index(line,'**MONITOR')
      nr=index(line,'REC=')
      nf=index(line,'FLAG =')
      
      read(line(ncr+7:nx-1),*)cr
      read(line(nx+6:ny-1),*)x
      read(line(ny+6:200),*)y
      read(line(nf+7:nf+7),*)iflag
      read(inpunit,5)ux,uy
c      modux(j)=ux
c      moduy(j)=uy
      
      if (allocated(state2)) then
        if (state2(is) .eq. '02') then
            igrid=2
        elseif (state2(is) .eq. '15') then
            igrid=3
        elseif (state2(is) .eq. '72' .or. state2(is) .eq. '78') then
          igrid=4
        else
           igrid=1
        endif
      else
         igrid=1
      endif
      if (nblk .gt. 0 .and. nmodb .gt. 0) then  !block
          recid=line(nblk+9:nblk+23)
          read(line(ncr+7:ncr+14),*)cr
          modblock(i)%recid=recid
          modblock(i)%colrow=cr
          modblock(i)%ux=ux
          modblock(i)%uy=uy
          modblock(i)%x=x
          modblock(i)%y=y
          modblock(i)%iflag=iflag
          modblock(i)%mastrec=findrec('BM',i)
c         version 20002
          if (cr .le. 0) then
              maxd=grddist
          else
              maxd=blkdist
          endif
          
c          if (.not. akhiprvi) then
c              maxd=blkdist
c          else
c              maxd=grddist
c          endif
          call recdist('mbl',recid,ux,uy,zerod,maxd,lwrite,recflag,
     +    dist)  !get distance
          modblock(i)%dist=dist
          if (.not. lwrite) modblock(i)%iflag=3  !fix for bug in aermod input program
c           write(108,*)modblock(i)%recid,modblock(i)%x,modblock(i)%y,
c     +modblock(i)%dist
          
      elseif (nzblk .gt. 0 .and. nmodz .gt. 0) then  !zero pop block
          recid=line(nzblk+13:nzblk+27)
          read(line(ncr+7:ncr+14),*)cr
          zmodblock(i)%recid=recid
          zmodblock(i)%colrow=cr
          zmodblock(i)%ux=ux
          zmodblock(i)%uy=uy
          zmodblock(i)%x=x
          zmodblock(i)%y=y
          zmodblock(i)%iflag=iflag
          zmodblock(i)%mastrec=findrec('ZB',i)
          if (.not. akhiprvi) then
              maxd=blkdist
          else
              maxd=grddist
          endif
          call recdist('zbl',recid,ux,uy,zerod,maxd,lwrite,recflag,
     +    dist)  !get distance
          zmodblock(i)%dist=dist
          if (.not. lwrite) zmodblock(i)%iflag=3  !fix for bug in aermod input program
          
      elseif (nmon .gt. 0 .and. nmodm .gt. 0) then !monitor
c         recid=line(nmon+11:nmon+19)
          recid=line(nmon+11:ncr-1)
          read(line(ncr+7:ncr+14),*)cr
          modmon(i)%recid=recid
          modmon(i)%colrow=cr
          modmon(i)%ux=ux
          modmon(i)%uy=uy
          modmon(i)%x=x
          modmon(i)%y=y
          modmon(i)%iflag=iflag
          modmon(i)%mastrec=findrec('MM',i)
c         version 20002
          maxd=blkdist
          
c          if (.not. akhiprvi) then
c              maxd=blkdist
c          else
c              maxd=grddist
c          endif
          call recdist('mmm',recid,ux,uy,zerod,maxd,lwrite,recflag,
     +    dist)  !get distance
          modmon(i)%dist=dist
          if (.not. lwrite) modmon(i)%iflag=3 !fix for bug in aermod input program
c          write(110,*)modmon(i)%recid,modmon(i)%x,modmon(i)%y,
c     +modmon(i)%dist
          
      else
          read(line(ncr+7:ncr+13),*)cr
          read(line(nr+4:nr+7),*)irec
          write(recid,'(a4,i6.6,i4.4)')gridst(igrid),cr,irec
          modgrid(i)%recid=recid
          modgrid(i)%colrow=cr
          modgrid(i)%ux=ux
          modgrid(i)%uy=uy
          modgrid(i)%x=x
          modgrid(i)%y=y
          modgrid(i)%iflag=iflag
          modgrid(i)%mastrec=findrec('GM',i)
          call recdist('mgr',recid,ux,uy,zerod,grddist,lwrite,recflag,
     +    dist)  !get distance
          modgrid(i)%dist=dist
          if (.not. lwrite) then
              modgrid(i)%iflag=3
          else !account for being close to one source but not another
              if(modgrid(i)%iflag .eq. 3) modgrid(i)%iflag=0
          endif
c         write(109,*)i,modgrid(i)%recid,modgrid(i)%x,modgrid(i)%y,
c    +modgrid(i)%dist,modgrid(i)%iflag,modgrid(i)%ux,modgrid(i)%uy,
c    +lwrite
      endif

5     format(11x,f13.3,1x,f13.3)
      return
      end          
c##########################################################################
c     subset master receptor list to within plus/minus 5 cells for faster processing
c     when getting master receptor numbers
      subroutine subsetrecs(ns) 
      use main1
      implicit none
      integer*8 i,c,r,col1,col2,row1,row2,modcol,modrow,ig,ns,get_state,
     +ig1
      character astate*2,st*2
      if (akhiprvi) then
          if (noncontr) then
              i=get_state(ncst)
              astate=states1(i)
          else
             if (state2(ns) .eq. 'AK' .or. state2(ns) .eq. '02') then
                astate='02'
               
             elseif (state2(ns) .eq. 'HI' .or. state2(ns) .eq. '15') 
     +       then
                astate='15'
               
             elseif (state2(ns) .eq. 'PR' .or. state2(ns) .eq. '78') 
     +       then
                astate='72'
                
             else
                astate='78'
               
              endif
          endif
          if (astate .eq. '02') then
              modcol=7
              modrow=7
              ig1=2
          elseif (astate .eq. '15') THEN
              modcol=20
              modrow=20
              ig1=3
          else
              modcol=20
              modrow=20
              ig1=4
          endif
          
          col1=col-modcol
          col2=col+modcol
          row1=row-modrow
          row2=row+modrow
      else
          astate='00'
          ig1=1
          if (.not. runcmv) then
c              modcol=5
c              modrow=5
              modcol=8
              modrow=8
          else
              modcol=40
              modrow=40
          endif
      endif

      col1=col-modcol
      col2=col+modcol
      row1=row-modrow
      row2=row+modrow
      
      ngrids1=0
c     gridded
      do i=1,ngrids
          if (gridrec(i)%col .ge. col1 .and. gridrec(i)%col .le. 
     +    col2 .and. gridrec(i)%row .ge. row1 .and. 
     +    gridrec(i)%row .le. row2 .and. gridrec(i)%recid(1:4) .eq. 
     +    gridst(ig1)) ngrids1=ngrids1+1
      enddo
      if (ngrids1 .gt. 0) then
          allocate(gridrec1(ngrids1))
          ig=0
          do i=1,ngrids
              if (gridrec(i)%col .ge. col1 .and. gridrec(i)%col .le. 
     +        col2 .and. gridrec(i)%row .ge. row1 .and. 
     +        gridrec(i)%row .le. row2 .and. gridrec(i)%recid(1:4) 
     +        .eq. gridst(ig1)) then
                  ig=ig+1
                  gridrec1(ig)%recid=gridrec(i)%recid
                  gridrec1(ig)%irec=gridrec(i)%irec1
              endif
          enddo
      endif

c     blocks
      nblocks1=0
      do i=1,nblocks
          if (.not. akhiprvi) then
              if (blockrec(i)%col .ge. col1 .and. blockrec(i)%col .le. 
     +        col2 .and. blockrec(i)%row .ge. row1 .and. 
     +        blockrec(i)%row .le. row2) nblocks1=nblocks1+1
          else
              st=blockrec(i)%recid(1:2)
c             account for fact that some receptors in VI/PR will overlap in 50 km domain
              if (((st .eq. '72' .or. st .eq. '78') .and. 
     +        (astate .eq. '72' .or. astate .eq. '78')) .or. 
     +        ((st .eq. '02' .or. st .eq. '15') .and. 
     +        st .eq. astate) .and. ((blockrec(i)%col .ge. col1 .and.  
     +        blockrec(i)%col .le.col2 .and. blockrec(i)%row .ge. row1  
     +        .and.blockrec(i)%row .le. row2) .or. 
     +        (blockrec(i)%col .eq. 0 .and. blockrec(i)%row .eq. 0))) 
     +        then
                  nblocks1=nblocks1+1
              endif
c              if (st .eq. astate) nblocks1=nblocks1+1
          endif
      enddo
      if (nblocks1 .gt. 0) then
          allocate(blockrec1(nblocks1))
          ig=0
          do i=1,nblocks
              if (.not. akhiprvi) then
                  if (blockrec(i)%col .ge. col1 .and. blockrec(i)%col 
     +            .le. col2 .and. blockrec(i)%row .ge. row1 .and. 
     +            blockrec(i)%row .le. row2) then
                      ig=ig+1
                      blockrec1(ig)%recid=blockrec(i)%recid
                      blockrec1(ig)%irec=blockrec(i)%irec1
                  endif
              else
                  st=blockrec(i)%recid(1:2)
                  if (((st .eq. '72' .or. st .eq. '78') .and. 
     +            (astate .eq. '72' .or. astate .eq. '78')) .or. 
     +            ((st .eq. '02' .or. st .eq. '15') .and. 
     +            st .eq. astate) .and. ((blockrec(i)%col .ge. col1  
     +            .and. blockrec(i)%col .le.col2 .and. blockrec(i)%row   
     +            .ge. row1 .and.blockrec(i)%row .le. row2) .or. 
     +            (blockrec(i)%col .eq. 0 .and. blockrec(i)%row 
     +            .eq. 0))) then
                      ig=ig+1
                      blockrec1(ig)%recid=blockrec(i)%recid
                      blockrec1(ig)%irec=blockrec(i)%irec1
                  endif
c                  if (st .eq. astate) then
c                      ig=ig+1
c                      blockrec1(ig)%recid=blockrec(i)%recid
c                      blockrec1(ig)%irec=blockrec(i)%irec
c                  endif
              endif
          enddo
      endif   
c     nonpopulated blocks  
      nzblocks1=0
      do i=1,nzblocks
          if (.not. akhiprvi) then
              if (zblockrec(i)%col .ge. col1 .and. zblockrec(i)%col .le. 
     +        col2 .and. blockrec(i)%row .ge. row1 .and. 
     +        zblockrec(i)%row .le. row2) nzblocks1=nzblocks1+1
          else
              st=zblockrec(i)%recid(1:2)
c             account for fact that some receptors in VI/PR will overlap in 50 km domain
              if (st .eq. '02' .and. st .eq. astate) then
                  nzblocks1=nzblocks1+1
              endif
c              if (st .eq. astate) nzblocks1=nzblocks1+1
          endif
      enddo
      if (nzblocks1 .gt. 0) then
c        write(*,*)'nzblocks1 ',nzblocks1
          allocate(zblockrec1(nzblocks1))
          ig=0
          do i=1,nzblocks
              if (.not. akhiprvi) then
                  if (zblockrec(i)%col .ge. col1 .and. zblockrec(i)%col 
     +            .le. col2 .and. zblockrec(i)%row .ge. row1 .and. 
     +            zblockrec(i)%row .le. row2) then
                      ig=ig+1
                      zblockrec1(ig)%recid=zblockrec(i)%recid
                      zblockrec1(ig)%irec=zblockrec(i)%irec1
                  endif
              else
                  st=zblockrec(i)%recid(1:2)
                  if (st .eq. '02' .and. st .eq. astate) then
                      ig=ig+1
                      zblockrec1(ig)%recid=zblockrec(i)%recid
                      zblockrec1(ig)%irec=zblockrec(i)%irec1
                  endif
c                  if (st .eq. astate) then
c                      ig=ig+1
c                      blockrec1(ig)%recid=blockrec(i)%recid
c                      blockrec1(ig)%irec=blockrec(i)%irec
c                  endif
              endif
          enddo
      endif        
c     monitors
      nmons1=0
      do i=1,nmons
          if (.not. akhiprvi) then
              if (monrec(i)%col .ge. col1 .and. monrec(i)%col .le. 
     +        col2 .and. monrec(i)%row .ge. row1 .and. 
     +        monrec(i)%row .le. row2) nmons1=nmons1+1
          else
              st=monrec(i)%fip(1:2)
c             account for fact that some receptors in VI/PR will overlap in 50 km domain
              if (((st .eq. '72' .or. st .eq. '78') .and. 
     +        (astate .eq. '72' .or. astate .eq. '78')) .or. 
     +        ((st .eq. '02' .or. st .eq. '15') .and. 
     +        st .eq. astate) .and. ((monrec(i)%col .ge. col1 .and.  
     +        monrec(i)%col .le.col2 .and. monrec(i)%row .ge. row1  
     +        .and. monrec(i)%row .le. row2) .or. 
     +        (monrec(i)%col .eq. 0 .and. monrec(i)%row .eq. 0))) then
                  nmons1=nmons1+1
              endif
c             if (st .eq. astate) nmons1=nmons1+1
          endif
      enddo
      if (nmons1 .gt. 0) then
          allocate(monrec1(nmons1))
          ig=0
          do i=1,nmons
              if (.not. akhiprvi) then
                  if (monrec(i)%col .ge. col1 .and. monrec(i)%col .le. 
     +            col2 .and. monrec(i)%row .ge. row1 .and. 
     +            monrec(i)%row .le. row2) then
                      ig=ig+1
                      monrec1(ig)%recid=monrec(i)%recid
                      monrec1(ig)%irec=monrec(i)%irec1
                  endif
              else
                  st=monrec(i)%fip(1:2)
                  if (((st .eq. '72' .or. st .eq. '78') .and. 
     +            (astate .eq. '72' .or. astate .eq. '78')) .or. 
     +            ((st .eq. '02' .or. st .eq. '15') .and. 
     +            st .eq. astate) .and. ((monrec(i)%col .ge. col1 .and.  
     +            monrec(i)%col .le.col2 .and. monrec(i)%row .ge. row1  
     +            .and. monrec(i)%row .le. row2) .or. (monrec(i)%col
     +             .eq. 0 .and. monrec(i)%row .eq. 0))) then
                      ig=ig+1
                      monrec1(ig)%recid=monrec(i)%recid
                      monrec1(ig)%irec=monrec(i)%irec1
                  endif
c                  if (st .eq. astate) then
c                      ig=ig+1
c                      monrec1(ig)%recid=monrec(i)%recid
c                      monrec1(ig)%irec=monrec(i)%irec
c                  endif
              endif
          enddo
      endif          
      return
      end
c##########################################################################
c     subroutine to assign modeled concentrations to arrays
      subroutine assignmod(x,y,conc,irec,imonth,grp) 
      use main1
      implicit none
      integer*8 irec,irec1,igrp,findsrc,imonth
      real*8 x,y,conc,dist,distance
      character grp*12 !grp*8
     
      if (runap) then !.or. runcmv) then !.or. noncontr) then
          igrp=1
      else
          igrp=findsrc(grp)
      endif

c     blocks, gridded, and monitors for point, airports, cmv, and tract sources except onroad tracts
      if (runpt .or. runap .or. runcmv .or. 
     +(noncontr .and. .not. runon .and. .not. runnon)) then   
          if (irec .le. nmodb) then !blocks 
              irec1=irec
              dist=distance(x,y,modblock(irec)%ux,modblock(irec)%uy)
              if (dist .le. 1.0) modbc(igrp,irec1,1)=conc
          else
              if (irec .le. nmodb+nmodz) then !nonpopulated blocks
                  irec1=irec-nmodb
                  dist=distance(x,y,zmodblock(irec1)%ux,
     +            zmodblock(irec1)%uy)
                  if (dist .le. 1.0) modzbc(igrp,irec1,1)=conc
              elseif (irec .le. nmodb+nmodg+nmodz) then  !gridded
                  irec1=irec-nmodb
                  dist=distance(x,y,modgrid(irec1)%ux,modgrid(irec1)%uy)  !should be the same receptor
                  if (dist .le. 1.0) modgc(igrp,irec1,1)=conc
              else   !monitor
                  irec1=irec-nmodb-nmodg-nmodz
                  dist=distance(x,y,modmon(irec1)%ux,modmon(irec1)%uy)
                  modmc(igrp,irec1,1)=conc
              endif
          endif
      endif
     
c     blocks/monitors for onroad non-conus
      if ((runon .or. runnon) .and. noncontr) then  !grid source onroad/nonroad
c          irec1=irec-(nmodb+nmodm+nmodz)*((iseas-1)*24+(ihr-1))
c         irec1=irec-(nmodb+nmodm+nmodz)*(imonth-1)
          irec1=irec
          if (irec1 .le. nmodb) then   !block
              dist=distance(x,y,modblock(irec1)%ux,modblock(irec1)%uy)  !should be the same receptor
              if (dist .le. 1.0) then
                  modbc(igrp,irec1,imonth)=modbc(igrp,irec1,imonth)+
     +            conc*real(nhours(imonth))
c                  if (igrp .eq. 1 .and. irec1 .eq. nmodb) 
c     +                seashrs(iseas,ihr)=nhrs
              endif
          elseif (irec1 .le. nmodb+nmodz) then   !non-pop block
              irec1=irec1-nmodb
c             write(*,*)irec,irec1,imonth,nmodb,nmodz
              dist=distance(x,y,zmodblock(irec1)%ux,zmodblock(irec1)%uy)  !should be the same receptor
              if (dist .le. 1.0) then
                  modzbc(igrp,irec1,imonth)=modzbc(igrp,irec1,imonth)+
     +            conc*real(nhours(imonth))
c                  if (igrp .eq. 1 .and. irec1 .eq. nmodb) 
c     +                seashrs(iseas,ihr)=nhrs
              endif    
c          write(115,*)modblock(irec1)%ux,modblock(irec1)%uy,
c     +modbc(igrp,irec1,iseas),iseas
          else  !monitor
c              write(135,*)irec,irec1,nmodb,nmodm,x,y,iseas,ihr
              irec1=irec1-nmodb-nmodz
              dist=distance(x,y,modmon(irec1)%ux,modmon(irec1)%uy)  !should be the same receptor
              if (dist .le. 1.0) then
                  modmc(igrp,irec1,imonth)=modmc(igrp,irec1,imonth)+
     +            conc*real(nhours(imonth))
c                  if (igrp .eq. 1 .and. irec1 .eq. nmodm) 
c     +                seashrs(iseas,ihr)=nhrs
              endif
c              write(117,*)modmon(irec1)%ux,modmon(irec1)%uy,
c     +modmc(igrp,irec1,iseas),iseas,igrp,nhrs,seashrs(iseas,ihr)
          endif
      endif
     
c     gridded sources except onroad and nonroad, no blocks and no monitors
      if ((runon .or. runnon) .and. .not. noncontr) then
c          irec1=irec-(nmodg)*((iseas-1)*24+(ihr-1))
c         irec1=irec-nmodg*(imonth-1)
          irec1=irec
          dist=distance(x,y,modgrid(irec1)%ux,modgrid(irec1)%uy)  !should be the same receptor
          if (dist .le. 1.0) then
          modgc(igrp,irec1,imonth)=modgc(igrp,irec1,imonth)+
     +    conc*real(nhours(imonth))
c          if (igrp .eq. 1 .and. irec1 .eq. nmodg) 
c     +        seashrs(iseas,ihr)=nhrs
          endif
      endif 
      
c     all other gridded sources, no blocks and no monitors
      if (.not. noncontr .and. (runrwc .or. runnphi .or. runnplow .or. 
     +runog .or. runag)) then
          irec1=irec
          dist=distance(x,y,modgrid(irec1)%ux,modgrid(irec1)%uy)  !should be the same receptor
          if (dist .le. 1.0) modgc(igrp,irec1,1)=conc
      endif
      
      return
      end
      
c##########################################################################
c     subroutine to get apppropriate receptors for interpolation
      subroutine int_recs(ns) 
      use main1
      implicit none
      
c     eof:        end of file indicator
c     col1:       western most column to search for receptors
c     row1:       southern most row to search for receptors
c     col2:       eastern most column to search for receptors
c     row2:       northern most row to search for receptors   
c     nblock:     number of block receptors written to AERMOD.INP
c     ngrid:      number of gridded receptors written to AERMOD.INP
c     nmon:       number of monitor receptors written to AERMOD.INP
c     c:          column of receptor
c     r:          row of receptor
c     colrow:     column/row of receptor
c     modcol:     number of columns to search on either side of facility column
c     modrow:     number of rows to search on either side of facility row
c     recflag:    flag indicating if receptor is within 30 m of a source
c     irec:       gridded receptor number (1-9 for 4 km, 1-44 for 1 km)
      integer*8 eof,col1,col2,row1,row2,c,r,colrow,modcol,modrow,
     +recflag,ns

c     x:          receptor UTM-x coordinate
c     y:          receptor UTM-y coordinate
c     cmaqx:      receptor CMAQ x-coordinate
c     cmaqy:      receptor CMAQ-Y coordinate

c     relev:      receptor elevation
c     hill:       receptor hill height
c     mindist:      minimum distance to check
      real*8 x,y,cmaqx,cmaqy,relev,hill,dist,mindist
      
c     recfil:     receptor filename
c     recid:      census block ID
c     monid:      monitor ID       
      character recfil*250,recid*15,monid*15,tfil*250,ares*2,arun*10,
     +ahap*3,runtyp1*20
      
c     lexist:     logical variable denoting receptor file exists
c     lwrite:     logical varible denoting to write receptor to AERMOD.INP
      logical lexist,lwrite
      
c     version 20002, now interpolate for all areas
c     initialize variables
      nintb=0
      nintm=0
      if (.not. runcmv .and. (runpt .or. runap)) then
c          modcol=5
c          modrow=5
          if (state2(ns) .eq. '02') then !9km
              modcol=7 !63 km
              modrow=7
          elseif (state2(ns) .eq. '15' .or. state2(ns) .eq. '72' .or.  
     +    state2(ns) .eq. '78') then !3km
              modcol=20 !60 km
              modrow=20
          else
              modcol=8
              modrow=8
          endif
      elseif (runcmv) then !account for large shapes covering several grid cells.
          modcol=40
          modrow=40
      else
         if (ncst .eq. 'AK') then
           modcol=7
           modrow=7
         elseif(ncst .eq. 'NA') then
           modcol=8
           modrow=8
         else
           modcol=20
           modrow=20
         endif
      endif
      col1=col-modcol
      col2=col+modcol
      row1=row-modrow
      row2=row+modrow
      
      if (runpt .or. runap .or. runcmv) then
          mindist=blkdist
      else
          mindist=0.0
      endif
  
      write(arun,'(i10)')runno !run #
      if (nhap1 .gt. 0) then
          write(ahap,'(i3.3)')nhap1
      else
          ahap='001'
      endif
      
      if (lsingfac) then
          write(runtyp1,'(3(a))')trim(adjustl(runtyp)),'_',
     +    trim(adjustl(singlfac))
      elseif (lsingst) then
          write(runtyp1,'(3(a))')trim(adjustl(runtyp)),'_',
     +    trim(adjustl(singlst))
      else
          runtyp1=trim(adjustl(runtyp))
      endif
      
c     write blocks within specified distances (10-50 km) for point, airports, and CMV
c     for gridded, all receptors within 50 km.
c     flag receptors within 30 m of a facility release point, 30 m within runway (entire length),
c     or side of a CMV polygon)
      if (col .ge. 0 .and. row .ge. 0) then
c         temporary file for getting receptors to interpolate
c          write(tfil,'(3(a))')'temp_recs_',trim(adjustl(singlhap)),
c     +    '.txt'
          if (runcmv .or. runpt .or. runap) then
              write(tfil,'(8(a))')'temp_recs_',trim(adjustl(singlhap)),
     +        '_',trim(adjustl(runtyp1)),trim(adjustl(arun)),'_',
     +        trim(adjustl(ahap)),'.txt'
          else
              write(ares,'(i2)')gridres
              write(tfil,'(10(a))')'temp_recs_',trim(adjustl(singlhap)),
     +        '_',trim(adjustl(runtyp1)),'_',trim(adjustl(ares)),
     +        trim(adjustl(arun)),'_',trim(adjustl(ahap)),'.txt'
          endif
          open(unit=98,file=tfil,status='unknown')
          write(logunit,*)'reading blocks'
          if (utmzone .ge. 10) then
              write(recfil,'(2a,i2,a)')trim(adjustl(recdir)),
     +        'blocks_utm_',utmzone,'.dat'
          else
              write(recfil,'(2a,i1,a)')trim(adjustl(recdir)),
     +        'blocks_utm_',utmzone,'.dat'
          endif
          lexist=.false.
          inquire(file=recfil,exist=lexist)
          if (.not. lexist) then
              write(logunit,35)trim(adjustl(recfil))
              stop
          endif
          open(unit=recunit,file=recfil,status='old')
          read(recunit,*,iostat=eof)recid,colrow,x,y,cmaqx,cmaqy,relev,
     +    hill
    5     if (eof .eq. 0) then
              c=int(real(colrow)/1000.0) !get column
              r=colrow-c*1000 !get row
              if ((col1 .le. c .and. c .le. col2 .and. row1 .le. r .and. 
     +        r .le. row2) .or. colrow .le. 0) then  !if within the 5x5 block of cells
c                 get distance and determine whether to write out
                  call recdist('blk',recid,x,y,mindist,grddist,lwrite,
     +            recflag,dist)
                  if (lwrite) then
                      write(98,20)recid,colrow,cmaqx,cmaqy,recflag,
     +                x,y,relev,hill,dist
                      nintb=nintb+1
                  endif
              endif
              read(recunit,*,iostat=eof)recid,colrow,x,y,cmaqx,cmaqy,
     +        relev,hill
              goto 5
          endif
      
          close(recunit)
      
c     monitors
          write(logunit,*)'reading monitors'
          if (utmzone .ge. 10) then
              write(recfil,'(2a,i2,a)')trim(adjustl(recdir)),
     +        'monitors_utm_',utmzone,'.dat'
          else
              write(recfil,'(2a,i1,a)')trim(adjustl(recdir)),
     +        'monitors_utm_',utmzone,'.dat'
          endif
          lexist=.false.
          inquire(file=recfil,exist=lexist)
          if (.not. lexist) then
              write(logunit,35)trim(adjustl(recfil))
!                 stop
          endif
          if (lexist) then
              open(unit=recunit,file=recfil,status='old')
              read(recunit,*,iostat=eof)monid,colrow,x,y,cmaqx,cmaqy,
     +        relev,hill
 45           if (eof .eq. 0) then
                  c=int(real(colrow)/1000.0) !get column
                  r=colrow-c*1000 !get row
                  if ((col1 .le. c .and. c .le. col2 .and. row1 .le. r 
     +            .and. r .le. row2) .or. colrow .le. 0) then  !if within the 5x5 block of cells
c                 get distance and determine whether to write out
                      call recdist('mon',monid,x,y,mindist,grddist,
     +                lwrite,recflag,dist)
                      if (lwrite) then
                          write(98,25)monid,colrow,cmaqx,cmaqy,
     +                    recflag,x,y,relev,hill,dist
                          nintm=nintm+1
                      endif
                  endif
                  read(recunit,*,iostat=eof)monid,colrow,x,y,cmaqx,
     +            cmaqy,relev,hill
                  goto 45
              endif
          endif
          close(recunit)
          close(98)
      
c          write(*,*)nintb,nintm
c      pause
          
!      write(aerunit,40)nblock,ngrid,nmon
          if (nintb .eq. 0 .and. nintm .eq. 0) then
              write(logunit,50)
          else
              call int_recs2
          endif
      else
          write(logunit,51)
      endif
20    format('**BLOCK= ',a,1x,'COLROW= ',i6,1x,'CMAQX= ',f11.0,1x,
     +'CMAQY= ',f11.0,' FLAG= ',i1,1x,f12.3,1x,f13.3,
     +2(1x,f9.3),1x,f11.5) 
25    format('**MONITOR= ',a,1x,'COLROW= ',i6,1x,'CMAQX= ',f11.0,1x,
     +'CMAQY= ',f11.0,' FLAG= ',i1,1x,f12.3,1x,f13.3,
     +2(1x,f9.3),1x,f11.5) 
    
            
 35   format('receptor file ',a,' does not exist, stopping program')

 50   format('WARNING NO RECEPTORS FOR SOURCE, SKIP INTERPOLATION!!!')
 51   format('SOURCE IS IN AK, HI, PR, VI, SKIP INTERPOLATION')
      return
      end
c##########################################################################
      subroutine int_recs2
      use main1
      implicit none
      character j1*10,j1a*10,recid*15,rt*10,tfil*250,ares*2,
     +arun*10,ahap*3,runtyp1*20
      integer*8 eof,cr,iflag,b,m,findrec,ns
      real*8 x,y,ux,uy,relev,hill,dist
      
c      write(*,*)'enter int_recs2'
      b=0
      m=0
c      write(*,*)'you are here1'
      if (nintb .gt. 0) then
          allocate(intbc(ngrps1,nintb,nmonths))
          allocate(intblock(nintb))
          intblock%recid='NA'
          intblock%x=0.0
          intblock%y=0.0
          intblock%ux=0.0
          intblock%uy=0.0
          intblock%colrow=0
          intblock%mastrec=0
          intblock%dist=-9.0
          intbc=0.0
      endif
      if (nintm .gt. 0) then
          allocate(intmc(ngrps1,nintm,nmonths))
          allocate(intmon(nintm))
          intmon%recid='NA'
          intmon%x=0.0
          intmon%y=0.0
          intmon%ux=0.0
          intmon%uy=0.0
          intmon%colrow=0
          intmon%mastrec=0
          intmon%dist=-9.0
          intmc=0.0
      endif
      write(arun,'(i10)')runno !run #
      if (nhap1 .gt. 0) then
          write(ahap,'(i3.3)')nhap1
      else
          ahap='001'
      endif
      
      if (lsingfac) then
          write(runtyp1,'(3(a))')trim(adjustl(runtyp)),'_',
     +    trim(adjustl(singlfac))
      elseif (lsingst) then
          write(runtyp1,'(3(a))')trim(adjustl(runtyp)),'_',
     +    trim(adjustl(singlst))
      else
          runtyp1=trim(adjustl(runtyp))
      endif
c      write(*,*)'you are here2'
c      write(tfil,'(3(a))')'temp_recs_',trim(adjustl(singlhap)),'.txt'
      if (runcmv .or. runpt .or. runap) then
              write(tfil,'(8(a))')'temp_recs_',trim(adjustl(singlhap)),
     +        '_',trim(adjustl(runtyp1)),trim(adjustl(arun)),'_',
     +        trim(adjustl(ahap)),'.txt'
          else
              write(ares,'(i2)')gridres
              write(tfil,'(10(a))')'temp_recs_',trim(adjustl(singlhap)),
     +        '_',trim(adjustl(runtyp1)),'_',trim(adjustl(ares)),
     +        trim(adjustl(arun)),'_',trim(adjustl(ahap)),'.txt'
      endif
      
      
      open(unit=98,file=tfil,status='unknown')
      eof=0
      read(98,*,iostat=eof)rt,recid,j1,cr,j1a,x,j1a,y,j1a,iflag,ux,uy,
     +relev,hill,dist
 10   if (eof .eq. 0) then
          if (index(rt,'BLOCK') .gt. 0) then
              b=b+1
              intblock(b)%recid=recid
              intblock(b)%x=x
              intblock(b)%y=y
              intblock(b)%ux=ux
              intblock(b)%uy=uy
              intblock(b)%colrow=cr
              intblock(b)%dist=dist
c             write(*,*)'b ',b
              intblock(b)%mastrec=findrec('BI',b)
c              write(111,*)intblock(b)%recid,intblock(b)%x,intblock(b)%y
          endif
          if (index(rt,'MONITOR') .gt. 0) then
              m=m+1
              intmon(m)%recid=recid
              intmon(m)%x=x
              intmon(m)%y=y
              intmon(m)%ux=ux
              intmon(m)%uy=uy
              intmon(m)%colrow=cr
              intmon(m)%dist=dist
c             write(*,*)'m ',m
              intmon(m)%mastrec=findrec('MI',m)
c              write(112,*)intmon(m)%recid,intmon(m)%x,intmon(m)%y
          endif
          read(98,*,iostat=eof)rt,recid,j1,cr,j1a,x,j1a,y,j1a,iflag,ux,
     +    uy,relev,hill,dist
          goto 10
      endif
      close(98,status='delete')
c      close(98)
c      write(*,*)'exit int_recs2'
      return
c 5    format(3(a),1x,i6,2(1x,a,f11.0),1x,a,1x,i1,1x,f12.3,1x,f13.3)
      end
c##########################################################################
c     subroutine to interpolate to blocks/monitors
      subroutine interp(ns)
      use main1
      implicit none
c     grecs:      array indices of 4 closest gridded receptors
      integer*8 ir,ig,grecs(4),is,ns
      real*8 ic
      character ofile*250
      write(logunit,*)'interpolating'
      call datetime
      
      if (writefac) then
          write(ofile,'(5(a))')trim(adjustl(hapdir)),
     +    trim(adjustl(singlhap)),'_',trim(adjustl(facid(ns))),
     +    '_interp.txt'
          open(unit=intunit,file=ofile,status='replace')
      endif
      
c     blocks
      write(logunit,*)'number of interpolated blocks ',nintb
      if (nintb .gt. 0) then
          do ir=1,nintb
c             find 4 closest gridded receptors
              call find4(intblock(ir)%x,intblock(ir)%y,grecs)
c             calculate terms of interpolation for each group
              do ig=1,ngrps1
                  do is=1,nmonths
                      call calcint(intblock(ir)%x,intblock(ir)%y,
     +                intblock(ir)%recid,ig,is,grecs,ic)
                      intbc(ig,ir,is)=ic
                  enddo
              enddo
          enddo
      endif
          
c     monitors
      if (nintm .gt. 0) then
          do ir=1,nintm
c             find 4 closest gridded receptors
              call find4(intmon(ir)%x,intmon(ir)%y,grecs)
c             calculate terms of interpolation for each group
              do ig=1,ngrps1
                  do is=1,nmonths
                      call calcint(intmon(ir)%x,intmon(ir)%y,
     +                intmon(ir)%recid,ig,is,grecs,ic)
                      intmc(ig,ir,is)=ic
                  enddo
              enddo
          enddo
      endif  
      if (writefac) close(intunit)
          
      write(logunit,*)'interpolation done'
      call datetime
      return
      end
c##########################################################################
c     subroutine to find the closest 4 gridded receptors to a block or monitor
      subroutine find4(x,y,grecs)
      use main1
      implicit none
c      integer ir,grecs(4),gr,gr1,g,gg,grecs1(nmodg)
      integer*8 grecs(4),gr,gr1,g,gg,ng,ig
c      real*8 dx,dy,mindx,mindy,dist(nmodg),distance,min1,mind(4)
      real*8 dx,dy,distance,min1,mind(4),dist1,x,y
      real*8, allocatable, dimension(:) :: dist
      integer, allocatable, dimension(:) :: grecs1
      grecs=0
c      grecs1=0
c      dist=0.0
      
c     4 points are denoted as below where x is receptor to be interpolated

c     2               3

c             x

c     1               4

c     point 1, dx <= 0 and dy <=0
c     point 2, dx <=0 and dy >=0
c     point 3, dx >=0 and dy >=0
c     point 4, dx >=0 0 and dy <=0

c     first get distance between block and all gridded receptors
c     only keep those within 10 km
      ng=0
      do gr=1,nmodg
          dist1=distance(x,y,modgrid(gr)%x,modgrid(gr)%y)
          if (dist1 .le. 10000.0)ng=ng+1
c          dist(gr)=distance(intblock(ir)%x,intblock(ir)%y,
c     +    modgrid(gr)%x,modgrid(gr)%y)
c          grecs1(gr)=gr
c          write(121,*)gr,dist(gr),grecs1(gr),modgrid(gr)%x,modgrid(gr)%y
      enddo
      
      allocate(dist(ng))
      allocate(grecs1(ng))
      dist=0.0
      grecs1=0.0
      ig=0
      do gr=1,nmodg
          dist1=distance(x,y,modgrid(gr)%x,modgrid(gr)%y)
          if (dist1 .le. 10000.0) then
              ig=ig+1
              dist(ig)=dist1
              grecs1(ig)=gr
          endif
c          dist(gr)=distance(intblock(ir)%x,intblock(ir)%y,
c     +    modgrid(gr)%x,modgrid(gr)%y)
c          grecs1(gr)=gr
c          write(121,*)gr,dist(gr),grecs1(gr),modgrid(gr)%x,modgrid(gr)%y
      enddo     

      
c     now sort based on increasing distance, i.e. first obs of new array is closest
c      do gr=1,nmodg-1
      do gr=1,ng-1
          min1=dist(gr)
          gg=grecs1(gr)
          g=gr
c          do gr1=gr+1,nmodg
          do gr1=gr+1,ng
              if (dist(gr1) .lt. min1) then
                  min1=dist(gr1)
                  gg=grecs1(gr1)
                  g=gr1
              endif
          enddo
          dist(g)=dist(gr)
          grecs1(g)=grecs1(gr)
          dist(gr)=min1
          grecs1(gr)=gg
      enddo
      
c      do gr=1,nmodg
c          write(122,*)gr,dist(gr),grecs1(gr)
c      enddo
      
c     pick 4 points surrounding block
c     in some cases, there may not be 4 due to
c     being on coast or near 50 km.
      mind=100000.
c      do gr=1,nmodg !4
      do gr=1,ng
          gg=grecs1(gr)
          dx=modgrid(gg)%x-x
          dy=modgrid(gg)%y-y
          if (dx .le. 0.0 .and. dy .le. 0.0) gr1=1  !point 1
          if (dx .le. 0.0 .and. dy .ge. 0.0) gr1=2  !point 2
          if (dx .ge. 0.0 .and. dy .ge. 0.0) gr1=3  !point 3
          if (dx .ge. 0.0 .and. dy .le. 0.0) gr1=4  !point 4
          if (dist(gr) .lt. mind(gr1)) then
              mind(gr1)=dist(gr)
              grecs(gr1)=gg
          endif
          
c          if (grecs(gr1) .eq. 0) grecs(gr1)=gg
      enddo
     
      deallocate(dist)
      deallocate(grecs1)
      return
      end
c##########################################################################
c     subroutine to calculate terms of interpolation
c     based on HEM but linear instead of logarithmic
      subroutine calcint(rx,ry,rid,ig,is,grecs,termc)
      use main1
      implicit none
      integer*8 ig,is,grecs(4),nzeros,n,g,gg
      real*8 terma,termb,termc,c(4),x(4),y(4),rx,ry,dx_4_1,dx_3_2,
     +calcterm
      character rid*15,sgrp*12 !sgrp*8

      c=0.0
      x=0.0
      y=0.0
      
c      if (rt .eq. 'blk') then
c          rx=intblock(ir)%x
c          ry=intblock(ir)%y
c          rid=intblock(ir)%recid
c      else
c          rx=intmon(ir)%x
c          ry=intmon(ir)%y
c          rid=intmon(ir)%recid
c      endif
      
      nzeros=0
      do g=1,4
          gg=grecs(g)
          if (gg .eq. 0) then
              nzeros=nzeros+1
          else
              c(g)=modgc(ig,gg,is)
              x(g)=modgrid(gg)%x
              y(g)=modgrid(gg)%y
          endif
      enddo
      
      terma=0.0
      termb=0.0
      termc=0.0
      if (nzeros .eq. 0) then ! all 4 gridded receptors are present
c         check dx and dy spacing between 4 terms
          dx_3_2=abs(x(3)-x(2))
          dx_4_1=abs(x(4)-x(1))
          if (dx_3_2 .ne. dx_4_1) then
              terma=calcterm(c(1),c(4),rx,x(1),x(4))
              termb=calcterm(c(2),c(3),rx,x(2),x(3))
              termc=calcterm(terma,termb,ry,y(1),y(2))
          else
              terma=calcterm(c(1),c(2),ry,y(1),y(2))
              termb=calcterm(c(4),c(3),ry,y(4),y(3))
              termc=calcterm(terma,termb,rx,x(1),x(4))
          endif
c          terma=c(1)+(c(2)-c(1))*(dy/dy1)
c          termb=c(4)+(c(3)-c(4))*(dy/dy1)
c          termc=terma+(termb-terma)*(dx/dx1)
      else !just take average
          n=0
          do g=1,4
              if (grecs(g) .ne. 0) then
                  termc=termc+c(g)
                  n=n+1
              endif
          enddo
          termc=termc/real(n)
      endif
      if (writefac)then
          if (runap) then
              sgrp='ALL'
          elseif (runpt) then
              sgrp=ptsrcs(ig)%srcid
          elseif (runcmv) then
 !             sgrp='ALL'  !fix later
              sgrp=cmvsrcs(ig)%srcid
          else
              sgrp=gridsrcs(ig)%srcid !fix later for gridded
          endif
          write(intunit,5)rid,sgrp,rx,ry,x(1),y(1),c(1),x(2),y(2),
     +    c(2),x(3),y(3),c(3),x(4),y(4),c(4),terma,termb,termc
      endif
      
c      write(131,'(a,2(1x,f11.0),4(2(1x,f11.0),1x,e13.6),3(1x,e13.6))')
c     +rid,rx,ry,x(1),y(1),c(1),x(2),y(2),c(2),x(3),y(3),
c     +c(3),x(4),y(4),c(4),terma,termb,termc
     
 5    format(a,1x,a,2(1x,f11.0),4(2(1x,f11.0),1x,e13.6),3(1x,e13.6))
      return
      end
c##########################################################################
c     calculate terms a and b of interpolation
      function calcterm(c1,c2,xr,x1,x2)
      use main1
      implicit none
      real*8 c1,c2,xr,x1,x2,calcterm
      
      calcterm=c1+(c2-c1)*((xr-x1)/(x2-x1))
      return
      end
c##########################################################################
c     subroutine to calculate HAP-specific concentrations
      subroutine hapconc(ns)
      use main1
      implicit none
      
      integer*8 ns,ihap,ilines,i,nlines1,ncas,r,e,ism,icat,
     +findsrc,findcat,irec,findhap,iskip,rt,eind
c      character cas*15,srcgrp*8,cats1*100,ofile*250,fips1*5
      character cas*15,srcgrp*12,cats1*100,ofile*250,fips1*5,ofile1*250,
     +zipline*300
      real*8 conc,hconc,monemis(12),sx,sy !emis1000(5),
      logical lexist
      
      integer*8, allocatable, dimension(:) :: prid
      integer*8, allocatable, dimension(:) :: uid
      integer*8, allocatable, dimension(:) :: rid
      real*8, allocatable, dimension(:,:) :: emis
*      integer*8, allocatable, dimension(:) :: hunit
      
      allocatable :: cas(:)
      allocatable :: cats1(:)
      allocatable :: srcgrp(:)
      allocatable :: fips1(:)
      
c     now using group specific seasonal emissions, version 17243      
c      emis1000(1)=8571.56256  !winter
c      emis1000(2)=8762.041728 !spring
c      emis1000(3)=8762.041728 !summer
c      emis1000(4)=8666.802144 !spring
c      emis1000(5)=34762.4872  !annual
      
   
c     write facility filename !13352
      if (writefac) then
          write(ofile,'(3(a))')trim(adjustl(hapdir)),
     +    trim(adjustl(facid(ns))),'_concs.txt'
c          hunit(ilines)=hapunit+ilines
           open(hapunit,file=ofile,status='replace')
      endif
c     write facility filename for HEID !13352
      if (writheid) then
          write(ofile1,'(6(a))')trim(adjustl(hapdir)),'facilities/',
     +    trim(adjustl(state2(ns))),'/',trim(adjustl(facid(ns))),
     +    '_block_concs.txt'
c          hunit(ilines)=hapunit+ilines
           open(hapunit+1,file=ofile1,status='replace')
      endif
c     get entries in facility emissions array and write to local array
      
      if (runap) then
          eind=1
          ilines=0
          do i=1,nlines
              if (airport(i)%state .eq. state2(ns) .and. 
     +        airport(i)%facid .eq. facid(ns)) ilines=ilines+1
          enddo
          nlines1=ilines
          ncas=nlines
          allocate(emis(nlines1,1))
          allocate(cas(nlines1))
          allocate(cats1(nlines1))
          allocate(srcgrp(nlines1))
c          allocate(hunit(nlines1))
          allocate(fips1(nlines1))
          
          fips1=''
          ilines=0
          do i=1,nlines
              if (airport(i)%state .eq. state2(ns) .and. 
     +        airport(i)%facid .eq. facid(ns)) then
                  ilines=ilines+1
                  emis(ilines,1)=airport(i)%emis(1)
                  cas(ilines)=airport(i)%cas
                  cats1(ilines)=airport(i)%cat
                  srcgrp(ilines)=airport(i)%srcgrp
c                  if (writefac) then
c                      write(ofile,'(5(a))')trim(adjustl(hapdir)),
c     +                trim(adjustl(cas(ilines))),'_',
c     +                trim(adjustl(facid(ns))),'_concs.txt'
c                      hunit(ilines)=hapunit+ilines
c                      inquire(file=ofile,exist=lexist)
c                      if (.not. lexist) then
c                          open(hunit(ilines),file=ofile,
c     +                    status='replace')
c                      endif
c                  endif
              endif
          enddo
c          do i=1,nlines1
c              write(136,*)i,srcgrp(i),emis(i,1),cas(i),cats1(i)
c          enddo
          
      elseif (runpt) then
          eind=1
          ilines=0
          do i=1,nlines
              if (facility(i)%state .eq. state2(ns) .and. 
     +        facility(i)%facid .eq. facid(ns)) ilines=ilines+1
          enddo
          nlines1=ilines
          ncas=nlines
          allocate(emis(nlines1,1))
          allocate(cas(nlines1))
          allocate(cats1(nlines1))
          allocate(srcgrp(nlines1))
c          allocate(hunit(nlines1))
          allocate(fips1(nlines1))
          allocate(prid(nlines1))
          allocate(uid(nlines1))
          allocate(rid(nlines1))
          fips1=''
          ilines=0
          do i=1,nlines
              if (facility(i)%state .eq. state2(ns) .and. 
     +        facility(i)%facid .eq. facid(ns)) then
                  ilines=ilines+1
                  emis(ilines,1)=facility(i)%emis(1)
                  cats1(ilines)=facility(i)%cat
                  srcgrp(ilines)=facility(i)%srcgrp
                  cas(ilines)=facility(i)%cas
                  ihap=findhap(cas(ilines))
                  prid(ilines)=facility(i)%procid
                  uid(ilines)=facility(i)%unitid
                  rid(ilines)=facility(i)%relid
c                  if (writefac) then
c                      write(ofile,'(5(a))')trim(adjustl(hapdir)),
c     +                trim(adjustl(cas(ilines))),'_',
c     +                trim(adjustl(facid(ns))),'_concs.txt'
c                      hunit(ilines)=hapunit+ihap
cc                     see if file already exists
c                      inquire(file=ofile,exist=lexist)
c                      if (.not. lexist) then
c                          open(hunit(ilines),file=ofile,
c     +                    status='replace')
c                      endif
c                  endif
                  
              endif
          enddo
      elseif (runcmv) then
          eind=1
          ilines=0
          do i=1,nlines
              if (cmv(i)%state .eq. state2(ns) .and. 
     +        cmv(i)%facid .eq. facid(ns)) ilines=ilines+1
          enddo
          nlines1=ilines
          ncas=nlines
          allocate(emis(nlines1,1))
          allocate(cas(nlines1))
          allocate(cats1(nlines1))
          allocate(srcgrp(nlines1))
c          allocate(hunit(nlines1))
          allocate(fips1(nlines1))
          fips1=''
          ilines=0
          do i=1,nlines
              if (cmv(i)%state .eq. state2(ns) .and. 
     +        cmv(i)%facid .eq. facid(ns)) then
                  ilines=ilines+1
                  emis(ilines,1)=cmv(i)%emis(1)
                  cats1(ilines)=cmv(i)%cat
                  srcgrp(ilines)=cmv(i)%srcgrp
                  cas(ilines)=cmv(i)%cas
                  ihap=findhap(cas(ilines))
c                  if (writefac) then
c                      write(ofile,'(5(a))')trim(adjustl(hapdir)),
c     +                trim(adjustl(cas(ilines))),'_',
c     +                trim(adjustl(facid(ns))),'_concs.txt'
c                      hunit(ilines)=hapunit+ihap
cc                     see if file already exists
c                      inquire(file=ofile,exist=lexist)
c                      if (.not. lexist) then
c                          open(hunit(ilines),file=ofile,
c     +                    status='replace')
c                      endif
c                  endif
                  
              endif
          enddo
c      elseif (noncontr) then
c          eind=5
c          ilines=0
c          do i=1,nlines
c              if (tracts(i)%facid .eq. facid(ns)) ilines=ilines+1
c          enddo
c          nlines1=ilines
c          ncas=nlines
c          allocate(emis(nlines1,13))
c          allocate(cas(nlines1))
c          allocate(cats1(nlines1))
c          allocate(srcgrp(nlines1))
c          allocate(fips1(nlines1))
c          allocate(hunit(nlines1))
c          fips1=''
c          ilines=0
c          do i=1,nlines
c              if (tracts(i)%facid .eq. facid(ns)) then
c                  ilines=ilines+1
c                  do e=1,13
c                      emis(ilines,e)=tracts(i)%emis(e)
c                  enddo
c                  cas(ilines)=tracts(i)%cas
c                  cats1(ilines)=tracts(i)%cat
c                  srcgrp(ilines)=tracts(i)%srcgrp
c                  ihap=findhap(cas(ilines))
cc                  if (writefac) then
cc                      write(ofile,'(5(a))')trim(adjustl(hapdir)),
cc     +                trim(adjustl(cas(ilines))),'_',
cc     +                trim(adjustl(facid(ns))),'_concs.txt'
cc                      hunit(ilines)=hapunit+ihap
ccc                     see if file already exists
cc                      inquire(file=ofile,exist=lexist)
cc                      if (.not. lexist) then
cc                          open(hunit(ilines),file=ofile,
cc     +                    status='replace')
cc                      endif
cc                  endif
c              endif
c          enddo
      else
          eind=13
          ilines=0
          do i=1,nlines
              if (gridded(i)%facid .eq. facid(ns)) ilines=ilines+1
          enddo
          nlines1=ilines
          ncas=nlines
          allocate(emis(nlines1,13))
          allocate(cas(nlines1))
          allocate(cats1(nlines1))
          allocate(srcgrp(nlines1))
          allocate(fips1(nlines1))
c          allocate(hunit(nlines1))
          ilines=0
          do i=1,nlines
              if (gridded(i)%facid .eq. facid(ns)) then
                  ilines=ilines+1
                  do e=1,13
                      emis(ilines,e)=gridded(i)%emis(e)
                  enddo
                  cas(ilines)=gridded(i)%cas
                  cats1(ilines)=gridded(i)%cat
                  srcgrp(ilines)=gridded(i)%srcgrp
                  fips1(ilines)=gridded(i)%fips
                  ihap=findhap(cas(ilines))
c                  if (writefac) then
c                      write(ofile,'(5(a))')trim(adjustl(hapdir)),
c     +                trim(adjustl(cas(ilines))),'_',
c     +                trim(adjustl(facid(ns))),'_concs.txt'
c                      hunit(ilines)=hapunit+ihap
cc                     see if file already exists
c                      inquire(file=ofile,exist=lexist)
c                      if (.not. lexist) then
c                          open(hunit(ilines),file=ofile,
c     +                    status='replace')
c                      endif
c                  endif
              endif
          enddo
      endif
      
c     get facility specific, HAP specific concentrations
      
      do i=1,nlines1
          monemis=0.0
          if (runap )then ! .or. runcmv) then
              ism=1
          else
              ism=findsrc(srcgrp(i))
          endif
          if (runpt) then
              sx=ptsrcs(ism)%utmx
              sy=ptsrcs(ism)%utmy
          elseif (runcmv) then
              sx=cmvsrcs(ism)%utmx1
              sy=cmvsrcs(ism)%utmy1
          elseif (runap) then
              sx=apsrcs(ism)%utmx1
              sy=apsrcs(ism)%utmy1
          else
              call calc_center(ism,sx,sy)
          endif
          ihap=findhap(cas(i))
          icat=findcat(cats1(i))
          if (nmonths .eq. 12) then
              do e=1,12
                  monemis(e)=emis(i,e)
              enddo
          endif
c          write(*,*)'s ',i, srcgrp(i),nlines1
c         gridded modeled receptors
          if (nmodg .gt. 0) then
              do r=1,nmodg
c                 check gridded sources to see if gridded receptor is > 50 km from that particular source, if so don't do calculations
                  if (.not. runap .and. .not. runpt .and. .not. 
     +            runcmv) then 
                      call check_rec_grid(modgrid(r)%ux,modgrid(r)%uy,
     +                srcgrp(i),r,iskip)
                  else
                      if (modgrid(r)%iflag .eq. 3) then
                          iskip=1
                      else
                          iskip=0
                      endif
                  endif  
                  if (iskip .eq. 0) then
                      if (nmonths .eq. 1) then
                          conc=modgc(ism,r,1)
                          hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5)
c                     write(121,*)modgc(ism,r,1),conc,emis(i,eind),hconc
c      write(137,*)r,conc,hconc,emis(i,1),emis1000(5),modgrid(r)%recid
                      else
                          call calcann('GM',hconc,ism,monemis,r)
c     +                    emis1000,r)  !calculate annual averages from seasonal averages
                      endif
                      rt=10

                      if (writefac) then
                          if (runpt) then
                              write(hapunit,15)modgrid(r)%recid,
     +                        rt,srcgrp(i),sx,sy,prid(i),uid(i),rid(i),
     +                        cas(i),cats1(i),hconc,
     +                        modgrid(r)%x,modgrid(r)%y,modgrid(r)%ux,
     +                        modgrid(r)%uy,modgrid(r)%dist,fips1(i)
                          else
                              write(hapunit,5)modgrid(r)%recid,
     +                        rt,srcgrp(i),sx,sy,cas(i),cats1(i),hconc,
     +                        modgrid(r)%x,modgrid(r)%y,modgrid(r)%ux,
     +                        modgrid(r)%uy,modgrid(r)%dist,fips1(i)
                          endif
                      endif
                      irec=modgrid(r)%mastrec
                      gridconc(ihap,irec,icat)=gridconc(ihap,irec,icat)+
     +                hconc
                  endif
                  
c      write(*,*)r,irec,gridrec(irec)%recid,gridconc(ihap,irec,icat),
c     +modgrid(r)%recid
c      pause
              enddo
          endif
c         block modeled receptors
          if (nmodb .gt. 0) then
              do r=1,nmodb
                  if (modblock(r)%iflag .ne. 3) then
                      if (nmonths .eq. 1) then
                          conc=modbc(ism,r,1)
                          hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5)
                      else
                          call calcann('BM',hconc,ism,monemis,r)
c     +                    emis1000,r)  !calculate annual averages from seasonal averages
                      endif
                      rt=20
                      if (writefac) then
                          if (runpt) then
                              write(hapunit,15)modblock(r)%recid,
     +                        rt,srcgrp(i),sx,sy,prid(i),uid(i),rid(i),
     +                        cas(i),cats1(i),hconc,
     +                        modblock(r)%x,modblock(r)%y,
     +                        modblock(r)%ux,modblock(r)%uy,
     +                        modblock(r)%dist,fips1(i)
                              if (writheid) write(hapunit+1,20)
     +                        modblock(r)%recid,srcgrp(i),sx,sy,prid(i),
     +                        uid(i),rid(i),cas(i),
     +                        hconc,modblock(r)%ux,modblock(r)%uy
                          else
                              write(hapunit,5)modblock(r)%recid,
     +                        rt,srcgrp(i),sx,sy,cas(i),cats1(i),hconc,
     +                        modblock(r)%x,modblock(r)%y,
     +                        modblock(r)%ux,modblock(r)%uy,
     +                        modblock(r)%dist,fips1(i)
                              if (writheid) write(hapunit+1,10)
     +                        modblock(r)%recid,srcgrp(i),sx,sy,cas(i),
     +                        hconc,modblock(r)%ux,modblock(r)%uy
                          endif
                      endif
                      irec=modblock(r)%mastrec
                      blokconc(ihap,irec,icat)=blokconc(ihap,irec,icat)+
     +                hconc
                  endif
              enddo
          endif
c         block interpolated receptors
          if (nintb .gt. 0) then
              do r=1,nintb
                  if (.not. runap .and. .not. runpt .and. .not. 
     +            runcmv) then 
                      call check_rec_grid(intblock(r)%ux,intblock(r)%uy,
     +                srcgrp(i),r,iskip)
                  else
                      iskip=0
                  endif
                  if (iskip .eq. 0) then
                      if (nmonths .eq. 1) then
                          conc=intbc(ism,r,1)
                          hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5)
                      else
                          call calcann('BI',hconc,ism,monemis,r)
c     +                    emis1000,r)  !calculate annual averages from seasonal averages
                      endif
                      rt=21
                      if (writefac) then
                          if (runpt) then
                              write(hapunit,15)intblock(r)%recid,
     +                        rt,srcgrp(i),sx,sy,prid(i),uid(i),rid(i),
     +                        cas(i),cats1(i),hconc,
     +                        intblock(r)%x,intblock(r)%y,
     +                        intblock(r)%ux,intblock(r)%uy,
     +                        intblock(r)%dist,fips1(i)
                          else
                              write(hapunit,5)intblock(r)%recid,
     +                        rt,srcgrp(i),sx,sy,cas(i),cats1(i),hconc,
     +                        intblock(r)%x,intblock(r)%y,
     +                        intblock(r)%ux,intblock(r)%uy,
     +                        intblock(r)%dist,fips1(i)
                          endif
                      endif
                      irec=intblock(r)%mastrec
                      blokconc(ihap,irec,icat)=blokconc(ihap,irec,icat)+
     +                hconc
                  endif
                  
              enddo
          endif
c         non-populated block modeled receptors
          if (nmodz .gt. 0) then
              do r=1,nmodz
                  if (zmodblock(r)%iflag .ne. 3) then
                      if (nmonths .eq. 1) then
                          conc=modzbc(ism,r,1)
                          hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5)
                      else
                          call calcann('ZB',hconc,ism,monemis,r)
c     +                    emis1000,r)  !calculate annual averages from seasonal averages
                      endif
                      rt=22
                      if (writefac) then
                          if (runpt) then
                              write(hapunit,15)zmodblock(r)%recid,
     +                        rt,srcgrp(i),sx,sy,prid(i),uid(i),rid(i),
     +                        cas(i),cats1(i),hconc,
     +                        zmodblock(r)%x,zmodblock(r)%y,
     +                        zmodblock(r)%ux,zmodblock(r)%uy,
     +                        zmodblock(r)%dist,fips1(i)
                          else
                              write(hapunit,5)zmodblock(r)%recid,
     +                        rt,srcgrp(i),sx,sy,cas(i),cats1(i),hconc,
     +                        zmodblock(r)%x,zmodblock(r)%y,
     +                        zmodblock(r)%ux,zmodblock(r)%uy,
     +                        zmodblock(r)%dist,fips1(i)
                          endif
                      endif
                      irec=zmodblock(r)%mastrec
                      zblokconc(ihap,irec,icat)=
     +                zblokconc(ihap,irec,icat)+hconc
                  endif
              enddo
          endif
c         monitor modeled receptors
          if (nmodm .gt. 0) then
              do r=1,nmodm
                  if (modmon(r)%iflag .ne. 3) then
                      if (nmonths .eq. 1) then
                          conc=modmc(ism,r,1)
                          hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5)
                      else
                          call calcann('MM',hconc,ism,monemis,r)
c     +                    emis1000,r)  !calculate annual averages from seasonal averages
                      endif
                      rt=30
                      if (writefac) then
                          if (runpt) then
                              write(hapunit,15)modmon(r)%recid,rt,
     +                        srcgrp(i),sx,sy,prid(i),uid(i),rid(i),
     +                        cas(i),cats1(i),hconc,
     +                        modmon(r)%x,modmon(r)%y,modmon(r)%ux,
     +                        modmon(r)%uy,modmon(r)%dist,fips1(i)
                          else
                              write(hapunit,5)modmon(r)%recid,rt,
     +                        srcgrp(i),sx,sy,cas(i),cats1(i),hconc,
     +                        modmon(r)%x,modmon(r)%y,modmon(r)%ux,
     +                        modmon(r)%uy,modmon(r)%dist,fips1(i)
                          endif
                      endif
                      
                      irec=modmon(r)%mastrec
                      monconc(ihap,irec,icat)=monconc(ihap,irec,icat)+
     +                hconc
                  endif
              enddo
              
              
          endif
c         monitor interpoloated receptors
          if (nintm .gt. 0) then
              do r=1,nintm
                  if (.not. runap .and. .not. runpt .and. .not. 
     +            runcmv) then 
                      call check_rec_grid(intmon(r)%ux,intmon(r)%uy,
     +                srcgrp(i),r,iskip)  !check to see if outside 50 km for gridded source
                  else
                      iskip=0
                  endif
                  if (iskip .eq. 0) then
                      if (nmonths .eq. 1) then
                          conc=intmc(ism,r,1)
                          hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5)
                      else
                          call calcann('MI',hconc,ism,monemis,r)
c     +                    emis1000,r)  !calculate annual averages from seasonal averages
                      endif
                      rt=31
                      if (writefac) then
                          if (runpt) then
                              write(hapunit,15)intmon(r)%recid,rt,
     +                        srcgrp(i),sx,sy,prid(i),uid(i),rid(i),
     +                        cas(i),cats1(i),hconc,
     +                        intmon(r)%x,intmon(r)%y,intmon(r)%ux,
     +                        intmon(r)%uy,intmon(r)%dist,fips1(i)
                          else
                              write(hapunit,5)intmon(r)%recid,rt,
     +                        srcgrp(i),sx,sy,cas(i),cats1(i),hconc,
     +                        intmon(r)%x,intmon(r)%y,intmon(r)%ux,
     +                        intmon(r)%uy,intmon(r)%dist,fips1(i)
                          endif
                      endif
                      irec=intmon(r)%mastrec
                      monconc(ihap,irec,icat)=monconc(ihap,irec,icat)+
     +                hconc
                  endif
              enddo
          endif
      enddo
      
c      do i=1,nlines1
c          close(hunit(i))
c      enddo
      if (writefac) then
          close(hapunit)
          write(zipline,'(a,1x,a)')'gzip -f',trim(adjustl(ofile))
          call system(zipline)
      endif
      
          
      if (writheid)close(hapunit+1)
      
c 5    format(a15,1x,i2,1x,a8,1x,a15,e13.6,2(1x,f11.0),2(1x,f13.3),1x,
c     +f11.5)
c 5    format(a15,1x,i2,1x,a8,1x,a45,e13.6,2(1x,f11.0),2(1x,f13.3),1x,
c     +f11.5,1x,a5)
 5    format(a15,1x,i2,1x,a8,2(1x,f13.3),1x,a15,1x,a45,e13.6,2(1x,f11.0)
     +,2(1x,f13.3),1x,f11.5,1x,a5)
 10   format(a15,1x,a8,2(1x,f13.3),1x,a15,e13.6,2(1x,f13.3))
 15   format(a15,1x,i2,1x,a8,2(1x,f13.3),3(1x,i15),1x,a15,1x,a45,e13.6,
     +2(1x,f11.0),2(1x,f13.3),1x,f11.5,1x,a5)

 20   format(a15,1x,a8,2(1x,f13.3),3(1x,i15),1x,a15,e13.6,2(1x,f13.3))
c     deallocate arrays
      deallocate(emis)
      deallocate(cas)
      deallocate(cats1)
      deallocate(srcgrp)
      deallocate(fips1)
      if (allocated(prid)) deallocate(prid)
      if (allocated(uid)) deallocate(uid)
      if (allocated(rid)) deallocate(rid)
      return
      end
c##########################################################################
c     subroutine to determine distance from receptor to source 
c     and determine if receptor within 30 m of source
      subroutine recdist(rt,rid,x,y,mindist,maxdist,lwrite,iflag,md)
      use main1
      implicit none

c     iflag:      flag denoting if source within 30 m of source
c     igrp:       source loop counter
      integer*8 iflag,igrp,i
      
c     x:          receptor UTM-x coordinate
c     y:          receptor UTM-y coordinate
c     dist:       distance from receptor to source
c     maxdist:    maximum distance to search
c     apdist:     distance from receptor to airport
c     ptdist:     distance from receptor to point source
      real*8 x,y,dist,mindist,maxdist,apdist,ptdist,md,avgx,avgy,
     +distance

c     lwrite:     logical variable denoting to write receptor to AERMOD.INP
      logical lwrite,ldup
      
      character rt*3,rid*15
      
      iflag=0
      lwrite=.false.
      md=maxdist+10000.0
c      if (trim(adjustl(rt)) .eq. 'mon') pause
      if (runap) then
          do igrp=1,ngrps
              dist=apdist(x,y,igrp)  !get distance from receptor to points along runways
              if (dist .gt. mindist .and. dist .le. maxdist) then
                  lwrite=.true.
                  if (dist .lt. md) md=dist
                  if (dist .le. 30.0) iflag=1
              endif
              
          enddo
      elseif (runpt) then
          do igrp=1,ngrps
              if (ptsrcs(igrp)%srctyp .eq. 'AREA') then
                  dist=apdist(x,y,igrp) !get distance from receptor to points along area source
              else
                  dist=ptdist(x,y,igrp) !get distance from receptor to stacks
              endif
              if (dist .gt. mindist .and. dist .le. maxdist) then 
                  lwrite=.true.
                  if (dist .lt. md) md=dist
                  if (dist .le. 30.0) iflag=1
              endif
          enddo 
      elseif (runcmv) then ! .or. noncontr) then 
          do igrp=1,ngrps
              call cmvcheck(x,y,igrp,dist)  !distance from receptor to CMV polygons
              if (dist .gt. mindist .and. dist .le. maxdist) then 
                  lwrite=.true.
                  if (dist .lt. md) md=dist
                  if (dist .le. 30.0) iflag=1
              else  
                  call polycenter(igrp,avgx,avgy)
                  dist=distance(x,y,avgx,avgy)
                  if (dist .le. maxdist) then
                      lwrite=.true.
                      if (dist .le. 30.0) iflag=1
                      if (dist .lt. md) md=dist
                  endif
              endif
          enddo
      else
          call recs4grid(x,y,grddist,md,lwrite)  !distance from receptor to center of gridded sources
      endif
      
c     if receptor is to be interpolated, make sure not a modeled receptor when near 10 km
c      if (lwrite .and. dist .le. 11000.0) then
c     if (lwrite .and. md .le. 11000.0) then
      if (lwrite) then !check all distances
          i=1
          ldup=.false.
          if (trim(adjustl(rt)) .eq. 'blk') then
              if (nmodb .gt. 0) then
                  ldup=.false.
                  do while (i .le. nmodb .and. .not. ldup)
                      if (trim(adjustl(rid)) .eq. modblock(i)%recid 
     +                .and. modblock(i)%iflag .ne. 3) 
     +                ldup=.true.
                      i=i+1
                  enddo
              endif
          endif
c         if (trim(adjustl(rt)) .eq. 'zbl') then
c             if (nmodz .gt. 0) then
c                 ldup=.false.
c                 do while (i .le. nmodz .and. .not. ldup)
c                     if (trim(adjustl(rid)) .eq. zmodblock(i)%recid 
c    +                .and. zmodblock(i)%iflag .ne. 3) 
c    +                ldup=.true.
c                     i=i+1
c                 enddo
c             endif
c         endif
          if (trim(adjustl(rt)) .eq. 'mon') then
              if (nmodm .gt. 0) then
                  do while (i .le. nmodm .and. .not. ldup)
                      if (trim(adjustl(rid)) .eq. modmon(i)%recid .and. 
     +                modmon(i)%iflag .ne. 3)
     +                ldup=.true.
                      i=i+1
                  enddo
              endif
          endif
          if (ldup)lwrite=.false.
      endif

      return
      end
c##########################################################################
c     subroutine to calculate polygon centroid
      subroutine polycenter(ig,avgx,avgy)
      use main1
      implicit none 
      integer*8 ig,j,nv1
      real*8 avgx,avgy,area
      
      if (runcmv) then
          nv1=cmvsrcs(ig)%nvert
      else
          nv1=tractsrcs(ig)%nvert
      endif
      
      avgx=0.0
      avgy=0.0
      area=0.0
      do j=1,nv1
          if (j .lt. nv1) then
              area=area+(polyx(ig,j)*polyy(ig,j+1)-polyx(ig,j+1)*
     +        polyy(ig,j))
              avgx=(polyx(ig,j)+polyx(ig,j+1))*(polyx(ig,j)*
     +        polyy(ig,j+1)-polyx(ig,j+1)*polyy(ig,j))+avgx
              avgy=(polyy(ig,j)+polyy(ig,j+1))*(polyx(ig,j)*
     +        polyy(ig,j+1)-polyx(ig,j+1)*polyy(ig,j))+avgy
          else
              area=area+(polyx(ig,j)*polyy(ig,1)-polyx(ig,1)*
     +        polyy(ig,j))
              avgx=(polyx(ig,j)+polyx(ig,1))*(polyx(ig,j)*polyy(ig,1)
     +        -polyx(ig,1)*polyy(ig,j))+avgx
              avgy=(polyy(ig,j)+polyy(ig,1))*(polyx(ig,j)*polyy(ig,1)-
     +        polyx(ig,1)*polyy(ig,j))+avgy
          endif
      enddo

      avgx=avgx/(3.0*area)
      avgy=avgy/(3.0*area)
      return
      end
c##########################################################################
c     calculate distance between receptor and point source
      function ptdist(x,y,ig)
      use main1
      implicit none

c     ig:         source index counter
      integer*8 ig
      
c     x:          receptor UTM-x coordinate
c     y:          receptor UTM-y coordinate
c     fx:         UTM-x coordinate of stack
c     fy:         UTM-y coordinate of stack
c     distance:   distance function
c     ptdist:     distance from receptor to stack
      real*8 x,y,fx,fy,distance,ptdist
      
      fx=ptsrcs(ig)%utmx
      fy=ptsrcs(ig)%utmy
      
      ptdist=distance(x,y,fx,fy)
      return
      end
c##########################################################################
      subroutine cmvcheck(x,y,ig,mindist)
      use main1
      implicit none

c     ig:         source index counter
c     eof:        end of file indicator
c     nv:         number of vertices for polygon (read from file)
c     nv1:        number of vertices for polygon being processed
c     j:          vertex counter
      integer*8 ig,eof,nv,nv1,j,i,inout1,inout
      
c     area1:      area of polygon
c     cmvdist:    distance of receptor from polygon edge function
c     dist:       distance of receptor from polygon edge
c     mindist:    minimum distance of receptor from all edges of all modeled polygons
c     frac1:      fraction polygon is of total port
c     h:          release height
c     sz:         sigma-z
c     ux:         array of UTM x-coordinates of vertices
c     uy:         array of UTM y-coordinates of vertices
      real*8 x,y,cmvdist,dist,mindist,xi,yi,xj,yj
      real*8, allocatable, dimension(:) :: ux
      real*8, allocatable, dimension(:) :: uy
      logical ix,iy,jx,jy,L_EOR
      
      eof=0
      mindist=900000.0
      
      if (runcmv) then
          nv=cmvsrcs(ig)%nvert+1 !+1 added 17243
      else
          nv=tractsrcs(ig)%nvert+1 !+1 added 17243
      endif
      
c     version 17234, first check to see if receptor inside the polygon before checking the distances
      allocate(ux(nv))
      allocate(uy(nv))
      
      do j=1,nv-1
          ux(j)=polyx(ig,j)
          uy(j)=polyy(ig,j)
      enddo
      ux(nv)=polyx(ig,1)
      uy(nv)=polyy(ig,1)
      
      INOUT1=-1

      DO I=1,NV
c          write(*,*)i
         XI=UX(I)-X
         YI=UY(I)-Y
C        CHECK WHETHER THE POINT IN QUESTION IS AT THIS VERTEX.
         IF (XI.EQ. 0.0 .AND. YI.EQ. 0.0) THEN
            INOUT1 = 0
            EXIT
         END IF
C        J IS NEXT VERTEX NUMBER OF POLYGON.
         J=1+MOD(I,NV)
         XJ=UX(J)-X
         YJ=UY(J)-Y
C        IS THIS LINE OF 0 LENGTH ?
         IF (XI.EQ.XJ .AND. YI.EQ.YJ) CYCLE
         IX=XI.GE.0.0D0
         IY=YI.GE.0.0D0
         JX=XJ.GE.0.0D0
         JY=YJ.GE.0.0D0

C        CHECK WHETHER (PX,PY) IS ON VERTICAL SIDE OF POLYGON.
         L_EOR = EOR(IY,JY)
         IF (XI.EQ. 0.0 .AND. XJ.EQ. 0.0 .AND. L_EOR) THEN
            INOUT1 = 0
            EXIT
         END IF
C        CHECK WHETHER (PX,PY) IS ON HORIZONTAL SIDE OF POLYGON.
         L_EOR = EOR(IX,JX)
         IF (YI.EQ. 0.0 .AND. YJ.EQ. 0.0 .AND. L_EOR) THEN
            INOUT = 0
            EXIT
         END IF
C        CHECK WHETHER BOTH ENDS OF THIS SIDE ARE COMPLETELY 1) TO RIGHT
C        OF, 2) TO LEFT OF, OR 3) BELOW (PX,PY).
         L_EOR = EOR(IX,JX)
         IF (.NOT.((IY.OR.JY).AND.L_EOR)) then
             CYCLE
          endif
C        DOES THIS SIDE OBVIOUSLY CROSS LINE RISING VERTICALLY FROM (PX,PY)
         L_EOR = EOR(IX,JX)
         IF (.NOT.(IY.AND.JY.AND.L_EOR)) THEN
            IF ((YI*XJ-XI*YJ)/(XJ-XI) .LT. 0.0) THEN
               CYCLE
            ELSE IF ((YI*XJ-XI*YJ)/(XJ-XI) .EQ. 0.0) THEN
               INOUT1 = 0
               EXIT
            ELSE
               INOUT1 = -INOUT1
            END IF
         ELSE
            INOUT1 = -INOUT1
         END IF

      END DO
      
       if (inout1 .ge. 0) then
          mindist=0.0
      else
          do j=1,nv-1
c         calculate distance of receptor to each edge (shortest distance from receptor to anywhere along edge)
c             dist=cmvdist(x,y,polyx(ig,j),polyy(ig,j),polyx(ig,j+1),
c    +        polyy(ig,j+1))
              dist=cmvdist(x,y,ux(j),uy(j),ux(j+1),uy(j+1))
c             get minimum distance of receptor to any polygon
              if (dist .le. mindist) mindist=dist
          enddo
      endif
          
      deallocate(ux)
      deallocate(uy)
      
c      rewind(parunit)
      return
      CONTAINS
        LOGICAL FUNCTION EOR(IX,IY)
          LOGICAL IX, IY,a,b
          a=(IX.OR.IY)
          b=.NOT.(IX.AND.IY)
          EOR = (IX.OR.IY) .AND. .NOT.(IX.AND.IY)
        END FUNCTION
      end
      
c##########################################################################
      function cmvdist(x,y,fx1,fy1,fx2,fy2)
c     calculate shortest distance between receptor and polygon edge
c     anglea is angle with vertex at x1,y1 and endpoints at x2,y2 and receptor x,y
c     a is distance between receptor and x2,y2
c     angleb is angle with vertex at x2,y2 and endpoints at x1,y1 and receptor x,y
c     b is distance between receptor and x1,y1
c     anglec is angle with vertex at receptor x,y and endpoints at x1,y1 and x2,y2
c     c is distance between x1,y1 and x2,y2
c      use law of cosines
      use main1
      implicit none
      
c     x:          receptor UTM-x coordinate
c     y:          receptor UTM-y coordinate
c     fx1:        UTM-x coordinate of vertex 1
c     fy1:        UTM-y coordinate of vertex 1 
c     fx2:        UTM-x coordinate of vertex 2
c     fy2:        UTM-y coordinate of vertex 2 
c     distance:   distance function
c     a:          distance between receptor and fx2,fy2
c     b:          distance between receptor and fx1,fy1
c     c:          distance between fx1,fy1 and fx2,fy2
c     anglea:     angle with vertex at fx1,fy1 and endpoints at fx2,fy2 and receptor x,y
c     angleb:     angle with vertex at fx2f,y2 and endpoints at fx1,fy1 and receptor x,y
c     anglec:     angle with vertex at receptor x,y and endpoints at fx1,fy1 and fx2,fy2
c     cmvdist:    distance between receptor and edge
c     angle:      angle calculation function based on law of cosines
      real*8 x,y,fx1,fy1,fx2,fy2,distance,a,b,c,anglea,angleb,anglec,
     +cmvdist,angle
      
c     calculate distances a, b, and c and anglea, angleb, and anglec
      a=distance(x,y,fx2,fy2)
      b=distance(x,y,fx1,fy1)
      c=distance(fx1,fy1,fx2,fy2)
c     calculate angles
      anglea=angle(b,c,a)
      angleb=angle(a,c,b)
      anglec=angle(a,b,c)
c     check angles to determine next step
      if (anglea .le. 90.0 .and. angleb .le. 90.0) then !receptor is alongside edge
c         calculate distance from receptor to centerline
          cmvdist=b*sin(anglea*d2r)
          
      else
          cmvdist=min(a,b)
      endif
c      write(104,*)x,y,fx1,fy1,fx2,fy2,a,b,c,anglea,angleb,anglec,cmvdist
      return
      end
c##########################################################################
      function apdist(x,y,ig)
c     calculate shortest distance between receptor and airport runway edge
c     or rectangular area source
c     anglea is angle with vertex at x1,y1 and endpoints at x2,y2 and receptor x,y
c     a is distance between receptor and x2,y2
c     angleb is angle with vertex at x2,y2 and endpoints at x1,y1 and receptor x,y
c     b is distance between receptor and x1,y1
c     anglec is angle with vertex at receptor x,y and endpoints at x1,y1 and x2,y2
c     c is distance between x1,y1 and x2,y2
c      use law of cosines
      use main1
      implicit none

c     ig:         source index counter
      integer*8 ig

c     x:          receptor UTM-x coordinate
c     y:          receptor UTM-y coordinate
c     fx1:        UTM-x coordinate of endpoint 1
c     fy1:        UTM-y coordinate of endpoint 1 
c     fx2:        UTM-x coordinate of endpoint 2
c     fy2:        UTM-y coordinate of endpoint 2 
c     distance:   distance function
c     a:          distance between receptor and fx2,fy2
c     b:          distance between receptor and fx1,fy1
c     c:          distance between fx1,fy1 and fx2,fy2
c     anglea:     angle with vertex at fx1,fy1 and endpoints at fx2,fy2 and receptor x,y
c     angleb:     angle with vertex at fx2f,y2 and endpoints at fx1,fy1 and receptor x,y
c     anglec:     angle with vertex at receptor x,y and endpoints at fx1,fy1 and fx2,fy2
c     cmvdist:    distance between receptor and edge
c     angle:      angle calculation function based on law of cosines
c     cprime:     1/2 width of runway for runways, length of area source side for
c                 other area sources
c     d1:         distance from receptor to centerline of runway or edge of area source
c     angled:     obtuse angle (anglea or angleb) -90 degrees
c     d:          distance to corner of runway or area source based on angle D
c     bprime:     distance of distance a or b within runway or area source
c     d2:         portion of b or a outside runway
c     d3:         distance based on right triangle
c     ddprime:    opposite of angle A based on distance c being the adjacent side of angle for right triangle
c     dx1:        = a if angleb > 90, b if anglea > 90
c     dx2:        =a if anglea > 90, b if angleb > 90
c     ang1:       =anglea if anglea > 90, angleb if angleb >90
c     ang2:       =angleb if anglea > 90, anglea if angleb > 90
c     apdist:     minimum of d,d2,dx1,or d3
c     d3a:        distance along runway for distance d3
      real*8 x,y,fx1,fy1,fx2,fy2,a,b,c,anglea,angleb,anglec,bprime,d1,
     +d,d2,d3,cprime,distance,angle,apdist,ang1,ang2,dx1,dx2,angled,
     +ddprime,d3a,aprime,bprime1
      
      apdist=-999.0
      if (runap) then 
          fx1=apsrcs(ig)%utmx1
          fy1=apsrcs(ig)%utmy1
c         if (runtyp .eq. 'APRUN') then
          if (apsrcs(ig)%utmx2 .ge. 0.0) then !18031
              fx2=apsrcs(ig)%utmx2
              fy2=apsrcs(ig)%utmy2
              cprime=apsrcs(ig)%xdim/2.0
          else
              fx2=apsrcs(ig)%utmx1+apsrcs(ig)%xdim/2.0
              fy2=apsrcs(ig)%utmy1+apsrcs(ig)%ydim
              cprime=apsrcs(ig)%xdim/2.0
          endif
      else
          cprime=ptsrcs(ig)%xdim/2.0
          a=sqrt((cprime*cprime)+(ptsrcs(ig)%ydim*ptsrcs(ig)%ydim))
          b=ptsrcs(ig)%ydim
          angleb=angle(a,cprime,b)
          anglec=angle(a,b,cprime)
          aprime=3.141592654-ptsrcs(ig)%angle-angleb-anglec  !remember angles are in radians, 3.14159 is 180
          fx1=ptsrcs(ig)%utmx+(cprime*sin(aprime))
          fy1=ptsrcs(ig)%utmy-(cprime*cos(aprime))
          bprime1=(3.141592654/2.0)-ptsrcs(ig)%angle-anglec !remember angles are in radians, 3.14159/2 is 90
          fx2=ptsrcs(ig)%utmx+(a*cos(bprime1))
          fy2=ptsrcs(ig)%utmy+(a*sin(bprime1))
          
c         fx1=ptsrcs(ig)%utmx+ptsrcs(ig)%xdim/2.0
c         fy1=ptsrcs(ig)%utmy
c         fx2=ptsrcs(ig)%utmx+ptsrcs(ig)%xdim/2.0
c         fy2=ptsrcs(ig)%utmy+ptsrcs(ig)%ydim
c         cprime=ptsrcs(ig)%xdim/2.0
      endif
      
!      write(*,*)'a ',fx1,fy1,fx2,fy2,cprime,x,y
!      pause
c     calculate distances
      a=distance(x,y,fx2,fy2)
      b=distance(x,y,fx1,fy1)
      c=distance(fx1,fy1,fx2,fy2)
c     calculate angles
      anglea=angle(b,c,a)
      angleb=angle(a,c,b)
      anglec=angle(a,b,c)
      
!      write(*,*)'b ',a,b,c,anglea,angleb,anglec
!      pause
c     check angles to determine next step
      if (anglea .le. 90.0 .and. angleb .le. 90.0) then !receptor is alongside runway or on runway
c         calculate distance from receptor to centerline
          d1=b*sin(anglea*d2r)
          apdist=d1-cprime  !distance to edge
!          write(*,*)'c ',d1
!          pause
      else
          if (anglea .gt. 90) then
              ang1=anglea
              dx1=b
              dx2=a
              ang2=angleb
          else
              ang1=angleb
              dx1=a
              dx2=b
              ang2=anglea
          endif
          angled=-90.0+ang1  !90-180+ang1
          d=sqrt((dx1*dx1)+(cprime*cprime)-
     +    (2*dx1*cprime*cos(angled*d2r)))  !distance to corner of runway based on angle D
          bprime=c/cos(ang2*d2r)
          ddprime=c*tan(ang2*d2r)  !opposite of angle A based on distance c being the adjacent side of angle for right triangle
!          write(*,*)'d ',d,ang1,dx1,dx2,ang2,angled,bprime,ddprime
!          pause
          if (ddprime .le. cprime) then
              d2=dx2-bprime
              d3=dx1*sin(angled*d2r)  !distance based on right triangle with dx1 as the hypotenuse
              d3a=dx1*cos(angled*d2r) !distance along runway edge to d3
              if (d3a .gt. cprime) then !d3a and associated d3 is past runway
                 apdist=min(dx2,d2,dx1)
              else
                 apdist=min(d,d2,dx1,d3)
              endif
!              write(*,*)'e ',d2,d3,apdist
!              pause
          else
              apdist=min(d,dx1)
!              write(*,*)'f ',dx1,apdist
!              pause
          endif
      endif
      
      
      return
      end
c##########################################################################
      function angle(d1,d2,d3)
c     calculate angles based on law of cosines
      
      use main1
      implicit none

c     angle:      angle returned
c     d1:         distance 1
c     d2:         distance 2
c     d3:         distance 3 (corresponds to opposite side of angle being determined)
c     xx:         ((d1*d1)+(d2*d2)-(d3*d3))/(2*d1*d2)
      real*8 angle,d1,d2,d3,xx
      
      
      xx=((d1*d1)+(d2*d2)-(d3*d3))/(2*d1*d2)
      angle=acos(xx)/d2r
      return
      end
c##########################################################################
      subroutine recs4grid(x,y,griddist,dist,lwrite)
      use main1
      implicit none
 
c     ig:         source index counter
      integer*8 ig
      
c     x:          receptor UTM-x coordinate
c     y:          receptor UTM-y coordinate
c     cenx:       UTM-x coordinate of grid cell center
c     cenx:       UTM-y coordinate of grid cell center
c     distance:   distance function
c     dist:       distance from receptor to center coordinates
c     griddist:   max distance to model
      real*8 x,y,dist,distance,cenx,ceny,griddist
      
c     logical variable denoting to write receptor to AERMOD.INP
      logical lwrite
      
      lwrite=.false.
      
c     only need to write the receptor once if within 50 km of multiple sources
c     if only 1 source then no problem
      ig=1
      do while(ig .le. ngrps .and. .not. lwrite)
          call calc_center(ig,cenx,ceny)
          dist=distance(x,y,cenx,ceny)
          if (dist .le. griddist) lwrite=.true.
          ig=ig+1
      enddo
      return
      end
c##########################################################################
      subroutine calc_center(ig,cx,cy)
      use main1
      implicit none
      
c     ig:         source index counter
      integer*8 ig
      
c     cx:         UTM-x coordinate of grid cell center
c     cx:         UTM-y coordinate of grid cell center
      real*8 cx,cy
      
cc      cx=0.0
cc      cy=0.0
      cx=(gridsrcs(ig)%utmx1+gridsrcs(ig)%utmx2+gridsrcs(ig)%utmx3+
     +gridsrcs(ig)%utmx4)/(real(gridsrcs(ig)%nvert))
      cy=(gridsrcs(ig)%utmy1+gridsrcs(ig)%utmy2+gridsrcs(ig)%utmy3+
     +gridsrcs(ig)%utmy4)/(real(gridsrcs(ig)%nvert))
      return
      end
c##########################################################################
      function distance(x,y,fx,fy)
      use main1
      implicit none
      
c     x:          x-coordinate of point 1
c     y:          y-coordinate of point 1
c     fx:         x-coordinate of point 2
c     fy:         y-coordinate of point 2
c     distance:   distance between point 1 and point 2
c     dx:         x-distance between point 1 and point 2
c     dy:         y-distance between point 1 and point 2     
      real*8 x,y,fx,fy,distance,dx,dy
      
      dx=x-fx
      dy=y-fy
      
      distance=sqrt(dx*dx+dy*dy)
      return
      end
      
c##########################################################################
      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
c##########################################################################
      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 (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

c##########################################################################
      function findsrc(src)
      use main1
      implicit none

c     findsrc:    source index in source array
      integer*8 findsrc

c     src:        AERMOD source id (8 characters)
c     src12:      AERMOD source ID (12 characters)
      character src*12 !,src12*12
      
c     lfound:     logical variable denoting source found in array
      logical lfound
      
c     src12=src
      findsrc=1
      lfound=.false.
      if (runpt) then
          do while(findsrc .le. ngrps .and. .not. lfound)
              if (trim(adjustl(src)) .eq. 
     +        trim(adjustl(ptsrcs(findsrc)%srcid)))then
                  lfound=.true.
              else
                  findsrc=findsrc+1
              endif
          enddo
      elseif (runap) then
          do while(findsrc .le. ngrps .and. .not. lfound)
              if (trim(adjustl(src)) .eq. 
     +        trim(adjustl(apsrcs(findsrc)%srcid)))then
                  lfound=.true.
              else
                  findsrc=findsrc+1
              endif
          enddo
      elseif (runcmv) then
          do while(findsrc .le. ngrps .and. .not. lfound)
              if (trim(adjustl(src)) .eq. cmvsrcs(findsrc)%srcid)then
                  lfound=.true.
             else
                  findsrc=findsrc+1
              endif
          enddo
c      elseif (noncontr) then
c          do while(findsrc .le. ngrps .and. .not. lfound)
c              if (trim(adjustl(src)) .eq. tractsrcs(findsrc)%srcid)then
c                  lfound=.true.
c             else
c                  findsrc=findsrc+1
c              endif
c          enddo
      else
          do while(findsrc .le. ngrps .and. .not. lfound)
c             write(*,*)src,gridsrcs(findsrc)%srcid
              if (trim(adjustl(src)) .eq. 
     +        trim(adjustl(gridsrcs(findsrc)%srcid)))then
                  lfound=.true.
              else
                  findsrc=findsrc+1
              endif
          enddo
      endif
      
      if (.not. lfound)findsrc=0
      return
      end
c##########################################################################
c     function to find facility index
      function findfac(fid)
      use main1
      implicit none
      integer findfac
      character fid*15
      logical lfound
      
      findfac=1
      lfound=.false.
      do while(findfac .le. nsrc .and. .not. lfound)
          if (trim(adjustl(fid)) .eq. trim(adjustl(facid(findfac))))then
              lfound=.true.
          else
              findfac=findfac+1
          endif
      enddo
     
      if (.not. lfound) findfac=0
      return
      end
c##########################################################################
c     subroutine to deallocate arrays
      subroutine cleanup
      use main1
      implicit none
      
      if (runap) then
          deallocate(apsrcs)
      elseif (runpt) then
          deallocate(ptsrcs)
      elseif (runcmv) then
          deallocate(cmvsrcs)
c      elseif (noncontr) then
c          deallocate(tractsrcs)
      else
          deallocate(gridsrcs)
      endif
      
      if (allocated(memis)) deallocate(memis)
      if (allocated(polyx)) deallocate(polyx)
      if (allocated(polyy)) deallocate(polyy)
      if (allocated(modgc)) deallocate(modgc)
      if (allocated(modgrid)) deallocate(modgrid)
      if (allocated(modbc)) deallocate(modbc)
      if (allocated(modblock)) deallocate(modblock)
      if (allocated(modmc)) deallocate(modmc)
      if (allocated(modmon)) deallocate(modmon)
      if (allocated(intbc)) deallocate(intbc)
      if (allocated(intblock)) deallocate(intblock)
      if (allocated(intmc)) deallocate(intmc)
      if (allocated(intmon)) deallocate(intmon)
      
      if (allocated(gridrec1)) deallocate(gridrec1)
      if (allocated(blockrec1)) deallocate(blockrec1)
      if (allocated(monrec1))deallocate(monrec1)
      
      if (allocated(modzbc)) deallocate(modzbc)
      if (allocated(zmodblock)) deallocate(zmodblock)
      return
      end
c##########################################################################
c     subroutine to calculate annual average concentration from seasonal
c     concentration
c      subroutine calcann(tag,conc,is,emis,emis1000,ir)
      subroutine calcann(tag,conc,is,emis,ir)
      use main1
      implicit none
      
      integer*8 is,ir,e,hrs,h
      real*8 conc,emis(12) !,emis1000(4)
      character tag*2
      
      conc=0.0
      hrs=0
      do e=1,12
          if (tag .eq. 'GM') then
              conc=conc+modgc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e)
          elseif (tag .eq. 'BM') then
              conc=conc+modbc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e)
          elseif (tag .eq. 'ZB') then
              conc=conc+modzbc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e)
          elseif (tag .eq. 'MM') then
              conc=conc+modmc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e)
          elseif (tag .eq. 'BI') then
              conc=conc+intbc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e)
          else
              conc=conc+intmc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e)
          endif
          hrs=hrs+nhours(e)
      enddo
c      if (tag .eq. 'MM') then
c          do e=1,4
c              write(113,*)ir,e,conc,modmc(is,ir,e),hrs
c              do h=1,24
c                  write(113,*)h,seashrs(e,h)
c              enddo
c          enddo
c      endif
      conc=conc/real(hrs)
      return
      end
c##########################################################################
c     function to find receptor index
      function findrec(tag,r)
      use main1
      implicit none
      integer*8 r,findrec,fr
      character tag*2
      real*8 x1,y1,x2,y2
      logical lfound
      fr=1
      lfound=.false.
      if (tag .eq. 'GM') then
          do while(fr .le. ngrids1 .and. .not. lfound)
              if (trim(adjustl(gridrec1(fr)%recid)) .eq. 
     +        trim(adjustl(modgrid(r)%recid))) then
                  lfound=.true.
              else
                  fr=fr+1
              endif
          enddo
          findrec=gridrec1(fr)%irec
      endif
      
      if (tag .eq. 'BM') then
          do while(fr .le. nblocks1 .and. .not. lfound)
              if (trim(adjustl(blockrec1(fr)%recid)) .eq. 
     +        trim(adjustl(modblock(r)%recid))) then
                  lfound=.true.
              else
                  fr=fr+1
              endif
          enddo
          findrec=blockrec1(fr)%irec
      endif
      
      if (tag .eq. 'ZB') then
          do while(fr .le. nzblocks1 .and. .not. lfound)
              if (trim(adjustl(zblockrec1(fr)%recid)) .eq. 
     +        trim(adjustl(zmodblock(r)%recid))) then
                  lfound=.true.
              else
                  fr=fr+1
              endif
          enddo
          findrec=zblockrec1(fr)%irec
      endif
      
      if (tag .eq. 'MM') then
          do while(fr .le. nmons1 .and. .not. lfound)
              if (trim(adjustl(monrec1(fr)%recid)) .eq. 
     +        trim(adjustl(modmon(r)%recid))) then
                  lfound=.true.
              else
                  fr=fr+1
              endif
          enddo
          findrec=monrec1(fr)%irec
      endif
      
      if (tag .eq. 'BI') then
          do while(fr .le. nblocks1 .and. .not. lfound)
              if (trim(adjustl(blockrec1(fr)%recid)) .eq. 
     +        trim(adjustl(intblock(r)%recid))) then
                  lfound=.true.
              else
                  fr=fr+1
              endif
          enddo
          findrec=blockrec1(fr)%irec
      endif
      
      if (tag .eq. 'MI') then
          do while(fr .le. nmons1 .and. .not. lfound)
              if (trim(adjustl(monrec1(fr)%recid)) .eq. 
     +        trim(adjustl(intmon(r)%recid))) then
                  lfound=.true.
              else
                  fr=fr+1
              endif
          enddo
          findrec=monrec1(fr)%irec
      endif
      return
      end
c##########################################################################
c     funtion to find HAP index
      function findhap(cas1)
      use main1
      implicit none
      integer*8 findhap
      character cas1*15
      logical lfound
      
      findhap=1
      lfound=.false.
      do while(findhap .le. nhaps .and. .not. lfound)
          if (trim(adjustl(cas1)) .eq. trim(adjustl(haps(findhap)))) 
     +    then
              lfound=.true.
          else
              findhap=findhap+1
          endif
      enddo
      if (.not. lfound) findhap=0
      return
      end
c##########################################################################
c     funtion to find category index
      function findcat(cat1)
      use main1
      implicit none
      integer*8 findcat
      character cat1*100
      logical lfound
      
      findcat=1
      lfound=.false.
      do while(findcat .le. ncats .and. .not. lfound)
          if (trim(adjustl(cat1)) .eq. trim(adjustl(cats(findcat)))) 
     +    then
              lfound=.true.
          else
              findcat=findcat+1
          endif
      enddo
      if (.not. lfound) findcat=0
      return
      end
c##########################################################################
c     subroutine to write output files
c     only write out when total concentration is > 0
      subroutine writeout
      use main1
      implicit none
      integer*8 h,r,c,rt
      character outfil*250,ares*2,form1*45,form2*25,acom*1,acom1,
     +rt1*30,arun*10,zipline*300
      real*8 totl,calctot
      logical lfill
      allocatable :: acom(:)
     
      write(logunit,*)'Writing output files..'
      call datetime
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),',ncats,
     +'(e12.6,a1))'
      write(form2,'(a,i3,a)')'(a,',ncats*2,'(a))'
      
      allocate(acom(ncats))
      acom1=','
      acom=','
      acom(ncats)=' '  !set last one to blank
      write(arun,'(i10)')runno !run #

c     if single facility or state, append that to filename
      if (lsingfac) then
          write(rt1,'(3(a))')trim(adjustl(runtyp)),'_',
     +    trim(adjustl(singlfac))
      elseif (lsingst) then
          write(rt1,'(3(a))')trim(adjustl(runtyp)),'_',
     +    trim(adjustl(singlst))
      else
          rt1=trim(adjustl(runtyp))
      endif

      do h=1,nhaps
          call datetime
          lfill=.false.
          write(logunit,'(2(a))')'Writing output for ',
     +    trim(adjustl(haps(h)))
          if (runpt .or. runap .or. runcmv) then
              write(outfil,'(7(a))')trim(adjustl(hapdir)),
     +        trim(adjustl(haps(h))),'_',trim(adjustl(rt1)),
     +        '_',trim(adjustl(arun)),'_concs.txt'
          else
              write(ares,'(i2.2)')gridres
               write(outfil,'(9(a))')trim(adjustl(hapdir)),
     +        trim(adjustl(haps(h))),'_',trim(adjustl(rt1)),'_',
     +        trim(adjustl(ares)),'_',trim(adjustl(arun)),'_concs.txt'
          endif
          open(hapunit,file=outfil,status='unknown')
c         write header record
          write(hapunit,form2)'receptor_id,recnum,rec_type,col,row,',
     +    (trim(adjustl(cats(c))),acom(c),c=1,ncats)
          do r=1,ngrids !gridded
              rt=1
              totl=calctot('G',h,r)
              if (totl .gt. 0.0) then
                  lfill=.true.
                  write(hapunit,form1)trim(adjustl(gridrec(r)%recid)),
     +            acom1,gridrec(r)%irec,acom1,rt,acom1,gridrec(r)%col,
     +            acom1,gridrec(r)%row,acom1,(gridconc(h,r,c),acom(c),
     +            c=1,ncats)
              endif
          enddo
          do r=1,nblocks !blocks
              rt=2
              totl=calctot('B',h,r)
              if (totl .gt. 0.0) then
                  if (.not. lfill) lfill=.true.
                  write(hapunit,form1)trim(adjustl(blockrec(r)%recid)),
     +            acom1,blockrec(r)%irec,acom1,rt,acom1,blockrec(r)%col,
     +            acom1,blockrec(r)%row,acom1,(blokconc(h,r,c),
     +            acom(c),c=1,ncats)
              endif
          enddo
          do r=1,nzblocks !blocks
              rt=4
              totl=calctot('Z',h,r)
              if (totl .gt. 0.0) then
                  if (.not. lfill) lfill=.true.
                  write(hapunit,form1)trim(adjustl(zblockrec(r)%recid)),
     +            acom1,zblockrec(r)%irec,acom1,rt,acom1,
     +            zblockrec(r)%col,acom1,zblockrec(r)%row,acom1,
     +            (zblokconc(h,r,c),acom(c),c=1,ncats)
              endif
          enddo
          do r=1,nmons !monitors
              rt=3
              totl=calctot('M',h,r)
              if (totl .gt. 0.0) then
                  if (.not. lfill) lfill=.true.
                  write(hapunit,form1)trim(adjustl(monrec(r)%recid)),
     +            acom1,monrec(r)%irec,acom1,rt,acom1,monrec(r)%col,
     +            acom1,monrec(r)%row,acom1,(monconc(h,r,c),acom(c),
     +            c=1,ncats)
              endif
          enddo
          write(hapunit,*)' '
          call datetime
          close(hapunit)
          write(zipline,'(a,1x,a)')'gzip -f',trim(adjustl(outfil))
          call system(zipline)  !zip file
          if (.not. lfill) write(logunit,'(2(a))')
     +    'No output written for ',trim(adjustl(haps(h)))
      enddo
      write(logunit,*)'Output files written'
      call datetime
      return
      end
c##########################################################################
c     fucntion to total up across the categories
      function calctot(rt,h,r)
      use main1
      implicit none
      integer*8 h,r,c
      real*8 calctot
      character rt*1
      calctot=0
      
      if (trim(adjustl(rt)) .eq. 'G') then
          do c=1,ncats
              calctot=gridconc(h,r,c)+calctot
          enddo
      elseif (trim(adjustl(rt)) .eq. 'B') then
          do c=1,ncats
              calctot=blokconc(h,r,c)+calctot
          enddo
      elseif (trim(adjustl(rt)) .eq. 'Z') then
          do c=1,ncats
              calctot=zblokconc(h,r,c)+calctot
          enddo
      else
          do c=1,ncats
              calctot=monconc(h,r,c)+calctot
          enddo
      endif
      return
      end
          
c##########################################################################
c     subroutine to check specific gridded source to receptor distance
      subroutine check_rec_grid(x,y,src,r,iskip)
      use main1
      implicit none
      character src*12
      integer*8 r,iskip,is,findsrc
      real*8 distance,d,cenx,ceny,x,y
      
c      x=modgrid(r)%ux
c      y=modgrid(r)%uy
c     find the source in the source list
      is=findsrc(src)
      call calc_center(is,cenx,ceny)
      d=distance(x,y,cenx,ceny)
      if (d .gt. grddist) then 
          iskip=1
      else
          iskip=0
      endif
c      write(137,'(a15,2(1x,f13.3),1x,a,3(1x,f13.3),1x,i1)')
c     +modgrid(r)%recid,x,y,src,cenx,ceny,d,iskip
      return
      end
c##########################################################################
c     function to get state index
      function get_state(sstate)
      use main1
      implicit none
      integer*8 i,get_state,ii
      character st1*2,sstate*2
      logical lfound,lnum
      lfound=.false.
      i=1
      ii=ichar(sstate(1:1))
      lnum=.false.
      if (ii .ge. 48 .and. ii .le. 57 .and. ii .ge. 48 .and. ii .le. 57)
     +lnum=.true.  !state is defined by fips code, not abbreviation
      do while (i .le. 54 .and. .not. lfound)
          if (lnum) then
              st1=states1(i)
          else
              st1=states(i)
          endif
          if (trim(adjustl(sstate)) .eq. trim(adjustl(st1))) then
              lfound=.true.
          else
              i=i+1
          endif
      enddo
      if (.not. lfound) then
          get_state=0
      else
          get_state=i
      endif
      
      return
      end
          
c##########################################################################
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.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
