c----------------------------------------------------------------------
c --- FGEMISS -- Emissions Processor for Flue Gas input
c             -- to CALPUFF via PTEMARB.DAT file
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 011107                   MAIN
c                D. Strimaitis, Earth Tech
c
c     Copyright (c) 1999 by TRC Environmental Corporation
c
c --- PURPOSE:  Prepare hourly water vapor (g/s) emission, exit velocity
c               and temperature (K) for a flug gas stream (point source)
c               whose operation varies with ambient temperature.
c               Pass all inputs required for CALPUFF variable emissions
c               file to output file.
c
c----------------------------------------------------------------------
c
c --- Include parameters
      include 'params.fge'
      include 'qa.fge'

c --- set version and level number of program
      ver='1.0'
      level='011107'

c --- SETUP PHASE -- initialization & program setup operations
      call SETUP

c --- COMPUTATIONAL PHASE -- basic time loop with output
      call COMP

c --- TERMINATION PHASE -- program termination functions
      call FIN

      stop
      end
c----------------------------------------------------------------------
      subroutine setup
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 011107                  SETUP
c                D. Strimaitis, Earth Tech
c
c --- SETUP PHASE of model run -- perform all initialization and
c     one-time setup operations
c
c --- INPUT:
c              none
c
c --- OUTPUT:
c
c     Common block /PT2/
c        TIEM2(7,mxpt2), BHT2(36,mxpt2), BWD2(36,mxpt2)
c
c --- Parameters used:
c        MXSPEC, IO6, MXPT2
c
c --- SETUP called by: MAIN
c --- SETUP calls:     DATETM,  COMLINE, READCF,
c                      OPENOT,  WRHDEM2, WRTIEM2
c                      
c----------------------------------------------------------------------
c
c --- Include parameters
      include 'params.fge'

c --- Include common blocks
      include 'filnam.fge'
      include 'pt2.fge'
      include 'qa.fge'

      character*80 aline

c --- Get date and time from system
      call DATETM(rdate,rtime,rcpu)

c --- Get the name of the control file from the command line
      call COMLINE(fgeinp)

c --- Open the control file
      open(io5,file=fgeinp,status='old')

      write(iomesg,*)'SETUP PHASE'

c --- Check that the version and level number in the parameter
c --- file matches those in the code itself
      if(ver.ne.mver.or.level.ne.mlevel)then
         write(iomesg,10)mmodel,ver,level,mver,mlevel
10       format(/1x,'ERROR in SUBR. SETUP -- The ',a12,' version ',
     1   'and level numbers do not match those in the parameter file'/
     2   5x,'    Model Code - Version: ',a12,' Level: ',a12/
     3   5x,'Parameter File - Version: ',a12,' Level: ',a12)
         stop
      endif

c --- Read control file inputs
      call READCF

c --- Open all other I/O files
      call OPENOT

c --- Report inputs to list file
      call LISTINP

c --- Read ISC header record (station,year for surface and upper)
      read(io7,'(a80)') aline
      write(*,*) 'ISC Header:'
      write(*,*) aline

c --- Write header records to the emissions file:
c     PTEMARB.DAT - io16 - Arbitrarily-varying pt. source emissions
      call WRHDEM2(io16)

c --- Write the time-invariant data to the PTEMARB.DAT file
c     (arbitrarily-varying point source emissions)
      call WRTIEM2(io16,lascii,npt2,cid,tiem2,bht2,bwd2)
c
      return
      end
c----------------------------------------------------------------------
      subroutine comp
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 011107                   COMP
c                D. Strimaitis,  Earth Tech
c
c --- PURPOSE:  Compute emissions for each hour in the period
c
c
c --- INPUTS:
c     Parameters:
c           MXSPEC, MXPT2, IO6, IO7, IO16 
c     Common block /PT2/:
c           IBDAT2, IBTIM2, IEDAT2, IETIM2
c
c --- OUTPUT:
c     Common block /QA/:
c           iyr, month, iday, ijul, ihr
c
c --- COMP called by:  MAIN
c --- COMP calls:      GRDAY, FGPROC, WREM2, INCR
c----------------------------------------------------------------------
c
c --- Include parameter statements
      include 'params.fge'

c --- Include common blocks
      include 'pt2.fge'
      include 'qa.fge'

      logical leof

c --- Set starting date-hour (assume emissions change hourly)
      idathr=ibdat2*100+ibtim2
c --- Set ending date-hour
      idathr0=iedat2*100+ietim2

c --- Set End-Of-File logical to .FALSE.
      leof=.FALSE.

      write(iomesg,*)'COMPUTATIONAL PHASE'

c --- Loop over hours
c -----------------------------------------------------------

c --- Read current temperature and relative humidity
10    call RDISC(idathr,tmpk,rhamb,kdathr,leof)
      if(LEOF) then
         write(io6,*)
         write(io6,*) 'EOF reached in MET data file'
         write(io6,*) 'Last period attempted YYYYJJJHH =',idathr
         write(io6,*)
         write(*,*)
         write(*,*) 'EOF reached in MET data file'
         write(*,*) 'Last period attempted YYYYJJJHH =',idathr
         write(*,*)
         return
      endif

c --- Set month, hour
      idate=idathr/100
      ihr=idathr-idate*100
      iyr =idate/1000
      ijul=idate-iyr*1000
      call GRDAY(iyr,ijul,month,iday)

c --- Compute flue gas properties for each source this hour
      call FGPROC(tmpk)

c --- Write record for this hour
      call WREM2(idate,ihr)

c --- Increment date-time for next hour
      nhrinc=1
      nyr=iyr
      njul=ijul
      nhr=ihr
      call INCR(nyr,njul,nhr,nhrinc)
      idathr=nyr*100000+njul*100+nhr

c --- Loop until last hour of period has been processed
      if(idathr.LE.idathr0) goto 10

      return
      end
c----------------------------------------------------------------------
      subroutine openot
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 011107                 OPENOT
c                D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Open all input/output files other than input control
c ---           file (IO5) and output list file (IO6)
c
c
c --- INPUTS:
c     Common block /FILNAM/ variables:
c           ISCDAT, PTDAT
c     Common block /PT2/ variables:
c           LASCII
c     Parameters: IO7, IO16
c
c --- OUTPUT:
c
c --- OPENOT called by:  SETUP
c --- OPENOT calls:      none
c----------------------------------------------------------------------
c
c --- Include parameter statements
      include 'params.fge'
c
c --- Include common blocks
      include 'filnam.fge'
      include 'pt2.fge'

c --- Open ISC meteorological data file (ISCMET.DAT)
      open(io7,file=iscdat,status='old')

c --- Open POINT SOURCE emissions file (PTEMARB.DAT)
c --- (stationary point sources with arbitrary variation in emissions)
      if(LASCII) then
         open(io16,file=ptdat,status='unknown')
      else
         open(io16,file=ptdat,status='unknown',form='unformatted')
      endif

      return
      end
c----------------------------------------------------------------------
      subroutine readcf
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 011107                 READCF
c                D. Strimaitis,  Earth Tech
c
c --- PURPOSE:  Read the structured control file (IO5)
c
c
c --- INPUTS:
c     Common block /QA/:
c           VER, LEVEL
c     Parameters:
c           MXSPEC, IO5, IO6,          
c           MMODEL
c
c --- OUTPUT:
c
c     Common block /FILNAM/:
c           fgelst,iscdat,ptdat
c     Common block /PT2/:
c           NPT2,
c           FNAME2, NSE2, IUTMZ2, IBDAT2, IBTIM2,
c           IEDAT2, IETIM2, VRS2, LABEL2, CSLST2(mxspec),
c           XMWEM2(mxspec),
c           TIEM2(7,mxpt2), BHT2(36,mxpt2),BWD2(36,mxpt2),
c           cid(mxpt2),tbin(12),
c           wexit12(12,mxpt2),texit12(12,mxpt2),qexit12(12,mxpt2)
c           
c
c --- READCF called by:  SETUP
c --- READCF calls:      QAYR4, YR4
c----------------------------------------------------------------------
c
c --- Include parameter statements
      include 'params.fge'

c --- Include common blocks
      include 'filnam.fge'
      include 'pt2.fge'
      include 'qa.fge'

c --- CID matching variable
      character*16 acid
      character*1 achar

c --- Section character field
      character*16 asection,asectdb(6)
      data asectdb/'--- Section 1 --','--- Section 2 --',
     &             '--- Section 3 --','--- Section 4 --',
     &             '--- Section 5 --','--- Section 6 --'/

c --- bdf additions:
      parameter(ndir=36)
      character*25 bpipstart
      character*70 bpipinp
      character*8 stkn(mxpt2),stkntemp
      integer npt2bpip
      real mxpbh(mxpt2,ndir),mxpbw(mxpt2,ndir)

c --- Skip to Section 1
1     read(io5,'(a16)',end=999) asection
      isect=1
      if(asection.NE.asectdb(isect)) goto 1
      read(io5,*)
c
c ----------------------------------------------
c --- General run control parameters -- Sect # 1
c ----------------------------------------------
c
c --- Read file names from control file
      read(io5,*) iscdat
      read(io5,*) fgelst
      read(io5,*) ptdat
      
c --- Open the output LIST file
      open(io6,file=fgelst,status='unknown')

c --- Write page header
      write(io6,1400)mmodel,ver,level
1400  format(40x,a12,3x,'Version: ',a12,3x,'Level: ',a12/
     1 1x,13('**********')/)
c
c --- Write clock date/time of run
      write(io6,1402)rtime,rdate
1402  format(//2x,'Clock time: ',a8/
     1         2x,'      Date: ',a10)

c --- Set format of variable emissions file
      read(io5,*) achar
      if(achar.EQ.'A') then
         lascii=.TRUE.
      elseif(achar.EQ.'B') then
         lascii=.FALSE.
      else
         write(*,*)'READCF:  PTEMARB format must be A,B : ',achar
         write(io5,*)'READCF:  PTEMARB format must be A,B : ',achar
         stop
      endif

c --- Set number of sources
      read(io5,*) npt2

c --- UTM Zone
      read(io5,*) iutmz2

c --- Processing period 
      read(io5,*) ibdat2
      read(io5,*) ibtim2
      read(io5,*) iedat2
      read(io5,*) ietim2

c --- Special QA on starting year of simulation
      ibyr=ibdat2/1000
      call QAYR4(io6,ibyr,0,ierr)
      if(ierr.NE.0) stop 'CORRECT THE CONTROL FILE'

c --- Do Y2K checks
      ieyr=iedat2/1000
      call YR4(io6,ieyr,ierr)
      if(ierr.NE.0) stop 'Halted in READCF'

c --- 12 Temperature range markers (K)
      read(io5,*)
      read(io5,*) tbin
c --- Check for increasing values
      do i=2,12
         if(tbin(i).LE.tbin(i-1)) then
            write(*,*)'READCF:  Temperatures (K) must increase!'
            write(*,'(10x,f10.2)') tbin
            stop
         endif
      enddo

c --- Skip to Section 2
2     read(io5,'(a16)',end=999) asection
      isect=2
      if(asection.NE.asectdb(isect)) goto 2
      read(io5,*)
c
c ----------------------------------------------
c --- Constant Source Data -- Sect # 2
c ----------------------------------------------
c
c --- Loop over sources
      do ic=1,npt2
         read(io5,*)cid(ic)
         read(io5,*)tiem2(1,ic)          ! x(km)
         read(io5,*)tiem2(2,ic)          ! y(km)
         read(io5,*)tiem2(3,ic)          ! Ht(m)
         read(io5,*)tiem2(4,ic)          ! Diam(m)
         read(io5,*)tiem2(5,ic)          ! Elev(m MSL)
         read(io5,*)tiem2(6,ic)          ! Building DW
         tiem2(7,ic)=0.                  ! User
      enddo

c --- Skip to Section 3
3     read(io5,'(a16)',end=999) asection
      isect=3
      if(asection.NE.asectdb(isect)) goto 3
      read(io5,*)
c
c ----------------------------------------------
c --- 12 Exit Velocities (m/s) -- Sect # 3
c ----------------------------------------------
c
c --- Loop over sources
      do i=1,npt2
         read(io5,*)acid
c ---    Match with CID array
         ic=0
         do j=1,npt2
            if(acid.EQ.cid(j)) ic=j
         enddo
         if(ic.EQ.0) stop 'No Match on source name'
         read(io5,*) (wexit12(k,ic),k=1,12)
      enddo

c --- Skip to Section 4
4     read(io5,'(a16)',end=999) asection
      isect=4
      if(asection.NE.asectdb(isect)) goto 4
      read(io5,*)
c
c ----------------------------------------------
c --- 12 Exit Temperatures (K) -- Sect # 4
c ----------------------------------------------
c
c --- Loop over sources
      do i=1,npt2
         read(io5,*)acid
c ---    Match with CID array
         ic=0
         do j=1,npt2
            if(acid.EQ.cid(j)) ic=j
         enddo
         if(ic.EQ.0) stop 'No Match on source name'
         read(io5,*) (texit12(k,ic),k=1,12)
      enddo

c --- Skip to Section 5
5     read(io5,'(a16)',end=999) asection
      isect=5
      if(asection.NE.asectdb(isect)) goto 5
      read(io5,*)
c
c -------------------------------------------------------
c --- 12 Additional Water Emission Rates (K) -- Sect # 5
c -------------------------------------------------------
c
c --- Loop over sources
      do i=1,npt2
         read(io5,*)acid
c ---    Match with CID array
         ic=0
         do j=1,npt2
            if(acid.EQ.cid(j)) ic=j
         enddo
         if(ic.EQ.0) stop 'No Match on source name'
         read(io5,*) (qexit12(k,ic),k=1,12)
      enddo

c --- Skip to Section 6
6     read(io5,'(a16)',end=997) asection
      isect=6
      if(asection.NE.asectdb(isect)) goto 6
      read(io5,*)
c
c ----------------------------------------------
c --- Building Data -- Sect # 6
c ----------------------------------------------
c
c --- read all building information present
      do i=1,mxpt2
            read(io5,*,end=995) stkn(i)
		if (stkn(i).eq.'--------') goto 995
            read(io5,*,end=995)
c ---       Heights
            read(io5,*,end=995) (mxpbh(i,k),k=1,36)
            read(io5,*,end=995)
c ---       Widths
            read(io5,*,end=995) (mxpbw(i,k),k=1,36)
	enddo

995	npt2bpip = i - 1

c	write (*,*) 'Number of cells found: ',npt2bpip
c	write (*,*) 'Number of cells in input file: ',npt2

c map data to correct source and corresponding arrays
c --- Loop over sources
      do i=1,npt2
	  if (tiem2(6,i).eq.1) then
c ---     Match with CID array
          write(stkntemp,'(a8)') cid(i)
          ic=0
          do j=1,npt2bpip
            if(stkn(j).EQ.stkntemp) ic=j
          enddo
c if no source is found with matching name, ask user
          if (ic.EQ.0) then
            write(*,*) 'Please enter matching source number'
		write(*,*) 'Current source name: ',cid(i)
		do k = 1,npt2bpip
		  write(*,*) k,' = ',stkn(k)
		enddo
		write(*,*)
		read(*,*) ic
	    endif
c fill in arrays with correct data
	    do k = 1,36
		bht2(k,i) = mxpbh(ic,k)
		bwd2(k,i) = mxpbw(ic,k)
	    enddo
	  else
	    do k = 1,36
	      bht2(k,i) = 999.9
	      bwd2(k,i) = 999.9
	    enddo
	  endif
      enddo
c
c	write(*,*) 'Section 6 taken from input file'
	write(sec6header,288)

	goto 998
c
c if there is no section 6, get information from BPIP file

997	continue

c
c ----------------------------------------------
c --- Retrieving Data from BPIP file -- Sect # 6 -- added by BdF
c ----------------------------------------------
c
c If section 6 is missing from the input file, this asks the user
c for a BPIP file to retrieve the information from.
c It reads the BPIP file and tries to match the cell names.
c If this is not possible, the user is asked to match the cells interactively.
c Section 6 is then appended to the input file for future use.	
c   There may be extra gaps between section 5 & 6 depending on the original
c   input file, however this does not affect future runs with this program.
c
c open BPIP output file
	write(*,*) 'Building information not in ip file,',
     &	' enter BPIP filename'
	write(*,*)
	read (*,*) bpipinp
	write(*,*)
	write(*,*) bpipinp

	open(io3,file=bpipinp,status='old')
c
c --- Skip to beginning of data
7     read(io3,'(a25)',end=997) bpipstart
      isect=6
      if(bpipstart.NE.' BPIP output is in meters') goto 7

c start reading

          DO 510 S = 1, mxpt2

	      read(io3,*,end=996)
	      read(io3,*,end=996)

c read cell heights
            DO 512 I = 1, 6
              J = (I-1) * 6 + 1
              K = I * 6
             read (io3,293,end=996) STKN(S), (MXPBH(S,D) , D = J,K)
512         CONTINUE
c read cell widths
            DO 514 I = 1, 6
              J = (I-1) * 6 + 1
              K = I * 6
             read (io3,293,end=996) STKN(S), (MXPBW(S,D) , D = J,K)
514         CONTINUE
510      CONTINUE
	
293   FORMAT(17x, A8, 6F8.2)
294   FORMAT(25x,6F8.2)

c end of file reached, no more cells
996	continue

	close (io3)

	npt2bpip = s-1
	write (*,*) 'No. sources in BPIP file: ',npt2bpip,
     &	' No. sources in ip file: ',npt2

c map data to correct source and corresponding arrays
c --- Loop over sources
      do i=1,npt2
	  if (tiem2(6,i).eq.1) then
c ---     Match with CID array
          write(stkntemp,'(a8)') cid(i)
          ic=0
          do j=1,npt2bpip
            if(stkn(j).EQ.stkntemp) ic=j
          enddo
c if no cell is found with matching name, ask user
          if (ic.EQ.0) then
             do k = 1,npt2bpip
                write(*,*) k,' = ',stkn(k)
             enddo
             write(*,*)'Please enter matching # for source: ',cid(i)
             write(*,*)'from the list above'
             write(*,*)
             read(*,*) ic
	    else
               write(*,*) 'matching source ',ic,' = ',cid(i)
	    endif
c fill in arrays with correct data
	    do k = 1,36
		bht2(k,i) = mxpbh(ic,k)
		bwd2(k,i) = mxpbw(ic,k)
	    enddo
	  else
	    write(*,*) 'no downwash for building ',cid(i),' no: ',i
	    do k = 1,36
	      bht2(k,i) = 999.9
	      bwd2(k,i) = 999.9
	    enddo
	  endif
      enddo
c
c write in section 6 in original input file
c NB: there may be empty space between the end of section 5 and the beginning
c of section 6, although this won't matter for this program
	write(*,*) 'Appending BPIP data to input file'
c write section 6 header
	write(io5,*)
	write(io5,280)
	write(io5,281)
	write(io5,287) bpipinp
	write(sec6header,287) bpipinp

c write information
	do i = 1,npt2
	  if (tiem2(6,i).eq.1) then
	    write(io5,282) cid(i)
	    write(io5,283)
	    write(io5,284) (bht2(j,i),j=1,ndir)
	    write(io5,285)
	    write(io5,284) (bwd2(j,i),j=1,ndir)
	  endif
	enddo
c write file footer
	write(io5,280)
	write(io5,286)
	write(io5,280)

280	format('-----------------------------------------------',
     &	'-----------------------')
281   format('--- Section 6 ---    Building Data by flow ',
     &	'direction (Repeat for each source)')
287   format('--- taken from BPIP file:',a55)
288   format('------------------------------------',
     &	' taken from INPUT file ')
282  	format("'",a16,"'",'  : Cell ID (16 characters or less)')
283	format('--- Height(m)')
284	format (6F7.2)
285	format('--- Width(m)')
286	format('--- END ---   ')
c
998   continue

c --- Close control file
      close(io5)
c

      return

999   write(*,*)'READCF Error in Control File: '
      write(*,*)'     Section boundary not correct in Section ',isect
      stop

      end
c----------------------------------------------------------------------
      subroutine listinp
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 011107                LISTINP
c                D. Strimaitis,  Earth Tech
c
c --- PURPOSE:  Write inputs to list file
c
c
c --- INPUTS:
c     Parameters:
c           MXSPEC, IO6,          
c     Common block /FILNAM/:
c           fgelst,iscdat,ptdat
c     Common block /PT2/:
c           NPT2,
c           FNAME2, NSE2, IUTMZ2, IBDAT2, IBTIM2,
c           IEDAT2, IETIM2, VRS2, LABEL2, CSLST2(mxspec),
c           XMWEM2(mxspec),
c           TIEM2(7,mxpt2), BHT2(36,mxpt2),BWD2(36,mxpt2),
c           cid(mxpt2),tbin(12),
c           wexit12(12,mxpt2),texit12(12,mxpt2),qexit12(12,mxpt2)
c
c --- OUTPUT:
c
c --- LISTINP called by:  SETUP
c --- LISTINP calls:      none
c----------------------------------------------------------------------
c
c --- Include parameter statements
      include 'params.fge'

c --- Include common blocks
      include 'filnam.fge'
      include 'pt2.fge'
c
      character*6 cr

      write(io6,*)
      write(io6,*)' ----------------------------------------------'
      write(io6,*)' --- General run control parameters -- Sect # 1'
      write(io6,*)' ----------------------------------------------'
      write(io6,*)
      write(io6,'(a16,a24)')' Met Data File: ',iscdat
      write(io6,'(a16,a24)')'     List File: ',fgelst
      write(io6,'(a16,a24)')'Emissions File: ',ptdat
      if(LASCII) then
         write(io6,*)'                 (ASCII)'
      else
         write(io6,*)'                 (BINARY)'
      endif
      write(io6,*)
      write(io6,*)'      UTM Zone: ',iutmz2
      write(io6,*)'         Start: ',ibdat2,ibtim2
      write(io6,*)'           End: ',iedat2,ietim2
      write(io6,*)'   Temp(K) Pts: '
      do i=1,12
         write(io6,'(20x,i2,5x,f7.3)') i,tbin(i)
      enddo

      write(io6,*)
      write(io6,*)' ----------------------------------------------'
      write(io6,*)' --- Constant Source Data -- Sect # 2'
      write(io6,*)' ----------------------------------------------'
c --- Loop over sources
      do ic=1,npt2
         write(io6,*)
         write(io6,*)'            Source: ',cid(ic)
         write(io6,*)'             x(km): ',tiem2(1,ic)
         write(io6,*)'             y(km): ',tiem2(2,ic)
         write(io6,*)'             Ht(m): ',tiem2(3,ic)
         write(io6,*)'           Diam(m): ',tiem2(4,ic)
         write(io6,*)'       Elev(m MSL): ',tiem2(5,ic)
         ibdw=NINT(tiem2(6,ic))
         write(io6,*)'Building DW(1=yes): ',ibdw
      enddo

      write(io6,*)
      write(io6,*)' ----------------------------------------------'
      write(io6,*)' --- 12 Exit Velocities (m/s) -- Sect # 3'
      write(io6,*)' ----------------------------------------------'
c --- Loop over cells
      do ic=1,npt2
         cr='T pnts'
         write(io6,*)
         write(io6,*)'     Source: ',cid(ic)
         write(io6,'(a6,a7,6f7.3)')cr,' 1-6 : ',(wexit12(k,ic),k=1,6)
         write(io6,'(a6,a7,6f7.3)')cr,' 7-12: ',(wexit12(k,ic),k=7,12)
      enddo

      write(io6,*)
      write(io6,*)' ----------------------------------------------'
      write(io6,*)' --- 12 Exit Temperatures (K) -- Sect # 4'
      write(io6,*)' ----------------------------------------------'
c --- Loop over sources
      do ic=1,npt2
         cr='T pnts'
         write(io6,*)
         write(io6,*)'     Source: ',cid(ic)
         write(io6,'(a6,a7,6f7.2)')cr,' 1-6 : ',(texit12(k,ic),k=1,6)
         write(io6,'(a6,a7,6f7.2)')cr,' 7-12: ',(texit12(k,ic),k=7,12)
      enddo

      write(io6,*)
      write(io6,*)' ----------------------------------------------'
      write(io6,*)' --- 12 Water Emission Rates (g/s) -- Sect # 5'
      write(io6,*)' ----------------------------------------------'
c --- Loop over sources
      do ic=1,npt2
         cr='T pnts'
         write(io6,*)
         write(io6,*)'     Source: ',cid(ic)
         write(io6,'(a6,a7,6e10.2)')cr,' 1-6 : ',(qexit12(k,ic),k=1,6)
         write(io6,'(a6,a7,6e10.2)')cr,' 7-12: ',(qexit12(k,ic),k=7,12)
      enddo

      write(io6,*)
      write(io6,*)' ----------------------------------------------'
      write(io6,*)' --- Building Downwash Data -- Sect # 6'
      write(io6,'(a80)') sec6header
c --- Loop over cells
      icount = 0
      do ic=1,npt2
         if(tiem2(6,ic).eq.1.) then
           write(io6,*)
           write(io6,*)'     Source: ',cid(ic)
           write(io6,*)'  Heights(m) '
           write(io6,'(a13,6f7.2)')'Range  1-6 : ',(bht2(k,ic),k=1,6)
           write(io6,'(a13,6f7.2)')'Range  7-12: ',(bht2(k,ic),k=7,12)
           write(io6,'(a13,6f7.2)')'Range 13-18: ',(bht2(k,ic),k=13,18)
           write(io6,'(a13,6f7.2)')'Range 19-24: ',(bht2(k,ic),k=19,24)
           write(io6,'(a13,6f7.2)')'Range 25-30: ',(bht2(k,ic),k=25,30)
           write(io6,'(a13,6f7.2)')'Range 31-36: ',(bht2(k,ic),k=31,36)
           write(io6,*)'   Widths(m) '
           write(io6,'(a13,6f7.2)')'Range  1-6 : ',(bwd2(k,ic),k=1,6)
           write(io6,'(a13,6f7.2)')'Range  7-12: ',(bwd2(k,ic),k=7,12)
           write(io6,'(a13,6f7.2)')'Range 13-18: ',(bwd2(k,ic),k=13,18)
           write(io6,'(a13,6f7.2)')'Range 19-24: ',(bwd2(k,ic),k=19,24)
           write(io6,'(a13,6f7.2)')'Range 25-30: ',(bwd2(k,ic),k=25,30)
           write(io6,'(a13,6f7.2)')'Range 31-36: ',(bwd2(k,ic),k=31,36)
           icount = icount + 1
         else
           write(io6,*)
           write(io6,*)'     Source: ',cid(ic)
           write(io6,*)' No downwash for this source'
         endif
      enddo

      write(io6,*)
      write(io6,*) 'No. of sources with downwash   : ',icount
      write(io6,*) 'No. of sources without downwash: ',npt2-icount

      return
      end
c----------------------------------------------------------------------
      subroutine rdisc(idathr,tmpk,rh,kdathr,leof)
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 011107                  RDISC
c                D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Read ISCMET.DAT meteorological data for one hour
c
c --- INPUTS:
c        IDATHR - integer       - Date & hour of required data
c                                 (YYYYJJJHH)
c         LECHO - logical       - Flag controlling printing of input
c                                 meteorological data fields
c    Parameters:
c         IO6
c
c --- OUTPUT:
c                   TMPK - real    - Air temperature (deg. K)
c                     RH - real    - Relative humidity (percent)
c                 KDATHR - integer - Date of hour of data (YYYYJJJHH)
c                   LEOF - logical - End-Of-File flag
c
c --- RDISC called by:  COMP
c --- RDISC calls:      JULDAY, INCR, GRDAY, YR4
c----------------------------------------------------------------------
c
c --- Include parameter statements
      include 'params.fge'
c
      character plfield*5,swfield*9,rhfield*3
      logical leof

c -------------------------------------------------------------------
c  -- Read one record of met data in extended ISC format (ASCII)
1     read(io7,5,end=999) iy,im,id,ih,fvec,wspd,tmpk,kst,rmix,umix,
     &                    ustr,xmon,z0m,ipc,pmmhr,dthdz,
     &                    plfield,swfield,rhfield
5     format(4i2,2f9.4,f6.1,i2,2f7.1,f9.4,f10.1,f8.4,i4,f7.2,
     &       f10.5,a5,a9,a3)
c
c --- Skip blank lines
      if(im.EQ.0 .OR. id.EQ.0) goto 1
c -------------------------------------------------------------------

c --- Do Y2K check
      call YR4(io6,iy,ierr)
      if(ierr.NE.0) stop 'Halted in RDISC'

c --- ISCMET.DAT uses a 1-24 clock (hour ending), whereas CALPUFF uses
c --- a 00-23 (hour ending) convention.  Therefore, convert hour "24"
c --- to hour "00" of the next day
      if(ih.EQ.24) then
c ---    Roll back 1 hour, then use INCR routine to move forward 1 hour
         ih=ih-1
         call JULDAY(iy,im,id,ijul)
         call INCR(iy,ijul,ih,1)
         call GRDAY(iy,ijul,im,id)
      else
         call JULDAY(iy,im,id,ijul)
      endif
      kdathr=iy*100000+ijul*100+ih

c --- Check to see if required date-time has been read
      if(kdathr.lt.idathr)then
c ---    Obtain next time period in MET file
         go to 1
      elseif(kdathr.gt.idathr)then
         write(io6,*)'ERROR in subr. RDISC -- current hour not found ',
     1   'in the ISCMET data file -- Current date/hour = ',idathr,
     2   ' Last date/hour read = ',kdathr
         stop
      endif

c --- Get relative humidity
      if(rhfield.EQ.' ') then
c ---    Blank field; just set to 85%
         irh=85
      else
         read(rhfield,'(bn,i3)') irh
      endif

c --- Return RH as a real
      rh=irh

      return

999   write(io6,*)'RDISC: End-Of-File found in ISCMET file'
      leof=.TRUE.
      return

      end
c----------------------------------------------------------------------
      subroutine WRHDEM2(iunit)
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 011107                WRHDEM2
c                D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Write the header records for the PTEMARB point source
c               emissions file (arbitrarily-varying emissions)
c
c --- INPUTS:
c            IUNIT - integer    - Fortran unit number of PTEMARB file
c     Common Block /PT2/ variables:
c            NPT2, LASCII,
c            IUTMZ2, IBDAT2, IBTIM2,
c            IEDAT2, IETIM2
c     Parameters:
c            IO6, MXSPEC
c
c
c --- OUTPUT:
c            none
c
c --- WRHDEM2 called by: SETUP
c --- WRHDEM2 calls:     none
c----------------------------------------------------------------------
c
      include 'params.fge'
      include 'pt2.fge'

c --- First, set "standard" values
c ---------------------------------
c --- Set File type and label
      fname2='PTEMARB'
      label2='FGEMISS'
c --- Set valid version number to 5
      vrs2='5'
c --- Set species for FOG modeling
      nse2=2
      cslst2(1)='H2O'
      cslst2(2)='TXS'
      xmwem2(1)=1.0
      xmwem2(2)=1.0

      if(LASCII) then

c ---    Write records in ASCII form
c -----------------------------------
c ---    Header Record #1 - No. sources, emitted species, beginning
c                           & ending dates/times
c
         write(iunit,101)fname2,npt2,nse2,iutmz2,ibdat2,ibtim2,
     1                   iedat2,ietim2,vrs2,label2

c ---    Header Record #2 - Species list
         write(iunit,102)(cslst2(n),n=1,nse2)

c ---    Header Record #3 - Molecular weights of each emitted species
         write(iunit,103)(xmwem2(n),n=1,nse2)

      else

c ---    Write records in Binary form
c -----------------------------------
c ---    Header Record #1 - No. sources, emitted species, beginning
c                           & ending dates/times
c
         write(iunit)fname2,npt2,nse2,iutmz2,ibdat2,ibtim2,
     1               iedat2,ietim2,vrs2,label2

c ---    Header Record #2 - Species list
         write(iunit)(cslst2(n),n=1,nse2)

c ---    Header Record #3 - Molecular weights of each emitted species
         write(iunit)(xmwem2(n),n=1,nse2)

      endif

101   format("'",a12,"'",1x,3i6,4i8,2(2x,"'",a12,"'"))
102   format(6a12)
103   format(6f12.3)

      return
      end
c----------------------------------------------------------------------
      subroutine WRTIEM2(iunit,lascii,nsrc2,cid,tiem2,bht2,bwd2)
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 011107                WRTIEM2
c                D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Write the time-invariant records to the PTEMARB point
c               source emissions file (arbitrarily-varying emissions)
c
c --- INPUTS:
c            IUNIT - integer    - Fortran unit number of PTEMARB file
c           LASCII - logical    - ASCII format for PTEMARB file?
c            NSRC2 - integer    - Number of sources in the PTEMARB file
c       CID(nsrc2) - C*16 array - Names for PTEMARB sources
c   TIEM2(7,nsrc2) - real array - Time-invariant PTEMARB source data
c   BHT2(36,nsrc2) - real array - Building heights every 10 deg from N
c   BWD2(36,nsrc2) - real array - Building widths every 10 deg from N
c
c --- OUTPUT:
c             none
c
c --- WRTIEM2 called by: SETUP
c --- WRTIEM2 calls:     none
c----------------------------------------------------------------------
c
      real tiem2(7,nsrc2),bht2(36,nsrc2),bwd2(36,nsrc2)
      character*16 cid(nsrc2)
      logical lascii


      if(LASCII) then

c ---    Write records in ASCII form
c -----------------------------------
c ---    Loop over sources
         do i=1,nsrc2
            write(iunit,101)cid(i),(tiem2(n,i),n=1,7)
            if(NINT(tiem2(6,i)).EQ.1) then
               write(iunit,102)(bht2(n,i),n=1,36)
               write(iunit,102)(bwd2(n,i),n=1,36)
            endif
         enddo

      else

c ---    Write records in Binary form
c -----------------------------------
c ---    Loop over sources
         do i=1,nsrc2
            write(iunit)cid(i),(tiem2(n,i),n=1,7)
            if(NINT(tiem2(6,i)).EQ.1) then
               write(iunit)(bht2(n,i),n=1,36)
               write(iunit)(bwd2(n,i),n=1,36)
            endif
         enddo

      endif

101   format("'",a16,"'",2x,2f11.5,2f8.3,3f8.2)
102   format(6f9.3)

      return
      end
c----------------------------------------------------------------------
      subroutine WREM2(ibdat,ibtim)
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 011107                  WREM2
c                D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Write the time-varying records of the PTEMARB point
c               source emissions file (arbitrarily-varying emissions)
c
c --- INPUTS:
c     Common block /PT2/:
c           NPT2,CID(mxpt2),LASCII,
c           TSTAK2(mxpt2),EXITW2(mxpt2),QSTAK2(mxspec,mxpt2)
c     Parameters:
c           MXPT2, MXSPEC, IO16
c
c --- WREM2 called by: COMP
c --- WREM2 calls:     none
c----------------------------------------------------------------------
c
c --- Include parameter statements and commons
      include 'params.fge'
      include 'pt2.fge'

      if(LASCII) then
      
c ---    Write records in ASCII form
c -----------------------------------
c ---    Write starting, ending date and time (begin=end for 1 hour)
         write(io16,101)ibdat,ibtim,ibdat,ibtim
c ---    Loop over sources
         do i=1,npt2
            write(io16,102)cid(i),tstak2(i),exitw2(i),
     &                     (qstak2(k,i),k=1,2)
         enddo

      else
      
c ---    Write records in Binary form
c -----------------------------------
c ---    Write starting, ending date and time (begin=end for 1 hour)
         write(io16)ibdat,ibtim,ibdat,ibtim
c ---    Loop over sources
         do i=1,npt2
            write(io16)cid(i),tstak2(i),exitw2(i),(qstak2(k,i),k=1,2)
         enddo

      endif

101   format(4i8)
102   format("'",a16,"'",2f10.2,2e14.6)

      return
      end
c----------------------------------------------------------------------
      subroutine FGPROC(tamb)
c----------------------------------------------------------------------
c --- FGEMISS    Version: 1.0      Level: 011107                 FGPROC
c                D. Strimaitis,  Earth Tech
c
c --- PURPOSE:  Flue Gas Emissions Processor
c               Computes the water vapor emission rate (in addition to
c               ambient), exit temperature, and exit velocity for a
c               given ambient temperature.  These properties are
c               assumed to vary linearly between values provided at
c               specific ambient temperatures (in control file).
c
c
c --- INPUTS:
c     Common block /PT2/:
c           NPT2, TIEM2(7,mxpt2),
c           wexit12(12,mxpt2),texit12(12,mxpt2),qexit12(12,mxpt2),
c           tbin(12)
c     Parameters:
c           MXPT2, MXSPEC
c
c --- OUTPUT:
c     Common block /PT2/:
c           TSTAK2(mxpt2),EXITW2(mxpt2),QSTAK2(mxspec,mxpt2)
c
c --- FGPROC called by:  COMP
c --- FGPROC calls:      none
c----------------------------------------------------------------------
c
c --- Include parameter statements and commons
      include 'params.fge'
      include 'pt2.fge'

      logical ldb

      data pi/3.1415927/

c     ldb=.TRUE.
      ldb=.FALSE.

c --- Compute temperature category (1-13) for variability
      itcat=13
      do i=12,1,-1
         if(tamb.LE.tbin(i)) itcat=i
      enddo
      itcatm1=itcat-1
c --- Fraction of interval from tbin(i-1) to tamb
      if(itcat.GT.1 .AND. itcat.LT.13) then
         fac=(tamb-tbin(itcatm1))/(tbin(itcat)-tbin(itcatm1))
      else
         fac=1.0
      endif

c --- DEBUG SECTION
      if(LDB) then
         write(io6,*)
         write(io6,*)'FGPROC:    itcat  = ',itcat
         write(io6,*)'            tamb  = ',tamb
         write(io6,*)
      endif

c --- Loop over sources
c ----------------------
      do is=1,npt2

c ---    Select data from arrays
         r0=0.5*tiem2(4,is)

c ---    Get exit velocity
         if(itcat.EQ.1) then               
            wexit=wexit12(1,is)
         elseif(itcat.EQ.13) then
            wexit=wexit12(12,is)
         else
            wexit=wexit12(itcatm1,is)+
     &             fac*(wexit12(itcat,is)-wexit12(itcatm1,is))
         endif

c ---    Get exit temperature
         if(itcat.EQ.1) then               
            texit=texit12(1,is)
         elseif(itcat.EQ.13) then
            texit=texit12(12,is)
         else
            texit=texit12(itcatm1,is)+
     &             fac*(texit12(itcat,is)-texit12(itcatm1,is))
         endif

c ---    Get emission rate
         if(itcat.EQ.1) then               
            qexit=qexit12(1,is)
         elseif(itcat.EQ.13) then
            qexit=qexit12(12,is)
         else
            qexit=qexit12(itcatm1,is)+
     &             fac*(qexit12(itcat,is)-qexit12(itcatm1,is))
         endif

c ---    DEBUG SECTION
         if(LDB) then
            write(io6,*)'Source:           = ',is
            write(io6,*)'wexit,texit,qexit = ',wexit,texit,qexit
            write(io6,*)
         endif

c ---    Temperature excess (min value is .01K)
         dtemp=AMAX1(0.01,texit-tamb)
         texit=tamb+dtemp

c ---    Area of release (m^2)
         a0=pi*r0**2

c ---    Volume flux (m^3/s)
         vf0=a0*wexit

c ---    Pass results to output arrays
         tstak2(is)=texit
         exitw2(is)=wexit
         qstak2(1,is)=qexit
         qstak2(2,is)=dtemp*vf0

50       continue

      enddo
c --- End Loop over sources

      return
      end
c----------------------------------------------------------------------
      subroutine julday(iyr,imo,iday,ijuldy)
c----------------------------------------------------------------------
c
c --- FGEMIS    Version: 1.0      Level: 981116                 JULDAY
c ---          J. Scire, SRC
c
c --- PURPOSE:  Compute the Julian day number from the Gregorian
c               date (month, day)
c
c --- INPUTS:
c           IYR - integer      - Year
c           IMO - integer      - Month
c          IDAY - integer      - Day
c        Parameters: IO6
c
c --- OUTPUT:
c          IJUL - integer      - Julian day
c
c --- JULDAY called by:  FIN, RDISC
c --- JULDAY calls:      none
c----------------------------------------------------------------------
c
c --- include parameters
      include 'params.fge'
c
      integer kday(12)
      data kday/0,31,59,90,120,151,181,212,243,273,304,334/
c
c --- Check for valid input data
      ierr=0
c --- Check for valid month
      if(imo.lt.1.or.imo.gt.12)ierr=1
c --- Check for valid day in 30-day months
      if(imo.eq.4.or.imo.eq.6.or.imo.eq.9.or.imo.eq.11)then
         if(iday.gt.30)ierr=1
      else if(imo.eq.2)then
         if(mod(iyr,4).eq.0)then
c ---       February in a leap year
            if(iday.gt.29)ierr=1
         else
c ---       February in a non-leap year
            if(iday.gt.28)ierr=1
         endif
      else
c ---    Check for valid day in 31-day months
         if(iday.gt.31)ierr=1
      endif
c
      if(ierr.eq.1)then
         write(io6,*)
         write(io6,*)'ERROR in SUBR. JULDAY'
         write(io6,*)'Invalid date - IYR = ',iyr,' IMO = ',
     1    imo,' IDAY = ',iday
         stop
      endif
c
c --- Compute the Julian day
      ijuldy=kday(imo)+iday
      if(imo.le.2)return
      if(mod(iyr,4).EQ.0)ijuldy=ijuldy+1
c
      return
      end
c----------------------------------------------------------------------
      subroutine grday(iyr,ijul,imo,iday)
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 940430                  GRDAY
c                J. Scire, SRC
c
c --- PURPOSE:  Compute the Gregorian date (month, day) from the
c               Julian day
c
c --- INPUTS:
c           IYR - integer      - Year
c          IJUL - integer      - Julian day
c     Parameters: IO6
c
c --- OUTPUT:
c           IMO - integer      - Month
c          IDAY - integer      - Day
c
c --- GRDAY called by:  COMP, RDISC
c --- GRDAY calls:      none
c----------------------------------------------------------------------
c
c --- Include parameter statements
      include 'params.fge'
c
      integer kday(12,2)
      data kday/31,59,90,120,151,181,212,243,273,304,334,365,
     1          31,60,91,121,152,182,213,244,274,305,335,366/
c
c
      ileap=1
      if(mod(iyr,4).eq.0)ileap=2
      if(ijul.lt.1.or.ijul.gt.kday(12,ileap))go to 11
c
      do 10 i=1,12
      if(ijul.gt.kday(i,ileap))go to 10
      imo=i
      iday=ijul
      if(imo.ne.1)iday=ijul-kday(imo-1,ileap)
      return
10    continue
c
11    continue
      write(io6,12)iyr,ijul
12    format(//2x,'ERROR in SUBR. GRDAY -- invalid Julian day '//2x,
     1 'iyr = ',i5,3x,'ijul = ',i5)
      stop
      end
c----------------------------------------------------------------------
      subroutine incr(iyr,ijul,ihr,nhrinc)
c----------------------------------------------------------------------
c
c --- FGEMISS    Version:  1.0     Level:  980304                  INCR
c                J. Scire, SRC
c
c --- PURPOSE:  Increment the time and date by "NHRINC" hours
c
c --- UPDATE
c --- V5.0-V5.0     980304  (DGS): Allow for a negative "increment" of
c                                  up to 24 hours
c                   980304  (DGS): Allow for arbitrarily large nhrinc
c
c --- INPUTS:
c       IYR    - integer - Current year
c       IJUL   - integer - Current Julian day
c       IHR    - integer - Current hour (00-23)
c       NHRINC - integer - Time increment (hours)
c
c               NOTE: "NHRINC" must >= -24
c                      Hour is between 00-23
c
c       Parameters: IO6
c
c --- OUTPUT:
c       IYR    - integer - Updated year
c       IJUL   - integer - Updated Julian day
c       IHR    - integer - Updated hour (00-23)
c
c --- INCR called by: COMP, RDISC
c --- INCR calls:     none
c----------------------------------------------------------------------
c
c --- Include parameter statements
      include 'params.fge'
c
c --- Check nhrinc
      if(nhrinc.lt.-24) then
         write(io6,*)'ERROR IN SUBR. INCR -- Invalid value of NHRINC ',
     1   '-- NHRINC = ',nhrinc
         stop
      endif

c --- Save increment remaining (needed if nhrinc > 8760)
      nleft=nhrinc
c
c --- Process change in hour
      if(nhrinc.gt.0)then
c
10       ninc=MIN0(nleft,8760)
         nleft=nleft-ninc
c
c ---    Increment time
         ihr=ihr+ninc
         if(ihr.le.23)return
c
c ---    Increment day
         ijul=ijul+ihr/24
         ihr=mod(ihr,24)
c
c ---    ILEAP = 0 (non-leap year) or 1 (leap year)
         if(mod(iyr,4).eq.0)then
            ileap=1
         else
            ileap=0
         endif
c
         if(ijul.gt.365+ileap) then
c ---       Update year
            iyr=iyr+1
            ijul=ijul-(365+ileap)
         endif
c
c ---    Repeat if more hours need to be added
         if(nleft.GT.0) goto 10
c
      elseif(nhrinc.lt.0)then
c ---    Decrement time
         ihr=ihr+nhrinc
         if(ihr.lt.0)then
            ihr=ihr+24
            ijul=ijul-1
            if(ijul.lt.1)then
               iyr=iyr-1
               if(mod(iyr,4).eq.0)then
                  ijul=366
               else
                  ijul=365
               endif
            endif
         endif
      endif
c
      return
      end
c----------------------------------------------------------------------
      subroutine fin
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 011107                    FIN
c ---            J. Scire, SRC
c
c --- PURPOSE:  Run termination routine -- compute runtime,
c               write last day processed
c
c --- INPUTS:
c       Common block /QA/
c          rdate, rtime, rcpu, iyr, month, iday, ijul, ihr
c       Parameters: IO6, IOMESG
c
c --- OUTPUT:  none
c
c --- FIN called by:  MAIN
c --- FIN calls:      DATETM, JULDAY, DELTT
c----------------------------------------------------------------------
c
c --- Include parameters
      include 'params.fge'
      include 'qa.fge'

      character*8 rtime2
      character*10 rdate2

      write(iomesg,*)'TERMINATION PHASE'

c --- Write last day/hour processed
      write(io6,5)iyr,month,iday,ijul,ihr
5     format(//2x,'LAST DAY/HOUR PROCESSED:'/5x,'Year: ',i4,2x,
     1   'Month: ',i2,3x,'Day: ',i2,3x,'Julian day: ',i3,3x,'Hour: ',
     2   i2)

c --- Get system date & time at end of run
      call datetm(rdate2,rtime2,rcpu)

c --- Compute runtime
      read(rtime(1:2),10)ihr1
      read(rtime(4:5),10)imin1
      read(rtime(7:8),10)isec1
10    format(i2)
      t1=ihr1*3600.+imin1*60.+isec1

      read(rtime2(1:2),10)ihr2
      read(rtime2(4:5),10)imin2
      read(rtime2(7:8),10)isec2
      t2=ihr2*3600.+imin2*60.+isec2

      if(rdate.eq.rdate2)then
         delt=t2-t1
      else
         read(rdate(1:2),10)imo1
         read(rdate(4:5),10)iday1
         read(rdate(7:10),'(i4)')iyr1
         call julday(iyr1,imo1,iday1,ijul1)

         read(rdate2(1:2),10)imo2
         read(rdate2(4:5),10)iday2
         read(rdate2(7:10),'(i4)')iyr2
         call julday(iyr2,imo2,iday2,ijul2)

c ---    Compute no. hours from beg. of first hour of run to
c ---    ending hour of ending day of the run
         call deltt(iyr1,ijul1,ihr1,iyr2,ijul2,ihr2,idelhr)

c ---    Adjust for minutes and seconds
         delt=idelhr*3600.-imin1*60.-isec1+imin2*60.+isec2
      endif

c --- On the PC, the runtime and CPU time are the same
c --- (DATETM provides RCPU = 0.0 on the PC)
      if(rcpu.EQ.0.0)rcpu=delt

      write(io6,1402)rtime2,rdate2,delt,rcpu
1402  format(//2x,'End of run -- Clock time: ',a8/
     1         2x,'                    Date: ',a10//
     2         2x,'      Elapsed Clock Time: ',f10.1,' (seconds)'//
     3         2x,'                CPU Time: ',f10.1,' (seconds)')
c
      return
      end
c----------------------------------------------------------------------
      subroutine deltt(j1yr,j1jul,j1hr,j2yr,j2jul,j2hr,jleng)
c----------------------------------------------------------------------
c
c --- FGEMISS    Version: 1.0      Level: 901030                  DELTT
c ---            J. Scire, SRC
c
c --- Compute the difference (in hours) between two dates & times
c ---    (time #2 - time #1)
c
c --- INPUTS:
c              J1YR - integer    - Year of date/time #1
c             J1JUL - integer    - Julian day of date/time #1
c              J1HR - integer    - Hour of date/time #1
c              J2YR - integer    - Year of date/time #2
c             J2JUL - integer    - Julian day of date/time #2
c              J2HR - integer    - Hour of date/time #2
c
c --- OUTPUT:
c             JLENG - integer    - Difference (#2 - #1) in hours
c
c --- DELTT called by:  FIN
c --- DELTT calls:      none
c----------------------------------------------------------------------
c
      jmin=min0(j1yr,j2yr)
c
c --- find the number of hours between Jan. 1 of the "base" year and
c --- the first date/hour
      if(j1yr.eq.jmin)then
         j1=0
      else
         j1=0
         j1yrm1=j1yr-1
         do 10 i=jmin,j1yrm1
         if(mod(i,4).eq.0)then
            j1=j1+8784
         else
            j1=j1+8760
         endif
10       continue
      endif
      j1=j1+(j1jul-1)*24+j1hr
c
c --- find the number of hours between Jan. 1 of the "base" year and
c --- the second date/hour
      if(j2yr.eq.jmin)then
         j2=0
      else
         j2=0
         j2yrm1=j2yr-1
         do 20 i=jmin,j2yrm1
         if(mod(i,4).eq.0)then
            j2=j2+8784
         else
            j2=j2+8760
         endif
20       continue
      endif
      j2=j2+(j2jul-1)*24+j2hr
c
c --- compute the time difference (in hours)
      jleng=j2-j1
c
      return
      end
c----------------------------------------------------------------------
      subroutine yr4(io,iyr,ierr)
c----------------------------------------------------------------------
c
c --- FGEMISS  Version: 1.0       Level: 011107                     YR4
c ---          D. Strimaitis    Earth Tech, Inc
c
c --- PURPOSE: Checks/converts 2-digit year to 4-digit year
c
c --- INPUTS:
c                IO - integer    - Unit for list file output
c               IYR - integer    - Year (YYYY or YY)
c
c     Common block /Y2K/:
c           IYYLO, ICCLO
c
c --- OUTPUT:
c               IYR - integer    - Year (YYYY)
c              IERR - integer    - Error code: 0=OK, 1=FATAL
c
c --- YR4 called by:  Input routines reading 'year' data
c --- YR4 calls:      none
c----------------------------------------------------------------------
      common /y2k/ iyylo,icclo

      ierr=0

c --- Test for 4-digit year (must exceed 1000)
      if(iyr.GT.1000) then
c ---    Passes 11th Century test (large year not trapped)
         return
      elseif(iyr.LT.100 .AND. iyr.GE.0) then
c ---    2-digit year
c ---    Construct 4-digit year
         if(iyr.LT.iyylo) then
            iyr=(icclo+1)*100+iyr
         else
            iyr=icclo*100+iyr
         endif
      else
c ---    Year not recognized
         ierr=1
         write(*,*)'ERROR in YR4 --- Year not recognized: ',iyr
         write(io,*)'ERROR in YR4 --- Year not recognized: ',iyr
      endif

      return
      end
c----------------------------------------------------------------------
      subroutine yr4c(iyr)
c----------------------------------------------------------------------
c
c --- FGEMISS  Version: 1.0       Level: 011107                    YR4C
c ---          D. Strimaitis    Earth Tech, Inc
c
c --- PURPOSE: Checks/converts 2-digit year to 4-digit year (CURRENT)
c
c --- INPUTS:
c               IYR - integer    - Year (YYYY or YY)
c
c --- OUTPUT:
c               IYR - integer    - Year (YYYY)
c
c --- YR4C called by:  DATETM
c --- YR4C calls:      none
c----------------------------------------------------------------------
c --- Set parameters for converting a current year (1999 - 2098)
c --- Use KCCLO as century digits for years GE KYYLO
      data kyylo/99/, kcclo/19/

c --- Test for 4-digit year (must exceed 1000)
      if(iyr.GT.1000) then
c ---    Passes 11th Century test (large year not trapped)
         return
      elseif(iyr.LT.100 .AND. iyr.GE.0) then
c ---    2-digit year
c ---    Construct 4-digit year
         if(iyr.LT.kyylo) then
            iyr=(kcclo+1)*100+iyr
         else
            iyr=kcclo*100+iyr
         endif
      else
c ---    Year not recognized
         write(*,*)'ERROR in YR4C --- Year not recognized: ',iyr
      endif

      return
      end
c----------------------------------------------------------------------
      subroutine qayr4(io,iyr,metrun,ierr)
c----------------------------------------------------------------------
c
c --- FGEMISS  Version: 1.0       Level: 011107                   QAYR4
c ---          D. Strimaitis    Earth Tech, Inc
c
c --- PURPOSE: Defines century and year markers to use in converting
c ---          2-digit year to 4-digit year
c ---          The IBYR (YYYY) must be provided in the control file
c
c --- INPUTS:
c                IO - integer    - Unit number for list file output
c               IYR - integer    - Year provided for start of run
c            METRUN - integer    - Flag to run period in met file
c                                  0 = do not run period
c                                  1 = run period
c
c --- OUTPUT:
c              IERR - integer    - Error code: 0=OK, 1=FATAL
c
c     Common block /Y2K/:
c           IYYLO, ICCLO
c
c --- QAYR4 called by:  READCF
c --- QAYR4 calls:      none
c----------------------------------------------------------------------
      common /y2k/ iyylo,icclo

c --- Sets parameters for the starting century marker (CC) and the
c --- 2-digit year (YY) used as the marker between the starting century
c --- and the next century.  For example, if CC=19 and YY=30, then a
c --- year less than 30 (say 15) is assumed to be 2015.  Any year
c --- greater than or equal to 30 (say 56) is assumed to be 1956.

c --- Set number of years prior to start of simulation that must not
c --- be placed in the next century
      data ibackyr/20/

      ierr=0

c --- Expect explicit starting year (YYYY)
c --- Test for 4-digit year (must exceed 1000)
      if(iyr.GT.1000) then
c ---    Passes 11th Century test (large year not trapped)
c ---    Back up IBACKYR years to set IYYLO
         kyr=iyr-ibackyr
c ---    Extract starting 2-digit century and 2-digit year
         icclo=kyr/100
         iyylo=kyr-icclo*100

c ---    Warn user that control file input is used to convert to YYYY
         iyr1=icclo*100+iyylo
         iyr2=(icclo+1)*100+iyylo-1
c
         write(io,*)
         write(io,*)'-------------------------------------------------'
         write(io,*)'NOTICE: Starting year in control file sets the'
         write(io,*)'        expected century for the simulation.  All'
         write(io,*)'        YY years are converted to YYYY years in'
         write(io,*)'        the range: ',iyr1,iyr2
         write(io,*)'-------------------------------------------------'
         write(io,*)
      else
         ierr=1
         write(*,*)
         write(*,*)'--------------------------------------------'
         write(*,*)'QAYR4 -- Start year must be 4-digits!: ',iyr
         if(metrun.EQ.1) then
            write(*,*)'         and must always be provided'
         endif
         write(*,*)'--------------------------------------------'
         write(*,*)
         write(io,*)
         write(io,*)'-------------------------------------------'
         write(io,*)'QAYR4 -- Start year must be 4-digits!: ',iyr
         if(metrun.EQ.1) then
            write(io,*)'         and must always be provided'
         endif
         write(io,*)'-------------------------------------------'
         write(io,*)
      endif

      return
      end
c----------------------------------------------------------------------
      subroutine datetm(rdate,rtime,rcpu)
c----------------------------------------------------------------------
c
c --- CALMET     Version: 5.3      Level: 010731!a               DATETM
c ---            J. Scire, SRC
c
c --- PURPOSE:  Get system date and time from system clock
c
c --- UPDATE
c --- V5.2-V5.3     010731!a (DGS): rdate changed to include YYYY
c                                  format for year, and YR4C moved in
c                                  to convert YY (if old format found)
c
c --- INPUTS:  none
c
c --- OUTPUT:  rdate  - C*10 - Current system date (MM-DD-YYYY)
c              rtime  - C*8  - Current system time (HH:MM:SS)
c               rcpu  - real - CPU time (sec) from system utility
c
c --- DATETM called by:  SETUP
c --- DATETM calls:      DATE, TIME (Lahey compiler system utilities
c                        providing date and time)
c                        ETIME (SUN CPU utility program)
c                        YR4C
c----------------------------------------------------------------------
      character*8  rtime
      character*10 rdate

c --- Local store
      character*11 stime
      character*8 sdate

c --- Get system date and time [COMPILER-SPECIFIC!]

c --- Lahey F77L Compiler (begin)
c -------------------------------
c --- System date in MM-DD-YY
c --- System clock in HH:MM:SS.ss, where ss = hundredths of seconds
      call DATE(sdate)
      call TIME(stime)
c --- Pass to output formats (MM-DD-YYYY) and (HH:MM:SS)
      rdate='  -  -    '
      rdate(1:2)=sdate(1:2)
      rdate(4:5)=sdate(4:5)
      rdate(7:8)='00'
      rdate(9:10)=sdate(7:8)
      do i=1,8
         rtime(i:i)=stime(i:i)
      enddo
c --- Lahey F77L Compiler (end)

cc --- COMPAQ DF90/95 Compiler (begin)
cc -----------------------------------
cc --- System date in CCYYMMDD
cc --- System clock in HHMMSS.sss, where sss = thousandths of seconds
c      call DATE_AND_TIME(sdate,stime)
cc --- Pass to output formats (MMDDYYYY) and (HH:MM:SS)
c      rdate='  -  -    '
c      rdate(1:2)=sdate(5:6)
c      rdate(4:5)=sdate(7:8)
c      rdate(7:10)=sdate(1:4)
c      rtime='  :  :  '
c      rtime(1:2)=stime(1:2)
c      rtime(4:5)=stime(3:4)
c      rtime(7:8)=stime(5:6)
cc --- COMPAQ DF90/95 Compiler (end)

c --- Construct 4-digit year from current 2-digit year (if found)
      read(rdate(7:10),'(i4)') iyr
      call YR4C(iyr)
      write(rdate(7:10),'(i4)') iyr

c --- Get CPU time from SUN system utility (or PC dummy)
      call etime(rcpu)

      return
      end
c----------------------------------------------------------------------
      subroutine etime(rcpu)
c----------------------------------------------------------------------
c
c --- CALMET     Version: 5.3      Level: 941215                  ETIME
c ---            J. Scire, SRC
c
c --- PURPOSE:  Dummy system CPU time routine for PC
c
c --- INPUTS:  none
c
c --- OUTPUT:  RCPU  - real - CPU time (sec) -- set to zero for PC
c
c --- ETIME called by:  DATETM
c --- ETIME calls:      none
c----------------------------------------------------------------------
      rcpu=0.0
c
      return
      end
c----------------------------------------------------------------------
      subroutine comline(ctext)
c----------------------------------------------------------------------
c
c --- CALMET     Version: 5.3      Level: 010731!a              COMLINE
c                J. Scire, SRC
c
c --- PURPOSE:  Call the compiler-specific system routine that will
c               pass back the command line argument after the text
c               that executed the program
c
c               This routine contains calls for several different
c               compilers, but only one should be active at any one
c               time.
c
c --- INPUTS:   none
c
c
c --- OUTPUT:
c
c          CTEXT - character*70  - Command line argument #1
c                                  If command line argument is
c                                  missing, CTEXT is assigned a
c                                  default of 'calpuff.inp'
c
c --- COMLINE called by:  SETUP
c --- COMLINE calls:      GETCL - Lahey system routine
c                         NARGS, GETARG - Microsoft routines
c                         NARGS, GETARG - Compaq DF routines
c                         IARGC, GETARG - Sun routines
c                         IARGC, GETARG - HP routines
c----------------------------------------------------------------------
c
c ----------------------------------------
c --- COMPAQ DF compiler library directive
c      USE DFLIB
c ----------------------------------------
c
      character*70 ctext,cdeflt
c
c --- Microsoft variables
c *** integer*2 iarg,istat
c
c --- HP declaration
c *** external getarg
c *** external iargc
c
      data cdeflt/'calmet.inp'/
c
c --- The following is for any system without a command line routine
c --- and is also used as a default
      ctext=cdeflt
c
c ------------------
c --- Lahey compiler
c ------------------
      call getcl(ctext)
c
c ----------------------
c --- COMPAQ DF compiler
c ----------------------
c ***  numargs=NARGS()
c ***  if(numargs.ge.1)then
c ***     call GETARG(1,ctext)
c ***  endif
c
c ----------------------
c --- Microsoft compiler
c ----------------------
c *** numargs=NARGS()
c *** if(numargs.ge.1)then
c ***    iarg=1
c ***    call GETARG(iarg,ctext,istat)
c *** endif
c
c ----------------
c --- Sun compiler
c ----------------
c *** numargs=IARGC()
c *** if(numargs.ge.1)then
c ***    call GETARG(1,ctext)
c *** endif
c
c -----------------------------------------------------
c --- HP compiler -- NOTE: needs +U77 switch on compile
c -----------------------------------------------------
c *** numargs=IARGC()
c *** if(numargs.ge.1)then
c ***    call GETARG(1,ctext)
c *** endif
c

c --- If no command line arguments, use default
      if(ctext(1:1).eq.' ')ctext=cdeflt

      return
      end
