module main1 c character variables c runtyp: AERMOD runtype see subroutine readargs for values c facid: facility ID c state2: 2-character state abbreviation c srcnam: source name used for AERMOD.INP file title c ptloc: directory containing AERMOD include location files c ptparam: directory containing AERMOD include source parameter files c ptcmaq: c pttemp: directory containing AERMOD include temporal factor files c or hourly emission files c metdir: directory containing AERMOD meteorology files c recdir: directory containing receptor files. these are not c in AERMOD format c outdir: directory for AERMOD output c locsrc: array of source IDs from location files c parsc: array of source IDs from parameter files c temsrc: array of source IDs from temporal files c locfil: facility specific AERMOD source locations file c paramfil: facility specific AERMOD source parameter file(s) c variable is 2D. 1 dim is point sources filename c 2nd dim is area sources filename c if no point or no area sources, then paramfil c is set to NA c temsrc: facility specific AERMOD temporal factor file or c hourly emissions file c singlhap: specific hap to process c gridst: grids character runtyp*7,facid*15,state2*2,ptloc*250,ptparam*250, +pttemp*250,ptcmaq*250,metdir*250,recdir*250,outdir*250, +hapdir*250,haps*15,xdir*250,rtype*15,cats*100,singlhap*15, +haplist*15,hap1*15,fac1*15,st2*2,singlfac*15,singlst*2,xfile*250, +states(54)*2,states1(54)*2,versn*5,logfil*50,ncst*2,gridst(4)*4!,grecid*15, c +monid*15,blockid*15,mfips*5 c allocatable :: locsrc(:) c allocatable :: parsrc(:) c allocatable :: temsrc(:) allocatable :: haps(:) allocatable :: haplist(:) allocatable :: facid(:) allocatable :: state2(:) allocatable :: rtype(:) allocatable :: cats(:) allocatable :: hap1(:) allocatable :: fac1(:) allocatable :: st2(:) allocatable :: outdir(:) allocatable :: xfile(:) c allocatable :: grecid(:) c allocatable :: monid(:) c allocatable :: blockid(:) c allocatable :: mfips(:) c lparam: 2-D array denoting parameter files found c 1 dim is point parameter file c 2 dim is area parameter file for point sources c ltemp: 3-D array denoting temporal files found c 1 dim is emissions factor file c 2 dim is hourly emissions file c 3 dim is both have been found c lcnty: logical variable denoting county for grid cell has been c found in RWC CMAQ to county cross reference c runpt: point source being processed c runap: airport source being processed c runon: onroad sources being processed c runnon: nonroad sources being processed c runrwc: RWC sources being processed c runcmv: CMV (port or underway) sources being processed c runnphi: nonpoint 10 m sources being processed c runnplow: low-level nonpoint sources being processed c akhiprvi: variable denoting that AK, HI, PR, or VI being processed c lgrids: denotes which grid files to read logical runpt,runap,runon,runnon,runrwc,runcmv,runnphi,runnplow, +runog,runag,akhiprvi,lsinghap,writefac,noncontr,conus(2),lsingfac, +lsingst,writetot,lgrids(5),readzero,writheid !write file for HEID purposes logical, allocatable, dimension(:) :: lhaps !when processing a set number of HAPs, to see if HAP is in rungroup c d2r: degrees to radians conversion factor c blkdist: distance to search from source for block receptors c grdddist: distance to search from source for gridded receptors c profbase: profbase written to AERMOD.INP real*8 d2r,blkdist,grddist,profbase real*8, allocatable, dimension(:,:,:) :: gridconc !# of categories, # of receptors real*8, allocatable, dimension(:,:,:) :: monconc !# of categories, # of receptors real*8, allocatable, dimension(:,:,:) :: blokconc !# of categories, # of receptors real*8, allocatable, dimension(:,:,:) :: zblokconc !# of categories, # of receptors for nonpopulated blocks in non-CONUS version 17236 real*8, allocatable, dimension(:,:,:) :: modgc !modeled gridded concentrations real*8, allocatable, dimension(:,:,:) :: modbc !modeled block concentrations real*8, allocatable, dimension(:,:,:) :: modzbc !modeled nonpopulated non-CONUS block concentrations version 17236 real*8, allocatable, dimension(:,:,:) :: modmc !modeled monitor concentrations real*8, allocatable, dimension(:,:,:) :: intbc !modeled block concentrations real*8, allocatable, dimension(:,:,:) :: intmc !modeled monitor concentrations real*8, allocatable, dimension(:,:) :: polyx !cmvx real*8, allocatable, dimension(:,:) :: polyy !cmvy real*8, allocatable, dimension(:,:) :: memis !seasonal emissions by group c real*8, allocatable, dimension(:) :: gridx !cmaq grid coordinates c real*8, allocatable, dimension(:) :: gridy !cmaq grid coordinates c real*8, allocatable, dimension(:) :: monx !cmaq grid coordinates c real*8, allocatable, dimension(:) :: mony !cmaq grid coordinates c real*8, allocatable, dimension(:) :: blockx !cmaq grid coordinates c real*8, allocatable, dimension(:) :: blocky !cmaq grid coordinates c col: CMAQ column c row: CMAQ row c inpunit: input_directories.txt file unit, lists directories (11) c locunit: locations file unit (12) c parunit: parameter file unit (13) c temunit: temporal file unit (14) c recunit: receptor file unit (16) c aerunit: AERMOD.INP file unit (27) c metunit: file of AK and HI met stations (32) c ngrps: # of source groups, usually equal to number of sources c utmzone: UTM zone of source c hrunit: hourly file unit (28) c gridres: grid resolution in km c cntyunit: county file unit (29) c fips: County FIPS code c wrfunit: WRF elevations file unit (31) c days: array of number of days/month c ncats: number of categories output by source group integer*8 col,row,inpunit,recunit,days(12),modyear, +ngrps,utmzone,gridres,fips,nlines,nhaps,nsrc,ngrids,nhours(12), +nblocks,nmons,nmodg,nmodb,nmodm,nmodz,nzblocks,nintm, +nintb,pltunit,ncats,nmonths,ngrps1,nzblocks1,seashrs(4,24), +hapunit,intunit,ngrids1,nblocks1,nmons1,noutdir,nxfile, +nruns,runno,nhap1,nhap2,ncolrow,istate,errunit,logunit integer*8, allocatable, dimension(:) :: srccr !array of source column/rows used for subsetting receptor arrays integer*8, allocatable, dimension(:) :: srccr1 !array of source column/rows used for subsetting receptor arrays,includes duplicates integer*8, allocatable, dimension(:) :: stgrid !array of grid identifiers for srccr integer*8, allocatable, dimension(:) :: stgrid1 !array of grid identifiers for srccr1 integer*8, allocatable, dimension(:) :: nsrcs !number of sources at facility c integer, allocatable, dimension(:) :: blockcr c integer, allocatable, dimension(:) :: moncr c AERMOD point source inputs type facinfo1 character srcid*12,srctyp*8 real*8 utmx,utmy,xdim,ydim,angle end type facinfo1 c AERMOD point source emission info type facinfo c character cat*100,facid*12,cas*15,state*2,srcgrp*8 character cat*100,facid*15,cas*15,state*2,srcgrp*12 integer*8 factype,procid,unitid,relid !process id, unit id, release id real*8 emis(1) end type facinfo c AERMOD CMV source inputs type cmvinfo1 character srcid*12,srctyp*8 real*8 utmx1,utmy1 integer*8 nvert end type cmvinfo1 c AERMOD CMV source emission info type cmvinfo character cat*100,facid*15,cas*15,state*2,srcgrp*12 real*8 emis(1) end type cmvinfo c AERMOD tract source inputs type tractinfo1 character srcid*12,srctyp*8 real*8 utmx1,utmy1 integer*8 nvert end type tractinfo1 c AERMOD tract source emission info type tractinfo character cat*100,facid*15,cas*15,state*2,srcgrp*12 real*8 emis(13) end type tractinfo c AERMOD airport inputs type apinfo1 character srcid*12,srctyp*8 real*8 area,utmx1,utmx2,utmy1,utmy2,frac,xdim,ydim,angle end type apinfo1 c airport facility emissions info type apinfo c character facid*12,cas*15,state*2,cat*100,srcgrp*8 character facid*15,cas*15,state*2,cat*100,srcgrp*12 real*8 emis(1) end type apinfo type gridinfo1 character srcid*12 real*8 utmx1,utmx2,utmx3,utmx4,utmy1,utmy2,utmy3,utmy4 integer*8 nvert end type gridinfo1 c onroad, nonroad type gridinfo c character facid*12,cas*15,fips*5,cat*100,srcgrp*8 character facid*15,cas*15,fips*5,cat*100,srcgrp*12 real*8 emis(13) !monthly emissions and annual. annual may not be used end type gridinfo c all gridded receptors type gridrecs character recid*15,rtype*15 real*8 x,y !x,y are cmaqx,cmaqy integer*8 colrow,iflag,col,row,irec,irec1 end type gridrecs type gridrecs1 character recid*15 integer*8 irec end type gridrecs1 c modeled gridded receptors type modgrids character recid*15 real*8 ux,uy,x,y,dist !ux,uy are utm, x,y are cmaqx,cmaqy integer*8 colrow,iflag,mastrec !mastrec is the index of the receptor in national gridded network end type modgrids c all block receptors type blkrecs character recid*15,rtype*15 real*8 x,y !x,y are cmaqx,cmaqy integer*8 colrow,iflag,col,row,irec,irec1 end type blkrecs type blkrecs1 character recid*15 integer*8 irec end type blkrecs1 c modeled block receptors type modblks character recid*15 real*8 ux,uy,x,y,dist !ux,uy are utm, x,y are cmaqx,cmaqy integer*8 colrow,iflag,mastrec !mastrec is the index of the receptor in national block network end type modblks c all nonpopulated non-CONUS block receptors type zblkrecs character recid*15,rtype*15 real*8 x,y !x,y are cmaqx,cmaqy integer*8 colrow,iflag,col,row,irec,irec1 end type zblkrecs type zblkrecs1 character recid*15 integer*8 irec end type zblkrecs1 c modeled nonpopulated non-CONUSblock receptors type zmodblks character recid*15 real*8 ux,uy,x,y,dist !ux,uy are utm, x,y are cmaqx,cmaqy integer*8 colrow,iflag,mastrec !mastrec is the index of the receptor in national block network end type zmodblks c interpolated blocks type intblks character recid*15 real*8 ux,uy,x,y,dist !ux,uy are utm, x,y are cmaqx,cmaqy integer*8 colrow,mastrec !mastrec is the index of the receptor in national block network end type intblks c all monitors type monrecs character recid*15,fip*5,rtype*5 real*8 x,y !x,y are cmaqx,cmaqy integer*8 colrow,iflag,col,row,irec,irec1 end type monrecs type monrecs1 character recid*15 integer*8 irec end type monrecs1 c modeled monitor receptors type modmons character recid*15 real*8 ux,uy,x,y,dist !ux,uy are utm, x,y are cmaqx,cmaqy integer*8 colrow,iflag,mastrec !mastrec is the index of the receptor in national monitor network end type modmons c interpolated monitors type intmons character recid*15 real*8 ux,uy,x,y,conc,dist !ux,uy are utm, x,y are cmaqx,cmaqy integer*8 colrow,mastrec !mastrec is the index of the receptor in national monitor network end type intmons type(facinfo1), dimension(:),allocatable:: ptsrcs type(facinfo), dimension(:),allocatable:: facility type(cmvinfo1), dimension(:),allocatable:: cmvsrcs type(cmvinfo), dimension(:),allocatable:: cmv type(tractinfo1), dimension(:),allocatable:: tractsrcs type(tractinfo), dimension(:),allocatable:: tracts type(apinfo), dimension(:),allocatable:: airport type(gridinfo), dimension(:),allocatable:: gridded c type(mobinfo), dimension(:),allocatable:: mobile type(apinfo1), dimension(:),allocatable:: apsrcs type(gridinfo1), dimension(:),allocatable:: gridsrcs type(gridrecs), dimension(:),allocatable:: gridrec type(gridrecs1), dimension(:),allocatable:: gridrec1 type(modgrids), dimension(:),allocatable:: modgrid type(blkrecs), dimension(:),allocatable:: blockrec type(blkrecs1), dimension(:),allocatable:: blockrec1 type(modblks), dimension(:),allocatable:: modblock type(intblks), dimension(:),allocatable:: intblock type(zblkrecs), dimension(:),allocatable:: zblockrec type(zblkrecs1), dimension(:),allocatable:: zblockrec1 type(zmodblks), dimension(:),allocatable:: zmodblock type(monrecs), dimension(:),allocatable:: monrec type(monrecs1), dimension(:),allocatable:: monrec1 type(modmons), dimension(:),allocatable:: modmon type(intmons), dimension(:),allocatable:: intmon c parameters parameter(inpunit=11,recunit=16,xwalk=38,pltunit=39, +hapunit=41,intunit=26,errunit=27,d2r=3.141592654/180.0, +grddist=50000.0,blkdist=10000.0,logunit=37) c assign number of days per month. February set to 28 c J F M A M J J A S O N D data days /31,28,31,30,31,30,31,31,30,31,30,31/ data versn /'20016'/ data states /'AL','AK','AZ','AR','CA','CO','CT','DE','DC','FL', +'GA','HI','ID','IL','IN','IA','KS','KY','LA','ME','MD','MA','MI', +'MN','MS','MO','MT','NE','NV','NH','NJ','NM','NY','NC','ND','OH', +'OK','OR','PA','RI','SC','SD','TN','TX','UT','VT','VA','WA','WV', +'WI','WY','PR','VI','TB'/ data states1 /'01','02','04','05','06','08','09','10','11','12', +'13','15','16','17','18','19','20','21','22','23','24','25','26', +'27','28','29','30','31','32','33','34','35','36','37','38','39', +'40','41','42','44','45','46','47','48','59','50','51','53','54', +'55','56','72','78','88'/ data gridst /'12US','09AK','03HI','03PR'/ end module main1 c########################################################################## program postaermod use main1 implicit none c read command line arguments to get run type call readargs c get directories call getdir c get list of HAPs to process call gethaps c read crosswalk lgrids=.false. if (runap) then call read_air elseif (runpt) then call read_point elseif (runcmv) then call read_cmv c elseif (noncontr:q:q) then c call read_tract else call read_grid endif c determine if source array has mix of CONUS and non-CONUS sources and get column/rows of sources call src_cr c get grid cells that are covered by c sources being processed to help c cut down number of receptors in master array c call haps_sources !now called in read_ subroutines c create matrix of all blocks, gridded and monitors call rec_matrix c read in modeled sources and receptors call read_model c write output if (writetot) call writeout close(logunit) end c########################################################################## c subroutine to read command line arguments subroutine readargs use main1 implicit none c i: loop counter c ilen: length of argument c agrid: character string of grid resolution (km) c awrite: character string to denonte write facility specific output (yes or no) integer*8 i,ilen,istat,ndash,n1,n2,get_state character agrid*2,awrite*3,aruns*10,arun*10,noncon*1,ahap1*3, +ahap2*3,runtyp1*12,awrite1*3,awrite2*3 c write(*,'(2(a))')'Starting program version: ',versn c call datetime c initialize logical flags for run types runpt=.false. runap=.false. runcmv=.false. runon=.false. runnon=.false. runrwc=.false. runnphi=.false. runnplow=.false. runog=.false. runag=.false. c initialize variable denoting if AK, HI, PR, or VI c being processed akhiprvi=.false. noncontr=.false. c get the requested run type c call getarg(1,runtyp) call getarg(1,runtyp1) c these are subject to change c convert to upper case c ilen=len_trim(runtyp) ilen=len_trim(runtyp1) ndash=0 do i=1,ilen c call upcase(runtyp(i:i)) call upcase(runtyp1(i:i)) if (runtyp1(i:i) .eq. '_') ndash=ndash+1 enddo if (ndash .eq. 0) then runtyp=trim(adjustl(runtyp1)) elseif (ndash .eq. 1) then n1=index(runtyp1,'_') runtyp=trim(adjustl(runtyp1(1:n1-1))) c read(runtyp1(n1+1:ilen),*)gridres read(runtyp1(n1+1:ilen),*)agrid call checkgrid(agrid) if (gridres .ne. 4 .and. gridres .ne. 12) then write(*,5)agrid stop endif elseif (ndash .eq. 2) then n1=index(runtyp1,'_') n2=index(runtyp1(n1+1:ilen),'_')+n1 runtyp=trim(adjustl(runtyp1(1:n1-1))) c read(runtyp1(n1+1:n2-1),*)gridres read(runtyp1(n1+1:n2-1),*)agrid call checkgrid(agrid) if (gridres .ne. 4 .and. gridres .ne. 12 .and. gridres + .ne. 3 .and. gridres .ne. 9) then write(*,5)agrid stop endif read(runtyp1(n2+1:ilen),'(a)')noncon ilen=len_trim(noncon) do i=1,ilen call upcase(noncon(i:i)) enddo c allow first letter of state to be equivalent of yes if (noncon .eq. 'Y' .or. noncon .eq. 'A' .or. noncon .eq. 'H' + .or .noncon .eq. 'P' .or. noncon .eq. 'V') then noncontr=.true. if (noncon .eq. 'A') then ncst='AK' elseif (noncon .eq. 'H') then ncst='HI' elseif (noncon .eq. 'P') then ncst='PR' elseif (noncon .eq. 'VI') then ncst='VI' else ncst='NA' endif elseif (noncon .eq. 'N') then noncontr=.false. else write(*,5)noncon stop endif else write(*,5)runtyp1 stop endif if (runtyp .eq. 'PT') then runpt=.true. nmonths=1 elseif (runtyp .eq. 'AP') then runap=.true. nmonths=1 elseif (runtyp .eq. 'PORT' .or. runtyp .eq. 'UNDER' .or. runtyp +.eq. 'CMV') then runcmv=.true. nmonths=1 elseif (runtyp .eq. 'HDON' .or. runtyp .eq. 'HDOFF' .or. + runtyp .eq. 'LDON' .or. runtyp .eq. 'LDOFF' .or. + runtyp .eq. 'ONOFF' .or. runtyp .eq. 'HOTEL') then runon=.true. nmonths=12 elseif (runtyp .eq. 'NONRD') then runnon=.true. nmonths=12 c if (runnon .and. noncontr) nmonths=1 !reset nonroad to 1 season for tract sources elseif (runtyp .eq. 'NPHI') then runnphi=.true. nmonths=1 elseif (runtyp .eq. 'NPLO') then runnplow=.true. nmonths=1 elseif (runtyp .eq. 'RWC') then runrwc=.true. nmonths=1 elseif (runtyp .eq. 'OILGAS') then runog=.true. nmonths=1 elseif (runtyp .eq. 'AG') then runag=.true. nmonths=1 else write(*,5)runtyp stop endif c single or all haps call getarg(2,singlhap,istat) ilen=len_trim(singlhap) do i=1,ilen call upcase(singlhap(i:i)) enddo if (singlhap .eq. 'ALL') then lsinghap=.false. else lsinghap=.true. if (trim(adjustl(singlhap)) .eq. 'BENZN') singlhap='BENZENE' !make consistent across run groups endif c HAP #s call getarg(3,ahap1,istat) if (lsinghap) then nhap1=-1 nhap2=-1 else read(ahap1,*)nhap1 call getarg(4,ahap2) read(ahap2,*)nhap2 endif c single facility or state call getarg(5,singlfac,istat) if (istat .eq. -1) then lsingfac=.false. lsingst=.false. else ilen=len_trim(singlfac) do i=1,ilen call upcase(singlfac(i:i)) enddo if (ilen .eq. 2) then !state lsingst=.true. singlst=trim(adjustl(singlfac)) istate=get_state(singlst) !get index of state lsingfac=.false. else lsingst=.false. if (singlfac .eq. 'ALL') then lsingfac=.false. else lsingfac=.true. endif endif endif c write facility level output call getarg(6,awrite,istat) if (istat .eq. -1) then writefac=.false. else ilen=len_trim(awrite) do i=1,ilen call upcase(awrite(i:i)) enddo if (awrite .eq. 'YES') then writefac=.true. else writefac=.false. endif endif c # of runs call getarg(7,aruns,istat) if (istat .eq. -1 .or. lsingfac) then nruns=1 else read(aruns,*)nruns endif c run number call getarg(8,arun,istat) if (istat .eq. -1 .or. lsingfac) then runno=1 else read(arun,*)runno endif c write total output across all processed facilities c may not want to do when just one facility call getarg(9,awrite1,istat) if (istat .eq. -1) then writetot=.false. else ilen=len_trim(awrite1) do i=1,ilen call upcase(awrite1(i:i)) enddo if (awrite1 .eq. 'YES') then writetot=.true. else writetot=.false. endif endif c write facility level data in format for HEID purpsoes call getarg(10,awrite2,istat) if (istat .eq. -1) then writheid=.false. else ilen=len_trim(awrite2) do i=1,ilen call upcase(awrite2(i:i)) enddo if (awrite2 .eq. 'YES') then writheid=.true. else writheid=.false. endif endif if (runno .gt. nruns) then write(*,10)runno,nruns stop endif if ((nhap1 .gt. 0 .and. nhap1 .gt. nhap2) .or. nhap1 .eq. 0) then write(*,*)'Number of HAP range invalid' stop endif if (lsingst .and. istate .eq. 0) then write(*,15)trim(adjustl(singlst)) stop endif if (.not. lsingst) istate=1 !set to 1 for later logic c write filename for logfile based on arguments if (lsingfac) then if (lsinghap) then write(logfil,'(4(a))')trim(adjustl(singlfac)),'_', + trim(adjustl(singlhap)),'.log' else write(logfil,'(2(a),i3.3,a,i3.3,a)') + trim(adjustl(singlfac)),'_',nhap1,'_',nhap2,'.log' endif elseif (lsingst) then if (lsinghap) then write(logfil,'(4(a))')trim(adjustl(singlst)),'_', + trim(adjustl(singlhap)),'.log' else write(logfil,'(2(a),i3.3,a,i3.3,a)') + trim(adjustl(singlst)),'_',nhap1,'_',nhap2,'.log' endif else if (lsinghap) then write(logfil,'(4(a))')trim(adjustl(runtyp)),'_', + trim(adjustl(singlhap)),'.log' else write(logfil,'(2(a),i3.3,a,i3.3,a)') + trim(adjustl(runtyp)),'_',nhap1,'_',nhap2,'.log' endif endif open(unit=logunit,file=logfil,status='replace') write(logunit,'(2(a))')'Starting program version: ',versn call datetime 5 format(a,'is not valid, stopping program') 10 format('Run number: ',i10/'is greater than number of runs: ',i10/ +'stopping program') 15 format('State identifier ',a,' is not valid, stopping program') return end c########################################################################## c subroutine to get pathnames from input_directories.dat c edit on sol once moved over subroutine getdir use main1 implicit none c eof: end of file indicator c n: index variable c i: loop counter integer*8 eof,n,nout,nx c line: line of text read from input_directories.dat c tag: source parameter file tag, point or area c or airport runway or non-runway c or port or underway for CMV c fname: source parameter filename written to paramfil array c acol: character string of CMAQ column c arow: character string of CMAQ row c ares: character string of grid resolution in km character line*300,ares*2,infil*250,outfil(2)*250,xfile1(4)*250 c lloc: denotes if locations folder keyword found or locations file exists c lparam: 2D array denotes if parameter file(s) exist c 1D = point, 2D=area c lcmaq: denotes if facility grid-cell xref folder keyword found c lmet: denotes if met directory keyword found c lrec: denotes if receptor directory keyword found c lout: denotes if output directory keyword found c lparam1: denotes if parameter folder keyword found c ltemp1: denotes if temporal folder keyword found c lout: denotes if output folder keyword found c lexist: denotes that input_directories.dat is in directory logical lloc,lparam2(2),lcmaq,lmet,lrec,lparam1, +ltemp1,lout,lexist,lxwalk,lhap,lxfile c initialize variables eof=0 lloc=.false. c lparam=.false. c ltemp=.false. lcmaq=.false. lmet=.false. lrec=.false. lparam1=.false. ltemp1=.false. lout=.false. c lcnty=.false. lxwalk=.false. lhap=.false. outfil='NA' nout=0 lxfile=.false. xfile1='NA' nx=0 nxfile=0 c paramfil='NA' c temfil='NA' c locfil='NA' c create facility ID for gridded sources based on column and row ! if (runon .or. runnon .or. runnphi .or. runnplow .or. runrwc) then ! write(acol,'(i3)')col ! write(arow,'(i3)')row ! write(facid,'(4a)')'G',trim(adjustl(acol)),'R', ! + trim(adjustl(arow)) ! endif c check to see if input_directories.dat is in directory c version 17037, now file has rungroup in name if (runpt .or. runcmv .or. runap) then write(infil,'(3(a))')'input_',trim(adjustl(runtyp)), + '_directories_orig.dat' else write(ares,'(i2)')gridres if (.not. noncontr) then write(infil,'(5(a))')'input_',trim(adjustl(runtyp)), + '_',trim(adjustl(ares)),'_directories.dat' else write(infil,'(7(a))')'input_',trim(adjustl(runtyp)), + '_',trim(adjustl(ares)),'_directories_',ncst,'.dat' endif endif inquire(file=infil,exist=lexist) if (lexist) then open(unit=inpunit,file=infil,status='old') else write(logunit,*)trim(adjustl(infil)),' does not exist' write(logunit,50) stop endif read(11,5,iostat=eof)line 10 if (eof .eq. 0) then if (index(line,'PTLOC ') .gt. 0) then n=index(line,'PTLOC ') ptloc=line(n+6:300) lloc=.true. endif if (index(line,'PTPARAM ') .gt. 0) then n=index(line,'PTPARAM ') ptparam=line(n+8:300) lparam1=.true. endif if (index(line,'PTTEMP ') .gt. 0) then n=index(line,'PTTEMP ') pttemp=line(n+7:300) ltemp1=.true. endif if (index(line,'PTCMAQ ') .gt. 0) then n=index(line,'PTCMAQ ') ptcmaq=line(n+7:300) lcmaq=.true. endif if (index(line,'METDIR ') .gt. 0) then n=index(line,'METDIR ') metdir=line(n+7:300) lmet=.true. endif if (index(line,'RECDIR ') .gt. 0) then n=index(line,'RECDIR ') recdir=line(n+7:300) lrec=.true. endif if (index(line,'XWALK ') .gt. 0) then !emissions directory n=index(line,'XWALK ') xdir=line(n+6:300) lxwalk=.true. endif if (index(line,'HAPDIR ') .gt. 0) then n=index(line,'HAPDIR ') hapdir=line(n+7:300) lhap=.true. endif if (index(line,'OUTPUT ') .gt. 0) then n=index(line,'OUTPUT ') nout=nout+1 outfil(nout)=line(n+7:300) lout=.true. endif if (index(line,'XFILE ') .gt. 0) then n=index(line,'XFILE ') nx=nx+1 xfile1(nx)=line(n+6:300) lxfile=.true. endif read(11,5,iostat=eof)line goto 10 endif nxfile=nx allocate(xfile(nxfile)) do n=1,nxfile write(xfile(n),'(2(a))')trim(adjustl(xdir)), + trim(adjustl(xfile1(n))) enddo noutdir=nout allocate(outdir(noutdir)) do n=1,noutdir outdir(n)=outfil(n) enddo c if any keyword not found, stop program if (.not. lloc) write(logunit,15) if (.not. lparam1) write(logunit,20) if (.not. ltemp1) write(logunit,25) if (.not. lcmaq) write(logunit,30) if (.not. lmet) write(logunit,35) if (.not. lrec) write(logunit,40) if (.not. lout) write(logunit,45) if (.not. lxwalk) write(logunit,46) if (.not. lhap) write(logunit,47) if (.not. lloc .or. .not. lparam1 .or. .not. ltemp1 .or. +.not. lcmaq .or. .not. lmet .or. .not. lrec .or. .not. lxwalk +.or. .not. lout .or. .not. lhap) then write(logunit,50) stop endif 5 format(a300) 15 format('PTLOC keyword not found') 20 format('PTPARAM keyword not found') 25 format('PTEMP keyword not found') 30 format('PTCMAQ keyword not found') 35 format('METDIR keyword not found') 40 format('RECDIR keyword not found') 45 format('OUTDIR keyword not found') 46 format('XWALK keyword not found') 47 format('HAPDIR keyword not found') 50 format('Aborting run..') 55 format('location file not found for') 60 format('parameter file(s) not found for') 65 format('temporal file(s) not found for') 70 format('County file not found for') return end c########################################################################## c subroutine to get haps subroutine gethaps use main1 implicit none integer eof,i,j character ahap*15,hapfile*250 logical lexist if (nhap1 .gt. 0) then write(hapfile,'(2(a))')trim(adjustl(ptcmaq)), + 'hap_priority_list.csv' inquire(file=hapfile,exist=lexist) if (lexist) then nhaps=nhap2-nhap1+1 allocate(haps(nhaps)) allocate(lhaps(nhaps)) open(unit=hapunit,file=hapfile,status='old') eof=0 i=0 j=0 read(hapunit,*,iostat=eof)ahap read(hapunit,*,iostat=eof)ahap 5 if (eof .eq. 0) then if (ahap .eq. 'BENZN') ahap='BENZENE' i=i+1 if (i .ge. nhap1 .and. i .le. nhap2) then j=j+1 haps(j)=trim(adjustl(ahap)) lhaps(j)=.false. endif read(hapunit,*,iostat=eof)ahap goto 5 endif close(hapunit) else write(logunit,*)'HAP list does not exist' write(logunit,*)'use single hap or all haps' endif else !single hap or all HAPs if (lsinghap) then nhaps=1 allocate(haps(nhaps)) allocate(lhaps(nhaps)) haps(1)=trim(adjustl(singlhap)) lhaps(1)=.false. else write(hapfile,'(2(a))')trim(adjustl(ptcmaq)), + 'hap_priority_list.csv' inquire(file=hapfile,exist=lexist) if (lexist) then open(unit=hapunit,file=hapfile,status='old') eof=0 i=0 read(hapunit,*,iostat=eof)ahap read(hapunit,*,iostat=eof)ahap 10 if (eof .eq. 0) then if (ahap .eq. 'BENZN') ahap='BENZENE' i=i+1 read(hapunit,*,iostat=eof)ahap goto 10 endif rewind(hapunit) nhaps=i allocate(haps(nhaps)) allocate(lhaps(nhaps)) eof=0 i=0 read(hapunit,*,iostat=eof)ahap read(hapunit,*,iostat=eof)ahap 15 if (eof .eq. 0) then if (ahap .eq. 'BENZN') ahap='BENZENE' i=i+1 haps(i)=trim(adjustl(ahap)) lhaps(i)=.false. read(hapunit,*,iostat=eof)ahap goto 15 endif close(hapunit) else write(logunit,*)'HAP list does not exist' endif endif endif return end c########################################################################## c subroutine to read airport pollutant/category/AERMOD source cross-walk c for airports subroutine read_air use main1 implicit none integer*8 eof,ilines,i,ilen,ihap,findhap,nlines1,findfac,is,ityp character facid2*15,relid*8,scc1*8,cas1*15,src*12, +junk*1,st2a*2,junk2*20 real*8 ann logical lexist,lgo write(logunit,*)'Reading emissions file' c check for file existence inquire(file=xfile(1),exist=lexist) if (.not. lexist) then write(logunit,5)trim(adjustl(xfile(1))) stop endif c read in cross walk and assign variables to appropriate arrays lgo=.false. open(unit=xwalk,file=xfile,status='old') ilines=0 read(xwalk,*)junk ncats=1 read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann 10 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then if ((lsingst .and. (trim(adjustl(st2a)) .eq. + trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) + .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. + trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. + (.not. lsingst .and. .not. lsingfac)) then ilines=ilines+1 lhaps(ihap)=.true. lgo=.true. endif endif read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann goto 10 endif if (.not. lgo) then close(xwalk) write(logunit,*)'Requested not HAPs found' stop endif rewind(xwalk) c nlines=ilines nlines1=ilines c allocate(hap1(nlines1)) allocate(st2(nlines1)) allocate(fac1(nlines1)) fac1='NA' c allocate(airport(nlines)) c airport%srcgrp='ALL' c airport%cat='NR-Point-Airports' eof=0 ilines=0 read(xwalk,*)junk read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann 15 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then if ((lsingst .and. (trim(adjustl(st2a)) .eq. + trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) + .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. + trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. + (.not. lsingst .and. .not. lsingfac)) then ilines=ilines+1 c hap1(ilines)=cas1 st2(ilines)=st2a fac1(ilines)=facid2 endif endif read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann goto 15 endif c now get unique list of haps/sources. sources based on run# call haps_sources(nlines1) rewind(xwalk) eof=0 ilines=0 read(xwalk,*)junk read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann 20 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0)then is=findfac(facid2) if (is .gt. 0) ilines=ilines+1 endif read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann goto 20 endif nlines=ilines allocate(airport(nlines)) airport%srcgrp='ALL' airport%cat='NR-Point-Airports' c now assign to the array based on HAPs and sources rewind(xwalk) ilines=0 eof=0 read(xwalk,*)junk read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann 25 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0 )then is=findfac(facid2) if (is .gt. 0) then ilines=ilines+1 airport(ilines)%state=st2a if (st2a .eq. '02') then lgrids(2)=.true. elseif (st2a .eq. '15') then lgrids(3)=.true. elseif (st2a .eq. '72' .or. st2a .eq. '78') then lgrids(4)=.true. else lgrids(1)=.true. endif airport(ilines)%facid=facid2 airport(ilines)%cas=cas1 airport(ilines)%emis=ann endif endif read(xwalk,*,iostat=eof)st2a,facid2,junk2,ityp,cas1,ann goto 25 endif c write(*,*)nlines close(xwalk) c do ilines=1,nlines c write(107,'(3(a,1x),e13.6)') c + airport(ilines)%state,airport(ilines)%facid, c + airport(ilines)%cas,airport(ilines)%emis c enddo 5 format(a,/' cross-walk file does not exist...stopping program') return end c########################################################################## c subroutine to read pollutant/category/AERMOD source cross-walk c for point subroutine read_point use main1 implicit none integer*8 eof,ilines,ftype,i,ilen,ihap,findhap,nlines1,findfac,is, + prid,uid,rid,n1,n2,nn,nn1 character facid2*15,relid*8,scc1*8,cas1*15,src*12, +junk*1,st2a*2,junk2*20,ares*2,eline*200 real*8 ann logical lexist,lgo,lfix write(logunit,*)'Reading emissions file' c check for file existence inquire(file=xfile(1),exist=lexist) if (.not. lexist) then write(logunit,5)trim(adjustl(xfile(1))) stop endif c write the cross walk to a temporary cross walk and put in quotes for the c facility name c close the existing cross walk and open the new one using the same file unit open(unit=xwalk,file=xfile,status='old') open(unit=105,file='temp_emis.txt',status='replace') eof=0 do while (eof .eq. 0) read(xwalk,'(a200)',iostat=eof)eline if (eof .eq. 0) then nn=1 n1=0 n2=0 nn1=0 lfix=.false. if (index(eline,'"') .eq. 0) then c find 2nd and 3rd commas do while (nn .le. 200 .and. .not. lfix) if (eline(nn:nn) .eq. ',') then nn1=nn1+1 if (nn1 .eq. 2)n1=nn if (nn1 .eq. 3)n2=nn endif if (n1 .gt. 0 .and. n2 .gt. 0)lfix=.true. nn=nn+1 enddo if (eline(n1+1:n1+1) .ne. '"')eline(n1+1:n1+1)='"' if (eline(n2-1:n2-1) .ne. '"') eline(n2-1:n2-1)='"' endif write(105,'(a200)')eline endif enddo close(xwalk) close(105) stop lgo=.false. c read in cross walk and assign variables to appropriate arrays open(unit=xwalk,file=xfile,status='old') c open(unit=xwalk,file='temp_emis.txt',status='old') ilines=0 read(xwalk,*)junk c set initial value of ftype to 999 before reading in c because some records are missing the value ftype=999 read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,rid, +cas1,ann 10 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then if ((lsingst .and. (trim(adjustl(st2a)) .eq. + trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) + .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. + trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. + (.not. lsingst .and. .not. lsingfac)) then ilines=ilines+1 lhaps(ihap)=.true. lgo=.true. endif endif ftype=999 read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid, + rid,cas1,ann goto 10 endif if (.not. lgo) then write(logunit,*)'Requested HAPs not found' stop endif rewind(xwalk) c nlines=ilines nlines1=ilines c allocate(hap1(nlines1)) allocate(st2(nlines1)) allocate(fac1(nlines1)) c allocate(facility(nlines)) eof=0 ilines=0 read(xwalk,*)junk ftype=999 read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,rid, +cas1,ann 15 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then if ((lsingst .and. (trim(adjustl(st2a)) .eq. + trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) + .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. + trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. + (.not. lsingst .and. .not. lsingfac)) then ilines=ilines+1 c hap1(ilines)=cas1 st2(ilines)=st2a fac1(ilines)=facid2 endif endif ftype=999 read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid, + rid,cas1,ann goto 15 endif c now get unique list of haps/sources. sources based on run# call haps_sources(nlines1) rewind(xwalk) eof=0 ilines=0 ftype=999 read(xwalk,*)junk read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,rid, +cas1,ann 20 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then is=findfac(facid2) if (is .gt. 0) ilines=ilines+1 endif ftype=999 read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid, + rid,cas1,ann goto 20 endif rewind(xwalk) eof=0 nlines=ilines allocate(facility(nlines)) ilines=0 read(xwalk,*)junk ftype=999 read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid,rid, +cas1,ann 25 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then is=findfac(facid2) if (is .gt. 0) then ilines=ilines+1 facility(ilines)%state=st2a if (st2a .eq. '02') then lgrids(2)=.true. elseif (st2a .eq. '15') then lgrids(3)=.true. elseif (st2a .eq. '72' .or. st2a .eq. '78') then lgrids(4)=.true. else lgrids(1)=.true. endif facility(ilines)%facid=facid2 facility(ilines)%cas=cas1 facility(ilines)%emis=ann facility(ilines)%factype=ftype facility(ilines)%srcgrp=src facility(ilines)%procid=prid facility(ilines)%unitid=uid facility(ilines)%relid=rid if (ftype .eq. 151) then facility(ilines)%cat='NR-Point_Railyards' else facility(ilines)%cat='Point' endif endif endif ftype=999 read(xwalk,*,iostat=eof)st2a,facid2,junk2,ftype,src,prid,uid, + rid,cas1,ann goto 25 endif close(xwalk) do ilines=1,nlines c write(107,'(4(a,1x),i3,1x,e13.6)') c + facility(ilines)%state,facility(ilines)%facid, c + facility(ilines)%srcgrp,facility(ilines)%cas, c + facility(ilines)%factype,facility(ilines)%emis enddo 5 format(a,/' cross-walk file does not exist...stopping program') return end c########################################################################## c subroutine to read pollutant/category/AERMOD source cross-walk c for CMV sources subroutine read_cmv use main1 implicit none integer*8 eof,ilines,i,ilen,ihap,findhap,is,findfac,nlines1,nn character facid2*15,st2a*2,rungrp*3,src*12,cat1*100,cas1*15, +tag*8,junk*1 real*8 ann logical lexist,lgo write(logunit,*)'Reading emissions file' c if (trim(adjustl(runtyp)) .eq. 'PORT') then c tag='ports' c else c tag='underway' c endif c write(xfile,'(3(a))')trim(adjustl(xdir)),trim(adjustl(tag)), c +'_emis.csv' c check for file existence inquire(file=xfile(1),exist=lexist) if (.not. lexist) then write(logunit,5)trim(adjustl(xfile(1))) stop endif lgo=.false. c read in cross walk and assign variables to appropriate arrays open(unit=xwalk,file=xfile,status='old') ilines=0 read(xwalk,*)junk read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann 10 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then if ((lsingst .and. (trim(adjustl(st2a)) .eq. + trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) + .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. + trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. + (.not. lsingst .and. .not. lsingfac)) then ilines=ilines+1 lhaps(ihap)=.true. lgo=.true. endif endif read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann goto 10 endif if (.not. lgo) then write(logunit,*)'Requested HAPs not found' stop endif rewind(xwalk) nlines1=ilines c nlines=ilines c allocate(cmv(nlines)) c allocate(hap1(nlines1)) allocate(st2(nlines1)) allocate(fac1(nlines1)) eof=0 ilines=0 read(xwalk,*)junk read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann 15 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then if ((lsingst .and. (trim(adjustl(st2a)) .eq. + trim(adjustl(states(istate))) .or. trim(adjustl(st2a)) + .eq. trim(adjustl(states1(istate))))) .or. (lsingfac .and. + trim(adjustl(facid2)) .eq. trim(adjustl(singlfac))) .or. + (.not. lsingst .and. .not. lsingfac)) then ilines=ilines+1 c hap1(ilines)=cas1 st2(ilines)=st2a fac1(ilines)=facid2 endif endif read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann goto 15 endif c now get unique list of haps/sources. sources based on run# call haps_sources(nlines1) rewind(xwalk) eof=0 ilines=0 read(xwalk,*)junk read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann 20 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then is=findfac(facid2) if (is .gt. 0) ilines=ilines+1 endif read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann goto 20 endif rewind(xwalk) eof=0 nlines=ilines allocate(cmv(nlines)) ilines=0 read(xwalk,*)junk read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann 25 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then is=findfac(facid2) if (is .gt. 0) then ilines=ilines+1 cmv(ilines)%state=st2a if (st2a .eq. '02') then lgrids(2)=.true. elseif (st2a .eq. '15') then lgrids(3)=.true. elseif (st2a .eq. '72' .or. st2a .eq. '78') then lgrids(4)=.true. else lgrids(1)=.true. endif cmv(ilines)%facid=facid2 cmv(ilines)%cas=cas1 cmv(ilines)%emis=ann nn=index(src,'c') write(cmv(ilines)%srcgrp,'(2(a1))')'C',src(nn+1:nn+1) c cmv(ilines)%srcgrp=src if (runtyp .eq. 'PORT' .or. runtyp .eq. 'UNDER') then write(cmv(ilines)%cat,'(2(a))')'NR-CMV_', + trim(adjustl(cat1)) else cmv(ilines)%cat=cat1 endif endif endif read(xwalk,*,iostat=eof)st2a,rungrp,facid2,src,cat1,cas1,ann goto 25 endif close(xwalk) do ilines=1,nlines write(107,'(4(a,1x),1x,e13.6)') + cmv(ilines)%state,cmv(ilines)%facid, + cmv(ilines)%srcgrp,cmv(ilines)%cas, + cmv(ilines)%emis enddo 5 format(a,/' cross-walk file does not exist...stopping program') return end c########################################################################## c subroutine to read pollutant/category/AERMOD source cross-walk c for tract sources subroutine read_tract use main1 implicit none integer*8 eof,ilines,i,ilen,is,ihap,findhap,findfac,nlines1,is1,m character st2a*2,rungrp*3,src*12,cat1*100,cas1*15, +junk*10,stfips1(4)*2,ares*2,fips6*6,poly*4, +tract1*11,facid2*16 real*8 ann,memis1(12)!,win,spr,sum,fall logical lexist(4),lgo write(logunit,*)'Reading emissions file' stfips1(1)='02' stfips1(2)='15' stfips1(3)='72' stfips1(4)='78' write(ares,'(i2)')gridres lgo=.false. do is=1,4 c check for file existence inquire(file=xfile(is),exist=lexist(is)) enddo c check to see if any files exist if (.not. lexist(1) .and. .not. lexist(2) .and. .not. lexist(3) +.and. .not. lexist(4)) then write(logunit,*)' no x-walk files exist...stopping program' stop endif ilines=0 do is=1,4 if (lexist(is)) then c read in cross walk and assign variables to appropriate arrays open(unit=xwalk,file=xfile(is),status='old') read(xwalk,*)junk read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly, + facid2,src,cat1,cas1,ann 10 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then if ((lsingst .and. (trim(adjustl(st2a)) .eq. + trim(adjustl(states(istate))) .or. + trim(adjustl(st2a)) .eq. + trim(adjustl(states1(istate))))) .or. (lsingfac + .and. trim(adjustl(src)) .eq. + trim(adjustl(singlfac))) .or. (.not. lsingst .and. + .not. lsingfac)) then ilines=ilines+1 lhaps(ihap)=.true. lgo=.true. endif endif read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly, + facid2,src,cat1,cas1,ann goto 10 endif close(xwalk) endif enddo if (.not. lgo) then write(logunit,*)'Requested HAPs not found' stop endif nlines1=ilines c allocate(hap1(nlines1)) allocate(st2(nlines1)) allocate(fac1(nlines1)) c nlines=ilines c allocate(tracts(nlines)) ilines=0 do is=1,4 if (lexist(is)) then open(unit=xwalk,file=xfile(is),status='old') eof=0 read(xwalk,*)junk read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly, + facid2,src,cat1,cas1,ann 15 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then if ((lsingst .and. (trim(adjustl(st2a)) .eq. + trim(adjustl(states(istate))) .or. + trim(adjustl(st2a)).eq. + trim(adjustl(states1(istate))))) .or. (lsingfac + .and. trim(adjustl(src)) .eq. + trim(adjustl(singlfac))).or. (.not. lsingst .and. + .not. lsingfac)) then ilines=ilines+1 c hap1(ilines)=cas1 st2(ilines)=st2a fac1(ilines)=src endif endif read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly, + facid2,src,cat1,cas1,ann goto 15 endif close(xwalk) endif enddo c now get unique list of haps/sources. sources based on run# call haps_sources(nlines1) ilines=0 do is=1,4 if (lexist(is)) then open(unit=xwalk,file=xfile(is),status='old') eof=0 read(xwalk,*)junk read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly, + facid2,src,cat1,cas1,ann 20 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then is1=findfac(src) if (is1 .gt. 0) ilines=ilines+1 endif read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly, + facid2,src,cat1,cas1,ann goto 20 endif close(xwalk) endif enddo nlines=ilines ilines=0 allocate(tracts(nlines)) do is=1,4 if (lexist(is)) then open(unit=xwalk,file=xfile(is),status='old') eof=0 read(xwalk,*)junk if (runon) then read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly, + facid2,src,cat1,cas1,ann,(memis1(m),m=1,12) !,win,spr,sum,fall else read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1,poly, + facid2,src,cat1,cas1,ann endif 25 if (eof .eq. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then is1=findfac(src) if (is1 .gt. 0) then ilines=ilines+1 tracts(ilines)%state=st2a tracts(ilines)%facid=src !for tract, the AERMOD source ID is the facility ID tracts(ilines)%cas=cas1 if (runon .or. runnon) then do m=1,12 tracts(ilines)%emis(m)=memis1(m) enddo endif tracts(ilines)%emis(13)=ann tracts(ilines)%srcgrp=src tracts(ilines)%cat=cat1 endif endif if (runon) then read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1, + poly,facid2,src,cat1,cas1,ann,(memis1(m),m=1,12) !win,spr,sum,fall else read(xwalk,*,iostat=eof)rungrp,st2a,fips6,tract1, + poly,facid2,src,cat1,cas1,ann endif goto 25 endif close(xwalk) endif enddo do ilines=1,nlines write(107,'(4(a,1x),1x,e13.6)') + tracts(ilines)%state,tracts(ilines)%facid, + tracts(ilines)%srcgrp,tracts(ilines)%cas, + tracts(ilines)%emis(5) enddo c do ilines=1,nlines c write(105,*)tracts(ilines)%cas c enddo 5 format(a,/' cross-walk file does not exist...stopping program') return end c########################################################################## c subroutine to readpollutant/category/AERMOD source cross-walk c for gridded sources subroutine read_grid use main1 implicit none integer*8 eof,ilines,i,ilen,ihap,findhap,isrc,findfac,ngsrc,nrows, +ncols,r,c,nsrc1,j,jj,k,singc,singr,m integer*8, allocatable, dimension(:,:) :: igrps character facid2a*15,relid*8,scc1*8,cas1*15,src*12, +junk*1,st2a*2,junk2*20,ares*2,afips*5,cat1*100,facid2*15,fid*15 real*8 ann,memis1(12) !win,spr,sum,fall logical lexist,lf,lgo allocatable :: facid2(:) write(logunit,*)'Reading emissions file' write(ares,'(i2)')gridres c check for file existence inquire(file=xfile(1),exist=lexist) if (.not. lexist) then write(logunit,5)trim(adjustl(xfile(1))) stop endif c if single state if (lsingst) then open(unit=xwalk,file=xfile(1),status='old') ilines=0 read(xwalk,*)junk read(xwalk,*,iostat=eof)junk2,afips,facid2a 20 if (eof .eq. 0) then st2a=afips(1:2) !may change for v2, may be 5 characters instead of 6 if (trim(adjustl(st2a)) .eq. trim(adjustl(states(istate))) + .or. trim(adjustl(st2a)) .eq. + trim(adjustl(states1(istate)))) isrc=1 if (isrc .gt. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then read(facid2a,'(2(a1,i3))')junk,c,junk,r call gridsrc(r,c,fid,lf) if (lf) then ilines=ilines+1 lhaps(ihap)=.true. lgo=.true. endif endif endif read(xwalk,*,iostat=eof)junk2,afips,facid2a goto 20 endif if (.not. lgo) then write(logunit,40)trim(adjustl(states(istate))) stop endif nsrc1=ilines allocate(facid2(nsrc1)) rewind(xwalk) read(xwalk,*)junk read(xwalk,*,iostat=eof)junk2,afips,facid2a 25 if (eof .eq. 0) then st2a=afips(1:2) if (trim(adjustl(st2a)) .eq. trim(adjustl(states(istate))) + .or. trim(adjustl(st2a)) .eq. + trim(adjustl(states1(istate)))) isrc=1 if (isrc .gt. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then read(facid2a,'(2(a1,i3))')junk,c,junk,r call gridsrc(r,c,fid,lf) if (lf) then ilines=ilines+1 facid2(i)=fid endif endif endif read(xwalk,*,iostat=eof)junk2,afips,facid2a goto 25 endif close(xwalk) endif c get column/row of the single grid cell from command line if (lsingfac) then read(singlfac,'(2(a1,i3))')junk,singc,junk,singr call gridsrc(singr,singc,fid,lf) if (lf) then nsrc1=1 allocate(facid2(nsrc1)) facid2(1)=singlfac else write(logunit,45)trim(adjustl(singlfac)) stop endif endif c not single state or grid cell if (.not. lsingfac .or. lsingst) then ngsrc=0 nrows=246 ncols=396 do r=1,nrows do c=1,ncols call gridsrc(r,c,fid,lf) if (lf)ngsrc=ngsrc+1 enddo enddo nsrc1=ngsrc ngsrc=0 allocate(facid2(nsrc1)) do r=1,nrows do c=1,ncols call gridsrc(r,c,fid,lf) if (lf) then ngsrc=ngsrc+1 facid2(ngsrc)=fid endif enddo enddo endif c common to scenarios, single state, single facility, or multiple facilities nsrc=nsrc1/nruns if (mod(real(nsrc1),real(nruns)) .ne. 0.0) nsrc=nsrc+1 allocate(igrps(nruns,2)) c write(*,*)'a ',nsrc1,nruns,runno,nsrc c pause !delete later do j=1,nruns if (j .eq. 1) then igrps(j,1)=1 else igrps(j,1)=igrps(j-1,2)+1 endif igrps(j,2)=j*nsrc enddo igrps(nruns,2)=nsrc1 jj=igrps(runno,2)-igrps(runno,1)+1 allocate(facid(nsrc)) facid='NA' k=0 do i=1,nsrc1 if (i .ge. igrps(runno,1) .and. i .le. + igrps(runno,2)) then k=k+1 facid(k)=facid2(i) endif enddo write(logunit,'(a,1x,i7)')'Total sources:',nsrc1 write(logunit,'(a,1x,i7)')'Sources to process:',nsrc c read in cross walk and assign variables to appropriate arrays open(unit=xwalk,file=xfile,status='old') ilines=0 read(xwalk,*)junk read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1,cas1,ann 10 if (eof .eq. 0) then isrc=findfac(facid2a) if (isrc .gt. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups ihap=findhap(cas1) if (ihap .gt. 0) then ilines=ilines+1 lhaps(ihap)=.true. lgo=.true. endif endif read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1,cas1,ann goto 10 endif if (.not. lgo) then write(logunit,*)'Requested HAPs not found' close(xwalk) stop endif rewind(xwalk) nlines=ilines call haps_sources(nlines) allocate(gridded(nlines)) eof=0 ilines=0 read(xwalk,*)junk if (runnplow .or. runnphi .or. runrwc .or. runog .or. runag) then read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1,cas1,ann else read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1,cas1, + ann,(memis1(m),m=1,12) !win,spr,sum,fall endif 15 if (eof .eq. 0) then isrc=findfac(facid2a) if (isrc .gt. 0) then ilen=len_trim(cas1) do i=1,ilen call upcase(cas1(i:i)) enddo if (trim(adjustl(cas1)) .eq. 'BENZN') cas1='BENZENE' !consistent across run groups call fixcat(cat1) ihap=findhap(cas1) if (ihap .gt. 0) then call fixcat(cat1) ilines=ilines+1 gridded(ilines)%fips=afips if (afips(1:2) .eq. '02') then lgrids(2)=.true. elseif (afips(1:2) .eq. '15') then lgrids(3)=.true. elseif (afips(1:2) .eq. '72' .or. st2a .eq. '78') then lgrids(4)=.true. else lgrids(1)=.true. endif gridded(ilines)%facid=facid2a gridded(ilines)%cas=cas1 gridded(ilines)%emis(13)=ann gridded(ilines)%srcgrp=src gridded(ilines)%cat=cat1 if (runon .or. runnon) then do m=1,12 gridded(ilines)%emis(m)=memis1(m) enddo endif endif endif if (runnplow .or. runnphi .or. runrwc .or. runog + .or. runag) then read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1, + cas1,ann else read(xwalk,*,iostat=eof)junk2,afips,facid2a,src,cat1, + cas1,ann,(memis1(m),m=1,12) !win,spr,sum,fall endif goto 15 endif close(xwalk) deallocate(facid2) deallocate(igrps) 5 format(a,/' cross-walk file does not exist...stopping program') 40 format(a,/'AERMOD output not found or no sources for state ',a) 45 format(a,/'AERMOD output not found for ',a, +' or not listed in cross-walk') return end c########################################################################## c subroutine to replace blanks and commas with underscores in c category names subroutine fixcat(cat1) use main1 implicit none integer ilen,i character cat1*100 ilen=len_trim(cat1) do i=1,ilen if (cat1(i:i) .eq. ' ' .or. cat1(i:i) .eq. ',') cat1(i:i)='_' enddo return end c########################################################################## c subroutine to get unique list of HAPs,categories, and sources to be used when writing c HAP specific output files subroutine haps_sources(nlines1) use main1 implicit none integer*8 i,ihap,n,ihap1,nsrc1,j,k,jj,ihap2,nlines1 integer*8, allocatable, dimension(:,:) :: igrps character h*15, st*2,rt*15,facid2*15,st2a*2 allocatable :: h(:) allocatable :: facid2(:) allocatable :: st2a(:) logical ldup,lreset write(logunit,*)'Getting HAPs and sources' c NR-C1C2_ports c NR-C3_ports c NR-C1C2C3_underway c NP-AgriculturalLivestock c determine if any requsted HAPS need to be removed from list for rungroup c if (nhap1 .gt. 0) then lreset=.false. ihap1=nhaps ihap2=nhaps do ihap=1,nhaps if (.not. lhaps(ihap)) then ihap1=ihap1-1 lreset=.true. write(logunit,'(2(a))')'Removing ', + trim(adjustl(haps(ihap))) endif enddo if (lreset) then allocate(h(nhaps)) h=haps deallocate(haps) nhaps=ihap1 allocate(haps(nhaps)) ihap1=0 do ihap=1,ihap2 if (lhaps(ihap)) then ihap1=ihap1+1 haps(ihap1)=h(ihap) endif enddo deallocate(h) endif c endif c temporary array of haps, sources, duplicates included ihap1=nlines1 write(logunit,'(a,1x,i3)')'Number of HAPs to process:',nhaps c do ihap=1,nhaps c write(102,*)haps(ihap) c enddo c write(*,*)'haps done' c eliminate duplicate facid and state if applicable c write(*,*)'start sources' if (runcmv .or. runap .or. runpt) then ! .or. noncontr) then ihap=1 do i=2,ihap1 ldup=.false. c n=1 if (runcmv .or. runap .or. runpt) then n=i-1 else n=1 endif do while (n .lt. i .and. .not. ldup) if (fac1(i) .eq. fac1(n)) ldup=.true. c write(*,*)fac1(i),fac1(n),i,n,ldup,ihap c pause n=n+1 enddo if (.not. ldup) ihap=ihap+1 enddo nsrc=ihap nsrc1=ihap c write(*,*)'eliminated duplicates' c allocate(facid(nsrc)) c if (runpt .or. runap .or. runcmv) allocate(state2(nsrc)) allocate(facid2(nsrc1)) if (allocated(st2)) allocate(st2a(nsrc1)) c facid(1)=fac1(1) c if (runpt .or. runap .or. runcmv) state2(1)=st2(1) facid2(1)=fac1(1) if (allocated(st2a)) st2a(1)=st2(1) ihap=1 do i=2,ihap1 ldup=.false. c n=1 if (runcmv .or. runap .or. runpt) then n=i-1 else n=1 endif do while (n .lt. i .and. .not. ldup) if (fac1(i) .eq. fac1(n)) ldup=.true. n=n+1 enddo if (.not. ldup) then ihap=ihap+1 c facid(ihap)=fac1(i) c if (runpt .or. runap .or. runcmv) state2(ihap)=st2(i) facid2(ihap)=fac1(i) if (allocated(st2a)) st2a(ihap)=st2(i) endif enddo c break up nsrc=nsrc1/nruns if (nsrc1 .lt. nruns) then write(logunit,5)nruns,nsrc1 stop endif if (mod(real(nsrc1),real(nruns)) .ne. 0.0) nsrc=nsrc+1 allocate(igrps(nruns,2)) c write(*,*)'a ',nsrc1,nruns,runno,nsrc do j=1,nruns if (j .eq. 1) then igrps(j,1)=1 else igrps(j,1)=igrps(j-1,2)+1 endif igrps(j,2)=j*nsrc enddo igrps(nruns,2)=nsrc1 c do j=1,nruns c write(*,*)'b ',j,igrps(j,1),igrps(j,2) c enddo jj=igrps(runno,2)-igrps(runno,1)+1 c write(*,*)jj,igrps(runno,2),igrps(runno,2) allocate(facid(nsrc)) facid='NA' if (allocated(st2a)) then allocate(state2(nsrc)) state2='NA' endif k=0 do i=1,nsrc1 if (i .ge. igrps(runno,1) .and. i .le. + igrps(runno,2)) then k=k+1 facid(k)=facid2(i) if (allocated(st2a))state2(k)=st2a(i) endif enddo write(logunit,'(a,1x,i7)')'Total sources:',nsrc1 c write(*,'(a,1x,i7)')'Sources to process:',nsrc endif c write(*,'(a,1x,i7)')'Total sources:',nsrc1 write(logunit,'(a,1x,i7)')'Sources to process:',nsrc c do ihap=1,nsrc c write(104,*)facid(ihap) c enddo c write(*,*)'sources done' c duplicates for categories if applicable c version 17030, now assign categories, don't rely c on using xwalk because some HAPs don't have all c categories and want consistency between HAPs if (runap) then allocate(cats(1)) cats='NR-Point-Airports' ncats=1 elseif (runpt) then allocate(cats(2)) cats(1)='Point' cats(2)='NR-Point_Railyards' ncats=2 elseif (runcmv) then if (runtyp .eq. 'PORT') then allocate(cats(1)) cats(1)='NR-CMV_ports' ncats=1 elseif (runtyp .eq. 'UNDER') then allocate(cats(1)) cats(1)='NR-CMV_underway' ncats=1 else allocate(cats(3)) cats(1)='NR-C1C2_ports' cats(2)='NR-C3_ports' cats(3)='NR-C1C2C3_underway' ncats=3 endif elseif (runnphi) then allocate(cats(3)) cats(1)='NP-FuelCombustion_not_RWC' cats(2)='NP-Industrial' cats(3)='NP-WasteDisposal' ncats=3 elseif (runnplow) then allocate(cats(5)) cats(1)='NP-CommercialCooking' c cats(2)='NP-FuelCombustion_not_RWC' cats(2)='NP-MiscellaneousNonindustrial' cats(3)='NP-SolventsCoatings' cats(4)='NP-StorageTransfer_BulkTerminals_GasStage1' cats(5)='NR-Locomotives' ncats=5 elseif (runnon) then allocate(cats(7)) cats(1)='NR-Agriculture' cats(2)='NR-CommercialEquipment' cats(3)='NR-CommercialLawnGarden' cats(4)='NR-Construction' cats(5)='NR-Recreational-inc-PleasureCraft' cats(6)='NR-ResidentialLawnGarden' cats(7)='NR-AllOther' ncats=7 elseif (runog) then allocate(cats(1)) cats(1)='NP-OilGas' ncats=1 elseif (runon) then if (runtyp .eq. 'HDOFF') then allocate(cats(2)) cats(1)='OR-HeavyDuty-OffNetwork-Diesel' cats(2)='OR-HeavyDuty-OffNetwork-Gas' ncats=2 elseif (runtyp .eq. 'LDOFF') then allocate(cats(2)) cats(1)='OR-LightDuty-OffNetwork-Diesel' cats(2)='OR-LightDuty-OffNetwork-Gas' ncats=2 elseif (runtyp .eq. 'HDON') then allocate(cats(2)) cats(1)='OR-HeavyDuty-OnNetwork-Diesel' cats(2)='OR-HeavyDuty-OnNetwork-Gas' ncats=2 elseif (runtyp .eq. 'LDON') then allocate(cats(3)) cats(1)='OR-LightDuty-OnNetwork-Diesel' cats(2)='OR-LightDuty-OnNetwork-Gas' cats(3)='OR-Refueling' ncats=3 else !hotelling allocate(cats(1)) cats(1)='OR-HeavyDuty-Hoteling' ncats=1 endif elseif (runag) then allocate(cats(1)) cats(1)='NP-AgriculturalLivestock' ncats=1 else !rwc allocate(cats(1)) cats(1)='NP-ResidentialWoodCombustionRWC' ncats=1 endif c ihap=1 c do i=2,ihap1 c ldup=.false. c n=1 c do while (n .lt. i .and. .not. ldup) c if (cats1(i) .eq. cats1(n)) ldup=.true. c n=n+1 c enddo c if (.not. ldup) ihap=ihap+1 c enddo c ncats=ihap c allocate(cats(ncats)) c cats(1)=cats1(1) c ihap=1 c do i=2,ihap1 c ldup=.false. c n=1 c do while (n .lt. i .and. .not. ldup) c if (cats1(i) .eq. cats1(n)) ldup=.true. c n=n+1 c enddo c if (.not. ldup) then c ihap=ihap+1 c cats(ihap)=cats1(i) c endif c enddo c endif c do i=1,ncats c write(103,*)cats(i),ncats c enddo c write(*,*)'cats done' c deallocate(hap1) if (allocated(h)) deallocate(h) if (allocated(fac1)) deallocate(fac1) if (allocated(facid2)) deallocate(facid2) if (allocated(igrps)) deallocate(igrps) if (allocated(st2)) then deallocate(st2) deallocate(st2a) endif 5 format('# of runs: ',i10,/'exceeds # of sources: ',i10/ +'stopping program') return end c########################################################################## c subroutine to get gridded sources ids instead of checking emissions files subroutine gridsrc(r,c,fid,lf) use main1 implicit none integer*8 r,c logical lf character fid*15,ares*2,inpfil*300 write(fid,'(a1,i3.3,a1,i3.3)')'G',c,'R',r write(ares,'(i2.2)')gridres write(inpfil,'(8(a))')trim(adjustl(outdir(1))), +'aermod_',trim(adjustl(runtyp)),'_', +trim(adjustl(ares)),'_',trim(adjustl(fid)),'.out' inquire(file=inpfil,exist=lf) return end c########################################################################## subroutine src_cr use main1 implicit none integer*8 cr,ns,i,c,r,ig character fid*15,junk*1 c conus(1)=source is CONUS source c conus(2)=source is non-CONUS source conus=.false. readzero=.false. c determine if any sources are conus or nonconus, this is used in subroutine rec_matrix as well if (runcmv .or. runpt. or. runap) then ns=1 c see if any sources are non-CONUS sources do while (ns .le. nsrc .and. .not. conus(2)) if (state2(ns) .eq. 'AK' .or. state2(ns) .eq. 'HI' + .or. state2(ns) .eq. 'PR' .or. state2(ns) .eq. 'VI' + .or. state2(ns) .eq. '02' .or. state2(ns) .eq. '15' + .or. state2(ns) .eq. '72' .or. state2(ns) .eq. '78') + conus(2)=.true. ns=ns+1 enddo c see if any sources are CONUS sources ns=1 do while (ns .le. nsrc .and. .not. conus(1)) if (state2(ns) .ne. 'AK' .and. state2(ns) .ne. 'HI' + .and. state2(ns) .ne. 'PR' .and. state2(ns) .ne. 'VI' + .and. state2(ns) .ne. '02' .and. state2(ns) .ne. '15' + .and. state2(ns) .ne. '72' .and. state2(ns) .ne. '78') + conus(1)=.true. ns=ns+1 enddo else if (noncontr) then conus(2)=.true. !conus(1) is false else conus(1)=.true. !conus(2) is false endif c if conus(1) = true then get column/row of source c version 20002 now get column/row for all sources, not just CONUS if (runap .or. runcmv .or. runpt) then allocate(srccr1(nsrc)) allocate(stgrid1(nsrc)) srccr1=0 else allocate(srccr(nsrc)) allocate(stgrid(nsrc)) srccr=0 ncolrow=nsrc endif do ns=1,nsrc if (runap .or. runcmv .or. runpt) then call read_co(ns,cr) srccr1(ns)=cr if (srccr1(ns) .le. 0)readzero=.true. if (trim(adjustl(state2(ns))) .eq. '02') then ig=2 elseif (trim(adjustl(state2(ns))) .eq. '15') then ig=3 elseif (trim(adjustl(state2(ns))) .eq. '72' .or. + trim(adjustl(state2(ns))) .eq. '78') then ig=4 else ig=1 endif stgrid1(ns)=ig else fid=facid(ns) read(fid,'(a1,i3,a1,i3)')junk,c,junk,r if (c .le. 0 .or. r .le. 0) then readzero=.true. srccr(ns)=-1*(abs(c)*1000+abs(r)) else srccr(ns)=c*1000+r endif ! 20002, find grid for gridded source do ig=1,4 if (lgrids(ig)) stgrid(ns)=ig enddo endif enddo if (runap .or. runpt .or. runcmv) call sort_co return end c########################################################################## c subroutine to read CO pathway of AERMOD output file to get column and row subroutine read_co(ns,cr) use main1 implicit none integer*8 cr,nc,nr,i,eof,nu,ns character inpfil*300,line*250,rt(2)*10 logical lexist1,lgo,lexist2 c borrow from subroutine src_recs for non-gridded sources c for point check to see if non-egu or egu or airport is c runway or non-runway if (runap) then i=1 lexist1=.false. do while (i .le. 2 .and. .not. lexist1) !check for runway or non-runway write(inpfil,'(5(a))')trim(adjustl(outdir(i))), + trim(adjustl(state2(ns))),'/aermod_', + trim(adjustl(facid(ns))),'.out' inquire(file=inpfil,exist=lexist1) i=i+1 enddo if (.not. lexist1) then write(logunit,5)trim(adjustl(inpfil)) lgo=.false. goto 95 c stop endif elseif (runpt) then i=1 lexist1=.false. do while (i .le. 1 .and. .not. lexist1) !check for EGU or non-EGU write(inpfil,'(5(a))')trim(adjustl(outdir(i))), + trim(adjustl(state2(ns))),'/aermod_', + trim(adjustl(facid(ns))),'.out' inquire(file=inpfil,exist=lexist1) i=i+1 enddo if (.not. lexist1) then !need to see if source run was split i=1 lexist2=.false. do while (i .le. 2 .and. .not. lexist2) !check for EGU or non-EGU write(inpfil,'(5(a))')trim(adjustl(outdir(i))), + trim(adjustl(state2(ns))),'/aermod_', + trim(adjustl(facid(ns))),'_01.out' inquire(file=inpfil,exist=lexist2) i=i+1 enddo if (.not. lexist2) then write(logunit,5)trim(adjustl(inpfil)) lgo=.false. goto 95 endif c stop endif else write(inpfil,'(5(a))')trim(adjustl(outdir(1))), + trim(adjustl(state2(ns))),'/aermod_', + trim(adjustl(facid(ns))),'.out' inquire(file=inpfil,exist=lexist1) if (.not. lexist1) then write(logunit,5)trim(adjustl(inpfil)) lgo=.false. goto 95 endif endif eof=0 open(unit=inpunit,file=inpfil,status='old') read(inpunit,10,iostat=eof)line if (eof .eq. 0) then nc=index(line,'COL') nr=index(line,'ROW') nu=index(line,'UTM ZONE') if (index(line(nc+4:nr-1),'*') .gt. 0) then col=0 else read(line(nc+4:nr-1),*)col endif if (index(line(nr+4:nu-1),'*') .gt. 0) then row=0 else read(line(nr+4:nu-1),*)row endif cr=col*1000+row endif close(inpunit) 5 format('READ CO AERMOD file ',a,' not found'/'skip') 10 format(//a200) 95 return end c########################################################################## c subroutine to sort source column/row subroutine sort_co use main1 implicit none integer*8 i,n,ic logical ldup c 20002, now check state as well since using grids for all areas ic=1 do i=2,nsrc ldup=.false. n=1 do while (n .lt. i .and. .not. ldup) if (srccr1(i) .eq. srccr1(n) .and. stgrid1(i) .eq. + stgrid1(n)) ldup=.true. !duplicate if col/row and state match n=n+1 enddo if (.not. ldup) ic=ic+1 enddo ncolrow=ic ic=1 allocate(srccr(ncolrow)) allocate(stgrid(ncolrow)) srccr=0 stgrid=0 srccr(1)=srccr1(1) stgrid(1)=stgrid1(1) do i=2,nsrc ldup=.false. n=1 do while (n .lt. i .and. .not. ldup) if (srccr1(i) .eq. srccr1(n) .and. stgrid1(i) .eq. + stgrid1(n)) ldup=.true. n=n+1 enddo if (.not. ldup) then ic=ic+1 srccr(ic)=srccr1(i) stgrid(ic)=stgrid1(i) endif enddo return end c########################################################################## c create national matrix of receptors subroutine rec_matrix use main1 implicit none integer*8 eof,irec,iunit,c,r,gr,br,mr,cr,br1,mr1,gr1,i real*8 x,y,barea character rfile*300,rid*15,fcode*5 logical luse c version 17186 c allocate gridded network if in CONUS area c allocate only non-CONUS blocks/monitors when only non-CONUS sources c allocate only CONUS blocks/monitors when only CONUS areas c both if a mix of CONUS/non-CONUS write(logunit,*)'getting receptor arrays' iunit=60 c version 20002 c now have gridded receptors for all areas write(logunit,*)'Getting gridded receptors' call datetime c first read in gridded receptors file(s) c loop through the potential receptor files do i=1,4 if (lgrids(i)) then iunit=100*(60+i) write(rfile,'(3(a))')trim(adjustl(recdir)), + trim(adjustl(gridst(i))),'_gridded_receptors.csv' open(unit=iunit,file=rfile,status='old') eof=0 gr=0 read(iunit,*,iostat=eof)rid read(iunit,*,iostat=eof)c,r 5 if (eof .eq. 0) then call checkrec(c,r,i,luse) if (luse) gr=gr+1 read(iunit,*,iostat=eof)c,r goto 5 endif rewind(iunit) endif enddo ngrids=gr allocate(gridconc(nhaps,ngrids,ncats)) allocate(gridrec(ngrids)) gridrec%x=0.0 gridrec%y=0.0 gridrec%colrow=0 gridrec%col=0 gridrec%row=0 gridrec%irec=0 gridrec%irec1=0 gridrec%iflag=0 gridrec%recid='NA' gridrec%rtype='NA' gridconc=0.0 c now fill in array do i=1,4 if (lgrids(i)) then iunit=100*(60+i) eof=0 gr=0 gr1=0 read(iunit,*,iostat=eof)rid read(iunit,*,iostat=eof)c,r,irec,x,y 10 if (eof .eq. 0) then cr=c*1000+r gr1=gr1+1 call checkrec(c,r,i,luse) if (luse) then gr=gr+1 c version 20002, receptor ID now includes gridst write(rid,'(a4,i6.6,i4.4)') + trim(adjustl(gridst(i))),cr,irec c write(rid,'(i6,i4.4)')cr,irec gridrec(gr)%recid=rid gridrec(gr)%x=x gridrec(gr)%y=y gridrec(gr)%rtype='GRIDDED' gridrec(gr)%colrow=cr gridrec(gr)%col=c gridrec(gr)%row=r gridrec(gr)%irec=gr1 !master gridrec(gr)%irec1=gr endif read(iunit,*,iostat=eof)c,r,irec,x,y goto 10 endif close(iunit) endif enddo c now blocks c version 20002, not all blocks in non-CONUS modeled, some will be interpolated iunit=60 write(logunit,*)'Getting block receptors' call datetime write(rfile,'(2(a))')trim(adjustl(recdir)),'block_receptors.csv' open(unit=iunit,file=rfile,status='old') eof=0 br=0 read(iunit,*,iostat=eof)rid read(iunit,*,iostat=eof)rid,x,y,c,r,barea 15 if (eof .eq. 0) then if (rid(1:2) .eq. '02' .or. rid(1:2) .eq. '15' .or. + rid(1:2) .eq. '72' .or. rid(1:2) .eq. '78') then if (conus(2)) then if (rid(1:2) .eq. '02') then i=2 elseif (rid(1:2) .eq. '15') then i=3 else i=4 endif if (c .eq. 0 .and. r .eq. 0) then !keep blocks not in a grid luse=.true. else call checkrec(c,r,i,luse) endif if (luse) br=br+1 endif else if (conus(1)) then i=1 call checkrec(c,r,i,luse) if (luse) br=br+1 endif endif read(iunit,*,iostat=eof)rid,x,y,c,r,barea goto 15 endif rewind(iunit) nblocks=br br=0 br1=0 eof=0 ! allocate(blockid(nblocks)) if (nblocks .gt. 0) then allocate(blokconc(nhaps,nblocks,ncats)) allocate(blockrec(nblocks)) ! allocate(blockx(nblocks)) ! allocate(blocky(nblocks)) ! allocate(blockcr(nblocks)) blokconc=0.0 blockrec%x=0.0 blockrec%y=0.0 blockrec%colrow=0 blockrec%col=0 blockrec%row=0 blockrec%irec=0 blockrec%irec1=0 blockrec%iflag=0 blockrec%recid='NA' blockrec%rtype='NA' read(iunit,*,iostat=eof)rid read(iunit,*,iostat=eof)rid,x,y,c,r 20 if (eof .eq. 0) then cr=c*1000+r c br=br+1 br1=br1+1 if (rid(1:2) .eq. '02' .or. rid(1:2) .eq. '15' .or. + rid(1:2) .eq. '72' .or. rid(1:2) .eq. '78') then if (conus(2)) then if (rid(1:2) .eq. '02') then i=2 elseif (rid(1:2) .eq. '15') then i=3 else i=4 endif if (c .eq. 0 .and. r .eq. 0) then luse=.true. else call checkrec(c,r,i,luse) endif if (luse) br=br+1 else luse=.false. endif else if (conus(1)) then i=1 call checkrec(c,r,i,luse) if (luse) br=br+1 else luse=.false. endif endif if (luse) then blockrec(br)%recid=rid blockrec(br)%rtype='BLOCK' blockrec(br)%x=x blockrec(br)%y=y blockrec(br)%colrow=cr blockrec(br)%col=c blockrec(br)%row=r blockrec(br)%irec=br1 !master record # blockrec(br)%irec1=br endif ! blockid(br)=rid ! blockx(br)=x ! blocky(br)=y ! blockcr(br)=cr c write(108,*)blockrec(br)%recid,blockrec(br)%x,blockrec(br)%y, c + blockrec(br)%colrow read(iunit,*,iostat=eof)rid,x,y,c,r goto 20 endif endif close(iunit) c now nonpopulated blocks in non-conus areas c version 20002, now only if source is not in CMAQ domain c only in AK if (readzero) then write(logunit,*)'Getting nonpoulated block receptors' call datetime write(rfile,'(2(a))')trim(adjustl(recdir)), + 'zero_pop_block_receptors.csv' open(unit=iunit,file=rfile,status='old') eof=0 br=0 read(iunit,*,iostat=eof)rid read(iunit,*,iostat=eof)rid,x,y,c,r,barea 35 if (eof .eq. 0) then if (rid(1:2) .eq. '02' .and. c .eq. 0 .and. r .eq. 0 ) + br=br+1 read(iunit,*,iostat=eof)rid,x,y,c,r,barea goto 35 endif rewind(iunit) nzblocks=br br=0 br1=0 eof=0 ! allocate(blockid(nblocks)) if (nzblocks .gt. 0) then allocate(zblokconc(nhaps,nzblocks,ncats)) allocate(zblockrec(nzblocks)) ! allocate(blockx(nblocks)) ! allocate(blocky(nblocks)) ! allocate(blockcr(nblocks)) zblokconc=0.0 zblockrec%x=0.0 zblockrec%y=0.0 zblockrec%colrow=0 zblockrec%col=0 zblockrec%row=0 zblockrec%irec=0 zblockrec%irec1=0 zblockrec%iflag=0 zblockrec%recid='NA' zblockrec%rtype='NA' read(iunit,*,iostat=eof)rid read(iunit,*,iostat=eof)rid,x,y,c,r 40 if (eof .eq. 0) then cr=c*1000+r c br=br+1 br1=br1+1 if (rid(1:2) .eq. '02' .and. c .eq. 0 .and. + r .eq. 0) then br=br+1 luse=.true. else luse=.false. endif if (luse) then zblockrec(br)%recid=rid zblockrec(br)%rtype='BLOCK' zblockrec(br)%x=x zblockrec(br)%y=y zblockrec(br)%colrow=cr zblockrec(br)%col=c zblockrec(br)%row=r zblockrec(br)%irec=br1 !master record # zblockrec(br)%irec1=br endif ! blockid(br)=rid ! blockx(br)=x ! blocky(br)=y ! blockcr(br)=cr c write(108,*)blockrec(br)%recid,blockrec(br)%x,blockrec(br)%y, c + blockrec(br)%colrow read(iunit,*,iostat=eof)rid,x,y,c,r goto 40 endif endif close(iunit) endif c now monitors write(logunit,*)'Getting monitor receptors' call datetime write(rfile,'(2(a))')trim(adjustl(recdir)),'monitor_receptors.csv' open(unit=iunit,file=rfile,status='old') eof=0 mr=0 read(iunit,*,iostat=eof)rid read(iunit,*,iostat=eof)rid,fcode,x,y,c,r 25 if (eof .eq. 0) then if (fcode(1:2) .eq. '02' .or. fcode(1:2) .eq. '15' .or. + fcode(1:2) .eq. '72' .or. fcode(1:2) .eq. '78') then if (conus(2)) then if (fcode(1:2) .eq. '02') then i=2 elseif (fcode(1:2) .eq. '15') then i=3 else i=4 endif call checkrec(c,r,i,luse) if (luse) mr=mr+1 endif else if (conus(1)) then i=1 call checkrec(c,r,i,luse) if (luse) mr=mr+1 endif endif read(iunit,*,iostat=eof)rid,fcode,x,y,c,r goto 25 endif rewind(iunit) nmons=mr mr=0 mr1=0 eof=0 ! allocate(monid(nmons)) ! allocate(mfips(nmons)) if (nmons .gt. 0) then allocate(monconc(nhaps,nmons,ncats)) allocate(monrec(nmons)) ! allocate(monx(nmons)) ! allocate(mony(nmons)) ! allocate(moncr(nmons)) monconc=0.0 monrec%x=0.0 monrec%y=0.0 monrec%colrow=0 monrec%col=0 monrec%row=0 monrec%irec=0 monrec%irec1=0 monrec%iflag=0 monrec%recid='NA' monrec%rtype='NA' monrec%fip='NA' read(iunit,*,iostat=eof)rid read(iunit,*,iostat=eof)rid,fcode,x,y,c,r 30 if (eof .eq. 0) then cr=c*1000+r mr1=mr1+1 if (fcode(1:2) .eq. '02' .or. fcode(1:2) .eq. '15' .or. + fcode(1:2) .eq. '72' .or. fcode(1:2) .eq. '78') then if (conus(2)) then if (fcode(1:2) .eq. '02') then i=2 elseif (fcode(1:2) .eq. '15') then i=3 else i=4 endif call checkrec(c,r,i,luse) if (luse) mr=mr+1 else luse=.false. endif else if (conus(1)) then i=1 call checkrec(c,r,i,luse) if (luse) mr=mr+1 else luse=.false. endif endif if (luse) then monrec(mr)%x=x monrec(mr)%y=y monrec(mr)%colrow=cr monrec(mr)%col=c monrec(mr)%row=r monrec(mr)%irec=mr1 monrec(mr)%irec1=mr monrec(mr)%recid=rid monrec(mr)%fip=fcode monrec(mr)%rtype='MONITOR' endif ! monid(mr)=rid ! mfips(mr)=fcode ! monx(mr)=x ! mony(mr)=y ! moncr(mr)=cr c write(109,*)monrec(mr)%recid,monrec(mr)%x,monrec(mr)%y, c + monrec(mr)%fip,monrec(mr)%colrow read(iunit,*,iostat=eof)rid,fcode,x,y,c,r goto 30 endif endif write(logunit,*)'number of receptors ',ngrids,nblocks,nmons c write(*,'(a,3(1x,i8.0))')'Number of receptors', c + ngrids,nblocks,nmons close(iunit) return end c########################################################################## c subroutine to check if a receptor is within a specified block c of grid cells of source cells subroutine checkrec(c,r,ig,luse) use main1 implicit none integer*8 c,r,modcol,modrow,i,cr,col1,row1,col2,row2,scol,srow,ig logical luse if (.not. runcmv) then c modcol=5 c modrow=5 if (ig .eq. 1) then !US modcol=8 modrow=8 elseif (ig .eq. 2) then !AK modcol=7 modrow=7 else !HI OR PR/VI modcol=20 modrow=20 endif else modcol=40 modrow=40 endif cr=c*1000+r i=1 luse=.false. c version 20002, make sure looking at correct grid do while (i .le. ncolrow .and. .not. luse) scol=int(srccr(i)/1000) srow=srccr(i)-scol*1000 col1=scol-modcol col2=scol+modcol row1=srow-modrow row2=srow+modrow if (c .ge. col1 .and. c .le. col2 .and. r .ge. row1 .and. + r .le. row2 .and. ig .eq. stgrid(i)) luse=.true. i=i+1 enddo return end c########################################################################## c subroutine to loop through sources and read model output subroutine read_model use main1 implicit none integer*8 ns logical lgo do ns=1,nsrc if (trim(adjustl(facid(ns))) .ne. 'NA') then if (runpt .or. runap .or. runcmv) then write(logunit,'(2(a),1x,a)')'Start ', + trim(adjustl(state2(ns))),trim(adjustl(facid(ns))) else write(logunit,'(2(a))')'Start ', + trim(adjustl(facid(ns))) endif call datetime if (runcmv .or. runpt. or. runap) then ! .or. noncontr) then if (state2(ns) .eq. 'AK' + .or. state2(ns) .eq. 'HI' .or. state2(ns) .eq. 'PR' + .or. state2(ns) .eq. 'VI' .or. state2(ns) .eq. '02' + .or. state2(ns) .eq. '15' .or. state2(ns) .eq. '72' + .or. state2(ns) .eq. '78') then akhiprvi=.true. else akhiprvi=.false. endif else if (noncontr) then akhiprvi=.true. else akhiprvi=.false. endif endif if (writefac) call cleanfiles(ns) !check and delete HAP facility files before processing call src_recs(ns,lgo) !get source locations and modeled receptors if (lgo) then call int_recs(ns) !get blocks and monitors to interpolate c if (.not. akhiprvi) then call interp(ns) !interpolate blocks and monitors c endif call hapconc(ns) !calculate HAP specific concentrations and output results call cleanup !deallocate arrays specific to the source else write(logunit,*)'No processing for ', + trim(adjustl(facid(ns))) endif call datetime if (runpt .or. runcmv .or. runpt) then write(logunit,'(2(a),1x,a)')'End ', + trim(adjustl(state2(ns))),trim(adjustl(facid(ns))) else write(logunit,'(2(a))')'End ',trim(adjustl(facid(ns))) endif endif enddo return end c########################################################################## c subroutine to clean HAP-facility output files before new processing subroutine cleanfiles(ns) use main1 implicit none integer*8 ns,h character ofile*250 logical lexist write(logunit,*)'Deleting old HAP-facility output files' do h=1,nhaps write(ofile,'(5(a))')trim(adjustl(hapdir)), + trim(adjustl(haps(h))),'_', + trim(adjustl(facid(ns))),'_concs.txt' inquire(file=ofile,exist=lexist) if (lexist) then open(unit=99,file=ofile,status='old') close(99,status='delete') endif enddo return end c########################################################################## c subroutine to read in aermod output and plot file subroutine src_recs(ns,lgo) use main1 implicit none integer*8 eof,igrid,iblock,isrc,irec,nc,nr,nu,imon,n,ns,i, +imonth,nv,izblock,findsrc,m1,iyr,idir,ist1,get_state,imonth1 real*8 x,y,conc,z1,z2,z3,emis(12) c character block*15,pltfil*300,inpfil*300,line*250,ave*6,srcgrp*8, character pltfil*300,inpfil*300,line*250,ave*6, +srcgrp*12,rt(2)*10,ares*2,errfil*300,srcgrp1*12 logical lexist1,lrecstart,lrecend,lexist2,lexist3,lgo,lmonth, +lsplit,lexist,lexist1a,lexist2a write(logunit,*)'Reading source and modeled receptors' lsplit=.false. c for point check to see if non-egu or egu or airport is c runway or non-runway if (runap) then i=1 lexist1=.false. lexist2=.false. do while (i .le. 2 .and. .not. lexist1) write(inpfil,'(5(a))')trim(adjustl(outdir(i))), + trim(adjustl(state2(ns))), + '/aermod_',trim(adjustl(facid(ns))),'.out' inquire(file=inpfil,exist=lexist1) write(pltfil,'(5(a))')trim(adjustl(outdir(i))), + trim(adjustl(state2(ns))), + '/',trim(adjustl(facid(ns))),'.plt' inquire(file=pltfil,exist=lexist2) i=i+1 enddo elseif (runpt) then i=1 lexist1=.false. lexist2=.false. do while (i .le. 1 .and. .not. lexist1) write(inpfil,'(5(a))')trim(adjustl(outdir(i))), + trim(adjustl(state2(ns))), + '/aermod_',trim(adjustl(facid(ns))),'.out' inquire(file=inpfil,exist=lexist1) write(pltfil,'(5(a))')trim(adjustl(outdir(i))), + trim(adjustl(state2(ns))), + '/',trim(adjustl(facid(ns))),'.plt' inquire(file=pltfil,exist=lexist2) i=i+1 enddo if (.not. lexist1 .and. .not. lexist2) then i=1 lexist1a=.false. lexist2a=.false. do while (i .le. 2 .and. .not. lexist1a) write(inpfil,'(5(a))')trim(adjustl(outdir(i))), + trim(adjustl(state2(ns))), + '/aermod_',trim(adjustl(facid(ns))),'_01.out' inquire(file=inpfil,exist=lexist1a) write(pltfil,'(5(a))')trim(adjustl(outdir(i))), + trim(adjustl(state2(ns))), + '/',trim(adjustl(facid(ns))),'_01.plt' inquire(file=pltfil,exist=lexist2a) i=i+1 enddo if (lexist1a .or. lexist2a) lsplit=.true. idir=i-1 endif elseif (runcmv) then i=1 write(inpfil,'(6(a))')trim(adjustl(outdir(i))), + trim(adjustl(state2(ns))),'/aermod_', + trim(adjustl(facid(ns))),'.out' write(pltfil,'(6(a))')trim(adjustl(outdir(i))), + trim(adjustl(state2(ns))),'/', + trim(adjustl(facid(ns))),'.plt' c elseif (noncontr) then c ist1=get_state(ncst) c write(ares,'(i2.2)')gridres c write(inpfil,'(9(a))')trim(adjustl(outdir(i))), c + trim(adjustl(states1(ist1))),'/aermod_',trim(adjustl(runtyp)), c + '_',trim(adjustl(ares)),'_',trim(adjustl(facid(ns))),'.out' c if (runon .or. runnon) then c write(pltfil,'(9(a))')trim(adjustl(outdir(i))), c + trim(adjustl(states1(ist1))),'/',trim(adjustl(runtyp)),'_' c + ,trim(adjustl(ares)),'_',trim(adjustl(facid(ns))), c + '_month.pst' c write(errfil,'(5(a),i2.2,3(a))')trim(adjustl(outdir(i))), c + trim(adjustl(states1(ist1))),'/',trim(adjustl(runtyp)),'_' c + ,gridres,'_',trim(adjustl(facid(ns))),'_errors.out' c else c write(pltfil,'(9(a))')trim(adjustl(outdir(i))), c + trim(adjustl(states1(ist1))),'/',trim(adjustl(runtyp)),'_' c + ,trim(adjustl(ares)),'_',trim(adjustl(facid(ns))),'.plt' c endif else !gridded, check for seasonal files for mobile, otherwise, annual files i=1 write(ares,'(i2.2)')gridres write(inpfil,'(8(a))')trim(adjustl(outdir(i))), + 'aermod_',trim(adjustl(runtyp)),'_', + trim(adjustl(ares)),'_',trim(adjustl(facid(ns))),'.out' if (runon .or. runnon) then write(pltfil,'(7(a))')trim(adjustl(outdir(i))), + trim(adjustl(runtyp)),'_',trim(adjustl(ares)),'_', + trim(adjustl(facid(ns))),'_month.pst' write(errfil,'(3(a),i2.2,3(a))')trim(adjustl(outdir(i))), + trim(adjustl(runtyp)),'_',gridres,'_', + trim(adjustl(facid(ns))),'_errors.out' else write(pltfil,'(7(a))')trim(adjustl(outdir(i))), + trim(adjustl(runtyp)),'_',trim(adjustl(ares)),'_', + trim(adjustl(facid(ns))),'.plt' endif endif inquire(file=inpfil,exist=lexist1) inquire(file=pltfil,exist=lexist2) if (runon .or. runnon) then inquire(file=errfil,exist=lexist3) if (.not. lexist1 .or. .not. lexist2 .or. .not. lexist3) then if (.not. lexist1) write(logunit,5)trim(adjustl(inpfil)) if (.not. lexist2) write(logunit,5)trim(adjustl(pltfil)) if (.not. lexist3) write(logunit,5)trim(adjustl(errfil)) lgo=.false. goto 95 endif else if (.not. lexist1 .or. .not. lexist2) then if (.not. lexist1) write(logunit,5)trim(adjustl(inpfil)) if (.not. lexist2) write(logunit,5)trim(adjustl(pltfil)) lgo=.false. goto 95 endif endif c seasonal files: runtyp_res_facid_season.plt res is 04 or 12 c read aermod.out file to get # of sources, # of gridded receptors, # of block receptors, and # of monitors c that were explicitly modeled (this number for later writing of HAP specific concentrations) igrid=0 iblock=0 izblock=0 eof=0 isrc=0 imon=0 lrecstart=.false. lrecend=.false. lgo=.true. modyear=0 open(unit=inpunit,file=inpfil,status='old') read(inpunit,10,iostat=eof)line 15 if (eof .eq. 0) then if (index(line,'CO TITLEONE') .gt. 0) + read(line(13:16),'(i4)')modyear if (index(line,'CO TITLETWO') .gt. 0) then nc=index(line,'COL') nr=index(line,'ROW') nu=index(line,'UTM ZONE') if (index(line(nc+4:nr-1),'*') .gt. 0) then col=0 else read(line(nc+4:nr-1),*)col endif if (index(line(nr+4:nu-1),'*') .gt. 0) then row=0 else read(line(nr+4:nu-1),*)row endif c read(line(nc+4:nr-1),*)col c read(line(nr+4:nu-1),*)row read(line(nu+8:nu+18),*)utmzone endif if (index(line,'SO LOCATION') .gt. 0 .or. + index(line,'** LOCATION') .gt. 0)isrc=isrc+1 if (index(line,'** NUMBER OF BLOCKS') .gt. 0) then n=index(line,'** NUMBER OF BLOCKS') read(line(n+20:250),*)iblock endif if (index(line,'** NUMBER OF ZERO BLOCKS') .gt. 0) then n=index(line,'** NUMBER OF ZERO BLOCKS') read(line(n+25:250),*)izblock endif if (index(line,'** NUMBER OF GRIDDED RECEPTORS') .gt. 0) then n=index(line,'** NUMBER OF GRIDDED RECEPTORS') read(line(n+31:250),*)igrid endif if (index(line,'** NUMBER OF MONITORS') .gt. 0) then n=index(line,'** NUMBER OF MONITORS') read(line(n+22:250),*)imon endif if (index(line,'AERMOD Finishes UN-successfully') .gt. 0) + lgo=.false. read(inpunit,10,iostat=eof)line goto 15 endif c read errors file to get # of valid hours per month if (runon .or. runnon) call readerr(errfil) ngrps=isrc if (runap) then ngrps1=1 else ngrps1=ngrps endif c if (runpt) then c ngrps=isrc c else c ngrps=1 c endif nmodg=igrid nmodb=iblock nmodm=imon nmodz=izblock c write(*,*)nmodb,nmodg,nmodm,ngrps c pause if (lgo) then if (runap) then allocate(apsrcs(ngrps)) apsrcs%utmx1=-999.0 apsrcs%utmy1=-999.0 apsrcs%utmx2=-999.0 apsrcs%utmy2=-999.0 apsrcs%frac=-999.0 apsrcs%xdim=-999.0 apsrcs%ydim=-999.0 apsrcs%area=-999.0 apsrcs%angle=-999.0 apsrcs%srctyp='NA' apsrcs%srcid='NA' elseif (runpt) then allocate(ptsrcs(ngrps)) ptsrcs%srcid='NA' ptsrcs%srctyp='NA' ptsrcs%utmx=-999.0 ptsrcs%utmy=-999.0 ptsrcs%xdim=-999.0 ptsrcs%ydim=-999.0 ptsrcs%angle=-999.0 elseif (runcmv) then allocate(cmvsrcs(ngrps)) cmvsrcs%srcid='NA' cmvsrcs%srctyp='NA' cmvsrcs%utmx1=-999.0 cmvsrcs%utmy1=-999.0 cmvsrcs%nvert=-9 c elseif (noncontr) then c allocate(tractsrcs(ngrps)) c tractsrcs%srcid='NA' c tractsrcs%srctyp='NA' c tractsrcs%utmx1=-999.0 c tractsrcs%utmy1=-999.0 c tractsrcs%nvert=-9 else allocate(gridsrcs(ngrps)) gridsrcs%utmx1=-999.0 gridsrcs%utmy1=-999.0 gridsrcs%utmx2=-999.0 gridsrcs%utmy2=-999.0 gridsrcs%utmx3=-999.0 gridsrcs%utmy3=-999.0 gridsrcs%utmx4=-999.0 gridsrcs%utmy4=-999.0 gridsrcs%nvert=-9 gridsrcs%srcid='NA' endif else goto 95 endif allocate(memis(ngrps,13)) memis=0.0 emis=0.0 if (nmonths .eq. 1) then do isrc=1,ngrps memis(isrc,13)=10000. enddo endif if (nmodg .gt. 0) then allocate(modgc(ngrps1,nmodg,nmonths)) allocate(modgrid(nmodg)) modgc=0.0 modgrid%ux=0.0 modgrid%uy=0.0 modgrid%x=0.0 modgrid%y=0.0 modgrid%colrow=0 modgrid%iflag=0 modgrid%recid='NA' modgrid%mastrec=0 modgrid%dist=0.0 endif if (nmodb .gt. 0) then allocate(modbc(ngrps1,nmodb,nmonths)) allocate(modblock(nmodb)) modbc=0.0 modblock%ux=0.0 modblock%uy=0.0 modblock%x=0.0 modblock%y=0.0 modblock%recid='NA' modblock%iflag=0 modblock%mastrec=0 modblock%dist=-9.0 endif if (nmodz .gt. 0) then allocate(modzbc(ngrps1,nmodz,nmonths)) allocate(zmodblock(nmodz)) modzbc=0.0 zmodblock%ux=0.0 zmodblock%uy=0.0 zmodblock%x=0.0 zmodblock%y=0.0 zmodblock%recid='NA' zmodblock%iflag=0 zmodblock%mastrec=0 zmodblock%dist=-9.0 endif if (nmodm .gt. 0) then allocate(modmc(ngrps1,nmodm,nmonths)) allocate(modmon(nmodm)) modmc=0.0 modmon%ux=0.0 modmon%uy=0.0 modmon%x=0.0 modmon%y=0.0 modmon%iflag=0.0 modmon%recid='NA' modmon%mastrec=0 modmon%dist=-9.0 endif call subsetrecs(ns) rewind(inpunit) c get gridded receptors Lambert coordinates,column/row c get block ID and column/row eof=0 igrid=0 iblock=0 izblock=0 imon=0 irec=0 lrecstart=.false. lrecend=.false. isrc=0 lmonth=.false. read(inpunit,10,iostat=eof)line 20 if (eof .eq. 0) then if (index(line,'SO LOCATION') .gt. 0 .or. + index(line,'** LOCATION') .gt. 0) then isrc=isrc+1 call source_info(line,isrc) endif if (index(line,'SO SRCPARAM') .gt. 0 .or. + index(line,'** SRCPARAM')) then call source_info(line,isrc) if (runcmv) then ! .or. noncontr) then nv=cmvsrcs(isrc)%nvert c if (runcmv) then c nv=cmvsrcs(isrc)%nvert c else c nv=tractsrcs(isrc)%nvert c endif if (isrc .eq. 1) then allocate(polyx(ngrps,nv)) allocate(polyy(ngrps,nv)) polyx=-999.0 polyy=-999.0 endif endif endif if (index(line,'SO AREAVERT') .gt. 0) then call source_info(line,isrc) endif c seasonal emissions version 17243 if (nmonths .eq. 12) then if (index(line,'** MONTHEMIS') .gt. 0) then if (.not. lmonth) then lmonth=.true. m1=1 else lmonth=.false. m1=7 endif read(line(14:250),*)srcgrp,(emis(i),i=m1,m1+5) isrc=findsrc(srcgrp) if (isrc .gt. 0) then do i=1,13 if (i .le. 12) then memis(isrc,i)=emis(i) else memis(isrc,13)=10000. endif enddo endif endif endif if (index(line,'RE STARTING') .gt. 0) lrecstart=.true. if (index(line,'RE FINISHED') .gt. 0) lrecend=.true. if (lrecstart .and. .not. lrecend) then if (index(line,'**') .gt. 0) then irec=irec+1 if (index(line,'**COLROW=') .gt. 0) then igrid=igrid+1 !gridded receptor call recinfo(ns,line,igrid,irec) endif if (index(line,'**BLOCK=') .gt. 0) then iblock=iblock+1 !block receptor call recinfo(ns,line,iblock,irec) endif if (index(line,'**ZEROBLOCK=') .gt. 0) then izblock=izblock+1 !block receptor call recinfo(ns,line,izblock,irec) endif if (index(line,'**MONITOR=') .gt. 0) then imon=imon+1 !block receptor call recinfo(ns,line,imon,irec) endif endif endif read(inpunit,10,iostat=eof)line goto 20 endif close(inpunit) c do i=1,tractsrcs(isrc)%nvert !JAT c write(101,*)polyx(i),polyy(i) c enddo c now read the plt file c order of receptors by source group is block,gridded,monitor eof=0 irec=0 write(logunit,*)'Reading model concentrations' if (lsplit) then i=1 lexist=.true. do while (lexist) write(pltfil,'(5(a),i2.2,a)')trim(adjustl(outdir(idir))), + trim(adjustl(state2(ns))),'/', + trim(adjustl(facid(ns))),'_',i,'.plt' inquire(file=pltfil,exist=lexist) if (lexist) then write(logunit,*)'reading ',trim(adjustl(pltfil)) open(unit=pltunit,file=pltfil,status='old') srcgrp1='A' read(pltunit,10,iostat=eof)line 40 if (eof .eq. 0) then if (line(1:1) .eq. '*') then irec=0 if (nmonths .eq. 12) then imonth1=0 else imonth1=1 endif else if (nmonths .eq. 12) then read(line,30)x,y,conc,z1,z2,z3,ave,srcgrp, + iyr,imonth else read(line,25)x,y,conc,z1,z2,z3,ave,srcgrp imonth=1 endif if (trim(adjustl(srcgrp)) .ne. + trim(adjustl(srcgrp1))) then irec=0 srcgrp1=srcgrp endif if (imonth .ne. imonth1) then irec=0 imonth1=imonth endif irec=irec+1 call assignmod(x,y,conc,irec,imonth,srcgrp) endif read(pltunit,10,iostat=eof)line goto 40 endif close(pltunit) endif i=i+1 enddo else open(unit=pltunit,file=pltfil,status='old') read(pltunit,10,iostat=eof)line 35 if (eof .eq. 0) then if (line(1:1) .eq. '*') then irec=0 if (nmonths .eq. 12) then imonth1=0 else imonth1=1 endif else if (nmonths .eq. 12) then read(line,30)x,y,conc,z1,z2,z3,ave,srcgrp,iyr, + imonth else read(line,25)x,y,conc,z1,z2,z3,ave,srcgrp imonth=1 endif if (trim(adjustl(srcgrp)) .ne. + trim(adjustl(srcgrp1))) then irec=0 srcgrp1=srcgrp endif if (imonth .ne. imonth1) then irec=0 imonth1=imonth endif irec=irec+1 call assignmod(x,y,conc,irec,imonth,srcgrp) endif read(pltunit,10,iostat=eof)line goto 35 endif close(pltunit) endif 5 format('SRC_REC AERMOD file ',a,' not found'/'skip') 10 format(a200) 25 format(2(1X,F13.5),1X,E13.6,3(1X,F8.2),2X,A6,2X,A8) c30 format(2(1X,F13.5),1X,E13.6,3(1X,F7.2),2X,A8,2X,2(i2)) 30 format(2(1X,F13.5),1X,E13.6,3(1X,F8.2),2X,A6,2X,A8,2X,2(i2.2)) 95 return end c########################################################################## subroutine readerr(efile) use main1 implicit none character efile*300,line*100 integer*8 eof,m,nh logical lleap c determine if year is leap year call leapyr(lleap) c calculate max # of hours per month do m=1,12 if (m .eq. 2 .and. lleap) then nhours(m)=(days(m)+1)*24 else nhours(m)=days(m)*24 endif enddo nh=0 c read the errors file and subtract 1 from the hourly total for the appropriate month c when a calm or missing hour message is encountered open(unit=errunit,file=efile,status='old') read(errunit,5,iostat=eof)line 10 if (eof .eq. 0) then if (index(line,'MX I440') .gt. 0 .or. + index(line,'MX 1460') .gt. 0) then read(line(88:89),'(i2)')m nhours(m)=nhours(m)-1 endif read(errunit,5,iostat=eof)line goto 10 endif close(errunit) write(logunit,15) do m=1,12 nh=nh+nhours(m) write(logunit,20)m,nhours(m) enddo write(logunit,25)nh 5 format(a100) 15 format(/'Number of valid hours per month'/) 20 format(i2.2,1x,i3) 25 format(/'Total valid hours:',1x,i4/) return end c########################################################################## c subroutine leapyr(iyear,leap) subroutine leapyr(lleap) use main1 implicit none c iyear: integer 4-digit year of processed year c ileap4: integer result of MOD(iyear,4) c ileap100: integer result of MOD(iyear,100) c ileap400: integer result of MOD(iyear,400) integer ileap4,ileap100,ileap400 logical lleap ileap4=mod(modyear,4) ileap100=mod(modyear,100) ileap400=mod(modyear,400) c if year is not divisible by 4 or year is divisible by 100 but not divisible by 400 then year c is not a leap year c 2001 is not a leap year because it is not divisible by 4 c 2000 is is a leap year because it is divisible by 400 if (ileap4 .ne. 0 .or. + (ileap100 .eq. 0 .and. ileap400 .ne. 0)) then lleap=.false. else lleap=.true. endif return end c########################################################################## subroutine source_info(line,is) use main1 implicit none integer*8 is,nn real*8 e,ht,j3,x1,y1,x2,y2 character line*250,src*12,stype*8,j1*2,j2*8 if (index(line,'SO LOCATION') .gt. 0 .or. +index(line,'** LOCATION') .gt. 0) then if (runap) then read(line,*)j1,j2,src,stype if (stype .eq. 'LINE') then read(line,*)j1,j2,apsrcs(is)%srcid,apsrcs(is)%srctyp, + apsrcs(is)%utmx1,apsrcs(is)%utmy1,apsrcs(is)%utmx2, + apsrcs(is)%utmy2 else read(line,*)j1,j2,apsrcs(is)%srcid,apsrcs(is)%srctyp, + apsrcs(is)%utmx1,apsrcs(is)%utmy1 endif elseif (runpt) then read(line,*)j1,j2,ptsrcs(is)%srcid,ptsrcs(is)%srctyp, + ptsrcs(is)%utmx,ptsrcs(is)%utmy elseif (runcmv) then c read(line,*)j1,j2,cmvsrcs(is)%srcid,cmvsrcs(is)%srctyp, read(line,*)j1,j2,src,cmvsrcs(is)%srctyp, + cmvsrcs(is)%utmx1,cmvsrcs(is)%utmy1 nn=index(src,'c') write(cmvsrcs(is)%srcid,'(a1,a1)')'C',src(nn+1:nn+1) c elseif (noncontr) then c read(line,*)j1,j2,tractsrcs(is)%srcid, c + tractsrcs(is)%srctyp,tractsrcs(is)%utmx1, c + tractsrcs(is)%utmy1 else read(line,*)j1,j2,gridsrcs(is)%srcid,stype, + gridsrcs(is)%utmx1,gridsrcs(is)%utmy1 endif endif if (index(line,'SO SRCPARAM') .gt. 0 .or. +index(line,'** SRCPARAM') .gt. 0) then if (runap) then read(line,*)j1,j2,src,e,ht,apsrcs(is)%xdim,j3 if (apsrcs(is)%srctyp .eq. 'AREA') apsrcs(is)%ydim=j3 elseif (runpt) then if(ptsrcs(is)%srctyp .eq. 'AREA') + read(line,*)j1,j2,src,e,ht,ptsrcs(is)%xdim,ptsrcs(is)%ydim elseif (runcmv) then read(line,*)j1,j2,src,e,ht,cmvsrcs(is)%nvert c elseif (noncontr) then c read(line,*)j1,j2,src,e,ht,tractsrcs(is)%nvert else read(line,*)j1,j2,src,e,ht,gridsrcs(is)%nvert endif endif if (index(line,'SO AREAVERT') .gt. 0) then c if ((runon .or. runnon .or. runrwc .or. runnplow c + .or. runnphi .or. runog) .and. .not. noncontr) then if (.not. runcmv) then read(line,*)j1,j2,src,x1,y1,x2,y2 if (gridsrcs(is)%utmx2 .eq. -999.0) then c gridsrcs(is)%utmx1=x1 c gridsrcs(is)%utmy1=y1 gridsrcs(is)%utmx2=x2 gridsrcs(is)%utmy2=y2 else gridsrcs(is)%utmx3=x1 gridsrcs(is)%utmy3=y1 gridsrcs(is)%utmx4=x2 gridsrcs(is)%utmy4=y2 endif else call get_vertices(line,is) endif endif c write(107,*)ptsrcs(is)%srcid,ptsrcs(is)%srctyp,ptsrcs(is)%utmx, c + ptsrcs(is)%utmy return end c########################################################################## c subroutine to fill in CMV vertices array subroutine get_vertices(line,is) use main1 implicit none integer*8 iv,i,ilen,j,j1,j2,v,nv,is,nvmaster character line*250,line1*250 logical lmiss real*8, allocatable, dimension(:) :: x real*8, allocatable, dimension(:) :: y c check to see where the vertices array has been filled in v=1 lmiss=.false. if (runcmv) then nvmaster=cmvsrcs(is)%nvert else nvmaster=tractsrcs(is)%nvert endif do while (v .le. nvmaster .and. .not. lmiss) if (polyx(is,v) .eq. -999.0) then lmiss=.true. else v=v+1 endif enddo c subset line to part of line past source ID line1=line(25:250) ilen=len_trim(line1) c loop through line1 to see how many pairs of numbers there will be on the line c most is 3 pairs iv=0 do i=1,ilen j=ichar(line1(i:i)) if (j .eq. 32 .and. i .gt. 1) then !ignore if first character is blank j1=ichar(line1(i-1:i-1)) j2=ichar(line1(i+1:i+1)) if (j1 .ge. 48 .and. j1 .le. 57 .and. j2 .ge. 48 .and. + j2 .le. 57)iv=iv+1 endif enddo allocate(x(iv)) allocate(y(iv)) nv=iv c now fill in array read(line1,*)(x(iv),y(iv),iv=1,nv) do iv=1,nv polyx(is,v+iv-1)=x(iv) polyy(is,v+iv-1)=y(iv) enddo deallocate(x) deallocate(y) return end c########################################################################## c subroutine to fill in receptor arrays subroutine recinfo(is,line,i,j) use main1 implicit none integer*8 i,j,ncr,nx,ny,nblk,nmon,cr,nf,nr,iflag,irec, +findrec,recflag,nzblk,is,igrid real*8 x,y,ux,uy,dist,distance,zerod,maxd character line*200,recid*15 logical lwrite zerod=0.0 nx=index(line,'CMAQX=') ny=index(line,'CMAQY=') ncr=index(line,'COLROW=') nblk=index(line,'**BLOCK=') nzblk=index(line,'**ZEROBLOCK=') nmon=index(line,'**MONITOR') nr=index(line,'REC=') nf=index(line,'FLAG =') read(line(ncr+7:nx-1),*)cr read(line(nx+6:ny-1),*)x read(line(ny+6:200),*)y read(line(nf+7:nf+7),*)iflag read(inpunit,5)ux,uy c modux(j)=ux c moduy(j)=uy if (state2(is) .eq. '02') then igrid=2 elseif (state2(is) .eq. '15') then igrid=3 elseif (state2(is) .eq. '72' .or. state2(is) .eq. '78') then igrid=4 else igrid=1 endif if (nblk .gt. 0 .and. nmodb .gt. 0) then !block recid=line(nblk+9:nblk+23) read(line(ncr+7:ncr+14),*)cr modblock(i)%recid=recid modblock(i)%colrow=cr modblock(i)%ux=ux modblock(i)%uy=uy modblock(i)%x=x modblock(i)%y=y modblock(i)%iflag=iflag modblock(i)%mastrec=findrec('BM',i) c version 20002 if (cr .le. 0) then maxd=grddist else maxd=blkdist endif c if (.not. akhiprvi) then c maxd=blkdist c else c maxd=grddist c endif call recdist('mbl',recid,ux,uy,zerod,maxd,lwrite,recflag, + dist) !get distance modblock(i)%dist=dist if (.not. lwrite) modblock(i)%iflag=3 !fix for bug in aermod input program c write(108,*)modblock(i)%recid,modblock(i)%x,modblock(i)%y, c +modblock(i)%dist elseif (nzblk .gt. 0 .and. nmodz .gt. 0) then !zero pop block recid=line(nzblk+13:nzblk+27) read(line(ncr+7:ncr+14),*)cr zmodblock(i)%recid=recid zmodblock(i)%colrow=cr zmodblock(i)%ux=ux zmodblock(i)%uy=uy zmodblock(i)%x=x zmodblock(i)%y=y zmodblock(i)%iflag=iflag zmodblock(i)%mastrec=findrec('ZB',i) if (.not. akhiprvi) then maxd=blkdist else maxd=grddist endif call recdist('zbl',recid,ux,uy,zerod,maxd,lwrite,recflag, + dist) !get distance zmodblock(i)%dist=dist if (.not. lwrite) zmodblock(i)%iflag=3 !fix for bug in aermod input program elseif (nmon .gt. 0 .and. nmodm .gt. 0) then !monitor c recid=line(nmon+11:nmon+19) recid=line(nmon+11:ncr-1) read(line(ncr+7:ncr+14),*)cr modmon(i)%recid=recid modmon(i)%colrow=cr modmon(i)%ux=ux modmon(i)%uy=uy modmon(i)%x=x modmon(i)%y=y modmon(i)%iflag=iflag modmon(i)%mastrec=findrec('MM',i) c version 20002 maxd=blkdist c if (.not. akhiprvi) then c maxd=blkdist c else c maxd=grddist c endif call recdist('mmm',recid,ux,uy,zerod,maxd,lwrite,recflag, + dist) !get distance modmon(i)%dist=dist if (.not. lwrite) modmon(i)%iflag=3 !fix for bug in aermod input program c write(110,*)modmon(i)%recid,modmon(i)%x,modmon(i)%y, c +modmon(i)%dist else read(line(ncr+7:ncr+13),*)cr read(line(nr+4:nr+7),*)irec write(recid,'(a4,i6.6,i4.4)')gridst(igrid),cr,irec modgrid(i)%recid=recid modgrid(i)%colrow=cr modgrid(i)%ux=ux modgrid(i)%uy=uy modgrid(i)%x=x modgrid(i)%y=y modgrid(i)%iflag=iflag modgrid(i)%mastrec=findrec('GM',i) call recdist('mgr',recid,ux,uy,zerod,grddist,lwrite,recflag, + dist) !get distance modgrid(i)%dist=dist if (.not. lwrite) then modgrid(i)%iflag=3 else !account for being close to one source but not another if(modgrid(i)%iflag .eq. 3) modgrid(i)%iflag=0 endif c write(109,*)i,modgrid(i)%recid,modgrid(i)%x,modgrid(i)%y, c +modgrid(i)%dist,modgrid(i)%iflag,modgrid(i)%ux,modgrid(i)%uy, c +lwrite endif 5 format(11x,f13.3,1x,f13.3) return end c########################################################################## c subset master receptor list to within plus/minus 5 cells for faster processing c when getting master receptor numbers subroutine subsetrecs(ns) use main1 implicit none integer*8 i,c,r,col1,col2,row1,row2,modcol,modrow,ig,ns,get_state, +ig1 character astate*2,st*2 if (akhiprvi) then if (noncontr) then i=get_state(ncst) astate=states1(i) else if (state2(ns) .eq. 'AK' .or. state2(ns) .eq. '02') then astate='02' elseif (state2(ns) .eq. 'HI' .or. state2(ns) .eq. '15') + then astate='15' elseif (state2(ns) .eq. 'PR' .or. state2(ns) .eq. '78') + then astate='72' else astate='78' endif endif if (astate .eq. '02') then modcol=7 modrow=7 ig1=2 elseif (astate .eq. '15') THEN modcol=20 modrow=20 ig1=3 else modcol=20 modrow=20 ig1=4 endif col1=col-modcol col2=col+modcol row1=row-modrow row2=row+modrow else astate='00' ig1=1 if (.not. runcmv) then c modcol=5 c modrow=5 modcol=8 modrow=8 else modcol=40 modrow=40 endif endif col1=col-modcol col2=col+modcol row1=row-modrow row2=row+modrow ngrids1=0 c gridded do i=1,ngrids if (gridrec(i)%col .ge. col1 .and. gridrec(i)%col .le. + col2 .and. gridrec(i)%row .ge. row1 .and. + gridrec(i)%row .le. row2 .and. gridrec(i)%recid(1:4) .eq. + gridst(ig1)) ngrids1=ngrids1+1 enddo if (ngrids1 .gt. 0) then allocate(gridrec1(ngrids1)) ig=0 do i=1,ngrids if (gridrec(i)%col .ge. col1 .and. gridrec(i)%col .le. + col2 .and. gridrec(i)%row .ge. row1 .and. + gridrec(i)%row .le. row2 .and. gridrec(i)%recid(1:4) + .eq. gridst(ig1)) then ig=ig+1 gridrec1(ig)%recid=gridrec(i)%recid gridrec1(ig)%irec=gridrec(i)%irec1 endif enddo endif c blocks nblocks1=0 do i=1,nblocks if (.not. akhiprvi) then if (blockrec(i)%col .ge. col1 .and. blockrec(i)%col .le. + col2 .and. blockrec(i)%row .ge. row1 .and. + blockrec(i)%row .le. row2) nblocks1=nblocks1+1 else st=blockrec(i)%recid(1:2) c account for fact that some receptors in VI/PR will overlap in 50 km domain if (((st .eq. '72' .or. st .eq. '78') .and. + (astate .eq. '72' .or. astate .eq. '78')) .or. + ((st .eq. '02' .or. st .eq. '15') .and. + st .eq. astate) .and. ((blockrec(i)%col .ge. col1 .and. + blockrec(i)%col .le.col2 .and. blockrec(i)%row .ge. row1 + .and.blockrec(i)%row .le. row2) .or. + (blockrec(i)%col .eq. 0 .and. blockrec(i)%row .eq. 0))) + then nblocks1=nblocks1+1 endif c if (st .eq. astate) nblocks1=nblocks1+1 endif enddo if (nblocks1 .gt. 0) then allocate(blockrec1(nblocks1)) ig=0 do i=1,nblocks if (.not. akhiprvi) then if (blockrec(i)%col .ge. col1 .and. blockrec(i)%col + .le. col2 .and. blockrec(i)%row .ge. row1 .and. + blockrec(i)%row .le. row2) then ig=ig+1 blockrec1(ig)%recid=blockrec(i)%recid blockrec1(ig)%irec=blockrec(i)%irec1 endif else st=blockrec(i)%recid(1:2) if (((st .eq. '72' .or. st .eq. '78') .and. + (astate .eq. '72' .or. astate .eq. '78')) .or. + ((st .eq. '02' .or. st .eq. '15') .and. + st .eq. astate) .and. ((blockrec(i)%col .ge. col1 + .and. blockrec(i)%col .le.col2 .and. blockrec(i)%row + .ge. row1 .and.blockrec(i)%row .le. row2) .or. + (blockrec(i)%col .eq. 0 .and. blockrec(i)%row + .eq. 0))) then ig=ig+1 blockrec1(ig)%recid=blockrec(i)%recid blockrec1(ig)%irec=blockrec(i)%irec1 endif c if (st .eq. astate) then c ig=ig+1 c blockrec1(ig)%recid=blockrec(i)%recid c blockrec1(ig)%irec=blockrec(i)%irec c endif endif enddo endif c nonpopulated blocks nzblocks1=0 do i=1,nzblocks if (.not. akhiprvi) then if (zblockrec(i)%col .ge. col1 .and. zblockrec(i)%col .le. + col2 .and. blockrec(i)%row .ge. row1 .and. + zblockrec(i)%row .le. row2) nzblocks1=nzblocks1+1 else st=zblockrec(i)%recid(1:2) c account for fact that some receptors in VI/PR will overlap in 50 km domain if (st .eq. '02' .and. st .eq. astate) then nzblocks1=nzblocks1+1 endif c if (st .eq. astate) nzblocks1=nzblocks1+1 endif enddo if (nzblocks1 .gt. 0) then c write(*,*)'nzblocks1 ',nzblocks1 allocate(zblockrec1(nzblocks1)) ig=0 do i=1,nzblocks if (.not. akhiprvi) then if (zblockrec(i)%col .ge. col1 .and. zblockrec(i)%col + .le. col2 .and. zblockrec(i)%row .ge. row1 .and. + zblockrec(i)%row .le. row2) then ig=ig+1 zblockrec1(ig)%recid=zblockrec(i)%recid zblockrec1(ig)%irec=zblockrec(i)%irec1 endif else st=zblockrec(i)%recid(1:2) if (st .eq. '02' .and. st .eq. astate) then ig=ig+1 zblockrec1(ig)%recid=zblockrec(i)%recid zblockrec1(ig)%irec=zblockrec(i)%irec1 endif c if (st .eq. astate) then c ig=ig+1 c blockrec1(ig)%recid=blockrec(i)%recid c blockrec1(ig)%irec=blockrec(i)%irec c endif endif enddo endif c monitors nmons1=0 do i=1,nmons if (.not. akhiprvi) then if (monrec(i)%col .ge. col1 .and. monrec(i)%col .le. + col2 .and. monrec(i)%row .ge. row1 .and. + monrec(i)%row .le. row2) nmons1=nmons1+1 else st=monrec(i)%fip(1:2) c account for fact that some receptors in VI/PR will overlap in 50 km domain if (((st .eq. '72' .or. st .eq. '78') .and. + (astate .eq. '72' .or. astate .eq. '78')) .or. + ((st .eq. '02' .or. st .eq. '15') .and. + st .eq. astate) .and. ((monrec(i)%col .ge. col1 .and. + monrec(i)%col .le.col2 .and. monrec(i)%row .ge. row1 + .and. monrec(i)%row .le. row2) .or. + (monrec(i)%col .eq. 0 .and. monrec(i)%row .eq. 0))) then nmons1=nmons1+1 endif c if (st .eq. astate) nmons1=nmons1+1 endif enddo if (nmons1 .gt. 0) then allocate(monrec1(nmons1)) ig=0 do i=1,nmons if (.not. akhiprvi) then if (monrec(i)%col .ge. col1 .and. monrec(i)%col .le. + col2 .and. monrec(i)%row .ge. row1 .and. + monrec(i)%row .le. row2) then ig=ig+1 monrec1(ig)%recid=monrec(i)%recid monrec1(ig)%irec=monrec(i)%irec1 endif else st=monrec(i)%fip(1:2) if (((st .eq. '72' .or. st .eq. '78') .and. + (astate .eq. '72' .or. astate .eq. '78')) .or. + ((st .eq. '02' .or. st .eq. '15') .and. + st .eq. astate) .and. ((monrec(i)%col .ge. col1 .and. + monrec(i)%col .le.col2 .and. monrec(i)%row .ge. row1 + .and. monrec(i)%row .le. row2) .or. (monrec(i)%col + .eq. 0 .and. monrec(i)%row .eq. 0))) then ig=ig+1 monrec1(ig)%recid=monrec(i)%recid monrec1(ig)%irec=monrec(i)%irec1 endif c if (st .eq. astate) then c ig=ig+1 c monrec1(ig)%recid=monrec(i)%recid c monrec1(ig)%irec=monrec(i)%irec c endif endif enddo endif return end c########################################################################## c subroutine to assign modeled concentrations to arrays subroutine assignmod(x,y,conc,irec,imonth,grp) use main1 implicit none integer*8 irec,irec1,igrp,findsrc,imonth real*8 x,y,conc,dist,distance character grp*12 !grp*8 if (runap) then !.or. runcmv) then !.or. noncontr) then igrp=1 else igrp=findsrc(grp) endif c blocks, gridded, and monitors for point, airports, cmv, and tract sources except onroad tracts if (runpt .or. runap .or. runcmv .or. +(noncontr .and. .not. runon .and. .not. runnon)) then if (irec .le. nmodb) then !blocks irec1=irec dist=distance(x,y,modblock(irec)%ux,modblock(irec)%uy) if (dist .le. 1.0) modbc(igrp,irec1,1)=conc else if (irec .le. nmodb+nmodz) then !nonpopulated blocks irec1=irec-nmodb dist=distance(x,y,zmodblock(irec1)%ux, + zmodblock(irec1)%uy) if (dist .le. 1.0) modzbc(igrp,irec1,1)=conc elseif (irec .le. nmodb+nmodg+nmodz) then !gridded irec1=irec-nmodb dist=distance(x,y,modgrid(irec1)%ux,modgrid(irec1)%uy) !should be the same receptor if (dist .le. 1.0) modgc(igrp,irec1,1)=conc else !monitor irec1=irec-nmodb-nmodg-nmodz dist=distance(x,y,modmon(irec1)%ux,modmon(irec1)%uy) modmc(igrp,irec1,1)=conc endif endif endif c blocks/monitors for onroad non-conus if ((runon .or. runnon) .and. noncontr) then !grid source onroad/nonroad c irec1=irec-(nmodb+nmodm+nmodz)*((iseas-1)*24+(ihr-1)) c irec1=irec-(nmodb+nmodm+nmodz)*(imonth-1) irec1=irec if (irec1 .le. nmodb) then !block dist=distance(x,y,modblock(irec1)%ux,modblock(irec1)%uy) !should be the same receptor if (dist .le. 1.0) then modbc(igrp,irec1,imonth)=modbc(igrp,irec1,imonth)+ + conc*real(nhours(imonth)) c if (igrp .eq. 1 .and. irec1 .eq. nmodb) c + seashrs(iseas,ihr)=nhrs endif elseif (irec1 .le. nmodb+nmodz) then !non-pop block irec1=irec1-nmodb c write(*,*)irec,irec1,imonth,nmodb,nmodz dist=distance(x,y,zmodblock(irec1)%ux,zmodblock(irec1)%uy) !should be the same receptor if (dist .le. 1.0) then modzbc(igrp,irec1,imonth)=modzbc(igrp,irec1,imonth)+ + conc*real(nhours(imonth)) c if (igrp .eq. 1 .and. irec1 .eq. nmodb) c + seashrs(iseas,ihr)=nhrs endif c write(115,*)modblock(irec1)%ux,modblock(irec1)%uy, c +modbc(igrp,irec1,iseas),iseas else !monitor c write(135,*)irec,irec1,nmodb,nmodm,x,y,iseas,ihr irec1=irec1-nmodb-nmodz dist=distance(x,y,modmon(irec1)%ux,modmon(irec1)%uy) !should be the same receptor if (dist .le. 1.0) then modmc(igrp,irec1,imonth)=modmc(igrp,irec1,imonth)+ + conc*real(nhours(imonth)) c if (igrp .eq. 1 .and. irec1 .eq. nmodm) c + seashrs(iseas,ihr)=nhrs endif c write(117,*)modmon(irec1)%ux,modmon(irec1)%uy, c +modmc(igrp,irec1,iseas),iseas,igrp,nhrs,seashrs(iseas,ihr) endif endif c gridded sources except onroad and nonroad, no blocks and no monitors if ((runon .or. runnon) .and. .not. noncontr) then c irec1=irec-(nmodg)*((iseas-1)*24+(ihr-1)) c irec1=irec-nmodg*(imonth-1) irec1=irec dist=distance(x,y,modgrid(irec1)%ux,modgrid(irec1)%uy) !should be the same receptor if (dist .le. 1.0) then modgc(igrp,irec1,imonth)=modgc(igrp,irec1,imonth)+ + conc*real(nhours(imonth)) c if (igrp .eq. 1 .and. irec1 .eq. nmodg) c + seashrs(iseas,ihr)=nhrs endif endif c all other gridded sources, no blocks and no monitors if (.not. noncontr .and. (runrwc .or. runnphi .or. runnplow .or. +runog .or. runag)) then irec1=irec dist=distance(x,y,modgrid(irec1)%ux,modgrid(irec1)%uy) !should be the same receptor if (dist .le. 1.0) modgc(igrp,irec1,1)=conc endif return end c########################################################################## c subroutine to get apppropriate receptors for interpolation subroutine int_recs(ns) use main1 implicit none c eof: end of file indicator c col1: western most column to search for receptors c row1: southern most row to search for receptors c col2: eastern most column to search for receptors c row2: northern most row to search for receptors c nblock: number of block receptors written to AERMOD.INP c ngrid: number of gridded receptors written to AERMOD.INP c nmon: number of monitor receptors written to AERMOD.INP c c: column of receptor c r: row of receptor c colrow: column/row of receptor c modcol: number of columns to search on either side of facility column c modrow: number of rows to search on either side of facility row c recflag: flag indicating if receptor is within 30 m of a source c irec: gridded receptor number (1-9 for 4 km, 1-44 for 1 km) integer*8 eof,col1,col2,row1,row2,c,r,colrow,modcol,modrow, +recflag,ns c x: receptor UTM-x coordinate c y: receptor UTM-y coordinate c cmaqx: receptor CMAQ x-coordinate c cmaqy: receptor CMAQ-Y coordinate c relev: receptor elevation c hill: receptor hill height c mindist: minimum distance to check real*8 x,y,cmaqx,cmaqy,relev,hill,dist,mindist c recfil: receptor filename c recid: census block ID c monid: monitor ID character recfil*250,recid*15,monid*15,tfil*250,ares*2,arun*10, +ahap*3,runtyp1*20 c lexist: logical variable denoting receptor file exists c lwrite: logical varible denoting to write receptor to AERMOD.INP logical lexist,lwrite c version 20002, now interpolate for all areas c initialize variables nintb=0 nintm=0 if (.not. runcmv) then c modcol=5 c modrow=5 if (state2(ns) .eq. '02') then !9km modcol=7 !63 km modrow=7 elseif (state2(ns) .eq. '15' .or. state2(ns) .eq. '72' .or. + state2(ns) .eq. '78') then !3km modcol=20 !60 km modrow=20 else modcol=8 modrow=8 endif 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 if (runpt .or. runap .or. runcmv) then mindist=blkdist else mindist=0.0 endif write(arun,'(i10)')runno !run # if (nhap1 .gt. 0) then write(ahap,'(i3.3)')nhap1 else ahap='001' endif if (lsingfac) then write(runtyp1,'(3(a))')trim(adjustl(runtyp)),'_', + trim(adjustl(singlfac)) elseif (lsingst) then write(runtyp1,'(3(a))')trim(adjustl(runtyp)),'_', + trim(adjustl(singlst)) else runtyp1=trim(adjustl(runtyp)) endif c write blocks within specified distances (10-50 km) for point, airports, and CMV c for gridded, all receptors within 50 km. c flag receptors within 30 m of a facility release point, 30 m within runway (entire length), c or side of a CMV polygon) if (col .ge. 0 .and. row .ge. 0) then c temporary file for getting receptors to interpolate c write(tfil,'(3(a))')'temp_recs_',trim(adjustl(singlhap)), c + '.txt' if (runcmv .or. runpt .or. runap) then write(tfil,'(8(a))')'temp_recs_',trim(adjustl(singlhap)), + '_',trim(adjustl(runtyp1)),trim(adjustl(arun)),'_', + trim(adjustl(ahap)),'.txt' else write(ares,'(i2)')gridres write(tfil,'(10(a))')'temp_recs_',trim(adjustl(singlhap)), + '_',trim(adjustl(runtyp1)),'_',trim(adjustl(ares)), + trim(adjustl(arun)),'_',trim(adjustl(ahap)),'.txt' endif open(unit=98,file=tfil,status='unknown') write(logunit,*)'reading blocks' if (utmzone .ge. 10) then write(recfil,'(2a,i2,a)')trim(adjustl(recdir)), + 'blocks_utm_',utmzone,'.dat' else write(recfil,'(2a,i1,a)')trim(adjustl(recdir)), + 'blocks_utm_',utmzone,'.dat' endif lexist=.false. inquire(file=recfil,exist=lexist) if (.not. lexist) then write(logunit,35)trim(adjustl(recfil)) stop endif open(unit=recunit,file=recfil,status='old') read(recunit,*,iostat=eof)recid,colrow,x,y,cmaqx,cmaqy,relev, + hill 5 if (eof .eq. 0) then c=int(real(colrow)/1000.0) !get column r=colrow-c*1000 !get row if ((col1 .le. c .and. c .le. col2 .and. row1 .le. r .and. + r .le. row2) .or. colrow .le. 0) then !if within the 5x5 block of cells c get distance and determine whether to write out call recdist('blk',recid,x,y,mindist,grddist,lwrite, + recflag,dist) if (lwrite) then write(98,20)recid,colrow,cmaqx,cmaqy,recflag, + x,y,relev,hill,dist nintb=nintb+1 endif endif read(recunit,*,iostat=eof)recid,colrow,x,y,cmaqx,cmaqy, + relev,hill goto 5 endif close(recunit) c monitors write(logunit,*)'reading monitors' if (utmzone .ge. 10) then write(recfil,'(2a,i2,a)')trim(adjustl(recdir)), + 'monitors_utm_',utmzone,'.dat' else write(recfil,'(2a,i1,a)')trim(adjustl(recdir)), + 'monitors_utm_',utmzone,'.dat' endif lexist=.false. inquire(file=recfil,exist=lexist) if (.not. lexist) then write(logunit,35)trim(adjustl(recfil)) ! stop endif if (lexist) then open(unit=recunit,file=recfil,status='old') read(recunit,*,iostat=eof)monid,colrow,x,y,cmaqx,cmaqy, + relev,hill 45 if (eof .eq. 0) then c=int(real(colrow)/1000.0) !get column r=colrow-c*1000 !get row if ((col1 .le. c .and. c .le. col2 .and. row1 .le. r + .and. r .le. row2) .or. colrow .le. 0) then !if within the 5x5 block of cells c get distance and determine whether to write out call recdist('mon',monid,x,y,mindist,grddist, + lwrite,recflag,dist) if (lwrite) then write(98,25)monid,colrow,cmaqx,cmaqy, + recflag,x,y,relev,hill,dist nintm=nintm+1 endif endif read(recunit,*,iostat=eof)monid,colrow,x,y,cmaqx, + cmaqy,relev,hill goto 45 endif endif close(recunit) close(98) c write(*,*)nintb,nintm c pause ! write(aerunit,40)nblock,ngrid,nmon if (nintb .eq. 0 .and. nintm .eq. 0) then write(logunit,50) else call int_recs2 endif else write(logunit,51) endif 20 format('**BLOCK= ',a,1x,'COLROW= ',i6,1x,'CMAQX= ',f11.0,1x, +'CMAQY= ',f11.0,' FLAG= ',i1,1x,f12.3,1x,f13.3, +2(1x,f9.3),1x,f11.5) 25 format('**MONITOR= ',a,1x,'COLROW= ',i6,1x,'CMAQX= ',f11.0,1x, +'CMAQY= ',f11.0,' FLAG= ',i1,1x,f12.3,1x,f13.3, +2(1x,f9.3),1x,f11.5) 35 format('receptor file ',a,' does not exist, stopping program') 50 format('WARNING NO RECEPTORS FOR SOURCE, SKIP INTERPOLATION!!!') 51 format('SOURCE IS IN AK, HI, PR, VI, SKIP INTERPOLATION') return end c########################################################################## subroutine int_recs2 use main1 implicit none character j1*10,j1a*10,recid*15,rt*10,tfil*250,ares*2, +arun*10,ahap*3,runtyp1*20 integer*8 eof,cr,iflag,b,m,findrec,ns real*8 x,y,ux,uy,relev,hill,dist c write(*,*)'enter int_recs2' b=0 m=0 c write(*,*)'you are here1' if (nintb .gt. 0) then allocate(intbc(ngrps1,nintb,nmonths)) allocate(intblock(nintb)) intblock%recid='NA' intblock%x=0.0 intblock%y=0.0 intblock%ux=0.0 intblock%uy=0.0 intblock%colrow=0 intblock%mastrec=0 intblock%dist=-9.0 intbc=0.0 endif if (nintm .gt. 0) then allocate(intmc(ngrps1,nintm,nmonths)) allocate(intmon(nintm)) intmon%recid='NA' intmon%x=0.0 intmon%y=0.0 intmon%ux=0.0 intmon%uy=0.0 intmon%colrow=0 intmon%mastrec=0 intmon%dist=-9.0 intmc=0.0 endif write(arun,'(i10)')runno !run # if (nhap1 .gt. 0) then write(ahap,'(i3.3)')nhap1 else ahap='001' endif if (lsingfac) then write(runtyp1,'(3(a))')trim(adjustl(runtyp)),'_', + trim(adjustl(singlfac)) elseif (lsingst) then write(runtyp1,'(3(a))')trim(adjustl(runtyp)),'_', + trim(adjustl(singlst)) else runtyp1=trim(adjustl(runtyp)) endif c write(*,*)'you are here2' c write(tfil,'(3(a))')'temp_recs_',trim(adjustl(singlhap)),'.txt' if (runcmv .or. runpt .or. runap) then write(tfil,'(8(a))')'temp_recs_',trim(adjustl(singlhap)), + '_',trim(adjustl(runtyp1)),trim(adjustl(arun)),'_', + trim(adjustl(ahap)),'.txt' else write(ares,'(i2)')gridres write(tfil,'(10(a))')'temp_recs_',trim(adjustl(singlhap)), + '_',trim(adjustl(runtyp1)),'_',trim(adjustl(ares)), + trim(adjustl(arun)),'_',trim(adjustl(ahap)),'.txt' endif open(unit=98,file=tfil,status='unknown') eof=0 read(98,*,iostat=eof)rt,recid,j1,cr,j1a,x,j1a,y,j1a,iflag,ux,uy, +relev,hill,dist 10 if (eof .eq. 0) then if (index(rt,'BLOCK') .gt. 0) then b=b+1 intblock(b)%recid=recid intblock(b)%x=x intblock(b)%y=y intblock(b)%ux=ux intblock(b)%uy=uy intblock(b)%colrow=cr intblock(b)%dist=dist c write(*,*)'b ',b intblock(b)%mastrec=findrec('BI',b) c write(111,*)intblock(b)%recid,intblock(b)%x,intblock(b)%y endif if (index(rt,'MONITOR') .gt. 0) then m=m+1 intmon(m)%recid=recid intmon(m)%x=x intmon(m)%y=y intmon(m)%ux=ux intmon(m)%uy=uy intmon(m)%colrow=cr intmon(m)%dist=dist c write(*,*)'m ',m intmon(m)%mastrec=findrec('MI',m) c write(112,*)intmon(m)%recid,intmon(m)%x,intmon(m)%y endif read(98,*,iostat=eof)rt,recid,j1,cr,j1a,x,j1a,y,j1a,iflag,ux, + uy,relev,hill,dist goto 10 endif close(98,status='delete') c close(98) c write(*,*)'exit int_recs2' return c 5 format(3(a),1x,i6,2(1x,a,f11.0),1x,a,1x,i1,1x,f12.3,1x,f13.3) end c########################################################################## c subroutine to interpolate to blocks/monitors subroutine interp(ns) use main1 implicit none c grecs: array indices of 4 closest gridded receptors integer*8 ir,ig,grecs(4),is,ns real*8 ic character ofile*250 write(logunit,*)'interpolating' call datetime if (writefac) then write(ofile,'(5(a))')trim(adjustl(hapdir)), + trim(adjustl(singlhap)),'_',trim(adjustl(facid(ns))), + '_interp.txt' open(unit=intunit,file=ofile,status='replace') endif c blocks write(logunit,*)'number of interpolated blocks ',nintb if (nintb .gt. 0) then do ir=1,nintb c find 4 closest gridded receptors call find4(intblock(ir)%x,intblock(ir)%y,grecs) c calculate terms of interpolation for each group do ig=1,ngrps1 do is=1,nmonths call calcint(intblock(ir)%x,intblock(ir)%y, + intblock(ir)%recid,ig,is,grecs,ic) intbc(ig,ir,is)=ic enddo enddo enddo endif c monitors if (nintm .gt. 0) then do ir=1,nintm c find 4 closest gridded receptors call find4(intmon(ir)%x,intmon(ir)%y,grecs) c calculate terms of interpolation for each group do ig=1,ngrps1 do is=1,nmonths call calcint(intmon(ir)%x,intmon(ir)%y, + intmon(ir)%recid,ig,is,grecs,ic) intmc(ig,ir,is)=ic enddo enddo enddo endif if (writefac) close(intunit) write(logunit,*)'interpolation done' call datetime return end c########################################################################## c subroutine to find the closest 4 gridded receptors to a block or monitor subroutine find4(x,y,grecs) use main1 implicit none c integer ir,grecs(4),gr,gr1,g,gg,grecs1(nmodg) integer*8 grecs(4),gr,gr1,g,gg,ng,ig c real*8 dx,dy,mindx,mindy,dist(nmodg),distance,min1,mind(4) real*8 dx,dy,distance,min1,mind(4),dist1,x,y real*8, allocatable, dimension(:) :: dist integer, allocatable, dimension(:) :: grecs1 grecs=0 c grecs1=0 c dist=0.0 c 4 points are denoted as below where x is receptor to be interpolated c 2 3 c x c 1 4 c point 1, dx <= 0 and dy <=0 c point 2, dx <=0 and dy >=0 c point 3, dx >=0 and dy >=0 c point 4, dx >=0 0 and dy <=0 c first get distance between block and all gridded receptors c only keep those within 10 km ng=0 do gr=1,nmodg dist1=distance(x,y,modgrid(gr)%x,modgrid(gr)%y) if (dist1 .le. 10000.0)ng=ng+1 c dist(gr)=distance(intblock(ir)%x,intblock(ir)%y, c + modgrid(gr)%x,modgrid(gr)%y) c grecs1(gr)=gr c write(121,*)gr,dist(gr),grecs1(gr),modgrid(gr)%x,modgrid(gr)%y enddo allocate(dist(ng)) allocate(grecs1(ng)) dist=0.0 grecs1=0.0 ig=0 do gr=1,nmodg dist1=distance(x,y,modgrid(gr)%x,modgrid(gr)%y) if (dist1 .le. 10000.0) then ig=ig+1 dist(ig)=dist1 grecs1(ig)=gr endif c dist(gr)=distance(intblock(ir)%x,intblock(ir)%y, c + modgrid(gr)%x,modgrid(gr)%y) c grecs1(gr)=gr c write(121,*)gr,dist(gr),grecs1(gr),modgrid(gr)%x,modgrid(gr)%y enddo c now sort based on increasing distance, i.e. first obs of new array is closest c do gr=1,nmodg-1 do gr=1,ng-1 min1=dist(gr) gg=grecs1(gr) g=gr c do gr1=gr+1,nmodg do gr1=gr+1,ng if (dist(gr1) .lt. min1) then min1=dist(gr1) gg=grecs1(gr1) g=gr1 endif enddo dist(g)=dist(gr) grecs1(g)=grecs1(gr) dist(gr)=min1 grecs1(gr)=gg enddo c do gr=1,nmodg c write(122,*)gr,dist(gr),grecs1(gr) c enddo c pick 4 points surrounding block c in some cases, there may not be 4 due to c being on coast or near 50 km. mind=100000. c do gr=1,nmodg !4 do gr=1,ng gg=grecs1(gr) dx=modgrid(gg)%x-x dy=modgrid(gg)%y-y if (dx .le. 0.0 .and. dy .le. 0.0) gr1=1 !point 1 if (dx .le. 0.0 .and. dy .ge. 0.0) gr1=2 !point 2 if (dx .ge. 0.0 .and. dy .ge. 0.0) gr1=3 !point 3 if (dx .ge. 0.0 .and. dy .le. 0.0) gr1=4 !point 4 if (dist(gr) .lt. mind(gr1)) then mind(gr1)=dist(gr) grecs(gr1)=gg endif c if (grecs(gr1) .eq. 0) grecs(gr1)=gg enddo deallocate(dist) deallocate(grecs1) return end c########################################################################## c subroutine to calculate terms of interpolation c based on HEM but linear instead of logarithmic subroutine calcint(rx,ry,rid,ig,is,grecs,termc) use main1 implicit none integer*8 ig,is,grecs(4),nzeros,n,g,gg real*8 terma,termb,termc,c(4),x(4),y(4),rx,ry,dx_4_1,dx_3_2, +calcterm character rid*15,sgrp*12 !sgrp*8 c=0.0 x=0.0 y=0.0 c if (rt .eq. 'blk') then c rx=intblock(ir)%x c ry=intblock(ir)%y c rid=intblock(ir)%recid c else c rx=intmon(ir)%x c ry=intmon(ir)%y c rid=intmon(ir)%recid c endif nzeros=0 do g=1,4 gg=grecs(g) if (gg .eq. 0) then nzeros=nzeros+1 else c(g)=modgc(ig,gg,is) x(g)=modgrid(gg)%x y(g)=modgrid(gg)%y endif enddo terma=0.0 termb=0.0 termc=0.0 if (nzeros .eq. 0) then ! all 4 gridded receptors are present c check dx and dy spacing between 4 terms dx_3_2=abs(x(3)-x(2)) dx_4_1=abs(x(4)-x(1)) if (dx_3_2 .ne. dx_4_1) then terma=calcterm(c(1),c(4),rx,x(1),x(4)) termb=calcterm(c(2),c(3),rx,x(2),x(3)) termc=calcterm(terma,termb,ry,y(1),y(2)) else terma=calcterm(c(1),c(2),ry,y(1),y(2)) termb=calcterm(c(4),c(3),ry,y(4),y(3)) termc=calcterm(terma,termb,rx,x(1),x(4)) endif c terma=c(1)+(c(2)-c(1))*(dy/dy1) c termb=c(4)+(c(3)-c(4))*(dy/dy1) c termc=terma+(termb-terma)*(dx/dx1) else !just take average n=0 do g=1,4 if (grecs(g) .ne. 0) then termc=termc+c(g) n=n+1 endif enddo termc=termc/real(n) endif if (writefac)then if (runap) then sgrp='ALL' elseif (runpt) then sgrp=ptsrcs(ig)%srcid elseif (runcmv) then ! sgrp='ALL' !fix later sgrp=cmvsrcs(ig)%srcid else sgrp=gridsrcs(ig)%srcid !fix later for gridded endif write(intunit,5)rid,sgrp,rx,ry,x(1),y(1),c(1),x(2),y(2), + c(2),x(3),y(3),c(3),x(4),y(4),c(4),terma,termb,termc endif c write(131,'(a,2(1x,f11.0),4(2(1x,f11.0),1x,e13.6),3(1x,e13.6))') c +rid,rx,ry,x(1),y(1),c(1),x(2),y(2),c(2),x(3),y(3), c +c(3),x(4),y(4),c(4),terma,termb,termc 5 format(a,1x,a,2(1x,f11.0),4(2(1x,f11.0),1x,e13.6),3(1x,e13.6)) return end c########################################################################## c calculate terms a and b of interpolation function calcterm(c1,c2,xr,x1,x2) use main1 implicit none real*8 c1,c2,xr,x1,x2,calcterm calcterm=c1+(c2-c1)*((xr-x1)/(x2-x1)) return end c########################################################################## c subroutine to calculate HAP-specific concentrations subroutine hapconc(ns) use main1 implicit none integer*8 ns,ihap,ilines,i,nlines1,ncas,r,e,ism,icat, +findsrc,findcat,irec,findhap,iskip,rt,eind c character cas*15,srcgrp*8,cats1*100,ofile*250,fips1*5 character cas*15,srcgrp*12,cats1*100,ofile*250,fips1*5,ofile1*250, +zipline*300 real*8 conc,hconc,monemis(12),sx,sy !emis1000(5), logical lexist integer*8, allocatable, dimension(:) :: prid integer*8, allocatable, dimension(:) :: uid integer*8, allocatable, dimension(:) :: rid real*8, allocatable, dimension(:,:) :: emis * integer*8, allocatable, dimension(:) :: hunit allocatable :: cas(:) allocatable :: cats1(:) allocatable :: srcgrp(:) allocatable :: fips1(:) c now using group specific seasonal emissions, version 17243 c emis1000(1)=8571.56256 !winter c emis1000(2)=8762.041728 !spring c emis1000(3)=8762.041728 !summer c emis1000(4)=8666.802144 !spring c emis1000(5)=34762.4872 !annual c write facility filename !13352 if (writefac) then write(ofile,'(3(a))')trim(adjustl(hapdir)), + trim(adjustl(facid(ns))),'_concs.txt' c hunit(ilines)=hapunit+ilines open(hapunit,file=ofile,status='replace') endif c write facility filename for HEID !13352 if (writheid) then write(ofile1,'(6(a))')trim(adjustl(hapdir)),'facilities/', + trim(adjustl(state2(ns))),'/',trim(adjustl(facid(ns))), + '_block_concs.txt' c hunit(ilines)=hapunit+ilines open(hapunit+1,file=ofile1,status='replace') endif c get entries in facility emissions array and write to local array if (runap) then eind=1 ilines=0 do i=1,nlines if (airport(i)%state .eq. state2(ns) .and. + airport(i)%facid .eq. facid(ns)) ilines=ilines+1 enddo nlines1=ilines ncas=nlines allocate(emis(nlines1,1)) allocate(cas(nlines1)) allocate(cats1(nlines1)) allocate(srcgrp(nlines1)) c allocate(hunit(nlines1)) allocate(fips1(nlines1)) fips1='' ilines=0 do i=1,nlines if (airport(i)%state .eq. state2(ns) .and. + airport(i)%facid .eq. facid(ns)) then ilines=ilines+1 emis(ilines,1)=airport(i)%emis(1) cas(ilines)=airport(i)%cas cats1(ilines)=airport(i)%cat srcgrp(ilines)=airport(i)%srcgrp c if (writefac) then c write(ofile,'(5(a))')trim(adjustl(hapdir)), c + trim(adjustl(cas(ilines))),'_', c + trim(adjustl(facid(ns))),'_concs.txt' c hunit(ilines)=hapunit+ilines c inquire(file=ofile,exist=lexist) c if (.not. lexist) then c open(hunit(ilines),file=ofile, c + status='replace') c endif c endif endif enddo c do i=1,nlines1 c write(136,*)i,srcgrp(i),emis(i,1),cas(i),cats1(i) c enddo elseif (runpt) then eind=1 ilines=0 do i=1,nlines if (facility(i)%state .eq. state2(ns) .and. + facility(i)%facid .eq. facid(ns)) ilines=ilines+1 enddo nlines1=ilines ncas=nlines allocate(emis(nlines1,1)) allocate(cas(nlines1)) allocate(cats1(nlines1)) allocate(srcgrp(nlines1)) c allocate(hunit(nlines1)) allocate(fips1(nlines1)) allocate(prid(nlines1)) allocate(uid(nlines1)) allocate(rid(nlines1)) fips1='' ilines=0 do i=1,nlines if (facility(i)%state .eq. state2(ns) .and. + facility(i)%facid .eq. facid(ns)) then ilines=ilines+1 emis(ilines,1)=facility(i)%emis(1) cats1(ilines)=facility(i)%cat srcgrp(ilines)=facility(i)%srcgrp cas(ilines)=facility(i)%cas ihap=findhap(cas(ilines)) prid(ilines)=facility(i)%procid uid(ilines)=facility(i)%unitid rid(ilines)=facility(i)%relid c if (writefac) then c write(ofile,'(5(a))')trim(adjustl(hapdir)), c + trim(adjustl(cas(ilines))),'_', c + trim(adjustl(facid(ns))),'_concs.txt' c hunit(ilines)=hapunit+ihap cc see if file already exists c inquire(file=ofile,exist=lexist) c if (.not. lexist) then c open(hunit(ilines),file=ofile, c + status='replace') c endif c endif endif enddo elseif (runcmv) then eind=1 ilines=0 do i=1,nlines if (cmv(i)%state .eq. state2(ns) .and. + cmv(i)%facid .eq. facid(ns)) ilines=ilines+1 enddo nlines1=ilines ncas=nlines allocate(emis(nlines1,1)) allocate(cas(nlines1)) allocate(cats1(nlines1)) allocate(srcgrp(nlines1)) c allocate(hunit(nlines1)) allocate(fips1(nlines1)) fips1='' ilines=0 do i=1,nlines if (cmv(i)%state .eq. state2(ns) .and. + cmv(i)%facid .eq. facid(ns)) then ilines=ilines+1 emis(ilines,1)=cmv(i)%emis(1) cats1(ilines)=cmv(i)%cat srcgrp(ilines)=cmv(i)%srcgrp cas(ilines)=cmv(i)%cas ihap=findhap(cas(ilines)) c if (writefac) then c write(ofile,'(5(a))')trim(adjustl(hapdir)), c + trim(adjustl(cas(ilines))),'_', c + trim(adjustl(facid(ns))),'_concs.txt' c hunit(ilines)=hapunit+ihap cc see if file already exists c inquire(file=ofile,exist=lexist) c if (.not. lexist) then c open(hunit(ilines),file=ofile, c + status='replace') c endif c endif endif enddo c elseif (noncontr) then c eind=5 c ilines=0 c do i=1,nlines c if (tracts(i)%facid .eq. facid(ns)) ilines=ilines+1 c enddo c nlines1=ilines c ncas=nlines c allocate(emis(nlines1,13)) c allocate(cas(nlines1)) c allocate(cats1(nlines1)) c allocate(srcgrp(nlines1)) c allocate(fips1(nlines1)) c allocate(hunit(nlines1)) c fips1='' c ilines=0 c do i=1,nlines c if (tracts(i)%facid .eq. facid(ns)) then c ilines=ilines+1 c do e=1,13 c emis(ilines,e)=tracts(i)%emis(e) c enddo c cas(ilines)=tracts(i)%cas c cats1(ilines)=tracts(i)%cat c srcgrp(ilines)=tracts(i)%srcgrp c ihap=findhap(cas(ilines)) cc if (writefac) then cc write(ofile,'(5(a))')trim(adjustl(hapdir)), cc + trim(adjustl(cas(ilines))),'_', cc + trim(adjustl(facid(ns))),'_concs.txt' cc hunit(ilines)=hapunit+ihap ccc see if file already exists cc inquire(file=ofile,exist=lexist) cc if (.not. lexist) then cc open(hunit(ilines),file=ofile, cc + status='replace') cc endif cc endif c endif c enddo else eind=13 ilines=0 do i=1,nlines if (gridded(i)%facid .eq. facid(ns)) ilines=ilines+1 enddo nlines1=ilines ncas=nlines allocate(emis(nlines1,13)) allocate(cas(nlines1)) allocate(cats1(nlines1)) allocate(srcgrp(nlines1)) allocate(fips1(nlines1)) c allocate(hunit(nlines1)) ilines=0 do i=1,nlines if (gridded(i)%facid .eq. facid(ns)) then ilines=ilines+1 do e=1,13 emis(ilines,e)=gridded(i)%emis(e) enddo cas(ilines)=gridded(i)%cas cats1(ilines)=gridded(i)%cat srcgrp(ilines)=gridded(i)%srcgrp fips1(ilines)=gridded(i)%fips ihap=findhap(cas(ilines)) c if (writefac) then c write(ofile,'(5(a))')trim(adjustl(hapdir)), c + trim(adjustl(cas(ilines))),'_', c + trim(adjustl(facid(ns))),'_concs.txt' c hunit(ilines)=hapunit+ihap cc see if file already exists c inquire(file=ofile,exist=lexist) c if (.not. lexist) then c open(hunit(ilines),file=ofile, c + status='replace') c endif c endif endif enddo endif c get facility specific, HAP specific concentrations do i=1,nlines1 monemis=0.0 if (runap )then ! .or. runcmv) then ism=1 else ism=findsrc(srcgrp(i)) endif if (runpt) then sx=ptsrcs(ism)%utmx sy=ptsrcs(ism)%utmy elseif (runcmv) then sx=cmvsrcs(ism)%utmx1 sy=cmvsrcs(ism)%utmy1 elseif (runap) then sx=apsrcs(ism)%utmx1 sy=apsrcs(ism)%utmy1 else call calc_center(ism,sx,sy) endif ihap=findhap(cas(i)) icat=findcat(cats1(i)) if (nmonths .eq. 12) then do e=1,12 monemis(e)=emis(i,e) enddo endif c write(*,*)'s ',i, srcgrp(i),nlines1 c gridded modeled receptors if (nmodg .gt. 0) then do r=1,nmodg c check gridded sources to see if gridded receptor is > 50 km from that particular source, if so don't do calculations if (.not. runap .and. .not. runpt .and. .not. + runcmv) then call check_rec_grid(modgrid(r)%ux,modgrid(r)%uy, + srcgrp(i),r,iskip) else if (modgrid(r)%iflag .eq. 3) then iskip=1 else iskip=0 endif endif if (iskip .eq. 0) then if (nmonths .eq. 1) then conc=modgc(ism,r,1) hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5) c write(121,*)modgc(ism,r,1),conc,emis(i,eind),hconc c write(137,*)r,conc,hconc,emis(i,1),emis1000(5),modgrid(r)%recid else call calcann('GM',hconc,ism,monemis,r) c + emis1000,r) !calculate annual averages from seasonal averages endif rt=10 if (writefac) then if (runpt) then write(hapunit,15)modgrid(r)%recid, + rt,srcgrp(i),sx,sy,prid(i),uid(i),rid(i), + cas(i),cats1(i),hconc, + modgrid(r)%x,modgrid(r)%y,modgrid(r)%ux, + modgrid(r)%uy,modgrid(r)%dist,fips1(i) else write(hapunit,5)modgrid(r)%recid, + rt,srcgrp(i),sx,sy,cas(i),cats1(i),hconc, + modgrid(r)%x,modgrid(r)%y,modgrid(r)%ux, + modgrid(r)%uy,modgrid(r)%dist,fips1(i) endif endif irec=modgrid(r)%mastrec gridconc(ihap,irec,icat)=gridconc(ihap,irec,icat)+ + hconc endif c write(*,*)r,irec,gridrec(irec)%recid,gridconc(ihap,irec,icat), c +modgrid(r)%recid c pause enddo endif c block modeled receptors if (nmodb .gt. 0) then do r=1,nmodb if (modblock(r)%iflag .ne. 3) then if (nmonths .eq. 1) then conc=modbc(ism,r,1) hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5) else call calcann('BM',hconc,ism,monemis,r) c + emis1000,r) !calculate annual averages from seasonal averages endif rt=20 if (writefac) then if (runpt) then write(hapunit,15)modblock(r)%recid, + rt,srcgrp(i),sx,sy,prid(i),uid(i),rid(i), + cas(i),cats1(i),hconc, + modblock(r)%x,modblock(r)%y, + modblock(r)%ux,modblock(r)%uy, + modblock(r)%dist,fips1(i) if (writheid) write(hapunit+1,20) + modblock(r)%recid,srcgrp(i),sx,sy,prid(i), + uid(i),rid(i),cas(i), + hconc,modblock(r)%ux,modblock(r)%uy else write(hapunit,5)modblock(r)%recid, + rt,srcgrp(i),sx,sy,cas(i),cats1(i),hconc, + modblock(r)%x,modblock(r)%y, + modblock(r)%ux,modblock(r)%uy, + modblock(r)%dist,fips1(i) if (writheid) write(hapunit+1,10) + modblock(r)%recid,srcgrp(i),sx,sy,cas(i), + hconc,modblock(r)%ux,modblock(r)%uy endif endif irec=modblock(r)%mastrec blokconc(ihap,irec,icat)=blokconc(ihap,irec,icat)+ + hconc endif enddo endif c block interpolated receptors if (nintb .gt. 0) then do r=1,nintb if (.not. runap .and. .not. runpt .and. .not. + runcmv) then call check_rec_grid(intblock(r)%ux,intblock(r)%uy, + srcgrp(i),r,iskip) else iskip=0 endif if (iskip .eq. 0) then if (nmonths .eq. 1) then conc=intbc(ism,r,1) hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5) else call calcann('BI',hconc,ism,monemis,r) c + emis1000,r) !calculate annual averages from seasonal averages endif rt=21 if (writefac) then if (runpt) then write(hapunit,15)intblock(r)%recid, + rt,srcgrp(i),sx,sy,prid(i),uid(i),rid(i), + cas(i),cats1(i),hconc, + intblock(r)%x,intblock(r)%y, + intblock(r)%ux,intblock(r)%uy, + intblock(r)%dist,fips1(i) else write(hapunit,5)intblock(r)%recid, + rt,srcgrp(i),sx,sy,cas(i),cats1(i),hconc, + intblock(r)%x,intblock(r)%y, + intblock(r)%ux,intblock(r)%uy, + intblock(r)%dist,fips1(i) endif endif irec=intblock(r)%mastrec blokconc(ihap,irec,icat)=blokconc(ihap,irec,icat)+ + hconc endif enddo endif c non-populated block modeled receptors if (nmodz .gt. 0) then do r=1,nmodz if (zmodblock(r)%iflag .ne. 3) then if (nmonths .eq. 1) then conc=modzbc(ism,r,1) hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5) else call calcann('ZB',hconc,ism,monemis,r) c + emis1000,r) !calculate annual averages from seasonal averages endif rt=22 if (writefac) then if (runpt) then write(hapunit,15)zmodblock(r)%recid, + rt,srcgrp(i),sx,sy,prid(i),uid(i),rid(i), + cas(i),cats1(i),hconc, + zmodblock(r)%x,zmodblock(r)%y, + zmodblock(r)%ux,zmodblock(r)%uy, + zmodblock(r)%dist,fips1(i) else write(hapunit,5)zmodblock(r)%recid, + rt,srcgrp(i),sx,sy,cas(i),cats1(i),hconc, + zmodblock(r)%x,zmodblock(r)%y, + zmodblock(r)%ux,zmodblock(r)%uy, + zmodblock(r)%dist,fips1(i) endif endif irec=zmodblock(r)%mastrec zblokconc(ihap,irec,icat)= + zblokconc(ihap,irec,icat)+hconc endif enddo endif c monitor modeled receptors if (nmodm .gt. 0) then do r=1,nmodm if (modmon(r)%iflag .ne. 3) then if (nmonths .eq. 1) then conc=modmc(ism,r,1) hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5) else call calcann('MM',hconc,ism,monemis,r) c + emis1000,r) !calculate annual averages from seasonal averages endif rt=30 if (writefac) then if (runpt) then write(hapunit,15)modmon(r)%recid,rt, + srcgrp(i),sx,sy,prid(i),uid(i),rid(i), + cas(i),cats1(i),hconc, + modmon(r)%x,modmon(r)%y,modmon(r)%ux, + modmon(r)%uy,modmon(r)%dist,fips1(i) else write(hapunit,5)modmon(r)%recid,rt, + srcgrp(i),sx,sy,cas(i),cats1(i),hconc, + modmon(r)%x,modmon(r)%y,modmon(r)%ux, + modmon(r)%uy,modmon(r)%dist,fips1(i) endif endif irec=modmon(r)%mastrec monconc(ihap,irec,icat)=monconc(ihap,irec,icat)+ + hconc endif enddo endif c monitor interpoloated receptors if (nintm .gt. 0) then do r=1,nintm if (.not. runap .and. .not. runpt .and. .not. + runcmv) then call check_rec_grid(intmon(r)%ux,intmon(r)%uy, + srcgrp(i),r,iskip) !check to see if outside 50 km for gridded source else iskip=0 endif if (iskip .eq. 0) then if (nmonths .eq. 1) then conc=intmc(ism,r,1) hconc=conc*emis(i,eind)/memis(ism,13) !emis1000(5) else call calcann('MI',hconc,ism,monemis,r) c + emis1000,r) !calculate annual averages from seasonal averages endif rt=31 if (writefac) then if (runpt) then write(hapunit,15)intmon(r)%recid,rt, + srcgrp(i),sx,sy,prid(i),uid(i),rid(i), + cas(i),cats1(i),hconc, + intmon(r)%x,intmon(r)%y,intmon(r)%ux, + intmon(r)%uy,intmon(r)%dist,fips1(i) else write(hapunit,5)intmon(r)%recid,rt, + srcgrp(i),sx,sy,cas(i),cats1(i),hconc, + intmon(r)%x,intmon(r)%y,intmon(r)%ux, + intmon(r)%uy,intmon(r)%dist,fips1(i) endif endif irec=intmon(r)%mastrec monconc(ihap,irec,icat)=monconc(ihap,irec,icat)+ + hconc endif enddo endif enddo c do i=1,nlines1 c close(hunit(i)) c enddo if (writefac) then close(hapunit) write(zipline,'(a,1x,a)')'gzip -f',trim(adjustl(ofile)) call system(zipline) endif if (writheid)close(hapunit+1) c 5 format(a15,1x,i2,1x,a8,1x,a15,e13.6,2(1x,f11.0),2(1x,f13.3),1x, c +f11.5) c 5 format(a15,1x,i2,1x,a8,1x,a45,e13.6,2(1x,f11.0),2(1x,f13.3),1x, c +f11.5,1x,a5) 5 format(a15,1x,i2,1x,a8,2(1x,f13.3),1x,a15,1x,a45,e13.6,2(1x,f11.0) +,2(1x,f13.3),1x,f11.5,1x,a5) 10 format(a15,1x,a8,2(1x,f13.3),1x,a15,e13.6,2(1x,f13.3)) 15 format(a15,1x,i2,1x,a8,2(1x,f13.3),3(1x,i15),1x,a15,1x,a45,e13.6, +2(1x,f11.0),2(1x,f13.3),1x,f11.5,1x,a5) 20 format(a15,1x,a8,2(1x,f13.3),3(1x,i15),1x,a15,e13.6,2(1x,f13.3)) c deallocate arrays deallocate(emis) deallocate(cas) deallocate(cats1) deallocate(srcgrp) deallocate(fips1) if (allocated(prid)) deallocate(prid) if (allocated(uid)) deallocate(uid) if (allocated(rid)) deallocate(rid) return end c########################################################################## c subroutine to determine distance from receptor to source c and determine if receptor within 30 m of source subroutine recdist(rt,rid,x,y,mindist,maxdist,lwrite,iflag,md) use main1 implicit none c iflag: flag denoting if source within 30 m of source c igrp: source loop counter integer*8 iflag,igrp,i c x: receptor UTM-x coordinate c y: receptor UTM-y coordinate c dist: distance from receptor to source c maxdist: maximum distance to search c apdist: distance from receptor to airport c ptdist: distance from receptor to point source real*8 x,y,dist,mindist,maxdist,apdist,ptdist,md,avgx,avgy, +distance c lwrite: logical variable denoting to write receptor to AERMOD.INP logical lwrite,ldup character rt*3,rid*15 iflag=0 lwrite=.false. md=maxdist+10000.0 c if (trim(adjustl(rt)) .eq. 'mon') pause if (runap) then do igrp=1,ngrps dist=apdist(x,y,igrp) !get distance from receptor to points along runways if (dist .gt. mindist .and. dist .le. maxdist) then lwrite=.true. if (dist .lt. md) md=dist if (dist .le. 30.0) iflag=1 endif enddo elseif (runpt) then do igrp=1,ngrps if (ptsrcs(igrp)%srctyp .eq. 'AREA') then dist=apdist(x,y,igrp) !get distance from receptor to points along area source else dist=ptdist(x,y,igrp) !get distance from receptor to stacks endif if (dist .gt. mindist .and. dist .le. maxdist) then lwrite=.true. if (dist .lt. md) md=dist if (dist .le. 30.0) iflag=1 endif enddo elseif (runcmv) then ! .or. noncontr) then do igrp=1,ngrps call cmvcheck(x,y,igrp,dist) !distance from receptor to CMV polygons if (dist .gt. mindist .and. dist .le. maxdist) then lwrite=.true. if (dist .lt. md) md=dist if (dist .le. 30.0) iflag=1 else call polycenter(igrp,avgx,avgy) dist=distance(x,y,avgx,avgy) if (dist .le. maxdist) then lwrite=.true. if (dist .le. 30.0) iflag=1 if (dist .lt. md) md=dist endif endif enddo else call recs4grid(x,y,grddist,md,lwrite) !distance from receptor to center of gridded sources endif c if receptor is to be interpolated, make sure not a modeled receptor when near 10 km c if (lwrite .and. dist .le. 11000.0) then c if (lwrite .and. md .le. 11000.0) then if (lwrite) then !check all distances i=1 ldup=.false. if (trim(adjustl(rt)) .eq. 'blk') then if (nmodb .gt. 0) then ldup=.false. do while (i .le. nmodb .and. .not. ldup) if (trim(adjustl(rid)) .eq. modblock(i)%recid + .and. modblock(i)%iflag .ne. 3) + ldup=.true. i=i+1 enddo endif endif c if (trim(adjustl(rt)) .eq. 'zbl') then c if (nmodz .gt. 0) then c ldup=.false. c do while (i .le. nmodz .and. .not. ldup) c if (trim(adjustl(rid)) .eq. zmodblock(i)%recid c + .and. zmodblock(i)%iflag .ne. 3) c + ldup=.true. c i=i+1 c enddo c endif c endif if (trim(adjustl(rt)) .eq. 'mon') then if (nmodm .gt. 0) then do while (i .le. nmodm .and. .not. ldup) if (trim(adjustl(rid)) .eq. modmon(i)%recid .and. + modmon(i)%iflag .ne. 3) + ldup=.true. i=i+1 enddo endif endif if (ldup)lwrite=.false. endif return end c########################################################################## c subroutine to calculate polygon centroid subroutine polycenter(ig,avgx,avgy) use main1 implicit none integer*8 ig,j,nv1 real*8 avgx,avgy,area if (runcmv) then nv1=cmvsrcs(ig)%nvert else nv1=tractsrcs(ig)%nvert endif avgx=0.0 avgy=0.0 area=0.0 do j=1,nv1 if (j .lt. nv1) then area=area+(polyx(ig,j)*polyy(ig,j+1)-polyx(ig,j+1)* + polyy(ig,j)) avgx=(polyx(ig,j)+polyx(ig,j+1))*(polyx(ig,j)* + polyy(ig,j+1)-polyx(ig,j+1)*polyy(ig,j))+avgx avgy=(polyy(ig,j)+polyy(ig,j+1))*(polyx(ig,j)* + polyy(ig,j+1)-polyx(ig,j+1)*polyy(ig,j))+avgy else area=area+(polyx(ig,j)*polyy(ig,1)-polyx(ig,1)* + polyy(ig,j)) avgx=(polyx(ig,j)+polyx(ig,1))*(polyx(ig,j)*polyy(ig,1) + -polyx(ig,1)*polyy(ig,j))+avgx avgy=(polyy(ig,j)+polyy(ig,1))*(polyx(ig,j)*polyy(ig,1)- + polyx(ig,1)*polyy(ig,j))+avgy endif enddo avgx=avgx/(3.0*area) avgy=avgy/(3.0*area) return end c########################################################################## c calculate distance between receptor and point source function ptdist(x,y,ig) use main1 implicit none c ig: source index counter integer*8 ig c x: receptor UTM-x coordinate c y: receptor UTM-y coordinate c fx: UTM-x coordinate of stack c fy: UTM-y coordinate of stack c distance: distance function c ptdist: distance from receptor to stack real*8 x,y,fx,fy,distance,ptdist fx=ptsrcs(ig)%utmx fy=ptsrcs(ig)%utmy ptdist=distance(x,y,fx,fy) return end c########################################################################## subroutine cmvcheck(x,y,ig,mindist) use main1 implicit none c ig: source index counter c eof: end of file indicator c nv: number of vertices for polygon (read from file) c nv1: number of vertices for polygon being processed c j: vertex counter integer*8 ig,eof,nv,nv1,j,i,inout1,inout c area1: area of polygon c cmvdist: distance of receptor from polygon edge function c dist: distance of receptor from polygon edge c mindist: minimum distance of receptor from all edges of all modeled polygons c frac1: fraction polygon is of total port c h: release height c sz: sigma-z c ux: array of UTM x-coordinates of vertices c uy: array of UTM y-coordinates of vertices real*8 x,y,cmvdist,dist,mindist,xi,yi,xj,yj real*8, allocatable, dimension(:) :: ux real*8, allocatable, dimension(:) :: uy logical ix,iy,jx,jy,L_EOR eof=0 mindist=900000.0 if (runcmv) then nv=cmvsrcs(ig)%nvert+1 !+1 added 17243 else nv=tractsrcs(ig)%nvert+1 !+1 added 17243 endif c version 17234, first check to see if receptor inside the polygon before checking the distances allocate(ux(nv)) allocate(uy(nv)) do j=1,nv-1 ux(j)=polyx(ig,j) uy(j)=polyy(ig,j) enddo ux(nv)=polyx(ig,1) uy(nv)=polyy(ig,1) INOUT1=-1 DO I=1,NV c write(*,*)i XI=UX(I)-X YI=UY(I)-Y C CHECK WHETHER THE POINT IN QUESTION IS AT THIS VERTEX. IF (XI.EQ. 0.0 .AND. YI.EQ. 0.0) THEN INOUT1 = 0 EXIT END IF C J IS NEXT VERTEX NUMBER OF POLYGON. J=1+MOD(I,NV) XJ=UX(J)-X YJ=UY(J)-Y C IS THIS LINE OF 0 LENGTH ? IF (XI.EQ.XJ .AND. YI.EQ.YJ) CYCLE IX=XI.GE.0.0D0 IY=YI.GE.0.0D0 JX=XJ.GE.0.0D0 JY=YJ.GE.0.0D0 C CHECK WHETHER (PX,PY) IS ON VERTICAL SIDE OF POLYGON. L_EOR = EOR(IY,JY) IF (XI.EQ. 0.0 .AND. XJ.EQ. 0.0 .AND. L_EOR) THEN INOUT1 = 0 EXIT END IF C CHECK WHETHER (PX,PY) IS ON HORIZONTAL SIDE OF POLYGON. L_EOR = EOR(IX,JX) IF (YI.EQ. 0.0 .AND. YJ.EQ. 0.0 .AND. L_EOR) THEN INOUT = 0 EXIT END IF C CHECK WHETHER BOTH ENDS OF THIS SIDE ARE COMPLETELY 1) TO RIGHT C OF, 2) TO LEFT OF, OR 3) BELOW (PX,PY). L_EOR = EOR(IX,JX) IF (.NOT.((IY.OR.JY).AND.L_EOR)) then CYCLE endif C DOES THIS SIDE OBVIOUSLY CROSS LINE RISING VERTICALLY FROM (PX,PY) L_EOR = EOR(IX,JX) IF (.NOT.(IY.AND.JY.AND.L_EOR)) THEN IF ((YI*XJ-XI*YJ)/(XJ-XI) .LT. 0.0) THEN CYCLE ELSE IF ((YI*XJ-XI*YJ)/(XJ-XI) .EQ. 0.0) THEN INOUT1 = 0 EXIT ELSE INOUT1 = -INOUT1 END IF ELSE INOUT1 = -INOUT1 END IF END DO if (inout1 .ge. 0) then mindist=0.0 else do j=1,nv-1 c calculate distance of receptor to each edge (shortest distance from receptor to anywhere along edge) c dist=cmvdist(x,y,polyx(ig,j),polyy(ig,j),polyx(ig,j+1), c + polyy(ig,j+1)) dist=cmvdist(x,y,ux(j),uy(j),ux(j+1),uy(j+1)) c get minimum distance of receptor to any polygon if (dist .le. mindist) mindist=dist enddo endif deallocate(ux) deallocate(uy) c rewind(parunit) return CONTAINS LOGICAL FUNCTION EOR(IX,IY) LOGICAL IX, IY,a,b a=(IX.OR.IY) b=.NOT.(IX.AND.IY) EOR = (IX.OR.IY) .AND. .NOT.(IX.AND.IY) END FUNCTION end c########################################################################## function cmvdist(x,y,fx1,fy1,fx2,fy2) c calculate shortest distance between receptor and polygon edge c anglea is angle with vertex at x1,y1 and endpoints at x2,y2 and receptor x,y c a is distance between receptor and x2,y2 c angleb is angle with vertex at x2,y2 and endpoints at x1,y1 and receptor x,y c b is distance between receptor and x1,y1 c anglec is angle with vertex at receptor x,y and endpoints at x1,y1 and x2,y2 c c is distance between x1,y1 and x2,y2 c use law of cosines use main1 implicit none c x: receptor UTM-x coordinate c y: receptor UTM-y coordinate c fx1: UTM-x coordinate of vertex 1 c fy1: UTM-y coordinate of vertex 1 c fx2: UTM-x coordinate of vertex 2 c fy2: UTM-y coordinate of vertex 2 c distance: distance function c a: distance between receptor and fx2,fy2 c b: distance between receptor and fx1,fy1 c c: distance between fx1,fy1 and fx2,fy2 c anglea: angle with vertex at fx1,fy1 and endpoints at fx2,fy2 and receptor x,y c angleb: angle with vertex at fx2f,y2 and endpoints at fx1,fy1 and receptor x,y c anglec: angle with vertex at receptor x,y and endpoints at fx1,fy1 and fx2,fy2 c cmvdist: distance between receptor and edge c angle: angle calculation function based on law of cosines real*8 x,y,fx1,fy1,fx2,fy2,distance,a,b,c,anglea,angleb,anglec, +cmvdist,angle c calculate distances a, b, and c and anglea, angleb, and anglec a=distance(x,y,fx2,fy2) b=distance(x,y,fx1,fy1) c=distance(fx1,fy1,fx2,fy2) c calculate angles anglea=angle(b,c,a) angleb=angle(a,c,b) anglec=angle(a,b,c) c check angles to determine next step if (anglea .le. 90.0 .and. angleb .le. 90.0) then !receptor is alongside edge c calculate distance from receptor to centerline cmvdist=b*sin(anglea*d2r) else cmvdist=min(a,b) endif c write(104,*)x,y,fx1,fy1,fx2,fy2,a,b,c,anglea,angleb,anglec,cmvdist return end c########################################################################## function apdist(x,y,ig) c calculate shortest distance between receptor and airport runway edge c or rectangular area source c anglea is angle with vertex at x1,y1 and endpoints at x2,y2 and receptor x,y c a is distance between receptor and x2,y2 c angleb is angle with vertex at x2,y2 and endpoints at x1,y1 and receptor x,y c b is distance between receptor and x1,y1 c anglec is angle with vertex at receptor x,y and endpoints at x1,y1 and x2,y2 c c is distance between x1,y1 and x2,y2 c use law of cosines use main1 implicit none c ig: source index counter integer*8 ig c x: receptor UTM-x coordinate c y: receptor UTM-y coordinate c fx1: UTM-x coordinate of endpoint 1 c fy1: UTM-y coordinate of endpoint 1 c fx2: UTM-x coordinate of endpoint 2 c fy2: UTM-y coordinate of endpoint 2 c distance: distance function c a: distance between receptor and fx2,fy2 c b: distance between receptor and fx1,fy1 c c: distance between fx1,fy1 and fx2,fy2 c anglea: angle with vertex at fx1,fy1 and endpoints at fx2,fy2 and receptor x,y c angleb: angle with vertex at fx2f,y2 and endpoints at fx1,fy1 and receptor x,y c anglec: angle with vertex at receptor x,y and endpoints at fx1,fy1 and fx2,fy2 c cmvdist: distance between receptor and edge c angle: angle calculation function based on law of cosines c cprime: 1/2 width of runway for runways, length of area source side for c other area sources c d1: distance from receptor to centerline of runway or edge of area source c angled: obtuse angle (anglea or angleb) -90 degrees c d: distance to corner of runway or area source based on angle D c bprime: distance of distance a or b within runway or area source c d2: portion of b or a outside runway c d3: distance based on right triangle c ddprime: opposite of angle A based on distance c being the adjacent side of angle for right triangle c dx1: = a if angleb > 90, b if anglea > 90 c dx2: =a if anglea > 90, b if angleb > 90 c ang1: =anglea if anglea > 90, angleb if angleb >90 c ang2: =angleb if anglea > 90, anglea if angleb > 90 c apdist: minimum of d,d2,dx1,or d3 c d3a: distance along runway for distance d3 real*8 x,y,fx1,fy1,fx2,fy2,a,b,c,anglea,angleb,anglec,bprime,d1, +d,d2,d3,cprime,distance,angle,apdist,ang1,ang2,dx1,dx2,angled, +ddprime,d3a,aprime,bprime1 apdist=-999.0 if (runap) then fx1=apsrcs(ig)%utmx1 fy1=apsrcs(ig)%utmy1 c if (runtyp .eq. 'APRUN') then if (apsrcs(ig)%utmx2 .ge. 0.0) then !18031 fx2=apsrcs(ig)%utmx2 fy2=apsrcs(ig)%utmy2 cprime=apsrcs(ig)%xdim/2.0 else fx2=apsrcs(ig)%utmx1+apsrcs(ig)%xdim/2.0 fy2=apsrcs(ig)%utmy1+apsrcs(ig)%ydim cprime=apsrcs(ig)%xdim/2.0 endif else cprime=ptsrcs(ig)%xdim/2.0 a=sqrt((cprime*cprime)+(ptsrcs(ig)%ydim*ptsrcs(ig)%ydim)) b=ptsrcs(ig)%ydim angleb=angle(a,cprime,b) anglec=angle(a,b,cprime) aprime=3.141592654-ptsrcs(ig)%angle-angleb-anglec !remember angles are in radians, 3.14159 is 180 fx1=ptsrcs(ig)%utmx+(cprime*sin(aprime)) fy1=ptsrcs(ig)%utmy-(cprime*cos(aprime)) bprime1=(3.141592654/2.0)-ptsrcs(ig)%angle-anglec !remember angles are in radians, 3.14159/2 is 90 fx2=ptsrcs(ig)%utmx+(a*cos(bprime1)) fy2=ptsrcs(ig)%utmy+(a*sin(bprime1)) c fx1=ptsrcs(ig)%utmx+ptsrcs(ig)%xdim/2.0 c fy1=ptsrcs(ig)%utmy c fx2=ptsrcs(ig)%utmx+ptsrcs(ig)%xdim/2.0 c fy2=ptsrcs(ig)%utmy+ptsrcs(ig)%ydim c cprime=ptsrcs(ig)%xdim/2.0 endif ! write(*,*)'a ',fx1,fy1,fx2,fy2,cprime,x,y ! pause c calculate distances a=distance(x,y,fx2,fy2) b=distance(x,y,fx1,fy1) c=distance(fx1,fy1,fx2,fy2) c calculate angles anglea=angle(b,c,a) angleb=angle(a,c,b) anglec=angle(a,b,c) ! write(*,*)'b ',a,b,c,anglea,angleb,anglec ! pause c check angles to determine next step if (anglea .le. 90.0 .and. angleb .le. 90.0) then !receptor is alongside runway or on runway c calculate distance from receptor to centerline d1=b*sin(anglea*d2r) apdist=d1-cprime !distance to edge ! write(*,*)'c ',d1 ! pause else if (anglea .gt. 90) then ang1=anglea dx1=b dx2=a ang2=angleb else ang1=angleb dx1=a dx2=b ang2=anglea endif angled=-90.0+ang1 !90-180+ang1 d=sqrt((dx1*dx1)+(cprime*cprime)- + (2*dx1*cprime*cos(angled*d2r))) !distance to corner of runway based on angle D bprime=c/cos(ang2*d2r) ddprime=c*tan(ang2*d2r) !opposite of angle A based on distance c being the adjacent side of angle for right triangle ! write(*,*)'d ',d,ang1,dx1,dx2,ang2,angled,bprime,ddprime ! pause if (ddprime .le. cprime) then d2=dx2-bprime d3=dx1*sin(angled*d2r) !distance based on right triangle with dx1 as the hypotenuse d3a=dx1*cos(angled*d2r) !distance along runway edge to d3 if (d3a .gt. cprime) then !d3a and associated d3 is past runway apdist=min(dx2,d2,dx1) else apdist=min(d,d2,dx1,d3) endif ! write(*,*)'e ',d2,d3,apdist ! pause else apdist=min(d,dx1) ! write(*,*)'f ',dx1,apdist ! pause endif endif return end c########################################################################## function angle(d1,d2,d3) c calculate angles based on law of cosines use main1 implicit none c angle: angle returned c d1: distance 1 c d2: distance 2 c d3: distance 3 (corresponds to opposite side of angle being determined) c xx: ((d1*d1)+(d2*d2)-(d3*d3))/(2*d1*d2) real*8 angle,d1,d2,d3,xx xx=((d1*d1)+(d2*d2)-(d3*d3))/(2*d1*d2) angle=acos(xx)/d2r return end c########################################################################## subroutine recs4grid(x,y,griddist,dist,lwrite) use main1 implicit none c ig: source index counter integer*8 ig c x: receptor UTM-x coordinate c y: receptor UTM-y coordinate c cenx: UTM-x coordinate of grid cell center c cenx: UTM-y coordinate of grid cell center c distance: distance function c dist: distance from receptor to center coordinates c griddist: max distance to model real*8 x,y,dist,distance,cenx,ceny,griddist c logical variable denoting to write receptor to AERMOD.INP logical lwrite lwrite=.false. c only need to write the receptor once if within 50 km of multiple sources c if only 1 source then no problem ig=1 do while(ig .le. ngrps .and. .not. lwrite) call calc_center(ig,cenx,ceny) dist=distance(x,y,cenx,ceny) if (dist .le. griddist) lwrite=.true. ig=ig+1 enddo return end c########################################################################## subroutine calc_center(ig,cx,cy) use main1 implicit none c ig: source index counter integer*8 ig c cx: UTM-x coordinate of grid cell center c cx: UTM-y coordinate of grid cell center real*8 cx,cy cc cx=0.0 cc cy=0.0 cx=(gridsrcs(ig)%utmx1+gridsrcs(ig)%utmx2+gridsrcs(ig)%utmx3+ +gridsrcs(ig)%utmx4)/(real(gridsrcs(ig)%nvert)) cy=(gridsrcs(ig)%utmy1+gridsrcs(ig)%utmy2+gridsrcs(ig)%utmy3+ +gridsrcs(ig)%utmy4)/(real(gridsrcs(ig)%nvert)) return end c########################################################################## function distance(x,y,fx,fy) use main1 implicit none c x: x-coordinate of point 1 c y: y-coordinate of point 1 c fx: x-coordinate of point 2 c fy: y-coordinate of point 2 c distance: distance between point 1 and point 2 c dx: x-distance between point 1 and point 2 c dy: y-distance between point 1 and point 2 real*8 x,y,fx,fy,distance,dx,dy dx=x-fx dy=y-fy distance=sqrt(dx*dx+dy*dy) return end c########################################################################## subroutine upcase(flag) c convert lowercase variable to upper case variable c flag: variable to checked c lower: string of lowercase letters c upper: string of uppercase letters c i: index of flag in lower integer*8 i character flag*1 character lower*26,upper*26 i=0 lower='abcdefghijklmnopqrstuvwxyz' upper='ABCDEFGHIJKLMNOPQRSTUVWXYZ' i=index(lower,flag) if (i .gt. 0) flag=upper(i:i) return end c########################################################################## subroutine datetime c subroutine to write starting and ending times to screen and log file c subroutine calls the intrinsic FORTRAN subroutine date_and_time c inputs are iflag which tells the subroutine if the program starts (iflag=1) c or ends (iflag not equal to 1) c variables c idattim: array of date/time returned by date_and_time with 8 elements c elements are: c (1) 4-digit year c (2) month of year (1 to 12) c (3) day of month c (4) time difference with respect to GMT in minutes c (5) hour of day (0 to 23) in local time c (6) minute of time c (7) second of time c (8) milliseconds of time c hr: integer hour of day c cdate: character variable of date c ctime: character variable of time c czone: character variable of time zone from GMT c ampm: character variable with value of 'AM' or 'PM' use main1 implicit none integer*8 idattim(8),hr character cdate*8,ctime*10,czone*5,ampm*2,amonths(12)*10 data amonths /'January','February','March','April','May','June', & 'July','August','September','October','November','December'/ call date_and_time(cdate,ctime,czone,idattim) if (idattim(5) .le. 11) then hr=idattim(5) ampm='AM' else hr=12+(idattim(5)-12) ampm='PM' endif write(logunit,10)trim(adjustl(amonths(idattim(2)))),idattim(3), + idattim(1),hr,idattim(6),idattim(7),ampm 10 format(1x,a,1x,i2.2,', ',i4,2x,i2,':', + i2.2,':',i2.2,1x,a2) return end c########################################################################## function findsrc(src) use main1 implicit none c findsrc: source index in source array integer*8 findsrc c src: AERMOD source id (8 characters) c src12: AERMOD source ID (12 characters) character src*12 !,src12*12 c lfound: logical variable denoting source found in array logical lfound c src12=src findsrc=1 lfound=.false. if (runpt) then do while(findsrc .le. ngrps .and. .not. lfound) if (trim(adjustl(src)) .eq. + trim(adjustl(ptsrcs(findsrc)%srcid)))then lfound=.true. else findsrc=findsrc+1 endif enddo elseif (runap) then do while(findsrc .le. ngrps .and. .not. lfound) if (trim(adjustl(src)) .eq. + trim(adjustl(apsrcs(findsrc)%srcid)))then lfound=.true. else findsrc=findsrc+1 endif enddo elseif (runcmv) then do while(findsrc .le. ngrps .and. .not. lfound) if (trim(adjustl(src)) .eq. cmvsrcs(findsrc)%srcid)then lfound=.true. else findsrc=findsrc+1 endif enddo c elseif (noncontr) then c do while(findsrc .le. ngrps .and. .not. lfound) c if (trim(adjustl(src)) .eq. tractsrcs(findsrc)%srcid)then c lfound=.true. c else c findsrc=findsrc+1 c endif c enddo else do while(findsrc .le. ngrps .and. .not. lfound) c write(*,*)src,gridsrcs(findsrc)%srcid if (trim(adjustl(src)) .eq. + trim(adjustl(gridsrcs(findsrc)%srcid)))then lfound=.true. else findsrc=findsrc+1 endif enddo endif if (.not. lfound)findsrc=0 return end c########################################################################## c function to find facility index function findfac(fid) use main1 implicit none integer findfac character fid*15 logical lfound findfac=1 lfound=.false. do while(findfac .le. nsrc .and. .not. lfound) if (trim(adjustl(fid)) .eq. trim(adjustl(facid(findfac))))then lfound=.true. else findfac=findfac+1 endif enddo if (.not. lfound) findfac=0 return end c########################################################################## c subroutine to deallocate arrays subroutine cleanup use main1 implicit none if (runap) then deallocate(apsrcs) elseif (runpt) then deallocate(ptsrcs) elseif (runcmv) then deallocate(cmvsrcs) c elseif (noncontr) then c deallocate(tractsrcs) else deallocate(gridsrcs) endif if (allocated(memis)) deallocate(memis) if (allocated(polyx)) deallocate(polyx) if (allocated(polyy)) deallocate(polyy) if (allocated(modgc)) deallocate(modgc) if (allocated(modgrid)) deallocate(modgrid) if (allocated(modbc)) deallocate(modbc) if (allocated(modblock)) deallocate(modblock) if (allocated(modmc)) deallocate(modmc) if (allocated(modmon)) deallocate(modmon) if (allocated(intbc)) deallocate(intbc) if (allocated(intblock)) deallocate(intblock) if (allocated(intmc)) deallocate(intmc) if (allocated(intmon)) deallocate(intmon) if (allocated(gridrec1)) deallocate(gridrec1) if (allocated(blockrec1)) deallocate(blockrec1) if (allocated(monrec1))deallocate(monrec1) if (allocated(modzbc)) deallocate(modzbc) if (allocated(zmodblock)) deallocate(zmodblock) return end c########################################################################## c subroutine to calculate annual average concentration from seasonal c concentration c subroutine calcann(tag,conc,is,emis,emis1000,ir) subroutine calcann(tag,conc,is,emis,ir) use main1 implicit none integer*8 is,ir,e,hrs,h real*8 conc,emis(12) !,emis1000(4) character tag*2 conc=0.0 hrs=0 do e=1,12 if (tag .eq. 'GM') then conc=conc+modgc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e) elseif (tag .eq. 'BM') then conc=conc+modbc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e) elseif (tag .eq. 'ZB') then conc=conc+modzbc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e) elseif (tag .eq. 'MM') then conc=conc+modmc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e) elseif (tag .eq. 'BI') then conc=conc+intbc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e) else conc=conc+intmc(is,ir,e)*emis(e)/memis(is,e) !emis1000(e) endif hrs=hrs+nhours(e) enddo c if (tag .eq. 'MM') then c do e=1,4 c write(113,*)ir,e,conc,modmc(is,ir,e),hrs c do h=1,24 c write(113,*)h,seashrs(e,h) c enddo c enddo c endif conc=conc/real(hrs) return end c########################################################################## c function to find receptor index function findrec(tag,r) use main1 implicit none integer*8 r,findrec,fr character tag*2 real*8 x1,y1,x2,y2 logical lfound fr=1 lfound=.false. if (tag .eq. 'GM') then do while(fr .le. ngrids1 .and. .not. lfound) if (trim(adjustl(gridrec1(fr)%recid)) .eq. + trim(adjustl(modgrid(r)%recid))) then lfound=.true. else fr=fr+1 endif enddo findrec=gridrec1(fr)%irec endif if (tag .eq. 'BM') then do while(fr .le. nblocks1 .and. .not. lfound) if (trim(adjustl(blockrec1(fr)%recid)) .eq. + trim(adjustl(modblock(r)%recid))) then lfound=.true. else fr=fr+1 endif enddo findrec=blockrec1(fr)%irec endif if (tag .eq. 'ZB') then do while(fr .le. nzblocks1 .and. .not. lfound) if (trim(adjustl(zblockrec1(fr)%recid)) .eq. + trim(adjustl(zmodblock(r)%recid))) then lfound=.true. else fr=fr+1 endif enddo findrec=zblockrec1(fr)%irec endif if (tag .eq. 'MM') then do while(fr .le. nmons1 .and. .not. lfound) if (trim(adjustl(monrec1(fr)%recid)) .eq. + trim(adjustl(modmon(r)%recid))) then lfound=.true. else fr=fr+1 endif enddo findrec=monrec1(fr)%irec endif if (tag .eq. 'BI') then do while(fr .le. nblocks1 .and. .not. lfound) if (trim(adjustl(blockrec1(fr)%recid)) .eq. + trim(adjustl(intblock(r)%recid))) then lfound=.true. else fr=fr+1 endif enddo findrec=blockrec1(fr)%irec endif if (tag .eq. 'MI') then do while(fr .le. nmons1 .and. .not. lfound) if (trim(adjustl(monrec1(fr)%recid)) .eq. + trim(adjustl(intmon(r)%recid))) then lfound=.true. else fr=fr+1 endif enddo findrec=monrec1(fr)%irec endif return end c########################################################################## c funtion to find HAP index function findhap(cas1) use main1 implicit none integer*8 findhap character cas1*15 logical lfound findhap=1 lfound=.false. do while(findhap .le. nhaps .and. .not. lfound) if (trim(adjustl(cas1)) .eq. trim(adjustl(haps(findhap)))) + then lfound=.true. else findhap=findhap+1 endif enddo if (.not. lfound) findhap=0 return end c########################################################################## c funtion to find category index function findcat(cat1) use main1 implicit none integer*8 findcat character cat1*100 logical lfound findcat=1 lfound=.false. do while(findcat .le. ncats .and. .not. lfound) if (trim(adjustl(cat1)) .eq. trim(adjustl(cats(findcat)))) + then lfound=.true. else findcat=findcat+1 endif enddo if (.not. lfound) findcat=0 return end c########################################################################## c subroutine to write output files c only write out when total concentration is > 0 subroutine writeout use main1 implicit none integer*8 h,r,c,rt character outfil*250,ares*2,form1*45,form2*25,acom*1,acom1, +rt1*30,arun*10,zipline*300 real*8 totl,calctot logical lfill allocatable :: acom(:) write(logunit,*)'Writing output files..' call datetime c write(form1,'(a,i3,a)')'(a15,1x,i1,',ncats,'(1x,e13.6))' c write(form1,'(a,i3,a)')'(a,a1,i1,a1,',ncats,'(e12.6,a1))' c write(form1,'(a,i3,a)')'(a,a1,i1,a1,2(i3.3,a1),',ncats, c +'(e12.6,a1))' write(form1,'(a,i3,a)')'(a,a1,i7,a1,i1,a1,2(i3.3,a1),',ncats, +'(e12.6,a1))' write(form2,'(a,i3,a)')'(a,',ncats*2,'(a))' allocate(acom(ncats)) acom1=',' acom=',' acom(ncats)=' ' !set last one to blank write(arun,'(i10)')runno !run # c if single facility or state, append that to filename if (lsingfac) then write(rt1,'(3(a))')trim(adjustl(runtyp)),'_', + trim(adjustl(singlfac)) elseif (lsingst) then write(rt1,'(3(a))')trim(adjustl(runtyp)),'_', + trim(adjustl(singlst)) else rt1=trim(adjustl(runtyp)) endif do h=1,nhaps call datetime lfill=.false. write(logunit,'(2(a))')'Writing output for ', + trim(adjustl(haps(h))) if (runpt .or. runap .or. runcmv) then write(outfil,'(7(a))')trim(adjustl(hapdir)), + trim(adjustl(haps(h))),'_',trim(adjustl(rt1)), + '_',trim(adjustl(arun)),'_concs.txt' else write(ares,'(i2.2)')gridres write(outfil,'(9(a))')trim(adjustl(hapdir)), + trim(adjustl(haps(h))),'_',trim(adjustl(rt1)),'_', + trim(adjustl(ares)),'_',trim(adjustl(arun)),'_concs.txt' endif open(hapunit,file=outfil,status='unknown') c write header record write(hapunit,form2)'receptor_id,recnum,rec_type,col,row,', + (trim(adjustl(cats(c))),acom(c),c=1,ncats) do r=1,ngrids !gridded rt=1 totl=calctot('G',h,r) if (totl .gt. 0.0) then lfill=.true. write(hapunit,form1)trim(adjustl(gridrec(r)%recid)), + acom1,gridrec(r)%irec,acom1,rt,acom1,gridrec(r)%col, + acom1,gridrec(r)%row,acom1,(gridconc(h,r,c),acom(c), + c=1,ncats) endif enddo do r=1,nblocks !blocks rt=2 totl=calctot('B',h,r) if (totl .gt. 0.0) then if (.not. lfill) lfill=.true. write(hapunit,form1)trim(adjustl(blockrec(r)%recid)), + acom1,blockrec(r)%irec,acom1,rt,acom1,blockrec(r)%col, + acom1,blockrec(r)%row,acom1,(blokconc(h,r,c), + acom(c),c=1,ncats) endif enddo do r=1,nzblocks !blocks rt=4 totl=calctot('Z',h,r) if (totl .gt. 0.0) then if (.not. lfill) lfill=.true. write(hapunit,form1)trim(adjustl(zblockrec(r)%recid)), + acom1,zblockrec(r)%irec,acom1,rt,acom1, + zblockrec(r)%col,acom1,zblockrec(r)%row,acom1, + (zblokconc(h,r,c),acom(c),c=1,ncats) endif enddo do r=1,nmons !monitors rt=3 totl=calctot('M',h,r) if (totl .gt. 0.0) then if (.not. lfill) lfill=.true. write(hapunit,form1)trim(adjustl(monrec(r)%recid)), + acom1,monrec(r)%irec,acom1,rt,acom1,monrec(r)%col, + acom1,monrec(r)%row,acom1,(monconc(h,r,c),acom(c), + c=1,ncats) endif enddo write(hapunit,*)' ' call datetime close(hapunit) write(zipline,'(a,1x,a)')'gzip -f',trim(adjustl(outfil)) call system(zipline) !zip file if (.not. lfill) write(logunit,'(2(a))') + 'No output written for ',trim(adjustl(haps(h))) enddo write(logunit,*)'Output files written' call datetime return end c########################################################################## c fucntion to total up across the categories function calctot(rt,h,r) use main1 implicit none integer*8 h,r,c real*8 calctot character rt*1 calctot=0 if (trim(adjustl(rt)) .eq. 'G') then do c=1,ncats calctot=gridconc(h,r,c)+calctot enddo elseif (trim(adjustl(rt)) .eq. 'B') then do c=1,ncats calctot=blokconc(h,r,c)+calctot enddo elseif (trim(adjustl(rt)) .eq. 'Z') then do c=1,ncats calctot=zblokconc(h,r,c)+calctot enddo else do c=1,ncats calctot=monconc(h,r,c)+calctot enddo endif return end c########################################################################## c subroutine to check specific gridded source to receptor distance subroutine check_rec_grid(x,y,src,r,iskip) use main1 implicit none character src*12 integer*8 r,iskip,is,findsrc real*8 distance,d,cenx,ceny,x,y c x=modgrid(r)%ux c y=modgrid(r)%uy c find the source in the source list is=findsrc(src) call calc_center(is,cenx,ceny) d=distance(x,y,cenx,ceny) if (d .gt. grddist) then iskip=1 else iskip=0 endif c write(137,'(a15,2(1x,f13.3),1x,a,3(1x,f13.3),1x,i1)') c +modgrid(r)%recid,x,y,src,cenx,ceny,d,iskip return end c########################################################################## c function to get state index function get_state(sstate) use main1 implicit none integer*8 i,get_state,ii character st1*2,sstate*2 logical lfound,lnum lfound=.false. i=1 ii=ichar(sstate(1:1)) lnum=.false. if (ii .ge. 48 .and. ii .le. 57 .and. ii .ge. 48 .and. ii .le. 57) +lnum=.true. !state is defined by fips code, not abbreviation do while (i .le. 54 .and. .not. lfound) if (lnum) then st1=states1(i) else st1=states(i) endif if (trim(adjustl(sstate)) .eq. trim(adjustl(st1))) then lfound=.true. else i=i+1 endif enddo if (.not. lfound) then get_state=0 else get_state=i endif return end c########################################################################## c subroutine to make sure grid resolution from input argument is valid subroutine checkgrid(agrid) use main1 implicit none integer i,ii,ilen character agrid*2 logical lbad gridres=0.0 ilen=len_trim(agrid) i=1 lbad=.false. do while (i .le. ilen .and. .not. lbad) ii=ichar(agrid(i:i)) if (ii .lt. 48 .or. ii .gt. 57) lbad=.true. i=i+1 enddo if (.not. lbad) read(agrid,*)gridres return end