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