c     program to create AERMOD input files for NATA
c     command line argument if facility ID
c##########################################################################
      module main1
      implicit none

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     fips:       County FIPS code      
c     onsrc:      onroad or RWC source IDs in county-cell crosswalk for onroad 4 km or RWC
c     ncstfip:    state fips code for non-CONUS, needed for filenames
c     metdom:     met domain
      character runtyp*7,facid*20,state2*2,ptloc*250,ptparam*250,
     +srcnam*75,pttemp*250,ptcmaq*250,metdir*250,recdir*250,
     +outdir*250,locsrc*12,parsrc*12,temsrc*12,locfil*300,elevfil*300,
     +paramfil(2)*300,temfil(2)*300,cntyfil*300,xdir*250,fips*5,
     +stfips*2,onsrc*12,rwcsrc*12,ncstfip*2,versn*5,logfil*250,metdom*8 !akmon(10)*15, ,aktract(10)*15

      allocatable :: locsrc(:)
      allocatable :: parsrc(:)
      allocatable :: temsrc(:)
      allocatable :: fips(:)
      allocatable :: onsrc(:)
      allocatable :: stfips(:)
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     runog:      oil and gas nonpoint sources being processed
c     akhiprvi:   variable denoting that AK, HI, PR, or VI being processed
c     lleap:      variable denoting year is a leap year
c     noncontr:   variable denoting AK, HI, PR, or VI tract being processed      
      logical lparam(2),ltemp(3),lcnty,runpt,runap,runon,runnon,runrwc,
     +runcmv,runnphi,runnplow,runog,runag,akhiprvi,lleap,runport,
     +noncontr !(for tracts)
      
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
c     polyx:       CMV x-coordinates
c     polyy:       CMV y-coordinates
c     polycx:      centroid of tract or cmv polygon
c     polycy:      centroid of trat or cmv polygon
c     seasemis:    seasonal emissions for 10000 tons rate (1=winter,2=spring,3=summer,4=fall)
      real*8 d2r,blkdist,grddist,profbase,annemis !polycx,polycy,
      
      real*8, allocatable, dimension(:,:) :: polyx
      real*8, allocatable, dimension(:,:) :: polyy
      
c     emisfact:   AERMOD emissions factors
      real*8, allocatable, dimension(:) :: emisfact
c     emisfact1:  allocated emissions, match emissions factors
      real*8, allocatable, dimension(:) :: emission
      
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     wrfunit:    WRF elevations file unit (31)
c     days:       array of number of days/month
c     modyear:    year being modeled
c     nday3:      # of weekdays,saturday,sunday per month
c     nday7:      # of days per month (monday, tuesday, etc.)
      integer col,row,inpunit,locunit,parunit,temunit,recunit,aerunit,
     +ngrps,utmzone,hrunit,gridres,cntyunit,wrfunit,days(12),
     +metunit,nstate,modyear,seasons(12),hrfix,nrec,ilogunit,coffset,
     +roffset
      

c     point sources
c     srcid:      AERMOD source ID
c     srctype:    AERMOD source type
c     qflag:      AERMOD emissions factor emissions flag
c     utmx:       UTM x-coordinate of AERMOD source (SW corner for area source)
c     utmy:       UTM y-coordinate of AERMOD source (SW corner for area source)
c     elev:       source elevation (m)
c     ht:         stack or release height (m)
c     stktemp:    stack exit temperature (K)
c     stkvel:     stack exit velocity (m/s)
c     stkdiam:    stack diameter (m)
c     xdim:       area source x-dimension (m)
c     ydim:       area source y-dimension (m)
c     szinit:     initial sigma-z for area sources (m)
c     angle:      area source angle relative to north (degrees)
c     hrsrc:      integer flag indicating if AERMOD source is an 
c                 hourly emissions source
      type facinfo
          character srcid*12,srctyp*8,qflag*8
          real*8 utmx,utmy,elev,ht,stktemp,stkvel,stkdiam,xdim,ydim,
     +    szinit,angle,annual
          integer hrsrc
      end type facinfo
      
c     airport (runway and non-runway)
c     srcid:      AERMOD source ID
c     srctype:    AERMOD source type
c     qflag:      AERMOD emissions factor emissions flag
c     utmx1:      UTM x-coordinate of runway endpoint or SW corner of area source
c     utmy1:      UTM y-coordinate of runway endpoint or SW corner of area source
c     utmx2:      UTM x-coordinate of runway endpoint or SE corner of area source
c     utmy2:      UTM y-coordinate of runway endpoint or SE corner of area source
c     elev:       source elevation (m) 
c     frac:       fraction runway is of airport area 
c     area:       area of runway or area source (sq meters)
c     ht:         release height (m)
c     xdim:       runway width or length of area source side
c     ydim:       width of area source
c     szinit:     initial sigma-z for area sources (m)
c     angle:      area source angle relative to north (degrees)
      
      type apinfo
          character srcid*12,srctyp*8,qflag*8
          real*8 elev,area,utmx1,utmx2,utmy1,utmy2,frac,ht,xdim,ydim,
     +    szinit,angle,annual
      end type apinfo
      
c     gridded sources (nonroad, onroad, nonpoint, RWC)
c     srcid:      AERMOD source ID
c     srctype:    AERMOD source type
c     qflag:      AERMOD emissions factor emissions flag
c     utmx1:      UTM x-coordinate of SW corner of grid cell
c     utmy1:      UTM y-coordinate of SW corner of grid cell
c     utmx2:      UTM x-coordinate of grid cell vertex
c     utmy2:      UTM y-coordinate of grid cell vertex
c     utmx3:      UTM x-coordinate of grid cell vertex
c     utmy3:      UTM y-coordinate of grid cell vertex
c     utmx4:      UTM x-coordinate of grid cell vertex
c     utmy4:      UTM y-coordinate of grid cell vertex
c     ht:         release height (m)
c     szinit:     initial sigma-z for area sources (m)
c     nvert:      number of vertices for grid cell (4)
      type gridinfo
          character srcid*12,srctyp*8,qflag*8,state*2
          real*8 utmx1,utmx2,utmx3,utmx4,utmy1,utmy2,utmy3,utmy4,ht,
     +    szinit,annual,monemis(12)
          integer nvert
      end type gridinfo
      
c     CMV (port and underway)
c     srcid:      AERMOD source ID
c     srctype:    AERMOD source type
c     qflag:      AERMOD emissions factor emissions flag
c     utmx1:      UTM x-coordinate of SW corner of port polygon
c     utmy1:      UTM y-coordinate of SW corner of port polygon
c     frac:       fraction polygon is of port facility
c     area:       area of port polygon
c     ht:         release height (m)
c     szinit:     initial sigma-z for area sources (m)
c     nvert:      number of vertices for grid cell (4)
      type cmvinfo
          character srcid*12,srctyp*8,qflag*8
          real*8 utmx1,utmy1,ht,szinit,area,avgx,avgy,annual
          integer nvert
      end type cmvinfo
c     non-CONUS tract sources
      type tractinfo
          character srcid*12,srctyp*8,qflag*8,fips*5,tract*11
          real*8 utmx1,utmy1,ht,szinit,area,avgx,avgy,annual,monemis(12)
          integer nvert
      end type tractinfo
      
      type(facinfo), dimension(:),allocatable:: facility
      type(apinfo), dimension(:),allocatable:: airport
      type(gridinfo), dimension(:),allocatable:: gridded
      type(cmvinfo), dimension(:),allocatable:: cmv
      type(tractinfo), dimension(:),allocatable:: tract
      
c     parameters
      parameter(inpunit=11,locunit=12,parunit=13,temunit=14,recunit=16,
     +aerunit=27,hrunit=28,cntyunit=29,wrfunit=31,metunit=32,
     +d2r=3.141592654/180.0,grddist=50000.0,blkdist=10000.0,
     +annemis=10000.,ilogunit=18)
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/  
      
c     assign season to each month
c     1=winter, 2=spring, 3=summer, 4=autumn
c                 J  F  M  A  M  J  J  A S  O  N  D
      data seasons /1,1,2,2,2,3,3,3,4,4,4,1/
      data versn /'21202'/
      
c      data akmon /'020130002','020680003','020900010','020900034',
c     +'020900035','02090NOAA_CRV','021220009','021700011',
c     +'02185NOAA_BRW','022909000'/
      
c      data aktract /'02013000100','02068000100','02090000100',
c     +'02090000500','02090001500','02090001900','02122000100',
c     +'02170000101','02185000100','02290000200'/
      end module main1
c##########################################################################
      program makeinps
      use main1
      implicit none

c     read command line arguments
      call readargs
c     write AERMOD input file 
      call writeinp
      write(ilogunit,*)'Finished'
      write(*,*)'Finished'
      close(ilogunit)
      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     acol:       character string of CMAQ column
c     arow:       character string of CMAQ row
      
      integer i,ilen,ndash,n1,n2,istat,ii
      character agrid*2,acol*3,arow*3,runtyp1*12,noncon*1,ayear*4
      logical lbad
c     initialize logical flags for run types
      runpt=.false.
      runap=.false.
      runcmv=.false.
      runon=.false.
      runnon=.false.
      runrwc=.false.
      runnphi=.false.
      runnplow=.false.
      runport=.false.
      runog=.false.
      runag=.false. !new
c     initialize variable denoting if AK, HI, PR, or VI
c     being processed
      akhiprvi=.false.
      noncontr=.false.

 
c     get year
      call getarg(1,ayear,istat)
      if (istat .ne. -1) ilen=len_trim(ayear)
      do i=1,ilen
          call upcase(ayear(i:i))
      enddo
      if (istat .eq. -1 .or. index(ayear,'HEL') .gt. 0) then   !give user the order of command line arguments
          write(*,4)
          stop
      endif
      
c      ilen=len_trim(ayear)
      i=1
      lbad=.false.
      do while (i .le. ilen .and. .not. lbad)
          ii=ichar(ayear(i:i))
          if (ii .lt. 48 .or. ii .gt. 57) lbad=.true.
          i=i+1
      enddo
      if (.not. lbad .and. ilen .ne. 4) then
          write(ilogunit,5)trim(adjustl(ayear))
          stop
      else
          read(ayear,*)modyear
          call leapyr(modyear)
      endif
c     get the requested run type 
      call getarg(2,runtyp1)
      
c     these are subject to change
c     convert to upper case
      ndash=0
      ilen=len_trim(runtyp1)
      do i=1,ilen
          call upcase(runtyp1(i:i))
          if (runtyp1(i:i) .eq. '_') ndash=ndash+1
      enddo
      
      if (ndash .eq. 0) then
          runtyp=trim(adjustl(runtyp1))
      elseif (ndash .eq. 1) then  !CONUS
          n1=index(runtyp1,'_')
          runtyp=trim(adjustl(runtyp1(1:n1-1)))
          read(runtyp1(n1+1:ilen),*)agrid
          call checkgrid(agrid)
          if (gridres .ne. 4 .and. gridres .ne. 12 .and. gridres .ne. 9 
     +    .and. gridres .ne. 3) then
              write(ilogunit,5)agrid
              stop
          endif
      elseif (ndash .eq. 2) then !non-CONUS
          n1=index(runtyp1,'_')
          n2=index(runtyp1(n1+1:ilen),'_')+n1
          runtyp=trim(adjustl(runtyp1(1:n1-1)))
          read(runtyp1(n1+1:n2-1),*)agrid
          call checkgrid(agrid)
          if (gridres .ne. 4 .and. gridres .ne. 12 .and. gridres .ne. 9 
     +    .and. gridres .ne. 3) then
              write(ilogunit,5)agrid
              stop
          endif
          read(runtyp1(n2+1:ilen),'(a)')noncon
          ilen=len_trim(noncon)
          do i=1,ilen
              call upcase(noncon(i:i))
          enddo
          if (noncon .eq. 'Y') then
              noncontr=.true.
          elseif (noncon .eq. 'N') then
              noncontr=.false.
          else
              write(ilogunit,5)noncon
              stop
          endif
      else
          write(ilogunit,*)runtyp1
          stop
      endif
      
      if (runtyp .eq. 'PT') then
          runpt=.true.
      elseif (runtyp .eq. 'APRUN' .or. runtyp .eq. 'APNON') then
          runap=.true.
      elseif (runtyp .eq. 'CMVP') then
          runport=.true.
      elseif (runtyp .eq. 'UNDER' .or. runtyp 
     +    .eq. 'CMVU') then
          runcmv=.true.
      elseif (runtyp .eq. 'HDON' .or. runtyp .eq. 'HDOFF' .or. 
     +    runtyp .eq. 'LDON' .or. runtyp .eq. 'LDOFF' .or. 
     +    runtyp .eq. 'HDOFF' .or. runtyp .eq. 'HOTEL') then 
          runon=.true.
      elseif (runtyp .eq. 'NONRD') then
          runnon=.true.
      elseif (runtyp .eq. 'NPHI') then
          runnphi=.true.
      elseif (runtyp .eq. 'NPLO') then
          runnplow=.true.
      elseif (runtyp .eq. 'RWC') then
          runrwc=.true.
      elseif (runtyp .eq. 'OILGAS') then
          runog=.true.
      elseif (runtyp .eq. 'AG') then
          runag=.true.
      else
          write(ilogunit,5)runtyp
          stop
      endif
      
c     read other arguments
      if (runpt .or. runap .or. runport) then ! .or. noncontr) then
c         get facility and 2-letter state from the command line
          call getarg(3,facid)
          call getarg(4,state2)
          ilen=len_trim(facid)
          do i=1,ilen
              call upcase(facid(i:i))
          enddo
              
      else
c         get column, row, grid resolution, and state
c         note for non-CONUS tract sources, column and row will be 0
          call getarg(3,acol)
          call getarg(4,arow)
          if (noncontr) call getarg(5,state2)
c          call getarg(4,agrid)
          read(acol,*)col
          read(arow,*)row
c          if (col .gt. 0 .and. noncontr) then
c              write(*,7)'col'
c              col=0
c          endif
c          if (row .gt. 0 .and. noncontr) then
c              write(*,7)'row'
c              row=0
c          endif
c          read(agrid,*)gridres
c          call getarg(5,state2)  !used for non-CONUS (for CONUS set to US)
c          call getarg(6,facid)  !only non-CONUS
      endif
      
      do i=1,2
          call upcase(state2(i:i))
      enddo

c    change when get new cmv; may have 85 for CMV but 85 for point is in
c    in gulf
      if (state2 .eq. 'AK' .or. state2 .eq. 'HI' .or. state2 .eq. 'PR' 
     +.or. state2 .eq. 'VI' .or. state2 .eq. '02' .or. state2 .eq. '15'
     + .or. state2 .eq. '72' .or. state2 .eq. '78') akhiprvi=.true. ! .or. 
c    + state2 .eq. '85') akhiprvi=.true. !85 is underway near VI
      
C     column and row offsets for AK, everywhere else is 0
      if (state2 .eq. 'AK' .or. state2 .eq. '02') then
          coffset=96
          roffset=63
      else
          coffset=0
          roffset=0
      endif
      
c      if (akhiprvi .and. (runon .or. runnon .or. runnphi .or. runnplow 
c     +.or. runrwc .or. runog)) then
      if (noncontr) then
c          noncontr=.true.
          if (state2 .eq. 'AK') then
              ncstfip='02'
          elseif (state2 .eq. 'HI') then
              ncstfip='15'
          elseif (state2 .eq. 'PR') then
              ncstfip='72'
          else
              ncstfip='78'
          endif
      endif

c     create facility ID for gridded sources based on column and row
      if (runon .or. runnon .or. runnphi .or. runnplow .or. runrwc .or. 
     +runog .or. runag .or. runcmv) then ! .and. .not. noncontr) then
          write(acol,'(i3.3)')col
          write(arow,'(i3.3)')row
          write(facid,'(4a)')'G',trim(adjustl(acol)),'R',
     +    trim(adjustl(arow))
      endif

      write(logfil,'(2(a))')trim(adjustl(facid)),'_inp.txt'
      open(unit=ilogunit,file=logfil,status='replace')
      write(ilogunit,6)versn


 4    format(/'Order of arguments'//'4-digit year'/'Run type'//
     +'Depending on run type the next two arguments are a choice of:'//
     +'Facility ID (for point, airports, CMV, or tract sources)'/
     +'State abbreviation or FIPS code (for point, airports, CMV, or ',
     +'tract sources)'/
     +'CMAQ column (for gridded or tract sources, 0 for tracts)'/
     +'CMAQ row (for gridded or tract sources, 0 for tracts)'//
     +'Valid run types:'/'pt',13x,'point'/'aprun',10x,'airport runways'/
     +'apnon',10x,'airport area sources'/'port',11x,'CMV ports'/
     +'under',10x,'CMV underway'/'cmv',12x,'CMV'/'hdon_4_y',7x,
     +'HDON 4 km tract'/'hdon_4_n',7x,'HDON 4 km gridded'/
     +'ldon_4_y',7x,'LDON 4 km tract'/
     +'ldon_4_n',7x,'LDON 4 km gridded'/
     +'hotel_4_y',6x,'HOTEL 4 km tract'/'hotel_4_n',6x,
     +'HOTEL 4 km gridded'/'hdoff_12_y',5x,'HDOFF 12 km tract'/
     +'hdoff_12_n',5x,'HDOFF 12 km gridded'/
     +'ldoff_12_y',5x,'LDOFF 12 km tract'/
     +'ldoff_12_n',5x,'LDOFF 12 km gridded'/
     +'nonrd_12_y',5x,'NONRD 12 km tract'/'nonrd_12_n',5x,
     +'NONRD12 km gridded'/'nphi_12_y',6x,'NPHI 12 km tract'/
     +'nphi_12_n',6x,'NPHI 12 km gridded'/
     +'nplo_12_y',6x,'NPLO 12 km tract'/'nplo_12_n',6x,
     +'NPLO 12 km gridded'/'rwc_12_y',7x,'RWC 12 km tract'/'rwc_12_n',
     +7x,'RWC 12 km gridded'/'oilgas_4_y',5x,'OILGAS 4 km tract'/
     +'oilgas_4_N',5x,'OILGAS 4 km gridded'/'ag_12_n',4x,
     +'AG 12 km griddd')
      
 5    format(a,'is not valid, stopping program')
 6    format('Beginning aermod inputs program version: ',a5)
 7    format('Resetting ',a3,' to 0')
      return
      end
c##########################################################################
c     subroutine to write AERMOD.INP file for specific source
      subroutine writeinp
      use main1
      implicit none

c     ifile:      integer file unit for plotfile
c     igrp:       group loop counter
c     i:          loop counter



      integer ifile,igrp,i,is,m,ngrprec,ijunk
     

c     infil:      AERMOD input file
c     errfil:     AERMOD error file

c     outfil:     AERMOD plot file, currently set to plotfile.plt
c     monthfil:    AERMOD seasonal/hour output file for onroad sources
c     grp:        AERMOD source groups
c     src1:       AERMOD source ID for specific AERMOD run group (point, airport, etc.)
c     flat:       AERMOD keyword for flat terrain
      character infil*250,errfil*250,outfil*300,
     +grp*8,src1*12,monthfil*300,flat*4,avgtim*12,outag*2,hrfil*250
      
c     lurban:     logical variable denoting if source is an urban source
c                 true=urban,false=rural
      logical lurban,lexist

           
c     get directories and filenames of locations, parameters, temporal files,
c     meteorological files and receptors
      call getdir
 
c     input filename
c      if (runpt .or. runap .or. runcmv) then
      if (runpt .or. runap .or. runport) then
          write(infil,'(3(a))')'aermod_',trim(adjustl(facid)),
     +    '.inp'
      else
          write(infil,'(3(a),i2.2,3(a))')'aermod_',
     +    trim(adjustl(runtyp)),'_',gridres,'_',trim(adjustl(facid)),
     +    '.inp'
      endif
      
c     errors file
c      if (runpt .or. runap .or. runcmv) then
      if (runpt .or. runap .or. runport) then
          write(errfil,'(5(a))')trim(adjustl(outdir)),
     +    trim(adjustl(state2)),'/',trim(adjustl(facid)),'_errors.out'
      elseif (noncontr) then
          write(errfil,'(5(a),i2.2,3(a))')trim(adjustl(outdir)),
     +    trim(adjustl(state2)),'/',trim(adjustl(runtyp)),'_',gridres,
     +    '_',trim(adjustl(facid)),'_errors.out'
      else
          write(errfil,'(3(a),i2.2,3(a))')trim(adjustl(outdir)),
     +    trim(adjustl(runtyp)),'_',gridres,'_',trim(adjustl(facid)),
     +    '_errors.out'
      endif
      
c     set outag to output annual plot files or not
      if (runon .or. runnon) then
         outag='**'
      else
         outag='OU'
      endif
c     output and seasonal filenames
c     note, seasonal file only written to AERMOD.INP for onroad AERMOD runs
c      if (runpt .or. runap .or. runcmv) then
      if (runpt .or. runap .or. runport) then
          write(outfil,'(5(a))')trim(adjustl(outdir)),
     +    trim(adjustl(state2)),'/',trim(adjustl(facid)),'.plt'
          write(monthfil,'(5(a))')trim(adjustl(outdir)),
     +    trim(adjustl(state2)),'/',trim(adjustl(facid)),'_month.pst'
      elseif (noncontr) then
          write(outfil,'(5(a),i2.2,3(a))')trim(adjustl(outdir)),
     +    trim(adjustl(state2)),'/',trim(adjustl(runtyp)),'_',gridres,
     +    '_',trim(adjustl(facid)),'.plt'
          write(monthfil,'(5(a),i2.2,3(a))')trim(adjustl(outdir)),
     +    trim(adjustl(state2)),'/',trim(adjustl(runtyp)),'_',gridres,
     +    '_',trim(adjustl(facid)),'_month.pst'
      else
          write(outfil,'(3(a),i2.2,3(a))')trim(adjustl(outdir)),
     +    trim(adjustl(runtyp)),'_',gridres,'_',trim(adjustl(facid)),
     +    '.plt'
          write(monthfil,'(3(a),i2.2,3(a))')trim(adjustl(outdir)),
     +    trim(adjustl(runtyp)),'_',gridres,'_',trim(adjustl(facid)),
     +    '_month.pst'
      endif
       
c     make sure sources are same in location, parameter, and temporal files
      call checksrc
    
c     assign source parameters based on run type
      srcnam='NA'
      flat=''
      if (runpt) then
          call readpoint
      elseif (runap) then
          call readairport
!      elseif (runcmv) then
      elseif (runport) then
          call readport
          flat='FLAT'
      else
          call readgrid !gridded for NON-CONUS v2
c          if (noncontr) then
c              call readtract
c          else
c              call readgrid
c         endif
          flat='FLAT'
      endif

      if (runon .or. runnon) then
          avgtim='ANNUAL MONTH'
      else
          avgtim='ANNUAL' 
      endif
      
c     begin writing AERMOD.INP file
c     open AERMOD.INP file
c      open(unit=aerunit,file='aermod.inp',status='unknown')
      open(unit=aerunit,file=infil,status='unknown')
      ifile=35
      
c     control pathway
      write(aerunit,5)'CO'
c     write facility ID, facility name, column, row, UTM ZONE
c      if (akhiprvi .and. (runpt .or. runcmv .or. runap)) then
c          col=0
c          row=0
c      endif
      
      write(aerunit,10)modyear,trim(adjustl(facid)),
     +trim(adjustl(srcnam)),col,row,utmzone,flat,avgtim,
     +trim(adjustl(errfil)) !title line
      
c     determine if source is an urban source and write URBANOPT keyword if so
c     keyword written in subroutine urbrur
      
      call urbrur(lurban)
     
      write(aerunit,15)'CO'

c     source pathway
      write(aerunit,5)'SO'
      write(aerunit,20)  !source elevations in meters

c     write locations and parameters
      call locs_params
      
     
c     if AERMOD emissions factors
      if (ltemp(1) .and. .not. runrwc .and. .not. runag .and.  
     +((.not. runon .and. .not. noncontr) .or. (noncontr .and. .not. 
     +runrwc))) call emisfact1
c     if hourly files, write appropriate hourly files and HOUREMIS keyword
      if (ltemp(2)) call writehour
    
c     RWC or AG county hourly factors 
      if (runrwc .or. runag) call rwc2grid

c     onroad hourly emissions based on county
      if (runon) call onrd2grid
      
      if (runon .or. runnon) then !modified 17354 for no tract sources
          do igrp=1,ngrps
              write(aerunit,65)gridded(igrp)%srcid,
     +        (gridded(igrp)%monemis(is),is=1,6)
              write(aerunit,65)gridded(igrp)%srcid,
     +        (gridded(igrp)%monemis(is),is=7,12)
c              if (noncontr) then
c                  write(aerunit,65)tract(igrp)%srcid,
c     +            (tract(igrp)%monemis(is),is=1,6)
c                  write(aerunit,65)tract(igrp)%srcid,
c     +            (tract(igrp)%monemis(is),is=7,12)
c              else
c                  write(aerunit,65)gridded(igrp)%srcid,
c     +            (gridded(igrp)%monemis(is),is=1,6)
c                  write(aerunit,65)gridded(igrp)%srcid,
c     +            (gridded(igrp)%monemis(is),is=7,12)
c              endif
          enddo
      endif
      
      if (runport) then
       write(hrfil,'(2a)')trim(adjustl(facid)),'_hourly.dat'
        do igrp=1,ngrps
          write(aerunit,80)trim(adjustl(hrfil)),cmv(igrp)%srcid
        enddo
      endif

      if (runcmv) then
       write(hrfil,'(2a)')trim(adjustl(facid)),'_hourly.dat'
        do igrp=1,ngrps
          write(aerunit,80)trim(adjustl(hrfil)),gridded(igrp)%srcid
        enddo
      endif

c     write annual emissions to screen
      call writeann
c     if urban write that all sources are urban
      
      if (lurban) write(aerunit,30)

c     write source groups, source groups are same as source IDs
c     for point and gridded sources
c     for airports, tracts and CMV, source group ALL is used
      if (runap )then ! .or. runcmv) then ! .or. noncontr) then
          write(aerunit,55)
c      elseif (runcmv) then 
      elseif (runport) then
         do igrp=1,ngrps
             ijunk=index(cmv(igrp)%srcid,'c')
             src1=cmv(igrp)%srcid(ijunk:ijunk+1)
             do ijunk=1,2
                call upcase(src1(ijunk:ijunk))
             enddo
             write(aerunit,35)trim(adjustl(src1)),
     +       trim(adjustl(cmv(igrp)%srcid))
         enddo
      else
          do igrp=1,ngrps
              call get_source(igrp,src1)
              write(aerunit,35)trim(adjustl(src1)),
     +        trim(adjustl(src1))
          enddo
      endif
      
      write(aerunit,15)'SO'
 
c     receptor pathway
      write(aerunit,5)'RE'

      call receptors
      write(aerunit,15)'RE'
     
c     met pathway     
c     write filenames
c     NON-CONUS now has gridded met
c     writemet will write the entire met pathway
      call writemet
 

c     output pathway
      write(aerunit,5)'OU'
      write(aerunit,45)
c     write source group outputs, all going to same file
      if (runap ) then ! .or. runcmv) then ! .or. noncontr) then
          write(aerunit,50)outag,'ALL',trim(adjustl(outfil)),
     +    ifile
c         if (noncontr .and. (runon .or. runnon)) 
c    +    write(aerunit,60)'ALL',trim(adjustl(monthfil)),ifile+1
      else
          do igrp=1,ngrps
              if (runpt) then 
                  grp=facility(igrp)%srcid
c              elseif (runcmv) then
              elseif (runport) then
                  ijunk=index(cmv(igrp)%srcid,'c')
                  grp=cmv(igrp)%srcid(ijunk:ijunk+1)
                  do ijunk=1,2
                     call upcase(grp(ijunk:ijunk))
                  enddo
              else
                  grp=gridded(igrp)%srcid
              endif
              write(aerunit,50)outag,trim(adjustl(grp)),
     +        trim(adjustl(outfil)),ifile
              if (runon .or. runnon) write(aerunit,60)trim(adjustl(grp))
     +        ,trim(adjustl(monthfil)),ifile+1
          
          enddo
      endif
      
      write(aerunit,15)'OU'
      ngrprec=ngrps*nrec
      write(aerunit,70)trim(adjustl(facid)),ngrps,nrec,ngrprec
      if (ngrprec .gt. 1000000)write(aerunit,71)
      close(aerunit)
      
 5    format(a,' STARTING')
c version 17200, take beta out of MODELOPT since using 16216 and MMIF is not beta
 10   format('CO TITLEONE ',i4,1x,a,1x,a,/'CO TITLETWO ','COL ',i3,
     +' ROW ',i3,' UTM ZONE ',i2/'CO MODELOPT CONC FASTALL',1x,a/
     +'CO AVERTIME',1x,a/'CO POLLUTID OTHER'/'CO RUNORNOT RUN'/
     +'CO ERRORFIL'1x,a)
c 10   format('CO TITLEONE ',a,1x,a,/'CO TITLETWO ','COL ',i3,' ROW ',i3,
c     +' UTM ZONE ',i2/'CO MODELOPT CONC BETA FASTALL',1x,a/
c     +'CO AVERTIME ANNUAL'/'CO POLLUTID OTHER'/'CO RUNORNOT RUN'/
c     +'CO ERRORFIL'1x,a)
 15   format(a,' FINISHED'/)
 20   format('SO ELEVUNITS METERS')
 25   format('SO INCLUDED ',a)
 30   format('SO URBANSRC ALL')
 35   format('SO SRCGROUP',2(1x,a))
**NEED TO UPDATE TO 2014
 
 45   format('OU FILEFORM EXP')
 55   format('SO SRCGROUP ALL')
 50   format(a2,1x,'PLOTFILE ANNUAL ',a8,1x,a,1x,i2)
c 60   format('OU SEASONHR ',a8,1x,a,1x,i2)
 60   format('OU POSTFILE MONTH ',a8,1x,'PLOT',1x,a,1x,i2)
 65   format('** MONTHEMIS',1x,a12,6(1x,f10.4))
 70   format('** NGRPS*NSRC',1x,a,1x,i3,1x,i6,i9)
 71   format('** GROUP-RECEPTOR EXCEEDS 1 MILLION, SPLIT?')
 80   format('SO HOUREMIS ',a,1x,a12)
      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 eof,n,i,np

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,tag*20,fname*250,acol*3,arow*3,ares*2,
     +pfil(2)*300,lfil*300,efil*300,tfil*300,cntyfil1*300

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,llocfil,lelev,lpfil,ltemfil,lcnty1
      
c     initialize variables
      eof=0
      lloc=.false.
      lparam=.false.
      ltemp=.false.
      lcmaq=.false.
      lmet=.false.
      lrec=.false.
      lparam1=.false.
      ltemp1=.false.
      lout=.false.
      lcnty=.false.
      lxwalk=.false.
      llocfil=.false.
      lelev=.false.
      lpfil=.false.
      ltemfil=.false.
      lparam=.false.
      ltemp=.false.
      paramfil='NA'
      temfil='NA'
      locfil='NA'
      elevfil='NA'
      tfil='NA'
      pfil='NA'
      metdir='NA'
      np=0
      lcnty1=.false.
      
      

      
c     check to see if input_directories.dat is in directory
      inquire(file='input_directories.dat',exist=lexist)
      if (lexist) then
          open(unit=inpunit,file='input_directories.dat',status='old')
      else
          write(ilogunit,*)'input_directories.dat does not exist'
          write(ilogunit,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
              n=index(line,'XWALK ')
              xdir=line(n+6:300)
              lxwalk=.true.
          endif
          if (index(line,'OUTPUT ') .gt. 0) then
              n=index(line,'OUTPUT ')
              outdir=line(n+7:300)
              lout=.true.
          endif
          if (index(line,'ELEVFILE ') .gt. 0) then
              n=index(line,'ELEVFILE ')
              efil=line(n+9:300)
              lelev=.true.
              write(elevfil,'(2a)')trim(adjustl(ptcmaq)),
     +        trim(adjustl(efil))
          endif
          if (index(line,'LOCFILE ') .gt. 0) then
              n=index(line,'LOCFILE ')
              lfil=line(n+8:300)
              llocfil=.true.
          endif
          if (index(line,'PARAMFIL ') .gt. 0) then
              np=np+1
              n=index(line,'PARAMFIL ')
              pfil(np)=line(n+9:300)
              lpfil=.true.
          endif
          if (index(line,'TEMPFILE ') .gt. 0) then
              ltemfil=.true.
              n=index(line,'TEMPFILE ')
              tfil=line(n+9:300)
          endif
          if (index(line,'CNTYFILE ') .gt. 0) then
              lcnty1=.true.
              n=index(line,'CNTYFILE ')
              cntyfil1=line(n+9:300)
          endif
          read(11,5,iostat=eof)line
          goto 10
      endif
      
c     if any keyword not found, stop program
      if (.not. lloc) write(ilogunit,15)facid
      if (.not. lparam1) write(ilogunit,20)facid    
      if (.not. ltemp1) write(ilogunit,25)facid
      if (.not. lcmaq) write(ilogunit,30)facid
      if (.not. lmet) write(ilogunit,35)facid
      if (.not. lrec) write(ilogunit,40)facid
      if (.not. lout) write(ilogunit,45)facid
      if (.not. lxwalk) write(ilogunit,46)facid
      if (.not. lloc .or. .not. lparam1 .or. .not. ltemp1 .or. 
     +.not. lcmaq .or. .not. lmet .or. .not. lrec .or. .not. lout) then
          write(ilogunit,50)
          stop
      endif

      
c     get location, parameter, and temporal files
!      write(locfil,'(4a)')trim(adjustl(ptloc)),'point_',
!     +trim(adjustl(state2)),'_locations_elevations.csv'
      
      write(locfil,'(2a)')trim(adjustl(ptloc)),trim(adjustl(lfil))
      inquire(file=locfil,exist=lloc)

      
      if (.not. runpt) then
c         parameter file
          write(paramfil(1),'(2(a))')trim(adjustl(ptparam)),
     +    trim(adjustl(pfil(1)))
          inquire(file=paramfil(1),exist=lparam(1))
      endif
c     AIRPORTS or CMV temporal
      if (runap) then! .or. runcmv) then    
c         temporal
          write(temfil(1),'(2(a))')trim(adjustl(pttemp)),
     +    trim(adjustl(tfil))
          inquire(file=temfil(1),exist=ltemp(1))
      endif
c     POINT SOURCES
      if (runpt) then
          
c         need to determine how many parameter files there are
          do i=1,2
              np=0
              call checkparam(pfil(i),np)
              if (np .gt. 0) then
                  write(paramfil(np),'(2(a))')trim(adjustl(ptparam)),
     +            trim(adjustl(pfil(i)))
                  inquire(file=paramfil(np),exist=lparam(np))
              endif
c              if (i .eq. 1) then
c                  tag='point_point' !subject to change
c              else
c                  tag='point_fug'
c              endif
c              write(fname,'(4a)')trim(adjustl(ptparam)),
c     +        trim(adjustl(tag)),'_','srcparam.csv'
c              inquire(file=fname,exist=lparam(i))
c              if (lparam(i)) paramfil(i)=fname
          enddo
c         need to determine how many temporal files there are
c         an EGU may have both hourly and temporal factors.
c         non-EGU will only have non-hourly
c         first check for emissions factor file
          write(fname,'(2(a))')trim(adjustl(pttemp)),trim(adjustl(tfil))  
          inquire(file=fname,exist=ltemp(1))
          if (ltemp(1)) temfil(1)=fname   
c         now check for presence of hourly emissions file
          write(fname,'(5(a))')trim(adjustl(pttemp)),
     +    trim(adjustl(facid)),'_',trim(adjustl(state2)),'_hourly.csv'
          inquire(file=fname,exist=ltemp(2))
          if (ltemp(2)) temfil(2)=fname
      endif

c     GRIDDED (onroad, nonroad, nonpoint, RWC)
      if (runon .or. runnon .or. runnphi .or. runnplow .or. runrwc 
     +    .or. runog .or. runag) then
          write(ares,'(i2)')gridres
          if (.not. noncontr) then
c             temporal (county hourly file for onroad and RWC)
              if (runrwc .or. runon .or. runag) then
                  write(cntyfil,'(2(a))')trim(adjustl(xdir)),
     +            trim(adjustl(cntyfil1))
                  inquire(file=cntyfil,exist=lcnty)
                  call check_county(ares)
c              elseif (runon) then
c                  write(cntyfil,'(4a)')trim(adjustl(xdir)),
c     +            trim(adjustl(runtyp)),trim(adjustl(ares)),
c     +            '_county-to-gridcell.csv'
c                  inquire(file=cntyfil,exist=lcnty)
c                  call check_county(ares)
              else
                  write(temfil(1),'(2(a))')trim(adjustl(pttemp)),
     +            trim(adjustl(tfil))
              endif
          else
              lcnty=.true.
c             temporal (county hourly file for onroad and RWC)
              if (runrwc .or. runag .or. runon) then
                  write(cntyfil,'(2(a))')trim(adjustl(xdir)),
     +            trim(adjustl(cntyfil1))
c                  write(cntyfil,'(4a)')trim(adjustl(xdir)),
c     +            trim(adjustl(runtyp)),trim(adjustl(ares)),
c     +            '_county-to-gridcell.csv'
c                  inquire(file=temfil(1),exist=ltemp(1))
                  call check_county(ares)
c               elseif (runon) then
c                   write(cntyfil,'(2(a))')trim(adjustl(xdir)),
c     +            trim(adjustl(cntyfil1))
c                  write(cntyfil,'(4a)')trim(adjustl(xdir)),
c     +            trim(adjustl(runtyp)),trim(adjustl(ares)),
c     +            '_county-to-gridcell.csv'
c                  inquire(file=cntyfil,exist=lcnty)
c                  call check_county(ares)
              else
                  write(temfil(1),'(2(a))')trim(adjustl(pttemp)),
     +            trim(adjustl(tfil))
              endif
          endif
          inquire(file=temfil(1),exist=ltemp(1))
      endif

      
      
c     make sure files exist
c     reset lloc, lparam, and ltemp to false

      
c     if files don't exist, stop program
      if (runpt .and. (.not. lparam(1) .or. .not. lparam(2))) 
     +write(ilogunit,75)facid
      if (.not. lloc) write(ilogunit,55)facid
      if (.not. lparam(1) .and. .not. lparam(2)) write(ilogunit,60)facid    
      if (.not. runon .and. .not. runrwc .and. .not. runag .and. .not. 
     +runport .and. .not. runcmv .and. ltemp(1) .and. .not. 
     + ltemp(2)) write(ilogunit,65)facid
      if ((runrwc .or. runon .or. runag) .and. .not. lcnty) 
     +write(ilogunit,70)facid
      if (.not. lloc .or. (.not. lparam(1) .and. .not. lparam(2)) .or. 
     + (.not. runon .and. .not. runrwc .and. .not. runag .and.  
     + .not. runport .and. .not. runcmv .and. .not. ltemp(1) .and. 
     +.not. ltemp(2)) .or.  
     +((runrwc .or. runon .or. runag) .and. .not. lcnty)) then
          write(ilogunit,50)
          stop
      endif
      
 5    format(a300)
 15   format('PTLOC keyword not found for ',a)
 20   format('PTPARAM keyword not found for ',a)
 25   format('PTEMP keyword not found for ',a)
 30   format('PTCMAQ keyword not found for ',a)
 35   format('METDIR keyword not found for ',a)
 40   format('RECDIR keyword not found for ',a)
 45   format('OUTDIR keyword not found for ',a)
 46   format('XWALK keyword not found for ',a)
 50   format('Aborting run..')
 55   format('location file not found for ',a)
 60   format('parameter file(s) not found for ',a)
 65   format('temporal file(s) not found for ',a)   
 70   format('County file not found for ',a)
 75   format('One or more paramter files missing for ',a)
      return
      end

c##########################################################################
c     subroutine to check point parameter files
      subroutine checkparam(fname,np)
      use main1
      implicit none
      integer np,eof
      character fname*300,fname1*300,line*300
      logical lexist
      write(fname1,'(2(a))')trim(adjustl(ptparam)),trim(adjustl(fname))
      inquire(file=fname1,exist=lexist)
      if (.not. lexist) then
          write(ilogunit,5)trim(adjustl(fname1)),trim(adjustl(facid))
          stop
      endif
      open(unit=parunit,file=fname1,status='old')
c     read first line
      read(parunit,10,iostat=eof)line
      if (eof .eq. 0) then
          if (index(line,'szinit') .gt. 0) then
              np=2 !fugitives
          elseif (index(line,'velocity') .gt. 0) then
              np=1 !stack
          else
              np=0
          endif
      endif
      close(parunit)
 5    format(a,' does not exist for facility ',a)
 10   format(a300)
      return
      end
c##########################################################################
c     subroutine to check all sources are in all location, parameter, and
c     temporal files
      subroutine checksrc
      use main1
      implicit none

c     nlocs:      number of sources read from locations file
c     nparams:    number of sources read from parameters file(s)
c     ntemp:      number of unique temporal sources
c     ilocs:      locations loop counter
c     iparams:    parameters loop counter
c     itemp:      temporal loop counter
c     eof:        end of file indicator
c     i,j,k:      loop counters
      integer nlocs,nparams,ntemp,ilocs,iparams,itemp,eof,i,j,k

c     ldup:       logical variable denoting duplicates found in temporary
c                 temporal array, tsrc
c     lmatch:     logical array denoting that a source in one source ID
c                 array is found in other two source ID arrays
c     lloc:       logical variable denoting facility found in location file
c     lparam2:    logical variable denoting facility is in one or more
c                 parameter files
c     ltemp1      logical variable denoting facility is in temporal file
c     lparam1     placeholder variable
      logical ldup,lmatch(2),lloc,lparam2,ltemp2,lparam1,ltemp1
      
      nlocs=0
      nparams=0
      ntemp=0
      ilocs=0
      iparams=0
      itemp=0
      eof=0
      lparam1=.true. !not used, placeholder only
      
c     read location file and get # of sources
      call readfil(locunit,locfil,nlocs,1,lloc)
      if (.not. lloc) then
          write(ilogunit,5)trim(adjustl(facid)),trim(adjustl(locfil))
      else
          rewind(locunit)
          allocate(locsrc(nlocs))
          ilocs=0
          call readfil(locunit,locfil,ilocs,2,lloc)
          ilocs=0
      endif
      
      
c     read in parameter file(s)
c     note that lparam is being reset from the values returned in
c     subroutine getdir
  
      do i=1,2
          if (trim(adjustl(paramfil(i))) .ne. 'NA') then
              call readfil(parunit,paramfil(i),nparams,1,lparam(i))
              close(parunit)
          else
              lparam(i)=.false.
          endif
      enddo
c     check to see if # of parameter sources = # of location sources
c     if not, stop program
      if (nparams .ne. nlocs) then
          write(ilogunit,30)nlocs,trim(adjustl(locfil)),nparams
          stop
      endif
      
      if ((lparam(1) .and. lparam(2)) .or. (.not. lparam(1) .and. 
     +lparam(2)) .or. (lparam(1) .and. .not. lparam(2))) then
          lparam2=.true.
      else
          lparam2=.false.
      endif
      

c     assign sources to array
      if (lparam2) then
          allocate(parsrc(nparams))
!          allocate(srctyp(nparams))
          iparams=0
          do i=1,2
              if (lparam(i)) then
                  open(unit=parunit,file=paramfil(i),status='old')
                  call readfil(parunit,paramfil(i),iparams,2,lparam1)
                  close(parunit)
              endif
          enddo
      endif

c     read in temporal file(s) for every run group except onroad ag, and RWC which use a county cross reference file
c     note that ltemp is being reset from the values returned in
c     subroutine getdir      
c     emissions factor file always first for point
c     only file for other run groups, other than RWC and ag
c      if (runpt .or. runap .or. runcmv .or. runnplow .or. runnphi .or.
      if (runpt .or. runap .or. runnplow .or. runnphi .or.
     +runnon .or. runog) then ! .or. (runon .and. noncontr)) then
          call readfil(temunit,temfil(1),ntemp,1,ltemp(1))
          if (.not. ltemp(1)) write(ilogunit,5)trim(adjustl(facid)),
     +    trim(adjustl(temfil(1)))
          close(temunit)
      
c         now check hourly emissions file for # of sources for point sources, i.e. EGU's
          if (runpt) then
              if (ltemp(2)) then
                  call checkhourly(temfil(2),temunit,ntemp,1,ltemp(2))
                  if (.not. ltemp(2)) write(ilogunit,5)
     +            trim(adjustl(facid)),trim(adjustl(temfil(1)))
              endif
          endif
          allocate(temsrc(ntemp))
          ltemp(3)=.true.
          if (.not. ltemp(1) .and. .not. ltemp(2)) ltemp(3)=.false.
      
          if (ltemp(3)) then
              if (ltemp(2)) then  !only hourly data
                  itemp=0
!              allocate(hrsrc(nhrsrc))
                  call checkhourly(temfil(2),temunit,itemp,2,ltemp1)
              endif
              if (ltemp(1)) then
                  open(unit=temunit,file=temfil(1),status='old')
                  call readfil(temunit,temfil(1),itemp,2,ltemp1)
              endif
          else
              write(ilogunit,65)trim(adjustl(facid))
              stop
          endif
      endif    

c     check the county cross reference file for RWC and onroad to 
c     make sure grid cell is in file
c      if (runrwc) then
c          call check_county
c          if (.not. lcnty) then
c              write(*,70)trim(adjustl(facid)),trim(adjustl(cntyfil))
c              stop
c          endif
c      endif
      
c     now check to make sure same sources are in each array 
c     start with locations

      do i=1,nlocs 
          lmatch=.false.
          do j=1,nparams
              if (locsrc(i) .eq. parsrc(j))lmatch(1)=.true.
          enddo
c          if (runpt .or. runap .or. runcmv .or. runnplow .or. runnphi 
          if (runpt .or. runap .or. runport .or. runnplow .or. runnphi
     +    .or. runnon .or. runog ) then !.or. (runon .and. noncontr)) then
              do k=1,ntemp
                  if (locsrc(i) .eq. temsrc(k))lmatch(2)=.true.
              enddo
          else
              lmatch(2)=.true. !for if-then statement below
          endif
      enddo
      
      if (.not. lmatch(1) .and. .not. lmatch(2)) then
          write(ilogunit,60)trim(adjustl(facid))
          stop
      endif
      
      do i=1,nparams 
          lmatch=.false.
          do j=1,nlocs
              if (locsrc(j) .eq. parsrc(i))lmatch(1)=.true.
          enddo
c          if (runpt .or. runap .or. runcmv .or. runnplow .or. runnphi 
          if (runpt .or. runap .or. runport .or. runnplow .or. runnphi
     +    .or. runnon .or. runog) then ! .or. (runon .and. noncontr)) then
              do k=1,ntemp
                  if (parsrc(i) .eq. temsrc(k))lmatch(2)=.true.
              enddo
          else
              lmatch(2)=.true.
          endif
      enddo 
      
      if (.not. lmatch(1) .and. .not. lmatch(2)) then
          write(ilogunit,60)trim(adjustl(facid))
          stop
      endif
      
c      if (runpt .or. runap .or. runcmv .or. runnplow .or. 
      if (runpt .or. runap .or. runnplow .or.
     +runnphi .or. runnon .or. runog) then ! .or. (runon .and. noncontr)) then
          do i=1,ntemp
              lmatch=.false.
              do j=1,nlocs
                  if (locsrc(j) .eq. temsrc(i))lmatch(1)=.true.
              enddo
              do k=1,nparams
                  if (temsrc(i) .eq. parsrc(k))lmatch(2)=.true.
              enddo
          enddo 
      
          if (.not. lmatch(1) .and. .not. lmatch(2)) then
              write(ilogunit,60)trim(adjustl(facid))
              stop
          endif
      endif
      
      ngrps=nlocs
      
c     allocate the arrays
      if (runpt) then
          allocate(facility(ngrps))
      elseif (runap) then
          allocate(airport(ngrps))
c      elseif (runcmv) then
      elseif (runport) then
          allocate(cmv(ngrps))
c      else if (noncontr) then
c          allocate(tract(ngrps))
      else
          allocate(gridded(ngrps))
      endif
      
      
      
      
  5   format('Facility ',a,' not found in'/,a)
30    format('Number of sources ',i3,' in',a,' location file is not ',
     +'equal to'/'number of sources ',i3,' in parameter files'/
     +'stopping program')
40    format('Number of sources ',i3,' in',a,' location file is not ',
     +'equal to'/'number of sources ',i3,' in parameter files or is',
     + ' not  equal to number of sources',i3,' in temporal files'/
     +'stopping program')
 60   format('Sources do not match in files for ',a,' stopping program')
 65   format('No temporal files found for ',a/'Stopping program')
 70   format(a,' is not listed in'/,a,/'Stopping program')
      return
      end
c##########################################################################  
c     read locations, parameter, or temporal factor file to get # of sources
      subroutine readfil(iunit,infil,nval,iswitch,lfound)
      use main1
      implicit none
      
c     iunit:      file unit being read
c     eof:        end of file indicator
c     nval:       number of sources counter
c     iswitch:    switch determining if subroutine is getting
c                 # of sources or writing to arrays
      integer iunit,eof,nval,iswitch

c     infil:      filename being read
c     src:        AERMOD source id being read
c     j1:         junk variable being read (usually source name)
c     facid1:     facility ID being read from file
c     header:     header record of file (not used)
c     fips1:      5-digit FIPS code
      character infil*300,src*12,src1*12,j1*75,facid1*20,header*100,
     +st*2,fips1*5,fips6*6,poly*4,tract1*11,facid2*20
      
c     lfound:     variable denoting facility found in file
      logical lfound
      
      lfound=.false.
      eof=0
      src1='NA'
      if (iswitch .eq. 1) then
          open(unit=iunit,file=infil,status='old')
      else
          rewind(iunit)
      endif
      
      read(iunit,*,iostat=eof)header
      if (runpt .or. runap) then 
          if (iunit .eq. locunit) then
              read(iunit,*,iostat=eof)st,facid1,j1,src
          else
              read(iunit,*,iostat=eof)facid1,j1,src
          endif
      elseif (runcmv) then
          if (iunit .eq. locunit) then
              read(iunit,*,iostat=eof)j1,fips1,facid1,src 
          else
              read(iunit,*,iostat=eof)j1,facid1,src
          endif
      elseif (runport) then
          if (iunit .eq. parunit) then
              read(iunit,*,iostat=eof)st,facid1,src
          else
              read(iunit,*,iostat=eof)st,fips1,facid1,src
          endif
c      elseif (noncontr) then
c          read(iunit,*,iostat=eof)j1,st,fips6,tract1,poly,facid2,src
      else
          read(iunit,*,iostat=eof)j1,facid1,src
      endif
      
 10   if (eof .eq. 0) then
cc          if ((.not. runcmv .and. .not. noncontr .and. facid1 .eq. 
cc     +    facid) .or. ((runcmv .or. noncontr) .and. src .eq. facid)) 
cc     +    then
c         if ((.not. runcmv .and. facid1 .eq. 
c    +    facid) .or. (runcmv .and. src .eq. facid)) 
c         if (facid1 .eq. facid) then
c    +    then
          if (facid1 .eq. facid) then
              if (runport) then
                  if (trim(adjustl(src1)) .ne. trim(adjustl(src))) 
     +            nval=nval+1
              else
                  nval=nval+1
              endif
              if (iswitch .eq. 1) then
                  lfound=.true.  !found facility
              else
c                 assign to arrays
                  if (iunit .eq. locunit) locsrc(nval)=src
                  if (iunit .eq. parunit) parsrc(nval)=src  
                  if (iunit .eq. temunit) temsrc(nval)=src
              endif
              
          endif
          src1=src
          if (runpt .or. runap) then 
              if (iunit .eq. locunit) then
                  read(iunit,*,iostat=eof)st,facid1,j1,src
              else
                  read(iunit,*,iostat=eof)facid1,j1,src
              endif
          elseif (runcmv) then
              if (iunit .eq. locunit) then
                 read(iunit,*,iostat=eof)j1,fips1,facid1,src
              else
                 read(iunit,*,iostat=eof)j1,facid1,src
              endif

          elseif (runport) then
              if (iunit .eq. parunit) then
                  read(iunit,*,iostat=eof)st,facid1,src
              else
                  read(iunit,*,iostat=eof)st,fips1,facid1,src
              endif
c          elseif (noncontr) then
c              read(iunit,*,iostat=eof)j1,st,fips6,tract1,poly,facid2,src
          else
              read(iunit,*,iostat=eof)j1,facid1,src
          endif
          goto 10
      endif
      if (iswitch .ne. 1) close(iunit)
      return
      end
c##########################################################################  
c     check county to grid cross reference file for RWC and onroad and get
c     temporal file for onroad
      subroutine check_county(ares)
      use main1
      implicit none
      
c     eof:        end of file indicator
      integer eof,isrc,ist,n,i,isrc1

c     facid1:     facility ID being read from file
c     j1:         junk variable being read (not used)
c     line:       header record (not used)

      character facid1*20,j1*10,line*100,ares*2,afips*5,src*12,st*2

      logical ldup
      
      allocatable :: st(:)
      lcnty=.false.
      isrc=0
      open(unit=cntyunit,file=cntyfil,status='old')
      read(cntyunit,*,iostat=eof)line
      read(cntyunit,*,iostat=eof)j1,afips,facid1,src
 5    if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              isrc=isrc+1
              lcnty=.true.
          endif
          read(cntyunit,*,iostat=eof)j1,afips,facid1,src
          goto 5
      endif     
      rewind(cntyunit)
      isrc1=isrc
      allocate(fips(isrc)) !isrc should be equal to ngrps
c      if (runon) then
      allocate(onsrc(isrc))
      allocate(st(isrc))
c      endif
      
      isrc=0
      read(cntyunit,*,iostat=eof)line
c      if (runon .and. gridres .eq. 4) then
      read(cntyunit,*,iostat=eof)j1,afips,facid1,src
c      else
c         read(cntyunit,*,iostat=eof)j1,afips,facid1
c      endif
 10   if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              isrc=isrc+1
              fips(isrc)=afips
              st(isrc)=afips(1:2)
              onsrc(isrc)=src
c             if (runon .and. gridres .eq. 4) onsrc(isrc)=src
          endif
c          if (runon .and. gridres .eq. 4) then
          read(cntyunit,*,iostat=eof)j1,afips,facid1,src
c          else
c              read(cntyunit,*,iostat=eof)j1,afips,facid1
c          endif
          goto 10
      endif
      
      if (.not. lcnty) then
          write(ilogunit,70)trim(adjustl(facid)),trim(adjustl(cntyfil))
          stop
      endif
      
c      if (runon .or. runrwc) then
      ist=1  
      do i=2,isrc1
          ldup=.false.
          n=1
          do while (n .lt. i .and. .not. ldup)
              if (st(i) .eq. st(n)) ldup=.true.
              n=n+1
          enddo 
          if (.not. ldup) ist=ist+1
      enddo
      nstate=ist

      allocate(stfips(nstate))
      stfips(1)=st(1)
      ist=1
      do i=2,isrc1
          ldup=.false.
          n=1
          do while (n .lt. i .and. .not. ldup)
              if (st(i) .eq. st(n)) ldup=.true.
              n=n+1
          enddo 
          if (.not. ldup) then
              ist=ist+1
              stfips(ist)=st(i)
          endif
      enddo               
c          stfips=fips(1:2)
c          write(temfil(1),'(6a)')trim(adjustl(pttemp)),
c    +        trim(adjustl(runtyp)),trim(adjustl(ares)),'_',stfips,
c    +        '_hourly.csv'
c      endif
      
c      do i=1,isrc1
c      write(*,'(3(1x,a))')fips(i),onsrc(i),st(i)
c      enddo
c      do i=1,nstate
c      write(*,*)stfips(i)
c      enddo

      close(cntyunit)
70    format(a,' is not listed in'/,a,/'Stopping program')
      return
      end
c########################################################################## 
c     subroutine to check hourly file for existence of facility
      subroutine checkhourly(infil,iunit,nval,iswitch,lfound)
      use main1
      implicit none
      
c     iunit:      file unit being read
c     eof:        end of file indicator
c     nval:       number of sources counter
c     iyear:      2-digit year
c     imonth:     2-digit month
c     iday:       2-digit day
c     ihour:      2-digit hour
c     iswitch:    switch determining if subroutine is getting
c                 # of sources or writing to arrays
      
      integer eof,iunit,nval,iyear,imonth,iday,ihour,iswitch
      
c     infil:      filename being read
c     src:        AERMOD source id being read
c     facid1:     facility ID being read from file
c     header:     header record of file (not used)
      character infil*300,facid1*20,src*12,header*100
      
c     lfound:     variable denoting facility found in file
      logical lfound
      
      lfound=.false.
      eof=0
      
c     CHANGE YEAR TO 14 FROM 11!!!!
      open(unit=iunit,file=infil,status='old')
      read(iunit,*,iostat=eof)header
      read(iunit,*,iostat=eof)facid1,src,iyear,imonth,iday,ihour  !first line
      if (ihour .eq. 0) then
          hrfix=1
      else
          hrfix=0
      endif
 5    if (eof .eq. 0) then
          if (facid1 .eq. facid)lfound=.true.
          ihour=ihour+hrfix
          if (imonth .eq. 1 .and. iday .eq. 1 .and. 
     +    ihour .eq. 1) then
              nval=nval+1
          endif
     
          if (iswitch .eq. 2) then
              temsrc(nval)=src
          endif
          read(iunit,*,iostat=eof)facid1,src,iyear,imonth,iday,ihour
          goto 5
      endif
      close(temunit)
      return
      end
c##########################################################################
c     subroutine to write out the locations and parameters to AERMOD.INP
      subroutine locs_params
      use main1
      implicit none

c     igrp:       source loop counter
      integer igrp
      
c     emis:       emissions
c                 for stacks 1000 g/s
c                 for all other sources 1000 g/s divided by area
c     nrang:      non-runway airport angle (0)
c     version 17726:  emis is now just 1.0
c     area will be included in emissions factor calculations which will be the actual emissions
      real*8 emis,nrang
      
      emis=1.0 !17234
       do igrp=1,ngrps
          if (runpt) then
c             location
              write(aerunit,5)facility(igrp)%srcid,
     +        facility(igrp)%srctyp,facility(igrp)%utmx,
     +        facility(igrp)%utmy,facility(igrp)%elev
c             source parameters
              if (facility(igrp)%srctyp .eq. 'AREA') then
c                  emis=1000.0/(facility(igrp)%xdim*facility(igrp)%ydim)
                  write(aerunit,20)facility(igrp)%srcid,emis,
     +            facility(igrp)%ht,facility(igrp)%xdim,
     +            facility(igrp)%ydim,facility(igrp)%angle,
     +            facility(igrp)%szinit
              else
c                  emis=1000.0
                  write(aerunit,25)facility(igrp)%srcid,emis,
     +            facility(igrp)%ht,facility(igrp)%stktemp,
     +            facility(igrp)%stkvel,facility(igrp)%stkdiam
              endif
          elseif (runap) then
              if (runtyp .eq. 'APRUN') then
c                  emis=1000.0*airport(igrp)%frac/airport(igrp)%area !apportion 1000 g/s to runway based on fraction
c                 location, runways
                  write(aerunit,15)airport(igrp)%srcid,
     +            airport(igrp)%srctyp,airport(igrp)%utmx1,
     +            airport(igrp)%utmy1,airport(igrp)%utmx2,
     +            airport(igrp)%utmy2,airport(igrp)%elev
                  
c                 source parameters, runways
                  write(aerunit,30)airport(igrp)%srcid,emis,
     +            airport(igrp)%ht,airport(igrp)%xdim,
     +            airport(igrp)%szinit
              else
c                  emis=1000.0/(airport(igrp)%xdim*airport(igrp)%ydim)
                  nrang=0.0
c                 locations, nonrunways
                  write(aerunit,5)airport(igrp)%srcid,
     +            airport(igrp)%srctyp,airport(igrp)%utmx1,
     +            airport(igrp)%utmy1,airport(igrp)%elev
                  
c                 parameters, nonrunways
                  write(aerunit,20),airport(igrp)%srcid,emis,
     +            airport(igrp)%ht,airport(igrp)%xdim,
     +            airport(igrp)%ydim,nrang,airport(igrp)%szinit
              endif
c          elseif (runcmv) then
          elseif (runport) then
c             locations, CMV
              write(aerunit,5)cmv(igrp)%srcid,
     +        cmv(igrp)%srctyp,cmv(igrp)%utmx1,
     +        cmv(igrp)%utmy1
c              emis=1000.0/cmv(igrp)%area !apportion 1000 g/s to shape based on fraction
              write(aerunit,35)cmv(igrp)%srcid,emis,cmv(igrp)%ht,
     +        cmv(igrp)%nvert,cmv(igrp)%szinit
c             because # of vertices can change by source, call subroutine writepoly
              call writepoly(igrp)
c          elseif (noncontr) then
cc             locations, non-CONUS tracts
c              write(aerunit,5)tract(igrp)%srcid,
c     +        tract(igrp)%srctyp,tract(igrp)%utmx1,
c     +        tract(igrp)%utmy1
cc              emis=1000.0/tract(igrp)%area !apportion 1000 g/s to shape based on fraction
c              write(aerunit,35)tract(igrp)%srcid,emis,tract(igrp)%ht,
c     +        tract(igrp)%nvert,tract(igrp)%szinit
cc             because # of vertices can change by source, call subroutine writepoly
c              call writepoly(igrp)
          else
c             location, gridded                  
              write(aerunit,10)gridded(igrp)%srcid,
     +        gridded(igrp)%srctyp,gridded(igrp)%utmx1,
     +        gridded(igrp)%utmy1
c             source parameters, gridded
c              emis=1000.0/(gridres*1000.0*gridres*1000.0)
              write(aerunit,35)gridded(igrp)%srcid,emis,
     +        gridded(igrp)%ht,gridded(igrp)%nvert,gridded(igrp)%szinit
c             vertices
              write(aerunit,40)gridded(igrp)%srcid,gridded(igrp)%utmx1,
     +        gridded(igrp)%utmy1,gridded(igrp)%utmx2,
     +        gridded(igrp)%utmy2
              write(aerunit,40)gridded(igrp)%srcid,
     +        gridded(igrp)%utmx3,gridded(igrp)%utmy3,
     +        gridded(igrp)%utmx4,gridded(igrp)%utmy4
          endif
          
      enddo
      
 5    format('SO LOCATION',1x,a12,1x,a8,2(1x,f13.5),1x,f11.4)
 10   format('SO LOCATION',1x,a12,1x,a8,2(1x,f13.5))
 15   format('SO LOCATION',1x,a12,1x,a8,4(1x,f13.5),1x,f11.4)
 20   format('SO SRCPARAM ',a12,e12.4,5(1x,f9.3))
 25   format('SO SRCPARAM ',a12,e12.4,4(1x,f9.3))
 30   format('SO SRCPARAM ',a12,e12.4,3(1x,f9.3))
 35   format('SO SRCPARAM ',a12,e12.4,1x,f9.3,1x,i4,1x,f9.3)
c 40   format('SO AREAVERT ',a12,8(1x,f13.5))
 40   format('SO AREAVERT ',a12,4(1x,f13.5))
      return
      end
c##########################################################################
c     subroutine to write out vertices for CMV or tract polygons
       subroutine writepoly(i)
       use main1
       implicit none
c     i:          source ID counter
c     nv:         number of vertices for polygon (read from file)
c     eof:        end of file indicator
c     nlines:     number of lines to write in AERMOD.INP file
c                 based on 3 pairs of coordinates per line
c     nv1:        number of vertices for polygon being processed
c     nv2:        nlines*3
c     j:          vertex counter
c     k:          vertex counter for writing out coordinates
c      integer i,nv,eof,nlines,nv1,nv2,j,k
      integer i,nv1,nv2,j,k,nlines
      character src*12
c     mod1:       modules of nv1 and 3
c     area1:      area of polygon
c     frac1:      fraction polygon is of total port
c     h:          release height
c     sz:         sigma-z
c     rlines:     nv1/3
c     ux:         array of UTM x-coordinates of vertices
c     uy:         array of UTM y-coordinates of vertices
c     avgx:       sum of ux
c     avgy:       sum of uy
c      real*8 mod1,area1,frac1,h,sz,rlines,avgx,avgy
c     area:       area of polygon
      real*8 mod1,avgx,avgy,rlines,area
c      real*8, allocatable, dimension(:) :: ux
c      real*8, allocatable, dimension(:) :: uy
      
c     facid1:     facility ID being read from file
c     src:        AERMOD source id being read
c     stype:      AERMOD source type
c     line:       header record of file (not used)
c      character facid1*20,src*12,stype*8,line*100,st*2,fips1*5

c     if (runcmv) then 
      if (runport) then
          nv1=cmv(i)%nvert
          src=cmv(i)%srcid
      else
          nv1=tract(i)%nvert
          src=tract(i)%srcid
      endif
      
c      read(parunit,*,iostat=eof)line
c      read(parunit,*,iostat=eof)st,fips1,facid1,src,stype,area1,h,nv
c 5    if (eof .eq. 0) then
c          if (src .eq. facid .and. src .eq. cmv(i)%srcid) then
c              nv1=nv
c              allocate(ux(nv1))
c              allocate(uy(nv1))
c              backspace(parunit)
c              read(parunit,*,iostat=eof)st,fips1,facid1,src,stype,area1,
c     +        h,nv,sz,(ux(j),uy(j),j=1,nv)
c          endif
c          read(parunit,*,iostat=eof)st,fips1,facid1,src,stype,area1,h,nv
c          goto 5
c      endif
c      rewind(parunit)
c      close(parunit)
c     write out to AERMOD.INP file
c     write 3 coordinate pairs per line
c     do j=1,nv1
c         write(99,*)polyx(j),polyy(j)
c     enddo
      
      mod1=mod(real(nv1),3.0)
      rlines=0.0 
      if (mod1 .eq. 0) then
          nlines=nv1/3
      else
          rlines=real(nv1)/3.0
          nlines=int(rlines)
      endif
      nv2=nlines*3
      do j=1,nv2,3
          write(aerunit,10)src,(polyx(i,k),polyy(i,k),k=j,j+2)
      enddo
      if (nv1 .ne. nv2) then
         do j=nv2+1,nv1
             write(aerunit,15)src,polyx(i,j),polyy(i,j)
         enddo
      endif 
      
c     get average coordinates for met station selection and receptor for non-CONUS
c     use AERMOD approach

c      if (akhiprvi) then
          avgx=0.0
          avgy=0.0
          area=0.0
          do j=1,nv1
              if (j .lt. nv1) then
                  area=area+(polyx(i,j)*polyy(i,j+1)-polyx(i,j+1)*
     +            polyy(i,j))
                  avgx=(polyx(i,j)+polyx(i,j+1))*(polyx(i,j)*
     +            polyy(i,j+1)-polyx(i,j+1)*polyy(i,j))+avgx
                  avgy=(polyy(i,j)+polyy(i,j+1))*(polyx(i,j)*
     +            polyy(i,j+1)-polyx(i,j+1)*polyy(i,j))+avgy
              else
                  area=area+(polyx(i,j)*polyy(i,1)-polyx(i,1)*
     +            polyy(i,j))
                  avgx=(polyx(i,j)+polyx(i,1))*(polyx(i,j)*polyy(i,1)
     +            -polyx(i,1)*polyy(i,j))+avgx
                  avgy=(polyy(i,j)+polyy(i,1))*(polyx(i,j)*polyy(i,1)-
     +            polyx(i,1)*polyy(i,j))+avgy
              endif
c              avgx=polyx(j)+avgx
c              avgy=polyy(j)+avgy
          enddo
c          cmv(i)%avgx=avgx/real(nv1)
c          cmv(i)%avgy=avgy/real(nv1)
c          if (runcmv) then
          if (runport) then
              cmv(i)%avgx=avgx/(3.0*area)
              cmv(i)%avgy=avgy/(3.0*area)
          else
              tract(i)%avgx=avgx/(3.0*area)
              tract(i)%avgy=avgy/(3.0*area)
          endif
c      endif
      
c      deallocate(ux)
c      deallocate(uy)
      
 10   format('SO AREAVERT ',a12,6(1x,f13.5))
 15   format('SO AREAVERT ',a12,2(1x,f13.5))
      return
      end
       
c##########################################################################
c     subroutine to get point source parameters
      subroutine readpoint
      use main1
      implicit none

c     eof:        end of file indicator
c     i:          source counter
c     findsrc:    function to find source index in source array
c     c:          CMAQ column
c     r:          CMAQ row
c     uz:         UTM zone
c     iyear:      2-digit year
c     imonth:     2-digit month
c     iday:       2-digit day
c     ihour:      2-digit hour
      integer eof,i,findsrc,c,r,uz,iyear,imonth,iday,ihour
      
c     line:       header record of file (not used)
c     src:        AERMOD source id being read
c     facid1:     facility ID being read from file
c     j1:         junk variable (source name)
c     qf:         AERMOD emissions factor type
c     stype:      AERMOD source type
c     locfil1:    locations file with elevations
c     st:         2-letter state abbreviation
      character line*100,src*12,facid1*20,j1*75,qf*8,stype*8,
     +st*2
      
c     ux:         source UTM x-coordinate (stack location or SW corner of area source)
c     uy:         source UTM y-coordinate (stack location or SW corner of area source)
c     e:          source elevation (m)
c     h:          stack or release height (m)
c     temp1:      stack exit temperature (K)
c     vel:        stack exit velocity (m/s)
c     diam:       stack diameter (m)
c     xd:         area source x-dimension (m)
c     yd:         area source y-dimension (m)
c     sz:         initial sigma-z for area sources (m)
c     ang:        area source angle relative to north (degrees)    
      real*8 ux,uy,e,h,temp1,vel,diam,xd,yd,ang,sz
     
c     lstop:      logical variable to stop program if source temporal
c                 mismatch
c     lexist:     logical variable denoting if elevation file exists
c     lfound:     logical variable denoting facility in file
      logical lstop,lexist,lfound

c     initialize   
      facility%utmx=-999.0
      facility%utmy=-999.0
      facility%elev=-999.0
      facility%ht=-999.0
      facility%stktemp=-999.0
      facility%stkvel=-999.0
      facility%stkdiam=-999.0
      facility%xdim=-999.0
      facility%ydim=-999.0
      facility%szinit=-999.0
      facility%angle=-999.0
      facility%annual=0.0
      facility%hrsrc=0
      facility%srctyp='NA'
      facility%qflag='NA'
      facility%srcid='NA'
          
      lstop=.false.
      
c     get source IDs from locations file
c     this will be the order of the sources
c     use file that has elevations
c      write(locfil1,'(2a)')trim(adjustl(ptcmaq)),
c     +'point_locations_elevations.csv'
      i=0
      eof=0
      inquire(file=elevfil,exist=lexist)
      if (.not. lexist) then
          write(ilogunit,*)
     +    'location/elevation file does not exist...stop'
          stop
      endif
      
c     JAT version 17278, no change to reading locations file
c     file is not helper file
      lfound=.false.
      open(unit=locunit,file=elevfil,status='old')
      read(locunit,*,iostat=eof)line
      read(locunit,*,iostat=eof)st,facid1,src,ux,uy,uz,c,r,e
 5    if (eof .eq. 0) then
          if (trim(adjustl(facid1)) .eq. trim(adjustl(facid))) then
              lfound=.true.
              i=i+1
              facility(i)%srcid=src
              facility(i)%utmx=ux
              facility(i)%utmy=uy
              facility(i)%elev=e
              col=c-coffset
              row=r-roffset
              utmzone=uz
          endif
          read(locunit,*,iostat=eof)st,facid1,src,ux,uy,uz,c,r,e
          goto 5
      endif
      close(locunit)
      
      if (.not. lfound) then
          write(ilogunit,55)trim(adjustl(facid)),trim(adjustl(elevfil)) 
      endif
c     read parameter file(s)
c     point sources
c     version 17278, no changes to format for V2
      if (lparam(1)) then
          eof=0
          lfound=.false.
          open(unit=parunit,file=paramfil(1),status='old')
          read(parunit,*,iostat=eof)line
          read(parunit,*,iostat=eof)facid1,j1,src,stype,h,
     +    temp1,vel,diam
 10       if (eof .eq. 0) then
              if (h .lt. 0) h=0.0
              if (facid1 .eq. facid) then
                  lfound=.true.
                  i=findsrc(src)
                  if (i .gt. 0) then
                      if (srcnam .eq. 'NA') srcnam=j1
                      facility(i)%srctyp=stype
                      facility(i)%ht=h
                      facility(i)%stktemp=temp1
                      facility(i)%stkvel=vel
                      facility(i)%stkdiam=diam
                  else
                      write(ilogunit,15)src,trim(adjustl(paramfil(1)))
                  endif
              endif
              read(parunit,*,iostat=eof)facid1,j1,src,stype,h,
     +        temp1,vel,diam
              goto 10
          endif
          close(parunit)
          if (.not. lfound) then
              write(ilogunit,55)trim(adjustl(facid)),
     +        trim(adjustl(paramfil(1))) 
          endif
      endif

      
c     fugitives
c     no change 17278
      if (lparam(2)) then
          eof=0
          lfound=.false.
          open(unit=parunit,file=paramfil(2),status='old')
          read(parunit,*,iostat=eof)line
          read(parunit,*,iostat=eof)facid1,j1,src,stype,h,
     +    xd,yd,ang,sz
 20       if (eof .eq. 0) then
              if (facid1 .eq. facid) then
                  i=findsrc(src)
                  lfound=.true.
                  if (i .gt. 0) then
                      if (srcnam .eq. 'NA') srcnam=j1
                      facility(i)%srctyp=stype
                      facility(i)%ht=h
                      facility(i)%xdim=xd
                      facility(i)%ydim=yd
                      facility(i)%angle=ang
                      facility(i)%szinit=sz
                  else
                      write(ilogunit,15)src,trim(adjustl(paramfil(2)))
                  endif
              endif
              read(parunit,*,iostat=eof)facid1,j1,src,stype,h,
     +        xd,yd,ang,sz
              goto 20
          endif
          close(parunit)
          if (.not. lfound) then
              write(ilogunit,55)trim(adjustl(facid)),
     +        trim(adjustl(paramfil(2))) 
          endif
      endif      
 
      
c     now check to see which sources are hourly if applicable
      if (ltemp(2)) then
          eof=0
          lfound=.false.
          open(unit=temunit,file=temfil(2),status='old')
          read(temunit,*,iostat=eof)line
          read(temunit,*,iostat=eof)facid1,src,iyear,imonth,iday,ihour
 25       if (eof .eq. 0) then
              if (facid1 .eq. facid) then  !should always be the case
                  lfound=.true.
                  if (imonth .eq. 1 .and. iday .eq. 1 .and. 
     +            ihour .eq. 1) then
                      i=findsrc(src)
                      if (i .gt. 0) then
                          facility(i)%hrsrc=1
                      else
                          write(ilogunit,15)src,trim(adjustl(temfil(2))) 
                      endif
                  endif
              endif
              read(temunit,*,iostat=eof)facid1,src,iyear,imonth,iday,
     +        ihour
              goto 25
          endif
          close(temunit)
          if (.not. lfound) then
              write(ilogunit,55)trim(adjustl(facid)),
     +        trim(adjustl(temfil(2))) 
          endif
      endif
      
 
c     now check emissions factor file to get emission factor type if applicable
c     no change 17278
      if (ltemp(1)) then
          eof=0
          lfound=.false.
          open(unit=temunit,file=temfil(1),status='old')
          read(temunit,*,iostat=eof)line
          read(temunit,*,iostat=eof)facid1,j1,src,qf
 30       if (eof .eq. 0) then
              if (facid1 .eq. facid) then
                  i=findsrc(src)
                  lfound=.true.
                  if (i .gt. 0) then
                     facility(i)%qflag=qf
                  else
                      write(ilogunit,15)src,trim(adjustl(temfil(1))) 
                  endif
              endif
              read(temunit,*,iostat=eof)facid1,j1,src,qf
              goto 30
          endif
          close(temunit)
          if (.not. lfound) then
              write(ilogunit,55)trim(adjustl(facid)),
     +        trim(adjustl(temfil(1))) 
          endif
      endif
      
c     double check that a source doesn't have both hourly and emissions factors
c     or a source has not factors
c     also check to make sure source has appropriate release parameters for source type
      do i=1,ngrps
          if (facility(i)%hrsrc .eq. 1 .and. facility(i)%qflag 
     +    .ne. 'NA') then
              write(ilogunit,35)trim(adjustl(facility(i)%srcid))
              lstop=.true.
          endif
          if (facility(i)%hrsrc .eq. 0 .and. facility(i)%qflag 
     +    .eq. 'NA') then
              write(ilogunit,40)trim(adjustl(facility(i)%srcid))
              lstop=.true.
          endif
          if ((facility(i)%srctyp .ne. 'AREA' .and. 
     +        facility(i)%xdim .ge. 0) .or. 
     +        (facility(i)%srctyp .eq. 'AREA' .and. 
     +        facility(i)%stktemp .ge. 0)) then
              write(ilogunit,45)trim(adjustl(facility(i)%srcid))
              lstop=.true.
          endif
      enddo
      
      
      if (lstop) then
          write(ilogunit,50)
          stop
      endif
                  
 15   format('WARNING! source ',a,' in file '/a,/,
     +'not found in source array')       
 35   format('Hourly source ',a,' has emissions factors')
 40   format('Non-hourly source ',a,' has no emissions factors')   
 45   format('Source ',a,' has incorrect parameters for source type')
 50   format('Stopping program')
 55    format('WARNING! facility ',a,' not found in file ',a)
      return
      end     
c##########################################################################
c     subroutine to get airport parameters
      subroutine readairport
      use main1
      implicit none
      
c     eof:        end of file indicator
c     i:          source counter
c     findsrc:    function to find source index in source array
c     c:          CMAQ column
c     r:          CMAQ row
c     uz:         UTM zone
c     eunit:      elevation file unit
      integer eof,i,findsrc,c,r,uz,eunit
      
c     line:       header record of file (not used)
c     src:        AERMOD source id being read
c     facid1:     facility ID being read from file
c     j1:         junk variable (source name)
c     qf:         AERMOD emissions factor type
c     stype:      AERMOD source type
c     st:         State abbreviation
c     efile:      elevations file
      character line*100,src*12,facid1*20,j1*75,qf*8,stype*8,st,
     +efile*250

c     ux1:        source UTM x-coordinate of runway endpoint
c     uy1:        source UTM y-coordinate of runway endpoint
c     ux2:        source UTM x-coordinate of runway endpoint
c     uy2:        source UTM y-coordinate of runway endpoint
c     lat:        area source SW corner latitude
c     lon:        area source SW corner longitude
c     e:          source elevation (m)
c     h:          release height (m)
c     xd:         runway width or area source x-dimension (m)
c     yd:         area source y-dimension (m)
c     frac1:      fraction runway is of total airport runway area 
c     area1:      area of runway or area source
c     sz:         initial sigma-z for area sources (m)
c     ang:        area source angle relative to north (degrees)  
c     gx:         CMAQ x-coordinate
c     gy:         CMAQ y-coordinate
      real*8 ux1,uy1,ux2,uy2,e,h,lat,lon,xd,yd,frac1,area1,sz,ang,gx,gy

c     lfound:     logical variable denoting facility in file
      logical lfound
      
c     initialize         
      airport%utmx1=-999.0
      airport%utmy1=-999.0
      airport%utmx2=-999.0
      airport%utmy2=-999.0
      airport%frac=-999.0
      airport%elev=-999.0
      airport%ht=-999.0
      airport%xdim=-999.0
      airport%ydim=-999.0
      airport%szinit=-999.0
      airport%area=-999.0
      airport%angle=-999.0
      airport%annual=0.0
      airport%srctyp='NA'
      airport%qflag='NA'
      airport%srcid='NA'
      eunit=32
c     read locations file
c     no change 17278
      i=0
      eof=0
      lfound=.false.
      open(unit=locunit,file=locfil,status='old')
      read(locunit,*,iostat=eof)line
      if (runtyp .eq. 'APRUN') then
          read(locunit,*,iostat=eof)st,facid1,j1,src,ux1,uy1,ux2,uy2,uz,
     +    c,r
      else
          read(locunit,*,iostat=eof)st,facid1,j1,src,gx,gy,lon,lat,ux1,
     +    uy1,uz,c,r
      endif
      
 5    if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              lfound=.true.
              i=i+1
              if (trim(adjustl(facid1)) .eq. '10605711') then  !reset to UTM 59
                  ux1=646557.02
                  uy1=5855859.95
                  uz=59
              endif
              if (trim(adjustl(facid1)) .eq. '10604811') then !reset to UTM 60
                  ux1=305028.43
                  uy1=5844171.09
                  uz=60
              endif
              if (srcnam .eq. 'NA') srcnam=j1
              airport(i)%srcid=src
              airport(i)%utmx1=ux1
              airport(i)%utmy1=uy1
              if (runtyp .eq. 'APRUN') then
                  airport(i)%utmx2=ux2
                  airport(i)%utmy2=uy2
              endif
              col=c-coffset
              row=r-roffset
              utmzone=uz
          endif
          if (runtyp .eq. 'APRUN') then
              read(locunit,*,iostat=eof)st,facid1,j1,src,ux1,uy1,ux2,
     +        uy2,uz,c,r
          else
              read(locunit,*,iostat=eof)st,facid1,j1,src,gx,gy,lon,lat,
     +        ux1,uy1,uz,c,r
          endif
          goto 5
      endif
      close(locunit)

      if (.not. lfound) then
          write(ilogunit,15)trim(adjustl(facid)),trim(adjustl(locfil)) 
      endif
      
c     read elevation file
c     no change, not helper file from SMOKE
      lfound=.false.
      write(efile,'(2a)')trim(adjustl(ptcmaq)),'airport_elevations.csv'
      open(unit=eunit,file=efile,status='old')
      eof=0
      read(eunit,*)line
      read(eunit,*,iostat=eof)st,facid1,e
 25   if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              lfound=.true.
              do i=1,ngrps
                  airport(i)%elev=e
              enddo
!              i=findsrc(src)
!              if (i .gt. 0) then
!                  airport(i)%elev=e
!              else
!                  write(*,15)src,trim(adjustl(efile)) 
!              endif
          endif
          read(eunit,*,iostat=eof)st,facid1,e
          goto 25
      endif
      close(eunit)
      
      if (.not. lfound) then
          write(ilogunit,15)trim(adjustl(facid)),trim(adjustl(efile)) 
      endif
      
c     read parameter file
c     no change 17278
      eof=0
      lfound=.false.
      open(unit=parunit,file=paramfil(1),status='old')
      read(parunit,*,iostat=eof)line
      if (runtyp .eq. 'APRUN') then
          read(parunit,*,iostat=eof)facid1,j1,src,stype,area1,frac1,h,
     +    xd,sz
      else
          read(parunit,*,iostat=eof)facid1,j1,src,h,xd,yd,ang,sz
      endif
 10   if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              lfound=.true.
              i=findsrc(src)
              if (i .gt. 0) then
                 if (runtyp .eq. 'APRUN') then
                     airport(i)%srctyp=stype
                     airport(i)%area=area1
                     airport(i)%frac=frac1
                     
                 else
                     airport(i)%ydim=yd
                     airport(i)%angle=ang
                     airport(i)%area=xd*yd
                     airport(i)%srctyp='AREA'
                     airport(i)%frac=1.0
                 endif
                 airport(i)%ht=h
                 airport(i)%xdim=xd
                 airport(i)%szinit=sz
              else
                  write(ilogunit,30)src,trim(adjustl(paramfil(1))) 
              endif
          endif
          if (runtyp .eq. 'APRUN') then
              read(parunit,*,iostat=eof)facid1,j1,src,stype,area1,frac1,
     +        h,xd,sz
          else
              read(parunit,*,iostat=eof)facid1,j1,src,h,xd,yd,ang,sz
          endif
          goto 10
      endif
                  
      close(parunit)
      
      if (.not. lfound) then
          write(ilogunit,15)trim(adjustl(facid)),
     +    trim(adjustl(paramfil(1))) 
      endif
      
c     check emissions factor file to get emissions factor type
c     no change 17278
      eof=0
      lfound=.false.
      open(unit=temunit,file=temfil(1),status='old')
      read(temunit,*,iostat=eof)line
      read(temunit,*,iostat=eof)facid1,j1,src,qf
 20   if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              lfound=.true.
              i=findsrc(src)
              if (i .gt. 0) then
                  airport%qflag=qf
              else
                  write(ilogunit,30)src,trim(adjustl(temfil(1))) 
              endif
          endif
          read(temunit,*,iostat=eof)facid1,j1,src,qf
          goto 20
      endif
      
      close(temunit)
      if (.not. lfound) then
          write(ilogunit,15)trim(adjustl(facid)),
     +    trim(adjustl(temfil(1))) 
      endif
      
15    format('WARNING! airport ',a,' not found in file ')
30    format('WARNING! source ',a,' in file '/a,/,
     +'not found in source array')        
      return
      end
c########################################################################## 
c     subroutine to read CMV (port) parameters
      subroutine readport
      use main1
      implicit none
c     eof:        end of file indicator
c     i:          source counter
c     findsrc:    function to find source index in source array
c     c:          CMAQ column
c     r:          CMAQ row
c     uz:         UTM zone     
c     nv:         number of vertices for polygon
      integer eof,i,findsrc,c,r,uz,nv,j,iyear,imon,iday,ihr,iyear4,
     +iyear2,ig,ld,i1,iv
c     line:       header record of file (not used)
c     st:         2-digit state FIPS code
c     fips1:      5-digit FIPS code
c     src:        AERMOD source id being read
c     facid1:     facility ID being read from file
c     j1:         junk variable (source name)
c     qf:         AERMOD emissions factor type
c     stype:      AERMOD source type
      character line*100,st*2,fips1*5,src*12,facid1*20,j1*75,qf*8,
     +stype*8,locfil1*300
c     ux1:        source UTM x-coordinate of runway endpoint
c     uy1:        source UTM y-coordinate of runway endpoint
c     ux2:        source UTM x-coordinate of runway endpoint
c     uy2:        source UTM y-coordinate of runway endpoint
c     e:          source elevation (m)
c     h:          release height (m)
c     frac1:      fraction polygon is of total port area 
c     area1:      area of polygon
c     sz:         initial sigma-z for area sources (m)   
      real*8 ux1,uy1,ux2,uy2,e,h,frac1,area1,sz,efact
      real*8, allocatable, dimension(:) :: x
      real*8, allocatable, dimension(:) :: y
      real*8, allocatable, dimension(:,:,:,:) :: emis
      character hrfil*250
c     lfound:     logical variable denoting facility in airport elevation file
      logical lfound
      
c     initialize         
      cmv%utmx1=-999.0
      cmv%utmy1=-999.0
c      cmv%frac=-999.0
      cmv%nvert=-9
      cmv%ht=-999.0
      cmv%szinit=-999.0
      cmv%area=-999.0
      cmv%annual=0.0
      cmv%srctyp='NA'
      cmv%qflag='NA'
      cmv%srcid='NA'
      
      
c     read locations file
      i=0
      eof=0
      lfound=.false.
      open(unit=locunit,file=locfil,status='old')
      read(locunit,*,iostat=eof)line
      read(locunit,*,iostat=eof)st,fips1,facid1,src,stype,c,r,ux1,uy1,uz
 5    if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              lfound=.true.
              i=i+1
              if (srcnam .eq. 'NA') srcnam=facid1
              cmv(i)%srcid=src
              cmv(i)%utmx1=ux1
              cmv(i)%utmy1=uy1
              col=c-coffset
              row=r-roffset
c             reset column and row for AK P20408F02188 TO 87, 177;
c             column/row appear to be bad in location file
c             coordinates okay
              if (trim(adjustl(facid1)) .eq. 'P20408F02188') then
                col=87
                row=177
              endif
              utmzone=uz
          endif
          read(locunit,*,iostat=eof)st,fips1,facid1,src,stype,c,r,ux1,
     +    uy1,uz
          goto 5
      endif
      close(locunit)
      if (.not. lfound) then
          write(ilogunit,55)trim(adjustl(facid)),trim(adjustl(locfil)) 
      endif
c     read parameter file
c     don't read in vertices here
c     version 17010, now read vertices here
      eof=0
      lfound=.false.
      i1=0
      open(unit=parunit,file=paramfil(1),status='old')
      read(parunit,*,iostat=eof)line
      read(parunit,*,iostat=eof)st,facid1,src,stype,area1,h,nv,sz
 10   if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              i=findsrc(src)
              lfound=.true.
              if (i .gt. 0) then
                  if (i1 .ne. i) then
                      cmv(i)%srctyp=stype
                      cmv(i)%area=area1
c                  cmv(i)%frac=frac1
                      cmv(i)%ht=h
                      cmv(i)%szinit=sz
                      if(allocated(x)) deallocate(x)
                      if (allocated(y)) deallocate(y)
                      allocate(x(nv))
                      allocate(y(nv))
                      x=0.0
                      y=0.0
                      nv=nv-1
                      cmv(i)%nvert=nv
                      j=0
                  endif
                  backspace(parunit)
                  j=j+1
                  read(parunit,*,iostat=eof)st,facid1,src,stype,
     +            area1,h,nv,sz,x(j),y(j)
c    +            area1,h,nv,sz,(polyx(j),polyy(j),j=1,nv)
c                  if (x(1) .eq. x(nv) .and. y(1) .eq. y(nv)) nv=nv-1
                  if (.not. allocated(polyx)) then
                     allocate(polyx(ngrps,nv))
                     allocate(polyy(ngrps,nv))
                  endif
c                  do j=1,nv
                   if (j .le. cmv(i)%nvert) then
                      polyx(i,j)=x(j)
                      polyy(i,j)=y(j)
                   endif
c                  enddo

c                  deallocate(x)
c                  deallocate(y)
              else
                  write(ilogunit,15)src,trim(adjustl(paramfil(1))) 
              endif
              i1=i
          endif
          read(parunit,*,iostat=eof)st,facid1,src,stype,area1,h,
     +    nv,sz
          goto 10
      endif
      if (.not. lfound) then
          write(ilogunit,55)trim(adjustl(facid)),
     +    trim(adjustl(paramfil(1))) 
      endif
c     don't close file, it will be used in subroutine writepoly      
c     close version 17010
c      rewind(parunit) 
      close(parunit)
c     check emissions factor file to get emissions factor type
      eof=0
      lfound=.false.
!     for 2017, hourly emissions by state
      if (state2 .eq. 'AK' .or. state2 .eq. '02') then
          write(temfil(1),'(2(a))')trim(adjustl(pttemp)),
     +    'CMVP_9AK_02_hourly.csv'
      elseif (state2 .eq. 'HI' .or. state2 .eq. '15') then
          write(temfil(1),'(2(a))')trim(adjustl(pttemp)),
     +    'CMVP_3HI_15_hourly.csv'
      elseif (state2 .eq. 'PR' .or. state2 .eq. '72') then
          write(temfil(1),'(2(a))')trim(adjustl(pttemp)),
     +    'CMVP_3PR_72_hourly.csv'
      elseif (state2 .eq. 'VI' .or. state2 .eq. '78')then
          write(temfil(1),'(2(a))')trim(adjustl(pttemp)),
     +    'CMVP_3PR_78_hourly.csv'
      else
          write(temfil(1),'(4(a))')trim(adjustl(pttemp)),
     +    'CMVP_12_',trim(adjustl(state2)),'_hourly.csv'
      endif
      
      allocate(emis(ngrps,12,31,24))
      
      open(unit=temunit,file=temfil(1),status='old')
      read(temunit,*,iostat=eof)j1,facid1,src,iyear,imon,iday,ihr,efact
 20   if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              i=findsrc(src)
              lfound=.true.
              if (i .gt. 0) then
                  if (iyear .lt. 100) then
                      iyear4=iyear+2000
                      iyear2=iyear
                  else
                      iyear4=iyear
                      iyear2=iyear-2000
                  endif
                  if (iyear4 .eq. modyear) then
                      emis(i,imon,iday,ihr)=efact*annemis*251.9957778/
     +                cmv(i)%area
                      if (emis(i,imon,iday,ihr) .lt. 0.0)
     +                emis(i,imon,iday,ihr)=0.0
                      cmv(i)%annual=cmv(i)%annual+(efact*annemis)
                  endif
              else
                  write(ilogunit,15)src,trim(adjustl(temfil(1))) 
              endif
          endif
          read(temunit,*,iostat=eof)j1,facid1,src,iyear,imon,iday,ihr,
     +   efact
          goto 20
      endif
      
      close(temunit)
      
      write(hrfil,'(2a)')trim(adjustl(facid)),'_hourly.dat'
      open(unit=hrunit,file=hrfil,status='replace')
c     write emissions factors to AERMOD input file
      do imon=1,12
          if (lleap .and. imon .eq. 2) then
              ld=1
          else
              ld=0
          endif
          do iday=1,days(imon)+ld
              do ihr=1,24
                  do ig=1,ngrps
                     src=cmv(ig)%srcid
c                      src=gridded(ig)%srcid
                      write(hrunit,25)iyear2,imon,iday,ihr,src,
     +                emis(ig,imon,iday,ihr)
                  enddo
              enddo
          enddo
      enddo
     
c      deallocate(emisfact2)
      deallocate(emis)
      close(hrunit)
      if (.not. lfound) then
          write(ilogunit,55)trim(adjustl(facid)),
     +    trim(adjustl(temfil(1))) 
      endif
30    format('SO HOUREMIS ',a,1x,a12)
15    format('WARNING! source ',a,' in file '/a,/,
     +'not found in source array')  

 25   format('SO HOUREMIS',4(1x,i2),1x,a12,1x,e13.4)
55    format('WARNING! CMV shape ',a,' not found in file ')      
      return
      end
c##########################################################################
c     subroutine to read in tract emissions
      subroutine readtract
      use main1
      implicit none
      
      integer eof,i,findsrc,uz,nv,j,nva
      character j1*75,header*100,st*2,fips6*6,poly*4,tract1*11,src*12,
     +facid2*20,qf*6,line*100,stype*8
      real*8 ux1,uy1,ux2,uy2,e,h,area1,sz,tarea
      real*8, allocatable, dimension(:) :: x
      real*8, allocatable, dimension(:) :: y
c     lfound:     logical variable denoting facility in airport elevation file
      logical lfound
      
c     initialize         
      tract%utmx1=-999.0
      tract%utmy1=-999.0
c      cmv%frac=-999.0
      tract%nvert=-9
      tract%ht=-999.0
      tract%szinit=-999.0
      tract%area=-999.0
      tract%annual=0.0
c      tract%seasemis=0.0
      tract%srctyp='NA'
      tract%qflag='NA'
      tract%srcid='NA'
      tract%fips='NA'
      tract%tract='NA'
      
c     read locations file
      i=0
      eof=0
      lfound=.false.
      open(unit=locunit,file=locfil,status='old')
      read(locunit,*,iostat=eof)line
      read(locunit,*,iostat=eof)j1,st,fips6,tract1,poly,facid2,src,
     +area1,tarea,ux1,uy1,uz
 5    if (eof .eq. 0) then
          if (src .eq. facid) then
              lfound=.true.
              i=i+1
              if (srcnam .eq. 'NA') srcnam=src
              tract(i)%srcid=src
              tract(i)%utmx1=ux1
              tract(i)%utmy1=uy1
              tract(i)%area=area1
              tract(i)%fips=fips6(2:6)
              tract(i)%tract=tract1
              utmzone=uz
          endif
          read(locunit,*,iostat=eof)j1,st,fips6,tract1,poly,facid2,src,
     +    area1,tarea,ux1,uy1,uz
          goto 5
      endif
      close(locunit)
      if (.not. lfound) then
          write(ilogunit,55)trim(adjustl(facid)),trim(adjustl(locfil)) 
      endif
c     read parameter file
c     don't read in vertices here
c     version 17010, now read vertices here
      eof=0
      lfound=.false.
      open(unit=parunit,file=paramfil(1),status='old')
      read(parunit,*,iostat=eof)line
      read(parunit,*,iostat=eof)j1,st,fips6,tract1,poly,facid2,src,
     +stype,h,nv,sz
 10   if (eof .eq. 0) then
          if (src .eq. facid) then
              i=findsrc(src)
              lfound=.true.
              if (i .gt. 0) then
                  tract(i)%srctyp=stype
                  tract(i)%ht=h
                  tract(i)%szinit=sz
                  allocate(x(nv))
                  allocate(y(nv))
                  backspace(parunit)
                  read(parunit,*,iostat=eof)j1,st,fips6,tract1,poly,
     +            facid2,src,stype,h,nv,sz,(x(j),y(j),j=1,nv) 
c    +            area1,h,nv,sz,(polyx(j),polyy(j),j=1,nv)
c                 do j=1,nv
c                     write(98,*)x(j),y(j)
c             enddo
              
                  if (x(1) .eq. x(nv) .and. y(1) .eq. y(nv)) nv=nv-1
                  tract(i)%nvert=nv
                  if (i .eq. 1) then
                      allocate(polyx(ngrps,nv))
                      allocate(polyy(ngrps,nv))
                  endif
                  do j=1,nv
                      polyx(i,j)=x(j)
                      polyy(i,j)=y(j)
                  enddo
                  deallocate(x)
                  deallocate(y)
              else
                  write(ilogunit,15)src,trim(adjustl(paramfil(1))) 
              endif
          endif
          read(parunit,*,iostat=eof)j1,st,fips6,tract1,poly,facid2,src,
     +    stype,h,nv,sz
          goto 10
      endif
      if (.not. lfound) then
          write(ilogunit,55)trim(adjustl(facid)),
     +    trim(adjustl(paramfil(1))) 
      endif
c     don't close file, it will be used in subroutine writepoly      
c     close version 17010
c      rewind(parunit) 
      close(parunit)
c     check emissions factor file to get emissions factor type
      if (.not. runrwc) then
          eof=0
          lfound=.false.
          open(unit=temunit,file=temfil(1),status='old')
          read(temunit,*,iostat=eof)line
          read(temunit,*,iostat=eof)j1,st,fips6,tract1,poly,facid2,src,
     +    qf
 20       if (eof .eq. 0) then
              if (src .eq. facid) then
                  i=findsrc(src)
                  lfound=.true.
                  if (i .gt. 0) then
                      tract%qflag=qf
                  else
                      write(ilogunit,15)src,trim(adjustl(temfil(1))) 
                  endif
              endif
              read(temunit,*,iostat=eof)j1,st,fips6,tract1,poly,facid2,
     +        src,qf
              goto 20
          endif
      
          close(temunit)
          if (.not. lfound) then
              write(ilogunit,55)trim(adjustl(facid)),
     +        trim(adjustl(temfil(1))) 
          endif
      endif
      
15    format('WARNING! source ',a,' in file '/a,/,
     +'not found in source array')  
55    format('WARNING! tract shape ',a,' not found in file ')      
      return
      end
c##########################################################################
c     subroutine to read in gridded emissions
      subroutine readgrid
      use main1
      implicit none

c     eof:        end of file indicator
c     i:          source counter
c     findsrc:    function to find source index in source array
c     c:          CMAQ column
c     r:          CMAQ row
c     uz:         UTM zone      
      integer eof,i,findsrc,uz,nv,m,iyear,iyear4,iyear2,imon,iday,ihr,
     +ig,ld,istate1
      
c     line:       header record of file (not used)
c     src:        AERMOD source id being read
c     facid1:     facility ID being read from file
c     j1:         junk variable (source name)
c     qf:         AERMOD emissions factor type
c     stype:      AERMOD source type (new for V2, 17278)
      character line*100,src*12,facid1*20,j1*75,qf*8,stype*8,st1*2,
     +hrfil*250 
    
c     ux1:        UTM x-coordinate of vertex
c     uy1:        UTM y-coordinate of vertex
c     ux2:        UTM x-coordinate of vertex
c     uy2:        UTM y-coordinate of vertex
c     ux3:        UTM x-coordinate of vertex
c     uy3:        UTM y-coordinate of vertex
c     ux4:        UTM x-coordinate of vertex
c     uy4:        UTM y-coordinate of vertex
c     h:          release height (m)
c     sz:         initial sigma-z for area sources (m)
c     lo1-lo4,la1-lo4: latlons (not used)
      real*8 ux1,uy1,ux2,uy2,ux3,uy3,ux4,uy4,h,sz,
     +lo1,la1,lo2,la2,lo3,la3,lo4,la4,efact
      real*8, allocatable, dimension(:,:,:,:) :: emis
c     lfound:     logical variable denoting facility in file
      logical lfound
      
c     initialize      
      gridded%utmx1=-999.0
      gridded%utmy1=-999.0
      gridded%utmx2=-999.0
      gridded%utmy2=-999.0
      gridded%utmx3=-999.0
      gridded%utmy3=-999.0
      gridded%utmx4=-999.0
      gridded%utmy4=-999.0
      gridded%nvert=-9
      gridded%ht=-999.0
      gridded%szinit=-999.0
      gridded%annual=0.0
c     gridded%monemis=0.0
      gridded%srctyp='NA'
      gridded%qflag='NA'
      gridded%srcid='NA'
      gridded%state='NA'

c     offset column and row
      col=col-coffset
      row=row-roffset
      
c     read locations file
c     modified 17278 to include stype in location file
      i=0
      eof=0
      lfound=.false.
      open(unit=locunit,file=locfil,status='old')
      read(locunit,*,iostat=eof)line
      if (runon) then
          read(locunit,*,iostat=eof)j1,facid1,src,ux1,uy1,uz
      elseif (runcmv) then
          read(locunit,*,iostat=eof)j1,st1,facid1,src,ux1,uy1,uz
      else
          read(locunit,*,iostat=eof)j1,facid1,src,stype,ux1,uy1,uz
      endif
 5    if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              lfound=.true.
              i=i+1
              if (srcnam .eq. 'NA') srcnam=facid1
              gridded(i)%srcid=src
c             gridded(i)%utmx1=ux1
c             gridded(i)%utmy1=uy1
              if (runon .or. runcmv) stype='AREAPOLY'
              gridded(i)%srctyp=stype
              if (runcmv) then
                 read(st1,'(i2)')istate1
                 write(st1,'(i2.2)')istate1
                 gridded(i)%state=st1
              endif
              utmzone=uz
              do m=1,12
                gridded(i)%monemis(m)=0.0
              enddo
          endif
          if (runon) then
              read(locunit,*,iostat=eof)j1,facid1,src,ux1,uy1,uz
          elseif (runcmv) then
              read(locunit,*,iostat=eof)j1,st1,facid1,src,ux1,uy1,uz
          else
              read(locunit,*,iostat=eof)j1,facid1,src,stype,ux1,uy1,uz
          endif
          goto 5
      endif
      close(locunit)
      if (.not. lfound) then
          write(ilogunit,55)trim(adjustl(facid)),trim(adjustl(locfil)) 
      endif
c     read parameter file
c     modified 17278
      eof=0
      lfound=.false.
      open(unit=parunit,file=paramfil(1),status='old')
      read(parunit,*,iostat=eof)line
      read(parunit,*,iostat=eof)j1,facid1,src,h,nv,sz,ux1,uy1,ux2,uy2,
     +ux3,uy3,ux4,uy4,lo1,la1,lo2,la2,lo3,la3,lo4,la4
 10   if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              i=findsrc(src)
              lfound=.true.
              if (i .gt. 0) then
c                  gridded(i)%srctyp='AREAPOLY'
                  gridded(i)%ht=h
                  gridded(i)%szinit=sz
                  gridded(i)%nvert=nv
                  gridded(i)%utmx1=ux1
                  gridded(i)%utmy1=uy1
                  gridded(i)%utmx2=ux2
                  gridded(i)%utmy2=uy2
                  gridded(i)%utmx3=ux3
                  gridded(i)%utmy3=uy3
                  gridded(i)%utmx4=ux4
                  gridded(i)%utmy4=uy4
              else
                  write(ilogunit,15)src,trim(adjustl(paramfil(1))) 
              endif
          endif
          read(parunit,*,iostat=eof)j1,facid1,src,h,nv,sz,ux1,uy1,ux2,
     +    uy2,ux3,uy3,ux4,uy4,lo1,la1,lo2,la2,lo3,la3,lo4,la4
          goto 10
      endif
      
      close(parunit)
      if (.not. lfound) then
          write(ilogunit,55)trim(adjustl(facid)),
     +    trim(adjustl(paramfil(1))) 
      endif
c     check emissions factor file to get emissions factor type
c     for nonroad, and nonpoint
c     no change 17278
      if (runnon .or. runnphi .or. runnplow .or. runog) then
          eof=0
          lfound=.false.
          open(unit=temunit,file=temfil(1),status='old')
          read(temunit,*,iostat=eof)line
          read(temunit,*,iostat=eof)j1,facid1,src,qf
  20      if (eof .eq. 0) then
              if (facid1 .eq. facid) then
                  i=findsrc(src)
                  lfound=.true.
                  if (i .gt. 0) then
                      gridded%qflag=qf
                  else
                      write(ilogunit,15)src,trim(adjustl(temfil(1))) 
                  endif
              endif
              read(temunit,*,iostat=eof)j1,facid1,src,qf
              goto 20
          endif
      
          close(temunit)
          if (.not. lfound) then
              write(ilogunit,55)trim(adjustl(facid)),
     +        trim(adjustl(temfil(1))) 
          endif
      endif
      
!     for 2017, hourly emissions by state for CMV underway
      if (runcmv) then
          if (gridded(1)%state .eq. '02') then
              write(temfil(1),'(2(a))')trim(adjustl(pttemp)),
     +        'CMVU_9AK_02_hourly.csv'
          elseif (gridded(1)%state .eq. '15') then
              write(temfil(1),'(2(a))')trim(adjustl(pttemp)),
     +        'CMVU_3HI_15_hourly.csv'
          elseif (gridded(1)%state .eq. '72') then
              write(temfil(1),'(2(a))')trim(adjustl(pttemp)),
     +        'CMVU_3PR_72_hourly.csv'
          elseif (gridded(1)%state .eq. '78') then
              write(temfil(1),'(2(a))')trim(adjustl(pttemp)),
     +        'CMVU_3PR_78_hourly.csv'
          else
              write(temfil(1),'(4(a))')trim(adjustl(pttemp)),
     +        'CMVU_12_',gridded(1)%state,'_hourly.csv'
          endif
          allocate(emis(ngrps,12,31,24))
          
          open(unit=temunit,file=temfil(1),status='old')
          read(temunit,*,iostat=eof)j1,facid1,src,iyear,imon,iday,ihr,
     +    efact
 40       if (eof .eq. 0) then
              if (facid1 .eq. facid) then
                  i=findsrc(src)
                  lfound=.true.
                  if (i .gt. 0) then
                      if (iyear .lt. 100) then
                          iyear4=iyear+2000
                          iyear2=iyear
                      else
                          iyear4=iyear
                          iyear2=iyear-2000
                      endif
                      if (iyear4 .eq. modyear) then
                          emis(i,imon,iday,ihr)=efact*annemis*
     +                    251.9957778/(gridres*1000*gridres*1000)
                          if (emis(i,imon,iday,ihr) .lt. 0.0) 
     +                    emis(i,imon,iday,ihr)=0.0
                          gridded(i)%annual=gridded(i)%annual+
     +                    (efact*annemis)
                      endif
                  else
                      write(ilogunit,15)src,trim(adjustl(temfil(1))) 
                  endif
              endif
              read(temunit,*,iostat=eof)j1,facid1,src,iyear,imon,iday,
     +        ihr,efact
              goto 40
          endif
      
          close(temunit)
      
          write(hrfil,'(2a)')trim(adjustl(facid)),'_hourly.dat'
          open(unit=hrunit,file=hrfil,status='replace')
c         do ig=1,ngrps
c             write(aerunit,30)trim(adjustl(hrfil)),gridded(ig)%srcid
c         enddo
c         write emissions factors to AERMOD input file
          do imon=1,12
              if (lleap .and. imon .eq. 2) then
                  ld=1
              else
                  ld=0
              endif
              do iday=1,days(imon)+ld
                  do ihr=1,24
                      do ig=1,ngrps
                         src=gridded(ig)%srcid
                          write(hrunit,25)iyear2,imon,iday,ihr,src,
     +                    emis(ig,imon,iday,ihr)
                      enddo
              enddo
          enddo
     
c         deallocate(emisfact2)
          deallocate(emis)
          close(hrunit)
      endif    
      
15    format('WARNING! source ',a,' in file '/a,/,
     +'not found in source array')  
30    format('SO HOUREMIS ',a,1x,a12)
25    format('SO HOUREMIS',4(1x,i2),1x,a12,1x,e13.4)
55    format('WARNING! grid cell ',a,' not found in file ')     
      return
      end
c##########################################################################      
c     subroutine to get urban/rural classification of facility
      subroutine urbrur(lurban)
      use main1
      implicit none
      
c     facid1:     facility ID in urban source file
c     infil:      urban source filename
      character facid1*20,infil*300
      
c     eof:        end of file indicator
c     urbpop:     urban population for source
      integer eof,urbpop

c     lfound:     variable denoting urban file exists
c     lurban:     variable denoting source is urban (true) or rural (false)
      logical lfound,lurban
      
      lurban=.false.
      urbpop=0
      
c    version 18011, for gridded non-CONUS, use state specific file 
      if (noncontr) then
         write(infil,'(4(a))')trim(adjustl(ptcmaq)),'urban_sources_',
     +   trim(adjustl(state2)),'_grid.dat'
      else
         write(infil,'(2(a))')trim(adjustl(ptcmaq)),'urban_sources.dat'
      endif
      
      inquire(file=infil,exist=lfound)
      if (.not. lfound) then 
          write(ilogunit,5)trim(adjustl(infil)),facid
          stop
      endif
      open(unit=12,file=infil,status='old')
      eof=0
      read(12,*,iostat=eof)facid1,urbpop
      

15    if (eof .eq. 0) then
          if (facid1 .eq. facid) then
              lurban=.true.
               write(aerunit,10)real(urbpop)
          else
              read(12,*,iostat=eof)facid1,urbpop
              goto 15
          endif
      endif
      close(12)

 5    format(a,' not found for ',a,/'abort')
 10   format('CO URBANOPT ',f12.1)
      return
      end
c##########################################################################
      function findsrc(src)
      use main1
      implicit none

c     findsrc:    source index in source array
      integer findsrc

c     src:        AERMOD source id
      character src*12
      
c     lfound:     logical variable denoting source found in array
      logical lfound
      
      
      findsrc=1
      lfound=.false.
      if (runpt) then
          do while(findsrc .le. ngrps .and. .not. lfound)
              if (trim(adjustl(src)) .eq. 
     +        trim(adjustl(facility(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(airport(findsrc)%srcid)))then
                  lfound=.true.
              else
                  findsrc=findsrc+1
              endif
          enddo
c      elseif (runcmv) then
      elseif (runport) then
          do while(findsrc .le. ngrps .and. .not. lfound)
              if (trim(adjustl(src)) .eq. 
     +        trim(adjustl(cmv(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 (src .eq. tract(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)
              if (trim(adjustl(src)) .eq. 
     +        trim(adjustl(gridded(findsrc)%srcid)))then
                  lfound=.true.
              else
                  findsrc=findsrc+1
              endif
          enddo
      endif
      
      if (.not. lfound)findsrc=0
      return
      end
c##########################################################################      
c     subroutine to get apppropriate receptors
      subroutine receptors
      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)
c     nzero:      number of zero pop blocks written to AERMOD.INP for non-CONUS sources
      integer eof,col1,col2,row1,row2,nblock,ngrid,nmon,c,r,cr,colrow,
     +modcol,modrow,recflag,irec,nzero
      
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
      real*8 x,y,cmaqx,cmaqy,relev,hill
      
c     recfil:     receptor filename
c     recid:      census block ID
c     monid:      monitor ID       
      character recfil*250,recid*15,monid*13
      
c     lexist:     logical variable denoting receptor file exists
c     lwrite:     logical varible denoting to write receptor to AERMOD.INP
      logical lexist,lwrite
      
c     initialize variables
      nblock=0
      ngrid=0
      nmon=0
      nzero=0  !added 17234
c      if (.not. runcmv) then
      if (.not. runport) then
c         modcol=5
c         modrow=5
          if (state2 .eq. '02') then !9km
              modcol=7 !63 km
              modrow=7
          elseif (state2 .eq. '15' .or. state2 .eq. '72' .or.  
     +    state2 .eq. '78') then !3km
              modcol=20 !60 km
              modrow=20
          else
              modcol=8
              modrow=8
          endif
      else !account for large shapes covering several grid cells.
          modcol=40
          modrow=40
      endif
      
      col1=col-modcol
      col2=col+modcol
      row1=row-modrow
      row2=row+modrow

c     write blocks within specified distances for point, airports, and CMV
c     flag receptors within 30 m of a facility release point, 30 m within runway (entire length),
c     or side of a CMV polygon)
c     for non-conus gridded, now have gridded receptors for CMAQ except in UTM zones 1, 59,60
c     still model blocks in those zone within 50 km
c      if (runpt .or. runap .or. runcmv .or. (noncontr .and. 
      if (runpt .or. runap .or. runport .or. (noncontr .and. 
     +(utmzone .eq. 1 .or. utmzone .ge. 59))) then
          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(ilogunit,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. akhiprvi) then  !if within the 5x5 block of cells
c                 get distance and determine whether to write out
                  if (akhiprvi .and. colrow .eq. 0) then
c                     if colrow = 0, then outside grid check to see within 50 km
                      call recdist(x,y,grddist,lwrite,recflag)
                  else
                      call recdist(x,y,blkdist,lwrite,recflag)
                  endif
                  if (lwrite) then
                      write(aerunit,20)recid,colrow,cmaqx,cmaqy,recflag,
     +                x,y,relev,hill
                      nblock=nblock+1
                  endif
              endif
              read(recunit,*,iostat=eof)recid,colrow,x,y,cmaqx,cmaqy,
     +        relev,hill
              goto 5
          endif
      
          close(recunit)
      endif
      
c     added 17234, non-populated blocks for non-CONUS areas
c     concentrations will be used for tract averages
c     flag receptors within 30 m of a facility release point, 30 m within runway (entire length),
c     or side of a CMV polygon)
c     if (((runpt .or. runap .or. runcmv) .and. akhiprvi) .or. 
c    +noncontr) then
c     only model non-populated blocks in AK, HI in UTM zones 1-4 or 59-60 and outside CMAQ domains
      if (akhiprvi .and. (utmzone .eq. 1 .or. utmzone .ge. 59)) then
          if (utmzone .ge. 10) then
              write(recfil,'(2a,i2,a)')trim(adjustl(recdir)),
     +        'zero_blocks_utm_',utmzone,'.dat'
          else
              write(recfil,'(2a,i1,a)')trim(adjustl(recdir)),
     +        'zero_blocks_utm_',utmzone,'.dat'
          endif
          lexist=.false.
          inquire(file=recfil,exist=lexist)
          if (.not. lexist) then
              write(ilogunit,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
 55       if (eof .eq. 0) then
              c=int(real(colrow)/1000.0) !get column
              r=colrow-c*1000 !get row
c     start here JAT
c              if ((col1 .le. c .and. c .le. col2 .and. row1 .le. r .and. 
c    +        r .le. row2) .or. (akhiprvi .and. colrow .eq. 0)) then  !if within the 5x5 block of cells
              if (colrow .eq. 0 .and. akhiprvi) then
c                 get distance and determine whether to write out
                  if (akhiprvi) then
c                     if tract source check to see if receptor is same tract, if so, keep
c                     version 17208, now check for non-CONUS states
c                      if (noncontr .and. state2 .eq. 'AK') then
c                     if (noncontr) then
c                          call intract(recid,lwrite,recflag)
c                          if (.not. lwrite) !outside of tract
                          call recdist(x,y,grddist,lwrite,recflag)
c                     else
c                         call recdist(x,y,grddist,lwrite,recflag)
c                     endif
                  else
                      call recdist(x,y,blkdist,lwrite,recflag)
                  endif
                  if (lwrite) then
                      write(aerunit,60)recid,colrow,cmaqx,cmaqy,recflag,
     +                x,y,relev,hill
                      nzero=nzero+1
                  endif
              endif
              read(recunit,*,iostat=eof)recid,colrow,x,y,cmaqx,cmaqy,
     +        relev,hill
              goto 55
          endif
      
          close(recunit)
      endif

c     gridded receptors
c     version 19134, now will modeled gridded receptors for non-CONUS
      if (.not. akhiprvi) then
          write(recfil,'(2a,i2,a)')trim(adjustl(recdir)),'grid_utm_',
     +    utmzone,'.dat'
      else
          if (state2 .eq. '02' .or. state2 .eq. '15') then
              if (utmzone .lt. 10) write(recfil,'(3a,i1,a)')
     +        trim(adjustl(recdir)),trim(adjustl(state2)),'_grid_utm_',
     +        utmzone,'.dat'
          else
              write(recfil,'(2a,i2,a)')trim(adjustl(recdir)),
     +        '72_78_grid_utm_',utmzone,'.dat'
          endif
      endif
      lexist=.false.
      inquire(file=recfil,exist=lexist)
      if (.not. lexist) then
          write(ilogunit,35)trim(adjustl(recfil))
          stop
      endif
c     write gridded receptors within specified distances
      open(unit=recunit,file=recfil,status='old')
      read(recunit,*,iostat=eof)colrow,x,y,cmaqx,cmaqy,relev,hill,irec
 10   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) then !if within the 5x5 block of cells
              call recdist(x,y,grddist+10000.,lwrite,recflag)
              if (lwrite) then
                  write(aerunit,30)colrow,irec,cmaqx,cmaqy,recflag,
     +            x,y,relev,hill
                  ngrid=ngrid+1
              endif
c              if (runpt .or. runap .or. runcmv) then
c                 get distance and determine whether to write out
!                  call recdist(x,y,grddist,lwrite,recflag)
c                  if (lwrite) then
c                      write(aerunit,30)colrow,cmaqx,cmaqy,x,y,relev,hill
c                      ngrid=ngrid+1
c                  endif
c              else
c                 get distance and determine whether to write out
!                  call recs4grid(x,y,lwrite)
c                  if (lwrite) then
c                      write(aerunit,30)colrow,cmaqx,cmaqy,x,y,relev,hill
c                      ngrid=ngrid+1
c                  endif
c              endif
                  
          endif
          read(recunit,*,iostat=eof)colrow,x,y,cmaqx,cmaqy,relev,hill,
     +    irec
          goto 10
      endif
      
      close(recunit)
c      endif
c     monitors
c      if (runpt .or. runap .or. runcmv .or. (noncontr .and.
      if (runpt .or. runap .or. runport .or. (noncontr .and.
     +(utmzone .lt. 4 .or. utmzone .ge. 59))) then
          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(ilogunit,36)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. akhiprvi) then  !if within the 5x5 block of cells
c                 get distance and determine whether to write out
                      if (akhiprvi .and. colrow .eq. 0) then
                          call recdist(x,y,grddist,lwrite,recflag)
                      else
                          call recdist(x,y,blkdist,lwrite,recflag)
                      endif
                      if (lwrite) then
                          write(aerunit,25)monid,colrow,cmaqx,cmaqy,
     +                    recflag,x,y,relev,hill
                          nmon=nmon+1
                      endif
                  endif
                  read(recunit,*,iostat=eof)monid,colrow,x,y,cmaqx,
     +            cmaqy,relev,hill
                  goto 45
              endif
              close(recunit)
          endif
      endif
      
c     calculate total receptors
      nrec=nblock+ngrid+nmon+nzero
c      if (runcmv) close(parunit)    
      write(aerunit,40)nblock,ngrid,nmon,nzero
      if (nblock .eq. 0 .and. ngrid .eq. 0 .and. nmon .eq. 0 .and. 
     +nzero .eq. 0) write(ilogunit,50)
20    format('**BLOCK= ',a,1x,'COLROW= ',i6,1x,'CMAQX= ',f11.0,1x,
     +'CMAQY= ',f11.0,' FLAG = ',i1/'RE DISCCART ',f12.3,1x,f13.3,
     +2(1x,f9.3)) 
25    format('**MONITOR= ',a,1x,'COLROW= ',i6,1x,'CMAQX= ',f11.0,1x,
     +'CMAQY= ',f11.0,' FLAG = ',i1/'RE DISCCART ',f12.3,1x,f13.3,
     +2(1x,f9.3)) 
30    format('**COLROW= ',i6,1x,'REC= ',i3,' CMAQX= ',f11.0,1x,
     +'CMAQY= ',f11.0,' FLAG = ',i1/'RE DISCCART ',f12.3,1x,f13.3,
     +2(1x,f9.3))      
            
60    format('**ZEROBLOCK= ',a,1x,'COLROW= ',i6,1x,'CMAQX= ',f11.0,1x,
     +'CMAQY= ',f11.0,' FLAG = ',i1/'RE DISCCART ',f12.3,1x,f13.3,
     +2(1x,f9.3)) 
 35   format('receptor file ',a,' does not exist, stopping program')
 36   format('receptor file ',a,' does not exist')
 40   format('** NUMBER OF BLOCKS ',i6
     +/'** NUMBER OF GRIDDED RECEPTORS ',i6
     +/'** NUMBER OF MONITORS ',i6
     +/'** NUMBER OF ZERO BLOCKS ',i6)
 50   format('WARNING NO RECEPTORS FOR SOURCE, AERMOD WILL CRASH!!!')
      return
      end
c##########################################################################
c     subroutine to check to see if a receptor is inside a tract source
      subroutine intract(rid,lwrite,rflag)
      use main1
      implicit none
      integer i,rflag
      character rid*15,tract1*11
      logical lfound,lwrite
      
      
      lwrite=.false.
      tract1=trim(adjustl(rid(1:11)))
      if (trim(adjustl(tract1)) .eq. 
     +trim(adjustl(tract(1)%tract))) lwrite=.true.
      
      if (lwrite) then
         rflag=0
      else
         rflag=1
      endif
      return
      end
          
                  
      
c##########################################################################
c     subroutine to determine distance from receptor to source 
c     and determine if receptor within 30 m of source
      subroutine recdist(x,y,maxdist,lwrite,iflag)
      use main1
      implicit none

c     iflag:      flag denoting if source within 30 m of source
c     igrp:       source loop counter
      integer iflag,igrp
      
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,maxdist,apdist,ptdist,distance

c     lwrite:     logical variable denoting to write receptor to AERMOD.INP
      logical lwrite
      
      iflag=0
      lwrite=.false.

c     18036 added code to reset iflag each iteration of do loop
      if (runap) then
          do igrp=1,ngrps
              dist=apdist(x,y,igrp)  !get distance from receptor to points along runways
              if (dist .le. maxdist) then
                  lwrite=.true.
                  if (dist .le. 30.0) iflag=1
                  if (dist .gt. grddist) then 
                     iflag=3
                  else
                     iflag=0
                  endif
              endif
          enddo
      elseif (runpt) then
          do igrp=1,ngrps
              if (facility(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 .le. maxdist) then 
                  lwrite=.true.
                  if (dist .le. 30.0) iflag=1
                  if (dist .gt. grddist) then
                    iflag=3
                  else
                    iflag=0
                  endif
              endif
          enddo 
c      elseif (runcmv) then ! .or. noncontr) then 
      elseif (runport) then 
          do igrp=1,ngrps
              call cmvcheck(x,y,igrp,dist)  !distance from receptor to CMV polygons
              if (dist .le. maxdist) then 
                  lwrite=.true.
                  if (dist .le. 30.0) iflag=1
                  if (dist .gt. grddist) then
                     iflag=3
                  else
                     iflag=0 
                  endif
              endif

c             check distance to center of CMV shape or tract. maybe greater than maxdist from edge
c             but within distance of center, i.e. inside the shape
c              else  
c                  if (runcmv) then
c                      dist=distance(x,y,cmv(igrp)%avgx,cmv(igrp)%avgy)
c                  else
c                      dist=distance(x,y,tract(igrp)%avgx,
c     +                tract(igrp)%avgy)
c                  endif
c                  if (dist .le. maxdist) then
c                      lwrite=.true.
c                      if (dist .le. 30.0) iflag=1
c                      if (dist .gt. grddist) iflag=3
c                  endif
c              endif
          enddo
      elseif (noncontr) then
          call recs4grid(x,y,grddist,lwrite)
      else
          call recs4grid(x,y,grddist+10000.,lwrite)  !distance from receptor to center of gridded sources
      endif
      
      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 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=facility(ig)%utmx
      fy=facility(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 ig,eof,nv,nv1,j,i,inout1,inout
      
c     facid1:     facility ID being read from file
c     src:        AERMOD source id being read
c     stype:      AERMOD source type
c     line:       header record of file (not used)
c     st:         2-digit state FIPS code
c     fips1:      5-digit FIPS code
      character facid1*20,st,fips1*5,src*12,stype*8,line*100
      
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     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,area1,h,sz,xi,yi,xj,yj
      real*8, allocatable, dimension(:) :: ux !uncommented 17234
      real*8, allocatable, dimension(:) :: uy !uncommented 17234
      logical ix,iy,jx,jy,L_EOR
      
      eof=0
      mindist=900000.0
      
c      if (runcmv) then
      if (runport) then
          nv=cmv(ig)%nvert+1 !+1 added 17234
      else
          nv=tract(ig)%nvert+1 !+1 added 17234
      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)
              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     read in number of vertices and vertex coordinates
c      read(parunit,*,iostat=eof)line
c      read(parunit,*,iostat=eof)st,fips1,facid1,src,stype,area1,h,nv
c 5    if (eof .eq. 0) then
c          if (src .eq. facid .and. src .eq. cmv(ig)%srcid) then
c              nv1=nv
c              allocate(ux(nv1))
c              allocate(uy(nv1))
c              backspace(parunit)
c              read(parunit,*,iostat=eof)st,fips1,facid1,src,stype,area1,
c     +        h,nv,sz,(ux(j),uy(j),j=1,nv)
c              do j=1,nv
cc                 calculate distance of receptor to each edge (shortest distance from receptor to anywhere along edge)
c                  if (j .lt. nv) then
c                      dist=cmvdist(x,y,ux(j),uy(j),ux(j+1),uy(j+1))
c                  else
c                      dist=cmvdist(x,y,ux(j),uy(j),ux(1),uy(1))
c                  endif
c                 get minimum distance of receptor to any polygon
c                  if (dist .le. mindist) mindist=dist
c             enddo
c          endif
c          read(parunit,*,iostat=eof)st,fips1,facid1,src,stype,area1,h,nv
c          goto 5
c      endif
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      c=distance(fx1,fx2,fy1,fy2)
      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
      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 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
c     aprime:     angle used to determine fx1,fy1 for point area source
c     bprime1:    angle used to determine fx2,fy2 for point area source
      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=airport(ig)%utmx1
          fy1=airport(ig)%utmy1
          if (runtyp .eq. 'APRUN') then
              fx2=airport(ig)%utmx2
              fy2=airport(ig)%utmy2
              cprime=airport(ig)%xdim/2.0
          else
              fx2=airport(ig)%utmx1+airport(ig)%xdim/2.0
              fy2=airport(ig)%utmy1+airport(ig)%ydim
              cprime=airport(ig)%xdim/2.0
          endif
      else
          cprime=facility(ig)%xdim/2.0
          a=sqrt((cprime*cprime)+(facility(ig)%ydim*facility(ig)%ydim))
          b=facility(ig)%ydim
          angleb=angle(a,cprime,b)
          anglec=angle(a,b,cprime)
          aprime=3.141592654-facility(ig)%angle-angleb-anglec  !remember angles are in radians, 3.14159 is 180
          fx1=facility(ig)%utmx+(cprime*sin(aprime))
          fy1=facility(ig)%utmy-(cprime*cos(aprime))
          bprime1=(3.141592654/2.0)-facility(ig)%angle-anglec !remember angles are in radians, 3.14159/2 is 90
          fx2=facility(ig)%utmx+(a*cos(bprime1))
          fy2=facility(ig)%utmy+(a*sin(bprime1))
c         fx1=facility(ig)%utmx+facility(ig)%xdim/2.0
c         fy1=facility(ig)%utmy
c         fx2=facility(ig)%utmx+facility(ig)%xdim/2.0
c          fy2=facility(ig)%utmy+facility(ig)%ydim
      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,lwrite)
      use main1
      implicit none
 
c     ig:         source index counter
      integer 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 ig
      
c     cx:         UTM-x coordinate of grid cell center
c     cx:         UTM-y coordinate of grid cell center
      real*8 cx,cy
      
      cx=0.0
      cy=0.0
      cx=(gridded(ig)%utmx1+gridded(ig)%utmx2+gridded(ig)%utmx3+
     +gridded(ig)%utmx4)/(real(gridded(ig)%nvert))
      cy=(gridded(ig)%utmy1+gridded(ig)%utmy2+gridded(ig)%utmy3+
     +gridded(ig)%utmy4)/(real(gridded(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##########################################################################
c     subroutine to get met files for AK, HI, PR, or VI
      subroutine getmetfiles(sfcfil,pflfil,swban,uwban)
      use main1
      implicit none
c     swban:      surface station WBAN
c     uwban:      upper air station WBAN
c     ig:         number of source groups
c     eof:        end of file indicator
c     swban:      surface station WBAN
c     wban1:      surface WBAN read from file
c     uwban:      upper station WBAN
c     uwban1:     upper WBAN read from file   
c     izone:      UTM zone read for station
      integer swban,uwban,ig,eof,swban1,uwban1,izone,wban,wban1,m
      
c     sfcfil:     AERMOD surface meteorological file
c     pflfil:     AERMOD profile meteorological file
c     infil:      file of met stations for state
c     state2a:    state abbreviation
      character sfcfil*250,pflfil*250,infil*250,state2a*2
   
c     fx:         facility average UTM x-coordinate
c     fy:         facility average UTM y-coordinate   
c     mindist:    minimum distance between source and met station
c     dist:       distance between source and met station being read
c     distance:   distance function
c     x:          met station UTM x-coordinate
c     y:          met station UTM y-coordinate
c     elev:       met station elevation (m)
c     profbase:   met station elevation (m) used
      real*8 fx,fy,mindist,distance,dist,x,y,elev
      
c     lexist:     logical variable denoting infil exists
      logical lexist
     
      if (state2 .eq. '02') then
         state2a='AK'
      elseif (state2 .eq. '15') then
         state2a='HI'
      else
         state2a=state2
      endif
 
      fx=0.0
      fy=0.0
c     if point or airport get average coordinates of source
      if (runpt) then
          do ig=1,ngrps
              fx=fx+facility(ig)%utmx
              fy=fy+facility(ig)%utmy
          enddo
          fx=fx/real(ngrps)
          fy=fy/real(ngrps)
      endif
      if (runap .and. runtyp .eq. 'APRUN') then
          if (runtyp .eq. 'APRUN') then
              do ig=1,ngrps
                  fx=fx+airport(ig)%utmx1+airport(ig)%utmx2
                  fy=fx+airport(ig)%utmy1+airport(ig)%utmy2
              enddo
              fx=fx/real(ngrps)
              fy=fy/real(ngrps)
          else
              fx=airport(ig)%utmx1
              fy=airport(ig)%utmy1
          endif
      endif
c     if (runcmv) then
      if (runport) then
          fx=cmv(1)%avgx
          fy=cmv(1)%avgy
      endif
      if (noncontr) then
          call calc_center(1,fx,fy)
c          fx=tract(1)%avgx
c          fy=tract(1)%avgy
      endif
      
      
  
c     now read through the list of met stations for AK or HI and pick
c     the station 
      mindist=9999999.0
      eof=0
      write(infil,'(3(a))')trim(adjustl(ptcmaq)),trim(adjustl(state2a)),
     +'_met.dat'
      inquire (file=infil,exist=lexist)
      if (.not. lexist) then
          write(ilogunit,10)trim(adjustl(infil))
          stop
      endif
      
      open(unit=metunit,file=infil,status='old')
      read(metunit,*,iostat=eof)wban,uwban1,x,y,izone,elev
 5    if (eof .eq. 0) then
          if (izone .eq. utmzone) then
              dist=distance(x,y,fx,fy)
              if (dist .lt. mindist) then
                  mindist=dist
                  swban=wban
                  uwban=uwban1
                  profbase=elev
              endif
          endif
          read(metunit,*,iostat=eof)wban,uwban1,x,y,izone,elev
          goto 5
      endif

c     find the correct met directory to use
      lexist=.false.
      m=1
      do while (m .le. 2 .and. .not. lexist)
          write(sfcfil,'(2(a),i5,4a)')trim(adjustl(metdir)),
     +    trim(adjustl(state2a)),swban,'.sfc'
          inquire (file=sfcfil,exist=lexist)
          m=m+1
      enddo
      m=m-1
      write(pflfil,'(2(a),i5,4a)')trim(adjustl(metdir)),
     +trim(adjustl(state2a)),swban,'.pfl'
      
 10   format('File ',a,' does not exist, stopping program')
      return
      end
          
c##########################################################################
c     subroutine to write EMISFACT to AERMOD.INP file
      subroutine emisfact1
      use main1
      implicit none

c     version 17234, add seasonal emissions by source for onroad and nonroad

c     eof:        end of file indicator
c     nfact:      number of AERMOD emissions factors based on qflag
c     i:          source loop counter
c     i1:         emissions factor counter
c     j:          loop counter for writing emissions factors to AERMOD.INP
c     k:          loop counter for writing emissions factors to AERMOD.INP
c     esrc:       integer flag if source is hourly source(esrc=1) or not (esrc=0)
      integer eof,nfact,i,i1,assignfact,k,j,esrc

c     line:       header record of file (not used)
c     src:        AERMOD source id being read from file
c     facid1:     facility ID being read from file
c     j1:         junk variable (source name)
c     qf:         AERMOD emissions factor type
c     src1:       AERMOD source ID being processed

      character line*100,facid1*20,j1*75,src*12,qf*8,src1*12,st2*2,
     +fips6*6,tract1*11,poly*4,facid2*20
 
      real*8 ann_emis,monthwt,monwt,calcemis
c      real*8 calcemis,monthwt,monwt
c     emisfact:   AERMOD emissions factors
c      real*8, allocatable, dimension(:) :: emisfact
c     emisfact1:  allocated emissions, match emissions factors
c      real*8, allocatable, dimension(:) :: emission
c     lfound:     source found (to eliminate duplicates)
      logical, allocatable, dimension(:) ::  lfound
      
      allocate(lfound(ngrps))
      
c     first check to see if any of the qflags are MHRDOW or MHRDOW7
c      ldow3=.false.
c      ldow7=.false.
c      do i=1,ngrps
c          if (trim(adjustl(facility(i)%qflag)) .eq. 'MHRDOW') 
c     +    ldow3=.true.
c          if (trim(adjustl(facility(i)%qflag)) .eq. 'MHRDOW7') 
c     +    ldow7=.true.
c      enddo
      
c      if (ldow3 .or. ldow7) then
c          if (ldow3) then
c              allocate(df3(864))
c              allocate(nday3(12,3))
c          endif
c          if (ldow7) then
c              allocate(nday7(12,7))
c              allocate(df7(2016))
c          endif
c          call monthtot
c      endif
      
      open(unit=temunit,file=temfil(1),status='old')
      read(temunit,*,iostat=eof)line
      
      
      
      lfound=.false.
      do i=1,ngrps
          call get_source(i,src1) !get source
          if (runpt) then
              esrc=facility(i)%hrsrc !determine if hourly emissions source or not (1=hourly, 0=not hourly)
          else
              esrc=0
          endif
          if (esrc .eq. 0) then
              if (runpt .or. runap) then
                  read(temunit,*,iostat=eof)facid1,j1,src,qf
c              elseif (runcmv) then  !check input JAT
c                  read(temunit,*,iostat=eof)st2,facid1,src,qf
c              elseif (noncontr) then
c                  read(temunit,*,iostat=eof)j1,st2,fips6,tract1,poly,
c     +            facid2,src,qf
              else
                  read(temunit,*,iostat=eof)j1,facid1,src,qf
              endif
 5            if (eof .eq. 0) then
c                  if (noncontr) facid1=src1 !need to check JAT
                  if (facid1 .eq. facid .and. src .eq. src1 .and. 
     +            .not. lfound(i)) then
c    +            .not. runcmv) .or. ((runcmv .or. noncontr) .and. src 
c    +            .eq. src1)) .and. .not. lfound(i)) then
                      lfound(i)=.true.
                      nfact=assignfact(qf)
                      allocate(emisfact(nfact))
                      allocate(emission(nfact))
                      backspace(temunit)
                      if (runpt .or. runap) then
                          read(temunit,*,iostat=eof)facid1,j1,src,qf,
     +                    (emisfact(i1),i1=1,nfact)
c                      elseif (runcmv) then
c                          read(temunit,*,iostat=eof)st2,facid1,src,qf,
c     +                    (emisfact(i1),i1=1,nfact)
c                      elseif (noncontr) then
c                          read(temunit,*,iostat=eof)j1,st2,fips6,tract1,
c     +                    poly,facid2,src,qf,(emisfact(i1),i1=1,nfact)
                      else
                          read(temunit,*,iostat=eof)j1,facid1,src,qf,
     +                    (emisfact(i1),i1=1,nfact)
                      endif
                      if (eof .eq. 0) then
c                         calculate emissions
c                         calculate weighting factor for MONTH factors, consistent with SMOKE-CMAQ
                          if (trim(adjustl(qf)) .eq. 'MONTH') then
                              monthwt=monwt(qf)
                          else
                              monthwt=0.0
                          endif
                          do i1=1,nfact
                              emission(i1)=calcemis(i,i1,monthwt,qf)                       
c                              if (runon .or. runnon) 
c     +                        call calcseasons(emission(i1),i1,qf)
                          enddo
c                         calculate annual emissions, if onroad or nonroad calculate seasonal emissions
                          call calcann(i,qf)
                          
                          if (nfact .eq. 4) then
                              write(aerunit,15)src,qf,
     +                        (emission(i1),i1=1,nfact)
                          else
                              do j=1,nfact,6
                                  write(aerunit,20)src,qf,
     +                            (emission(k),k=j,j+5)
                              enddo
                              deallocate(emisfact)
                              deallocate(emission)
                          endif
                      endif
                  endif
                  if (runpt .or. runap) then
                      read(temunit,*,iostat=eof)facid1,j1,src,qf
c                  elseif (runcmv) then
c                      read(temunit,*,iostat=eof)st2,facid1,src,qf
c                  elseif (noncontr) then
c                      read(temunit,*,iostat=eof)j1,st2,fips6,tract1,
c     +            poly,facid2,src,qf
                  else
                      read(temunit,*,iostat=eof)j1,facid1,src,qf
                  endif
                  goto 5
              endif           
          endif
          rewind(temunit)
          read(temunit,*,iostat=eof)line
      enddo
      
      close(temunit)
 15   format('SO EMISFACT',1x,a12,1x,a8,4(1x,e13.6))
 20   format('SO EMISFACT',1x,a12,1x,a8,6(1x,e13.6))
      return
      end
c##########################################################################
c     get source ID for use in subroutine emisfact1
      subroutine get_source(i,srcid1)
      use main1
      implicit none

c     i:      source id index in source array
      integer i
      
c     srcid1:     source ID returned from appropriate array
      character srcid1*12
      
      if (runpt) then
          srcid1=facility(i)%srcid
      elseif (runap) then
          srcid1=airport(i)%srcid
c      elseif (runcmv) then
      elseif (runport) then
          srcid1=cmv(i)%srcid
c      elseif (noncontr) then
c          srcid1=tract(i)%srcid
      else
          srcid1=gridded(i)%srcid
      endif
      return
      end
      
c########################################################################## 
c     assign # of emissions factors based on qflag
      function assignfact(qf)
      use main1
      implicit none
      
c     assignfact:     number of emissions factors based on qf
      integer assignfact
      
c     qf:             AERMOD emissions factor type
      character qf*8
      
      if (qf .eq. 'SEASON') then
          assignfact=4
      elseif (qf .eq. 'MONTH') then
          assignfact=12
      elseif (qf .eq. 'HROFDY') then
          assignfact=24
      elseif (qf .eq. 'SEASHR') then
          assignfact=96
      elseif (qf .eq. 'HRDOW') then
          assignfact=72
      elseif (qf .eq. 'HRDOW7') then
          assignfact=168
      elseif (qf .eq. 'SHRDOW') then
          assignfact=288
      elseif (qf .eq. 'SHRDOW7') then
          assignfact=672
      elseif (qf .eq. 'MHRDOW') then
          assignfact=864
      elseif (qf .eq. 'MHRDOW7') then
          assignfact=2016
      else
           write(ilogunit,*)'invalid QFLAG, stopping'
          stop
      endif
      return
      end
c##########################################################################
c     function to calculate month weight
      function monwt(qf)
      use main1
      implicit none
      character qf*8
      real*8 mw,monwt
      integer i,d
      
      monwt=0.0

      do i=1,12
          if (i .eq. 2 .and. lleap) then
              d=days(i)+1
          else
              d=days(i)
          endif
          monwt=monwt+(real(d)*emisfact(i))
      enddo
c     write(*,*)'monwt ',monwt
      return
      end

c##########################################################################
c     function to temporally allocate emissions
      function calcemis(igrp,i1,mw,qf)
      use main1
      implicit none
      real*8 calcemis,em1,frac,area,mw
      integer i1,igrp,ndays,is,i
      character*8 qf
      
      if (runpt) then
          frac=1.0
          if (trim(adjustl(facility(igrp)%srctyp)) .eq. 'AREA') then
              area=facility(igrp)%xdim*facility(igrp)%ydim
          else
              area=1.0
          endif
      elseif (runap) then
          if (trim(adjustl(runtyp)) .eq. 'APRUN') then
              frac=airport(igrp)%frac
              area=airport(igrp)%area
          else
              frac=1.0
              area=airport(igrp)%xdim*airport(igrp)%ydim
          endif
c      elseif (runcmv) then
      elseif (runport) then  
          frac=1.0
          area=cmv(igrp)%area
c      elseif (noncontr) then
c          frac=1.0
c          area=tract(igrp)%area
      else
          frac=1.0
          area=real(gridres*1000*gridres*1000)
      endif
      
      if (trim(adjustl(qf)) .eq. 'MONTH') then !monthly
          if (lleap .and. i1 .eq. 2) then
              ndays=days(i1)+1
          else
              ndays=days(i1)
          endif
          if (i1 .le. 2 .or. i1 .eq. 12) then
              is=1
          elseif (i1 .le. 5) then
              is=2
          elseif (i1 .le. 8) then
              is=3
          else
              is=4
          endif
          em1=(annemis*(emisfact(i1)/mw))/24.0 !consistent with SMOKE-CMAQ
c        write(*,*)i1,emisfact(i1),mw,em1
c         g/s       tpy    1/# avg days/month 1/24  tons/hr to g/s
c          em1=10000.0*em*(12.0/365.0)*(1/24)*251.9957778   !converts tons/year to g/s for the month
c          em1=10000.0*em*(ndays)*(1/24)*251.9957778   !converts tons/year to g/s for the month
      elseif (trim(adjustl(qf)) .eq. 'SEASON' .or. trim(adjustl(qf)) 
     +.eq. 'SEASHR') then !seasonal or seasonal hour of day
          if (i1 .eq. 1) then
              ndays=days(1)+days(2)+days(12)
          elseif (i1 .eq. 2) then
              ndays=days(3)+days(4)+days(5)
          elseif (i1 .eq. 3) then
              ndays=days(6)+days(7)+days(8)
          else
              ndays=days(9)+days(10)+days(11)
          endif
          if (trim(adjustl(qf)) .eq. 'SEASON') then
              em1=annemis*emisfact(i1)*(1.0/real(ndays))*(1/24)   !converts tons/year to g/s for the season
          else
              em1=annemis*emisfact(i1)*(1.0/real(ndays))   !converts tons/year to g/s for the season-hour-of-day
          endif
      elseif (trim(adjustl(qf)) .eq. 'HROFDY') then   !converts tons/year to g/s for the season) then !seasonal hour of day
          em1=(annemis*emisfact(i1))/365.0
          do is=1,4
              if (is .eq. 1) then
                  ndays=days(1)+days(2)+days(12)
                  if (lleap)ndays=ndays+1
              elseif (is .eq. 2) then
                  ndays=days(3)+days(4)+days(5)
              elseif (is .eq. 3) then
                  ndays=days(6)+days(7)+days(8)
              else
                  ndays=days(9)+days(10)+days(11)
              endif
          enddo
      else 
          em1=annemis*emisfact(i1) !*7.0 !converts tons/year to tons/hr for the month-hour-day-of-week or day-of-week7
      endif
      calcemis=(em1*frac/area)*251.9957778 !convert to g/s and/or area rate
c     write(*,*)calcemis,em1,frac,area
      return
      end
c##########################################################################  
c     subroutine to calculate annual emissions from the hourly emissions rate
      subroutine calcann(igrp,qf)
      use main1
      implicit none
      character qf*8
      integer igrp,im,idy,ihr,ndays,ii
      real*8 area,em,ann_emis,memis(12)
      
      ann_emis=0.0
      memis=0.0
      
      if (runpt) then
          if (trim(adjustl(facility(igrp)%srctyp)) .eq. 'AREA') then
              area=facility(igrp)%xdim*facility(igrp)%ydim
          else
              area=1.0
          endif
      elseif (runap) then
          if (trim(adjustl(runtyp)) .eq. 'APRUN') then
              area=airport(igrp)%area
          else
              area=airport(igrp)%xdim*airport(igrp)%ydim
          endif
c      elseif (runcmv) then
      elseif (runport) then
          area=cmv(igrp)%area
c      elseif (noncontr) then
c          area=tract(igrp)%area
      else
          area=real(gridres*1000*gridres*1000)
      endif
      
      do im=1,12
          if (lleap .and. im .eq. 2) then
              ndays=days(im)+1
          else
              ndays=days(im)
          endif
c          iseas=seasons(im)
          do idy=1,ndays
              do ihr=1,24
                  if (trim(adjustl(qf)) .eq. 'HROFDY') then
                      ii=ihr
                  elseif (trim(adjustl(qf)) .eq. 'MONTH') then
                      ii=im
                  else
                      call weekdays(qf,im,idy,ihr,ii)
                  endif
                  em=(emission(ii)*area)/251.9957778  !convert back to ton/hr
                  memis(im)=memis(im)+em
                  ann_emis=ann_emis+em
              enddo
          enddo
      enddo
      if (runpt) then
         facility(igrp)%annual=ann_emis
      elseif (runap) then
          airport(igrp)%annual=ann_emis
c      elseif (runcmv) then
      elseif (runport) then
          cmv(igrp)%annual=ann_emis
c      elseif (noncontr) then
c          tract(igrp)%annual=ann_emis
c          do im=1,12
c              tract(igrp)%monemis(im)=memis(im)
c          enddo
      else
          gridded(igrp)%annual=ann_emis
          do im=1,12
              gridded(igrp)%monemis(im)=memis(im)
          enddo
      endif
      return
      end
c##########################################################################  
      subroutine writehour
      use main1
      implicit none
      
c     eof:        end of file indicator
c     iyear:      2-digit year
c     imonth:     2-digit month
c     iday:       2-digit day
c     ihour:      2-digit hour
c     i:          source index from array
c     findsrc:    function to find source index
c     ld:         leap year addition for Feb. 28 (non-leap=0,leap=1)
c                 set to 0 for other months
c     iyear4:     4-digit year
      integer eof,iyear,imonth,iday,ihour,i,findsrc,ld,iyear4,iyear2
      
c     facid1:     facility id read from file
c     src:        source ID read from file
c     hrfil:      hourly emissions file
      character facid1*20,src*12,hrfil*250
      
c     ef:         emissions factor
c     xd:         x-dimension of area source
c     yd:         y-dimension of area source
c     emis:       array of hourly emissions
c     ef2:        emissions multiplier (8760000)
      real*8 ef,xd,yd,ef2
      real*8, allocatable, dimension(:,:,:,:) :: emis
      
c     lleap:      logical variable denoting if year is leap year (true=leap year)
c      logical lleap
      
      write(hrfil,'(4(a))')trim(adjustl(facid)),'_',
     +trim(adjustl(state2)),'_hourly.dat'
c     write HOUREMIS keyword to aermod input file
      do i=1,ngrps
          if (facility(i)%hrsrc .eq. 1) 
     +    write(aerunit,10)trim(adjustl(hrfil)),facility(i)%srcid  
      enddo
      
      allocate(emis(ngrps,12,31,24))
      emis=-999.0
      ef2=annemis
      
      open(unit=temunit,file=temfil(2),status='old')
      read(temunit,*,iostat=eof)facid1
      read(temunit,*,iostat=eof)facid1,src,iyear,imonth,iday,ihour,ef
 5    if (eof .eq. 0) then
          i=findsrc(src)
          if (iyear .lt. 100) then
              iyear4=iyear+2000
              iyear2=iyear
          else
              iyear4=iyear
              iyear2=iyear-2000
          endif
          ihour=ihour+hrfix
c          call leapyr(iyear4,lleap)
          if (facility(i)%srctyp .eq. 'AREA') then
              xd=facility(i)%xdim
              yd=facility(i)%ydim
              emis(i,imonth,iday,ihour)=ef*ef2*251.9957778/(xd*yd)  
          else
              emis(i,imonth,iday,ihour)=ef*ef2*251.9957778  
          endif
          if (emis(i,imonth,iday,ihour) .lt. 0.0)
     +    emis(i,imonth,iday,ihour)=0.0
          facility(i)%annual=facility(i)%annual+(ef*ef2)
          read(temunit,*,iostat=eof)facid1,src,iyear,imonth,iday,ihour,
     +    ef
          goto 5
      endif
      close(temunit)
c     now write the hourly emissions file
      open(unit=hrunit,file=hrfil,status='unknown')
c     CHANGE IYEAR TO 14 WHEN GET REAL FILES!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c      iyear=11
      do imonth=1,12
          if (lleap .and. imonth .eq. 2) then
              ld=1
          else
              ld=0
          endif
          do iday=1,days(imonth)+ld
              do ihour=1,24
                  do i=1,ngrps
                      if (emis(i,imonth,iday,ihour) .ge. 0.0) then
                          if (facility(i)%srctyp .eq. 'AREA') then
                              write(hrunit,15)iyear2,imonth,iday,ihour,
     +                     facility(i)%srcid,emis(i,imonth,iday,ihour)
                          else
                              write(hrunit,20)iyear2,imonth,iday,ihour,
     +                     facility(i)%srcid,emis(i,imonth,iday,ihour),
     +                     facility(i)%stktemp,facility(i)%stkvel
                          endif
                      endif
                  enddo
              enddo
          enddo
      enddo
      
      deallocate(emis)
      close(hrunit)
 10   format('SO HOUREMIS ',a,1x,a12)
 15   format('SO HOUREMIS',4(1x,i2),1x,a12,1x,e13.4)
 20   format('SO HOUREMIS',4(1x,i2),1x,a12,1x,e13.4,2(1x,f10.5))      

      
      return
      end
c##########################################################################
c     subroutine to write out hourly RWC or AG file
      subroutine rwc2grid
      use main1
      implicit none
 
c     eof:        end of file indicator

c     iyear:      2-digit year
c     modyr2:     2-digit model year
c     imonth:     2-digit month
c     iday:       2-digit day
c     ihour:      2-digit hour
c     i:          source loop counter
c     ld:         leap year addition for Feb. 28 (non-leap=0,leap=1)
c                 set to 0 for other months
      integer iyear,imon,iday,ihr,eof,i,ld,ist,ig,modyr2
      
c     hfactor:    hourly emissions factor read from file
c     emis:       hourly emissions written to hourly emissions file
      real*8 hfactor
      real*8, allocatable, dimension(:,:,:,:) :: emis
c     line:       header line read from file (not used)
c     hrfil:      hourly emissions file
c     afips:      FIPS read from file
      character line*100,hrfil*250,afips*5,ares*2,j1*10,src*12,str*20,
     +astate*3
      
c     lleap:      logical variable denoting if year is leap year (true=leap year)
      logical lexist,lfound
      
      modyr2=modyear-(modyear/100)*100

      allocate(emis(ngrps,12,31,24))
      emis=-999.0
      
      if (runag .or. runrwc) then !(runrwc .and. noncontr)) then
        str='_hourly.csv'
      else
        str='_hourly_conus.csv'
      endif
      write(ares,'(i2)')gridres
c     if (.not. noncontr) then
          do ist=1,nstate
              if (noncontr) then
                if (stfips(ist) .eq. '02') then
                  astate='AK_'
                elseif (stfips(ist) .eq. '15') then
                  astate='HI_'
                else
                  astate='PR_'
                endif
              else
                astate='_'
              endif
              write(temfil(1),'(6a)')trim(adjustl(pttemp)),
     +        trim(adjustl(runtyp)),trim(adjustl(ares)),
     +        trim(adjustl(astate)),stfips(ist),trim(adjustl(str))
c             write(*,*)temfil(1)
              inquire(file=temfil(1),exist=lexist)
              if (.not. lexist) then
                  write(ilogunit,65)trim(adjustl(facid))
                  stop
              endif
              open(temunit,file=temfil(1),status='old')
              eof=0
              read(temunit,*,iostat=eof)line
              read(temunit,*,iostat=eof)j1,afips,iyear,imon,iday,ihr,
     +        hfactor
 5            if (eof .eq. 0) then
                  do ig=1,ngrps
                      if (afips .eq. fips(ig)) then
                        emis(ig,imon,iday,ihr)=hfactor*annemis*
     +                  251.9957778/real(gridres*1000*gridres*1000)
                        if (emis(ig,imon,iday,ihr) > 0.0) 
     +                  emis(ig,imon,iday,ihr)=0.0 
                        gridded(ig)%annual=gridded(ig)%annual+
     +                  (hfactor*annemis)
                     endif
                  enddo
                  read(temunit,*,iostat=eof)j1,afips,iyear,imon,iday,
     +            ihr,hfactor
                  goto 5
              endif
              close(temunit)
          enddo
c     else
c         lfound=.false.
c         open(temunit,file=temfil(1),status='old')
c         eof=0
c         read(temunit,*,iostat=eof)line
c         read(temunit,*,iostat=eof)j1,afips,iyear,imon,iday,ihr,
c    +    hfactor
c25      if (eof .eq. 0) then
c             do ig=1,ngrps
c                 if (afips .eq. tract(ig)%fips) then
c                     lfound=.true.
c                     emis(ig,imon,iday,ihr)=hfactor*annemis*
c    +                251.9957778/tract(ig)%area
c                     tract(ig)%annual=tract(ig)%annual+
c    +                (hfactor*annemis)
c                 endif
c             enddo
c             read(temunit,*,iostat=eof)j1,afips,iyear,imon,iday,
c    +        ihr,hfactor
c             goto 25
c         endif
c         if (.not. lfound) then
c             write(ilogunit,65)trim(adjustl(facid))
c             stop
c         endif
c         close(temunit)
c     endif
      
      write(hrfil,'(2a,i2.2,3a)')trim(adjustl(runtyp)),'_',gridres,'_',
     +trim(adjustl(facid)),'_hourly.dat'
      open(unit=hrunit,file=hrfil,status='unknown')
      do ig=1,ngrps
          src=gridded(ig)%srcid
c         if (runon) then
c             src=onsrc(ig)
c         elseif (runrwc .and. .not. noncontr) then
c             src=gridded(ig)%srcid
c         else
c             src=tract(ig)%srcid
c         endif
          if (iyear .ne. modyr2) iyear=modyr2 
          do imon=1,12
              do iday=1,31
                  do ihr=1,24
                      if (emis(ig,imon,iday,ihr) .ge. 0.0)
     +                write(hrunit,10)iyear,imon,iday,ihr,
     +                src,emis(ig,imon,iday,ihr)
c     +                gridded(ig)%srcid,emis(ig,imon,iday,ihr)
                  enddo
              enddo
          enddo
      enddo
      
c      open(temunit,file=temfil(1),status='old')
c      open(unit=hrunit,file=hrfil,status='unknown')
c      eof=0
c      read(temunit,*,iostat=eof)line
c      read(temunit,*,iostat=eof)afips,iyear,imonth,iday,ihour,hfactor
c 5    if (eof .eq. 0) then
c          if (afips .eq. fips(1)) then
c              call leapyr(iyear,lleap)
c              if (lleap .and. imonth .eq. 2) then
c                  ld=1
c              else
c                  ld=0
c              endif
c              if (iday .le. days(imonth)+ld) then
c                  do i=1,ngrps
c                      emis=1000.0*hfactor/(gridres*1000*gridres*1000)
c                      write(hrunit,10)iyear-2000,imonth,iday,ihour,
c     +                gridded(i)%srcid,emis
c                  enddo
c              endif
c          endif
c          read(temunit,*,iostat=eof)afips,iyear,imonth,iday,ihour,
c     +    hfactor
c          goto 5
c      endif
      
      do i=1,ngrps
c         if (.not. noncontr) then
              write(aerunit,15)trim(adjustl(hrfil)),gridded(i)%srcid  
c         else
c             write(aerunit,15)trim(adjustl(hrfil)),tract(i)%srcid
c         endif
      enddo
      
c      close(temunit)
      close(hrunit)
 10   format('SO HOUREMIS',4(1x,i2),1x,a12,1x,e13.4) 
 15   format('SO HOUREMIS',1x,a,1x,a12)
 65   format('RWC2GRID temporal file(s) not found for ',a) 
      return
      end
c##########################################################################
c     write onroad hourly files from county file
      subroutine onrd2grid
      use main1
      implicit none
      integer eof,iyear,imon,iday,ihr,ig,k,ist,is,iyear2,iyear4,ld,isrc,
     +findsrc
      character line*100,src*12,j1*10,afips*5,ares*2,hrfil*250,fips1*5,
     +onsrc1*12,astate*2,str*17
      real*8 efact
      logical lexist
      real*8, allocatable, dimension(:,:,:,:) :: emis
      allocatable :: fips1(:)
      allocatable :: onsrc1(:)
      
      allocate(fips1(ngrps))
      allocate(onsrc1(ngrps))
      

c      real*8, allocatable, dimension(:,:,:,:) :: emisfact2
      write(hrfil,'(2a)')trim(adjustl(facid)),'_hourly.dat'

      do ig=1,ngrps
          write(aerunit,10)trim(adjustl(hrfil)),gridded(ig)%srcid 
c         find the index of the sources in the onsrc array in the gridded array
c         and write to the onsrc1 and fips1 arrays
          isrc=findsrc(onsrc(ig))
          onsrc1(isrc)=onsrc(ig)   
          fips1(isrc)=fips(ig)
      enddo

c     re-order the fips and onsrc arrays to match the gridded srcid array
      fips=''
      onsrc=''
      fips=fips1
      onsrc=onsrc1
      deallocate(fips1)
      deallocate(onsrc1)
      
c      allocate(emisfact2(ngrps,3,12,24))
      allocate(emis(ngrps,12,31,24))
      emis=-999.0
c      emisfact=0.0
c      qf='MHRDOW'
      write(ares,'(i2)')gridres
      if (runrwc) then
         str='_hourly_conus.csv'
      else
         str='_hourly.csv'
      endif
      do ist=1,nstate
          if (noncontr) then
              if (stfips(ist) .eq. '02') then
                 astate='AK'
              elseif (stfips(ist) .eq. '15') then
                 astate='HI'
              else
                 astate='PR'
              endif
              write(temfil(1),'(7(a))')trim(adjustl(pttemp)),
     +        trim(adjustl(runtyp)),trim(adjustl(ares)),astate,
     +        '_',stfips(ist),trim(adjustl(str))
          else
             write(temfil(1),'(6(a))')trim(adjustl(pttemp)),
     +       trim(adjustl(runtyp)),trim(adjustl(ares)),'_',stfips(ist),
     +       trim(adjustl(str))
          endif
c          write(*,*)temfil(1)
          inquire(file=temfil(1),exist=lexist)
          if (.not. lexist) then
              write(ilogunit,65)trim(adjustl(facid))
              stop
          endif
          open(temunit,file=temfil(1),status='old')
          eof=0
          read(temunit,*,iostat=eof)line
          read(temunit,*,iostat=eof)j1,afips,iyear,imon,iday,ihr,efact
 5        if (eof .eq. 0) then
              do ig=1,ngrps
                  if (afips .eq. fips(ig)) then
                      if (iyear .lt. 100) then
                          iyear4=iyear+2000
                          iyear2=iyear
                      else
                          iyear4=iyear
                          iyear2=iyear-2000
                      endif
                      if (iyear4 .eq. modyear)
     +                emis(ig,imon,iday,ihr)=efact*annemis*251.9957778/
     +                real(gridres*1000*gridres*1000)
                      if (emis(ig,imon,iday,ihr) .lt. 0.0) 
     +                emis(ig,imon,iday,ihr)=0.0
                      gridded(ig)%annual=gridded(ig)%annual+
     +                (efact*annemis)
                      gridded(ig)%monemis(imon)=
     +                gridded(ig)%monemis(imon)+(efact*annemis)
                  endif
c                  if (imon .le. 2 .or. imon .eq. 12) then
c                      is=1
c                  elseif (imon .le. 5) then
c                      is=2
c                  elseif (imon .le. 8) then
c                      is=3
c                  else
c                      is=4
c                  endif
c                  seasemis(is)=seasemis(is)+(annemis*efact)
              enddo
              read(temunit,*,iostat=eof)j1,afips,iyear,imon,iday,ihr,
     +        efact
              goto 5
          endif
          close(temunit)
      enddo
      
c     write emissions factors to AERMOD input file
      open(unit=hrunit,file=hrfil,status='replace')
      do imon=1,12
          if (lleap .and. imon .eq. 2) then
              ld=1
          else
              ld=0
          endif
          do iday=1,days(imon)+ld
              do ihr=1,24
                  do ig=1,ngrps
                     src=onsrc(ig)
c                      src=gridded(ig)%srcid
                      write(hrunit,15)iyear2,imon,iday,ihr,src,
     +                emis(ig,imon,iday,ihr)
                  enddo
              enddo
          enddo
      enddo
     
c      deallocate(emisfact2)
      deallocate(emis)
      close(hrunit)
c 10   format('SO EMISFACT',1x,a12,1x,a8,6(1x,e13.6))   
 10   format('SO HOUREMIS ',a,1x,a12)
 15   format('SO HOUREMIS',4(1x,i2),1x,a12,1x,e13.4)
 65   format('temporal file(s) not found for ',a) 
      return
      end     
      
c##########################################################################
      subroutine checkmet(sfc,pfl)
      use main1
      implicit none

c     sfc:        AERMOD surface filename
c     pfl:        AERMOD profile filename
      character sfc*250,pfl*250
      
c     lsfc:       surface file exists
c     lpfl:       profile file exists
      logical lsfc,lpfl
      
      lsfc=.false.
      lpfl=.false.
      
      inquire(file=sfc,exist=lsfc)
      inquire(file=pfl,exist=lpfl)
      
      if (.not. lsfc) write(ilogunit,5)sfc
      if (.not. lpfl)write(ilogunit,5)pfl
      if (.not. lsfc .or. .not. lpfl)then
          write(ilogunit,10)
c          stop
      endif
 5    format(a,' does not exist')
 10   format('AERMOD will abort')
      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 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##########################################################################
c      subroutine leapyr(iyear,leap)
      subroutine leapyr
      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 iyear,ileap4,ileap100,ileap400
c      logical leap
      
      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##########################################################################
c     subroutine to make sure grid resolution from input argument is valid
      subroutine checkgrid(agrid)
      use main1
      implicit none
      integer i,ii,ilen
      character agrid*2
      logical lbad
      
      gridres=0
      
      
      ilen=len_trim(agrid)
      i=1
      lbad=.false.
      do while (i .le. ilen .and. .not. lbad)
          ii=ichar(agrid(i:i))
          if (ii .lt. 48 .or. ii .gt. 57) lbad=.true.
          i=i+1
      enddo
      if (.not. lbad) read(agrid,*)gridres
      return
      end
c##########################################################################
c     subroutine to get # of days per month
      subroutine weekdays(qf,imonth,iday,ih,ii)
      use main1
      implicit none
      integer ia,imonth,iy,im,id,iday,id3,id7,ih,ii
      character qf*8
      
     
      IA = (14-IMONTH)/12
      IY = modyear - IA
      IM = IMONTH + 12*IA - 2
      ID = MOD((IDAY+IY+IY/4-IY/100+IY/400+(31*IM)/12),7)
      if (id .ge. 1 .and. id .le. 5) then
          id3=1 !weekday
      elseif (id .eq. 6) then
          id3=2 !saturday
      elseif (id .eq. 0) then
          id3=3 !sunday
      endif
      if (id .eq. 0) then
           id7=7
      else
          id7=id
      endif
          
      if (trim(adjustl(qf)) .eq. 'MHRDOW') then
          ii=ih+(imonth-1)*24+(id3-1)*288
      else
          ii=ih+(imonth-1)*24+(id7-1)*288
      endif
c     if (ii .eq. 865) write(*,*)ii,imonth,iday,ih
      return
      end

c##########################################################################
c     subroutine to write out annual emissions
      subroutine writeann
      use main1
      implicit none
      integer igrp
      character str*100,src1*12
      real*8 diff,ann_emis
      
      write(ilogunit,4)
      if (runap) then !do facility level, since unit emissions divided among sources
          ann_emis=0.0
          do igrp=1,ngrps
              ann_emis=ann_emis+airport(igrp)%annual/airport(igrp)%frac
          enddo
          ann_emis=ann_emis/real(ngrps)
          diff=((ann_emis-annemis)/annemis)*100.0
          if (abs(diff) .ge. 0.0) then
              write(str,'(2(a))')'WARNING modeled emissions ',
     +        'differ from input emissions'
          else
              str='Modeled emissions = input emissions'
          endif
          write(ilogunit,5)trim(adjustl(facid)),ann_emis,annemis,diff,
     +    trim(adjustl(str))
      else
          do igrp=1,ngrps
              if (runpt) then
                  ann_emis=facility(igrp)%annual
                  src1=facility(igrp)%srcid
c              elseif (runcmv) then
              elseif (runport) then
                  ann_emis=cmv(igrp)%annual
                  src1=cmv(igrp)%srcid
c              elseif (noncontr) then
c                  ann_emis=tract(igrp)%annual
c                  src1=tract(igrp)%srcid
              else
                  ann_emis=gridded(igrp)%annual
                  src1=gridded(igrp)%srcid
              endif
              diff=((ann_emis-annemis)/annemis)*100.0
              if (abs(diff) .ge. 1.0e-06) then
                  write(str,'(2(a))')'WARNING modeled emissions ',
     +            'differ from input emissions'
              else
                  str='Modeled emissions = input emissions'
              endif
              write(ilogunit,5)trim(adjustl(src1)),ann_emis,annemis,
     +        diff,trim(adjustl(str))
          enddo
          
      endif
 4    format('Source    Modeled    Unit rate   Difference')     
 5    format(a,1x,2(f8.1,1x),e8.1,1x,a)
      return
      end
c##########################################################################
c     subroutine to write out met pathway
      subroutine writemet
c     wrfcol:     WRF grid cell column
c     wrfrow:     WRF grid cell row
c     swban:      surface station WBAN (99999 for WRF)
c     uwban:      upper air station WBAN (99999 for WRF)
c     ioffset:    WRF to CMAQ offset WRF=CMAQ+ioffset
c     sfcfil:     AERMOD surface meteorological file
c     pflfil:     AERMOD profile meteorological file
      use main1
      implicit none
      integer eof,icol,irow,wrfcol,wrfrow,ioffset(2),ilen,swban,uwban
      real*8 elev,metdist,mindist
      character sfcfil*250,pflfil*250,wrfelevs*250,dir*30,acol*3,arow*3
      logical lfound,lexist
      
c     WRF column is offset by 18 from CMAQ column
c     WRF row is offset by 15 from CMAQ row

      
c     for non-conus domains, read the elevation file to get the
c     cell's elevation, for conus use file to get filename tag
      if (state2 .eq. '02' .or. state2 .eq. 'AK') then
          metdom='9AK'
          ioffset=6
      elseif (state2 .eq. '15' .or. state2 .eq. 'HI') then
          metdom='3HI'
          ioffset=9
      elseif (state2 .eq. '72' .or. state2 .eq. 'PR' .or. state2 .eq. 
     +'78' .or. state2 .eq. 'VI') then
          metdom='3PR'
          ioffset=9
      else
          metdom='12_CONUS'
          ioffset(1)=18
          ioffset(2)=15
      endif
      
      wrfcol=col+ioffset(1)
      wrfrow=row+ioffset(2)
      write(acol,'(i3)')wrfcol
      write(arow,'(i3)')wrfrow
      
      lfound=.false.
      profbase=0.0
      swban=99999
      uwban=99999
      sfcfil='NA'
      pflfil='NA'
 
      write(wrfelevs,'(3(a))')trim(adjustl(ptcmaq)),
     +trim(adjustl(metdom)),'_wrf_elevations.csv'
      
      inquire(file=wrfelevs,exist=lexist)
      if (.not. lexist) then
           write(ilogunit,5)trim(adjustl(wrfelevs))
          stop
      endif
      open(unit=wrfunit,file=wrfelevs,status='old')
      eof=0 
      if (.not. noncontr .and. .not. akhiprvi) then
          read(wrfunit,*,iostat=eof)dir,icol,irow,elev
      else
          read(wrfunit,*,iostat=eof)icol,irow,elev
      endif
 10   if (eof .eq. 0) then
          if (icol .eq. wrfcol .and. irow .eq. wrfrow) then  !match
              lfound=.true.
              if (.not. noncontr .and. .not. akhiprvi) then !make filenames
                  ilen=len_trim(dir)
                  write(sfcfil,'(5(a))')trim(adjustl(metdir)),
     +            trim(adjustl(metdom)),'/MMIF_AERMET_2018_',
     +            trim(adjustl(dir(6:ilen))),'.SFC'
                  write(pflfil,'(5(a))')trim(adjustl(metdir)),
     +            trim(adjustl(metdom)),'/MMIF_AERMET_2018_',
     +            trim(adjustl(dir(6:ilen))),'.PFL'
              else
                  write(sfcfil,'(13(a))')trim(adjustl(metdir)),
     +            trim(adjustl(metdom)),'/MMIF_',trim(adjustl(metdom)),
     +            '_',trim(adjustl(acol)),'_',trim(adjustl(arow)),
     +            '/MMIF_AERMET_2018_',trim(adjustl(acol)),'_',
     +            trim(adjustl(arow)),'.SFC'
                  write(pflfil,'(13(a))')trim(adjustl(metdir)),
     +            trim(adjustl(metdom)),'/MMIF_',trim(adjustl(metdom)),
     +            '_',trim(adjustl(acol)),'_',trim(adjustl(arow)),
     +            '/MMIF_AERMET_2018_',trim(adjustl(acol)),'_',
     +            trim(adjustl(arow)),'.PFL'
              endif
               profbase=elev
          endif
          
          if (.not. noncontr .and. .not. akhiprvi) then
              read(wrfunit,*,iostat=eof)dir,icol,irow,elev
          else
              read(wrfunit,*,iostat=eof)icol,irow,elev
          endif
          goto 10
      endif
      
c     if a file not found, most likely outside domain, find closest one
      if (.not. lfound) then
          write(ilogunit,*)'Getting closest WRF cell'
          mindist=1000.0
          rewind(wrfunit)
          eof=0 
          if (.not. noncontr .and. .not. akhiprvi) then
              read(wrfunit,*,iostat=eof)dir,icol,irow,elev
          else
              read(wrfunit,*,iostat=eof)icol,irow,elev
          endif
20        if (eof .eq. 0) then
              
              metdist=sqrt(real((wrfcol-icol)*(wrfcol-icol)+
     +        (wrfrow-irow)*(wrfrow-irow)))
              if (metdist <= mindist) then  !match
                  lfound=.true.
                  mindist=metdist
                  write(acol,'(i3)')icol
                  write(arow,'(i3)')irow
                  if (.not. noncontr .and. .not. akhiprvi) then !make filenames
                      ilen=len_trim(dir)
                      write(sfcfil,'(5(a))')trim(adjustl(metdir)),
     +                trim(adjustl(metdom)),'/MMIF_AERMET_2018_',
     +                trim(adjustl(dir(6:ilen))),'.SFC'
                      write(pflfil,'(5(a))')trim(adjustl(metdir)),
     +                trim(adjustl(metdom)),'/MMIF_AERMET_2018_',
     +                trim(adjustl(dir(6:ilen))),'.PFL'
                  else
                      write(sfcfil,'(13(a))')trim(adjustl(metdir)),
     +                trim(adjustl(metdom)),'/MMIF_',
     +                trim(adjustl(metdom)),'_',trim(adjustl(acol)),'_',
     +                trim(adjustl(arow)),'/MMIF_AERMET_2018_',
     +                trim(adjustl(acol)),'_',trim(adjustl(arow)),'.SFC'
                      write(pflfil,'(13(a))')trim(adjustl(metdir)),
     +                trim(adjustl(metdom)),'/MMIF_',
     +                trim(adjustl(metdom)),'_',trim(adjustl(acol)),'_',
     +                trim(adjustl(arow)),'/MMIF_AERMET_2018_',
     +                trim(adjustl(acol)),'_',trim(adjustl(arow)),'.PFL'
                  endif
                  profbase=elev
              endif
          
              if (.not. noncontr .and. .not. akhiprvi) then
                  read(wrfunit,*,iostat=eof)dir,icol,irow,elev
              else
                  read(wrfunit,*,iostat=eof)icol,irow,elev
              endif
              goto 20
          endif    
      endif
      close(wrfunit)
          
c     make sure files exist
      call checkmet(sfcfil,pflfil)    
c     write ME pathway
      write(aerunit,15)trim(adjustl(sfcfil)),trim(adjustl(pflfil)),
     +swban,modyear,wrfcol,wrfrow,uwban,modyear,wrfcol,wrfrow,profbase
      
5     format(a,'does not exist, stop')
15    format('ME STARTING',/'ME SURFFILE ',a,/'ME PROFFILE ',a/
     +'ME SURFDATA ',i5,1x,i4,2(1x,i3)/'ME UAIRDATA ',i5,1x,i4,
     +2(1x,i3)/'ME PROFBASE ',f8.2/'ME FINISHED')
      return
      end
      
              
              
              
              
              

      
          
              
      
      
      
      
      
      
      
