c----------------------------------------------------------------------
c --- CTEMISS -- Emissions Processor for COOLING TOWER input
c             -- to CALPUFF via PTEMARB.DAT file
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 061020                   MAIN
c                D. Strimaitis
c
c     Copyright (c) 1999-2006 by TRC Environmental Corporation
c
c --- PURPOSE:  Compute hourly water vapor (g/s) and temperature (K)
c               at the exit of a mechanical-draft cooling tower as a
c               function of the hourly relative humidity and temperature
c               with options for 
c
c               1) Heat release (BTU/hr) varies by either:
c                    month
c                    temperature
c
c               2) % BTUs contributed by DRY section (abated units)
c                  varies by either:
c                    month
c                    temperature
c
c               3) Exit velocity (m/s) varies by either:
c                    month
c                    temperature
c
c               4) Monthly variability in
c                    % ambient air drawn into DRY section (abated units)
c
c               Pass all inputs required for CALPUFF variable emissions
c               file to output
c
c
c --- UPDATE:
c     5.2(031013) to 5.3(061020) --- DGS
c          (1)   Modify structure of control file to include coordinate
c                MAP, DATUM, and Time Zone, and specify simulation
c                period by true starting and ending times (HHSSSS)
c          (2)   Write PTEMARB file as dataset 2.1
c          (3)   Retain hourly time structure (require start/end seconds
c                to be ZERO or 3600)
c          (4)   Drop option for BINARY output file
c                Added:     MAP.CTE
c                           MIDNITE
c                Modified:  PT2.CTE
c                           READCF, LISTINP, WRHDEM2, WRTIEM2, WREM2,
c                           COMP, FIN
c
c     5.1(021031) to 5.2(031013) --- DGS
c          (1)   Add Y2K subroutines, and replace related date/time
c                routines.  These are taken from the latest:
c                CALUTILS   Version: 2.1      Level: 030402
c          (2)   Renamed subroutine DIAG to DIAGUP to remove conflict
c                with /DIAG/ reported by F95 compiler
c
c     5.1(020503) to 5.1(021031) --- DGS
c          (1)   Refine tolerance on iterations in CTFLUXA to avoid endless loop
c
c     5.0(000203) to 5.1(020503) --- DGS
c          (1)   Add recirculation feature, where R% is specified by wind sector
c                for each cell (like building info) in new Section 7
c                Building data are now in Section 8.
c          (2)   Recirculation increases temperature and emitted water flux due
c                to altered 'ambient' conditions
c          (3)   Increment in ambient water flux due to recirculation is added
c                to computed flux
c
c --- UPDATE:
c     5.0(991217) to 5.0(000203) --- DGS
c          (1)   Exit velocity varies by EITHER month or temperature
c
c --- UPDATE:
c     5.0(990604) to 5.0(991217) --- DGS
c          (1)   Establish ctemiss.dat as default input file name
c
c --- UPDATE:
c     5.0(990325) to 5.0(990604) --- BdF
c          (1)   If section 7 is missing from input file,
c                load data from BPIP file (interactive prompt)
c                and append new data to input file
c          (2)   Bug fixed in section 7 to correctly handle downwash flag
c                User specified cell is used if matching cell not found
c
c --- UPDATE:
c     5.0(990225) to 5.0(990325) --- DGS
c          (1)  Allow zero % ambient air drawn into DRY section to
c               override any non-zero % BTUs contributed by DRY section
c               so that tower becomes a pure wet tower if the air
c               intakes for the dry section are closed.
c----------------------------------------------------------------------
c
c --- Include parameters
      include 'params.cte'
      include 'qa.cte'

c --- set version and level number of program
      ver='5.3'
      level='061020'

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 --- CTEMISS    Version: 5.3      Level: 990204                  SETUP
c                D. Strimaitis
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(8,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,  DIAGINI, WRHDEM2, WRTIEM2,
c                      
c----------------------------------------------------------------------
c
c --- Include parameters
      include 'params.cte'

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

      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(cteinp)

c --- Open the control file
      open(io5,file=cteinp,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 --- Initialize diagnostic variables
      call DIAGINI

c --- Read ISC header record (station,year for surface and upper)
      read(io7,'(a80)') 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 --- CTEMISS    Version: 5.3      Level: 061020                   COMP
c                D. Strimaitis
c
c --- PURPOSE:  Compute emissions for each hour in the period
c
c --- UPDATE:
c     5.2(031013) to 5.3(061020) --- DGS
c          (1)   Dataset 2.1 format: True Start/End times as HH SSSS
c
c     5.1(020503) to 5.2(031013) --- DGS
c          (1)   Renamed subroutine DIAG to DIAGUP to remove conflict
c                with /DIAG/ reported by F95 compiler
c
c     5.0(990125) to 5.1(020503) --- DGS
c          (1)   Pass wind flow direction to CTFLUXA for recirculation
c                feature (direction-dependent)
c
c
c --- INPUTS:
c     Parameters:
c           MXSPEC, MXPT2, IO6, IO7, IO16 
c     Common block /PT2/:
c           NPT2,
c           IBDAT2, IBHR2, IBSEC2, IEDAT2, IEHR2, IESEC2
c           TIEM2(8,mxpt2),
c           cid(mxpt2),bturate(mxpt2),mbtu(mxpt2),btufac12(12,mxpt2),
c           wexit12(12,mxpt2),btudry12(12,mxpt2),airdry12(12,mxpt2),
c           rhmin
c
c --- OUTPUT:
c     Common block /PT2/:
c           TSTAK2(mxpt2),EXITW2(mxpt2),QSTAK2(mxspec,mxpt2)
c     Common block /QA/:
c           iyr, month, iday, ijul, ihr, isec
c
c --- COMP called by:  MAIN
c --- COMP calls:      GRDAY, CTFLUXA, WREM2, INCR, DIAGUP, DIAGOUT
c----------------------------------------------------------------------
c
c --- Include parameter statements
      include 'params.cte'

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

      logical leof,ldb
      data ldb/.FALSE./

c --- Set starting date-hour (time at START of first hour)
      ibdathr=ibdat2*100+ibhr2
c --- Set ending date-hour (time at END of last hour)
      idathr0=iedat2*100+iehr2

c --- Set the output seconds so that each hour goes from
c --- 0000 to 3600
      ibsec=0
      isec=3600

c --- Store starting period time
      idate=ibdathr/100
      ibhr=ibdathr-idate*100
      ibyr =idate/1000
      ibjul=idate-ibyr*1000
c --- Move time marker to the END of the first hour
c --- Increment date-time by 1 hour
      nhrinc=1
      iyr=ibyr
      ijul=ibjul
      ihr=ibhr
      call INCR(io6,iyr,ijul,ihr,nhrinc)
      idathr=iyr*100000+ijul*100+ihr

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

      write(iomesg,*)'COMPUTATIONAL PHASE'

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

c --- Read current temperature, relative humidity, and flow direction
10    call RDISC(idathr,tmpk,rhamb,flow,kdathr,leof)
c --- Set flow direction bin index (1-36)
      if(flow.LT.5.0) then
         iflow=36
      else
         iflow=INT(0.1*(flow-5.))+1
      endif
      if(iflow.GT.36) iflow=36
      
      rhamb=AMAX1(rhamb,rhmin)
      if(LEOF) then
         write(io6,*)
         write(io6,*) 'EOF reached in MET or CHEM data file'
         write(io6,*) 'Last period attempted YYJJJHH =',idathr
         write(io6,*)
         write(*,*)
         write(*,*) 'EOF reached in MET or CHEM data file'
         write(*,*) 'Last period attempted YYJJJHH =',idathr
         write(*,*)
         return
      endif

c --- Set month
      call GRDAY(io6,iyr,ijul,month,iday)

c --- Assume std surface pressure 
      pamb=1013.0

c --- Compute cooling tower flux for each cell this hour
      call CTFLUXA(iflow,tmpk,pamb,rhamb,month)

c --- Write record for this hour
      call WREM2(ibyr,ibjul,ibhr,ibsec,ibyr,ibjul,ibhr,isec)

c --- Update diagnostic arrays
      call DIAGUP(tmpk,pamb,rhamb,ihr,month,irhplm)

c --- Extra DEBUG output
      if(LDB) then
         if(irhplm.LT.100 .AND. btudry12(month,1).LE.0.1) then
            write(io6,'(a28,i10,3i4)')'DEBUG: dathr,mo,day,PlumeRH = ',
     &                                 idathr,month,iday,irhplm
         endif
      endif

c --- Save start of current hour for reporting at end
      kyr=ibyr
      kjul=ibjul
      khr=ibhr
      ksec=ibsec
      call GRDAY(io6,kyr,kjul,kmonth,kday)
c --- Assign start of next hour to end of current hour
      ibyr=iyr
      ibjul=ijul
      ibhr=ihr
c --- Increment date-time for next hour
      nhrinc=1
      call INCR(io6,iyr,ijul,ihr,nhrinc)
      idathr=iyr*100000+ijul*100+ihr

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

c --- DONE

c --- Output diagnostics to list file
      call DIAGOUT

      return
      end
c----------------------------------------------------------------------
      subroutine openot
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 990125                 OPENOT
c                D. Strimaitis
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.cte'
c
c --- Include common blocks
      include 'filnam.cte'
      include 'pt2.cte'

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 --- CTEMISS    Version: 5.3      Level: 061020                 READCF
c                D. Strimaitis
c
c --- PURPOSE:  Read the structured control file (IO5)
c
c --- UPDATE:
c     5.2(031013) to 5.3(061020) --- DGS
c               Add dataset version as first record
c               Add coordinate info and time zone
c               Start/End times as HHSSSS
c               Drop option for Binary output file
c     5.1(020503) to 5.2(031013) --- DGS
c               Add Y2K structures
c     5.0(000203) to 5.1(020503) --- DGS
c               Add recirculation % by flow direction as Section 7, and
c               move building data to Section 8
c     5.0(990223) to 5.0(000203) --- DGS
c               Exit velocity varies by EITHER month or temperature
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           ctelst,iscdat,ptdat
c     Common block /MAP/:
c           lutm,llcc,lps,lem,llaza,lttm,
c           iutmzn,feast,fnorth,
c           pmap,utmhem,datum,daten,clat0,clon0,clat1,clat2
c     Common block /PT2/:
c           NPT2,
c           FNAME2, NSE2, IUTMZ2, IBDAT2, IBHR2, IBSEC2,
c           IEDAT2, IEHR2, IESEC2, VRS2, LABEL2, CSLST2(mxspec),
c           XMWEM2(mxspec),
c           TIEM2(8,mxpt2), BHT2(36,mxpt2),BWD2(36,mxpt2),
c           cid(mxpt2),bturate(mxpt2),mbtu(mxpt2),mbtudry(mxpt2),
c           mwexit(mxpt2),
c           btufac12(12,mxpt2),wexit12(12,mxpt2),btudry12(12,mxpt2),
c           airdry12(12,mxpt2),tbin(12),recirc(36,mxpt2),rhmin,
c           ATZONE
c
c --- READCF called by:  SETUP
c --- READCF calls:      QAYR4, YR4
c----------------------------------------------------------------------
c
c --- Include parameter statements
      include 'params.cte'

c --- Include common blocks
      include 'filnam.cte'
      include 'map.cte'
      include 'pt2.cte'
      include 'qa.cte'

c --- Control file version
      character*16 dataset,dataver

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

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

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 --- Recirculation (local)
      integer nptcirc
      real recircin(mxpt2,ndir)

c --- Set coordinate map logicals
      lutm=.FALSE.
      llcc=.FALSE.
      lttm=.FALSE.
      lem=.FALSE.
      llaza=.FALSE.
      lps=.FALSE.

c --- Verify correct type of control file
      read(io5,'(2a16)') dataset,dataver
      if(dataset(1:11).NE.'CTEMISS.INP') then
         write(*,*)'Invalid control file type found'
         write(*,*)'Expected dataset = ','CTEMISS.INP     '
         write(*,*)'Found    dataset = ',dataset
         stop
      endif
      if(dataver(1:3).NE.'2.0') then
         write(*,*)'Invalid control file version found'
         write(*,*)'Expected dataver = ','2.0             '
         write(*,*)'Found    dataver = ',dataver
         stop
      endif

c --- Skip to Section 1
      isect=1
1     read(io5,'(a16)',end=999) asection
      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,*) ctelst
      read(io5,*) ptdat
      
c --- Open the output LIST file
      open(io6,file=ctelst,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 --- Remove BINARY option
      lascii=.TRUE.
cc --- Set format of variable emissions file
c      read(io5,*) achar
c      if(achar.EQ.'A') then
c         lascii=.TRUE.
c      elseif(achar.EQ.'B') then
c         lascii=.FALSE.
c      else
c         write(*,*)'READCF:  PTEMARB format must be A,B : ',achar
c         write(io5,*)'READCF:  PTEMARB format must be A,B : ',achar
c         stop
c      endif

c --- Minimum RH (%) (avoid zeros in met file)
      read(io5,*) rhmin

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

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

c --- Special QA on starting year of simulation
      metrun=0
      ibyr=ibdat2/1000
      ibjul=ibdat2-ibyr*1000
      call QAYR4(io6,ibyr,metrun,ierrb)
      if(ierrb.NE.0) then
         write(*,*)
         stop 'CORRECT THE CONTROL FILE'
      endif
      ibdat2=ibyr*1000+ibjul
c --- Enforce YYYY format for end-year
      ieyr=iedat2/1000
      iejul=iedat2-ieyr*1000
      call YR4(io6,ieyr,ierr)
      if(ierr.NE.0) then
         write(*,*)
         stop 'Halted in READCF'
      endif
      iedat2=ieyr*1000+iejul

c --- Condition times
      ibhr2=ibtim2/10000
      ibsec2=ibtim2-ibhr2*10000
      iehr2=ietim2/10000
      iesec2=ietim2-iehr2*10000
c --- Check for hour=24
      nhrinc=1
      if(ibhr2.EQ.24) then
         ibhr2=23
         call INCR(io6,ibyr,ibjul,ibhr2,nhrinc)
      endif
      if(iehr2.EQ.24) then
         iehr2=23
         call INCR(io6,ieyr,iejul,iehr2,nhrinc)
      endif
c --- Check for sec=3600
      nhrinc=1
      if(ibsec2.EQ.3600) then
         call INCR(io6,ibyr,ibjul,ibhr2,nhrinc)
         ibsec2=0
      endif
      if(iesec2.EQ.3600) then
         call INCR(io6,ieyr,iejul,iehr2,nhrinc)
         iesec2=0
      endif
c --- Reset dates (may have changed)
      ibdat2=ibyr*1000+ibjul
      iedat2=ieyr*1000+iejul

c --- Times must be hourly (seconds are ZERO)
      if(ibsec2.NE.0 .OR. iesec2.NE.0) then
         write(*,*)'READCF:  Times must be hourly (ZERO seconds)'
         write(*,*)'IBSEC,IESEC = ',ibsec2,iesec2
         stop
      endif

c --- Coordinate information
      read(io5,'(a8)') pmap
      pmap=ADJUSTL(pmap)
      read(io5,'(a8)') datum
      datum=ADJUSTL(datum)
      read(io5,'(i4)') iutmzn
      read(io5,'(a4)') utmhem
      utmhem=ADJUSTL(utmhem)
      read(io5,'(a16)') clat0
      clat0=ADJUSTL(clat0)
      read(io5,'(a16)') clon0
      clon0=ADJUSTL(clon0)
      read(io5,'(a16)') clat1
      clat1=ADJUSTL(clat1)
      read(io5,'(a16)') clat2
      clat2=ADJUSTL(clat2)
      read(io5,*) feast
      read(io5,*) fnorth

      if(pmap.EQ.'UTM     ') then
         lutm=.TRUE.
      elseif(pmap.EQ.'TTM     ') then
         lttm=.TRUE.
      elseif(pmap.EQ.'LCC     ') then
         llcc=.TRUE.
      elseif(pmap.EQ.'EM      ') then
         lem=.TRUE.
      elseif(pmap.EQ.'LAZA    ') then
         llaza=.TRUE.
      elseif(pmap.EQ.'PS      ') then
         lps=.TRUE.
      else
         write(*,*)'READCF:  Invalid map projection: ',pmap
         stop
      endif

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
      isect=2
2     read(io5,'(a16)',end=999) asection
      if(asection.NE.asectdb(isect)) goto 2
      read(io5,*)
c
c ----------------------------------------------
c --- Constant Source Data -- Sect # 2
c ----------------------------------------------
c
c --- Loop over cells
      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,*)bturate(ic)          ! Heat(BTU/hr)
         read(io5,*)mbtu(ic)             ! BTU/hr Mode
         read(io5,*)mbtudry(ic)          ! %Dry BTU/hr Mode
         read(io5,*)mwexit(ic)           ! Exit Velocity Mode
         read(io5,*)tiem2(6,ic)          ! Building DW
         tiem2(7,ic)=1.                  ! Vertical momentum flux factor
         tiem2(8,ic)=0.                  ! User
      enddo
c --- Check BTU/hr, %DRY BTU/hr, and wexit Modes
      do ic=1,npt2
         ibtu=mbtu(ic)
         if(ibtu.LT.1 .OR. ibtu.GT.2) then
            write(*,*)'READCF:  Invalid MBTU (must be 1,2) = ',ibtu
            stop
         endif
         idry=mbtudry(ic)
         if(idry.LT.1 .OR. idry.GT.2) then
            write(*,*)'READCF:  Invalid MBTUDRY (must be 1,2) = ',idry
            stop
         endif
         iw=mwexit(ic)
         if(iw.LT.1 .OR. iw.GT.2) then
            write(*,*)'READCF:  Invalid MWEXIT (must be 1,2) = ',iw
            stop
         endif
      enddo

c --- Skip to Section 3
      isect=3
3     read(io5,'(a16)',end=999) asection
      if(asection.NE.asectdb(isect)) goto 3
      read(io5,*)
c
c ----------------------------------------------
c --- 12 Exit Velocities -- Sect # 3
c ----------------------------------------------
c
c --- Loop over cells
      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 cell name'
         read(io5,*) (wexit12(k,ic),k=1,12)
      enddo

c --- Skip to Section 4
      isect=4
4     read(io5,'(a16)',end=999) asection
      if(asection.NE.asectdb(isect)) goto 4
      read(io5,*)
c
c ----------------------------------------------
c --- 12 DRY BTU % Values -- Sect # 4
c ----------------------------------------------
c
c --- Loop over cells
      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 cell name'
         read(io5,*) (btudry12(k,ic),k=1,12)
      enddo

c --- Skip to Section 5
      isect=5
5     read(io5,'(a16)',end=999) asection
      if(asection.NE.asectdb(isect)) goto 5
      read(io5,*)
c
c ----------------------------------------------
c --- Monthly DRY Air Flow % -- Sect # 5
c ----------------------------------------------
c
c --- Loop over cells
      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 cell name'
         read(io5,*) (airdry12(k,ic),k=1,12)
      enddo

c --- Skip to Section 6
      isect=6
6     read(io5,'(a16)',end=999) asection
      if(asection.NE.asectdb(isect)) goto 6
      read(io5,*)
c
c ----------------------------------------------------
c --- BTU/hr Heat Rejection Rate Factors -- Sect # 6
c ----------------------------------------------------
c
c --- Loop over cells
      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 cell name'
         read(io5,*) (btufac12(k,ic),k=1,12)
      enddo

c --- Skip to Section 7
      isect=7
7     read(io5,'(a16)',end=999) asection
      if(asection.NE.asectdb(isect)) goto 7
      read(io5,*)
c
c ----------------------------------------------
c --- Recirculation Data -- Sect # 7
c ----------------------------------------------
c
c --- Read all recirculation information present
      nptcirc=0
      do i=1,mxpt2
         read(io5,*,end=994) stkn(i)
         if (stkn(i).eq.'--------') goto 994
         read(io5,*,end=994) (recircin(i,k),k=1,36)
         nptcirc=nptcirc+1
      enddo

994   continue

c --- QA and map data to correct cells
c --- Loop over cells
      do i=1,npt2
c ---    Match with CID array
         write(stkntemp,'(a8)') cid(i)
         ic=0
         do j=1,nptcirc
            if(stkn(j).EQ.stkntemp) ic=j
         enddo
         if(ic.GT.0) then
            do k = 1,36
               if(recircin(ic,k).GT.100.0) then
                  write(*,*) 'Recirculation exceeds 100%'
                  write(*,*) 'for cell name: ',cid(i)
                  write(*,*) '36 values read are:'
                  do j = 1,nptcirc
                     write(*,*) recircin(ic,j)
                  enddo
                  stop
               endif
               recirc(k,i) = 0.01*recircin(ic,k)
            enddo
         else
            do k = 1,36
               recirc(k,i) = 0.0
            enddo
         endif
      enddo

c --- Skip to Section 8
      isect=8
8     read(io5,'(a16)',end=997) asection
      if(asection.NE.asectdb(isect)) goto 8
      read(io5,*)
c
c ----------------------------------------------
c --- Building Data -- Sect # 8
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 cells and corresponding arrays
c --- Loop over cells
      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
            write(*,*) 'Please enter matching cell number'
		write(*,*) 'Current cell 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 8 taken from input file'
	write(sec8header,288)

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

997	continue

c
c ----------------------------------------------
c --- Retrieving Data from BPIP file -- Sect # 8 -- added by BdF 960604
c ----------------------------------------------
c
c If section 8 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 8 is then appended to the input file for future use.	
c   There may be extra gaps between section 7 & 8 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
9     read(io3,'(a25)',end=997) bpipstart
      isect=8
      if(bpipstart.NE.' BPIP output is in meters') goto 9

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. cells in BPIP file: ',npt2bpip,
     &	' No. cells in ip file: ',npt2

c map data to correct cells and corresponding arrays
c --- Loop over cells
      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 cell # for cell: ',cid(i)
             write(*,*)'from the list above'
             write(*,*)
             read(*,*) ic
	    else
               write(*,*) 'matching cell ',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 8 in original input file
c NB: there may be empty space between the end of section 7 and the beginning
c of section 8, although this won't matter for this program
	write(*,*) 'Appending BPIP data to input file'
c write section 8 header
	write(io5,*)
	write(io5,280)
	write(io5,281)
	write(io5,287) bpipinp
	write(sec8header,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 8 ---    Building Data by flow ',
     &	'direction (Repeat for each cell)')
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 --- CTEMISS    Version: 5.3      Level: 061020                LISTINP
c                D. Strimaitis
c
c --- PURPOSE:  Write inputs to list file
c
c --- UPDATE:
c     5.1(020503) to 5.3(061020) --- DGS
c               Add coordinate info and time zone
c               Start/End times as HHSSSS
c     5.0(000203) to 5.1(020503) --- DGS
c               Add recirculation array
c     5.0(990125) to 5.0(000203) --- DGS
c               Exit velocity varies by EITHER month or temperature
c
c
c --- INPUTS:
c     Parameters:
c           MXSPEC, IO6,          
c
c     Common block /FILNAM/:
c           ctelst,iscdat,ptdat
c     Common block /MAP/:
c           lutm,llcc,lps,lem,llaza,lttm,
c           iutmzn,feast,fnorth,
c           pmap,utmhem,datum,daten,clat0,clon0,clat1,clat2
c     Common block /PT2/:
c           NPT2,
c           FNAME2, NSE2, IUTMZ2, IBDAT2, IBHR2, IBSEC,
c           IEDAT2, IEHR2, IESEC, VRS2, LABEL2, CSLST2(mxspec),
c           XMWEM2(mxspec),
c           TIEM2(8,mxpt2), BHT2(36,mxpt2),BWD2(36,mxpt2),
c           cid(mxpt2),bturate(mxpt2),mbtu(mxpt2),btufac12(12,mxpt2),
c           wexit12(12,mxpt2),mwexit(mxpt2),
c           btudry12(12,mxpt2),airdry12(12,mxpt2),recirc(36,mxpt2),
c           rhmin, ATZONE
c
c --- OUTPUT:
c
c --- LISTINP called by:  SETUP
c --- LISTINP calls:      none
c----------------------------------------------------------------------
c
c --- Include parameter statements
      include 'params.cte'

c --- Include common blocks
      include 'filnam.cte'
      include 'map.cte'
      include 'pt2.cte'
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: ',ctelst
      write(io6,'(a16,a24)')'Emissions File: ',ptdat
      if(LASCII) then
         write(io6,*)'                 (ASCII)'
      else
         write(io6,*)'                 (BINARY)'
      endif
      write(io6,*)'Minimum RH (%): ',rhmin
      write(io6,*)
      write(io6,*)'         Start: ',ibdat2,ibhr2,ibsec2
      write(io6,*)'           End: ',iedat2,iehr2,iesec2
      write(io6,*)'     Time Zone: ',atzone
      write(io6,*)
      write(io6,*)'Coordinate Map Projection: ',pmap
      write(io6,*)'         Datum: ',datum
      if(LUTM) then
         write(io6,*)'      UTM zone: ',iutmzn,utmhem
      else
         write(io6,*)'Origin Lat,Lon: ',clat0,clon0
         write(io6,*)'     Lat1,Lat2: ',clat1,clat2
         write(io6,*)' False Easting: ',feast
         write(io6,*)'False Northing: ',fnorth
      endif
      write(io6,*)
      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 cells
      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)
         write(io6,*)' Base Heat(BTU/hr): ',bturate(ic)
         if(mbtu(ic).EQ.1) then
            write(io6,*)'       BTU/hr Mode: By Month'
         else
            write(io6,*)'       BTU/hr Mode: By Temperature'
         endif
         if(mbtudry(ic).EQ.1) then
            write(io6,*)'  %Dry BTU/hr Mode: By Month'
         else
            write(io6,*)'  %Dry BTU/hr Mode: By Temperature'
         endif
         if(mwexit(ic).EQ.1) then
            write(io6,*)'Exit Velocity Mode: By Month'
         else
            write(io6,*)'Exit Velocity Mode: By Temperature'
         endif
         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
         if(mwexit(ic).EQ.1) then
            cr='Months'
         else
            cr='T pnts'
         endif
         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,*)' --- DRY BTU (%) Variation -- Sect # 4'
      write(io6,*)' ----------------------------------------------'
c --- Loop over cells
      do ic=1,npt2
         if(mbtudry(ic).EQ.1) then
            cr='Months'
         else
            cr='T pnts'
         endif
         write(io6,*)
         write(io6,*)'     Source: ',cid(ic)
         write(io6,'(a6,a7,6f7.2)')cr,' 1-6 : ',(btudry12(k,ic),k=1,6)
         write(io6,'(a6,a7,6f7.2)')cr,' 7-12: ',(btudry12(k,ic),k=7,12)
      enddo

      write(io6,*)
      write(io6,*)' ----------------------------------------------'
      write(io6,*)' --- Monthly  DRY  Air Flow (%) -- Sect # 5'
      write(io6,*)' ----------------------------------------------'
c --- Loop over cells
      do ic=1,npt2
         write(io6,*)
         write(io6,*)'     Source: ',cid(ic)
         write(io6,'(a13,6f7.2)')'Months 1-6 : ',(airdry12(k,ic),k=1,6)
         write(io6,'(a13,6f7.2)')'Months 7-12: ',(airdry12(k,ic),k=7,12)
      enddo

      write(io6,*)
      write(io6,*)' ----------------------------------------------'
      write(io6,*)' --- BTU/hr Heat Rejection Factors -- Sect # 6'
      write(io6,*)' ----------------------------------------------'
c --- Loop over cells
      do ic=1,npt2
         if(mbtu(ic).EQ.1) then
            cr='Months'
         else
            cr='T pnts'
         endif
         write(io6,*)
         write(io6,*)'     Source: ',cid(ic)
         write(io6,'(a6,a7,6f7.4)')cr,' 1-6 : ',(btufac12(k,ic),k=1,6)
         write(io6,'(a6,a7,6f7.4)')cr,' 7-12: ',(btufac12(k,ic),k=7,12)
      enddo

      write(io6,*)
      write(io6,*)' ----------------------------------------------'
      write(io6,*)' --- Recirculation Data -- Sect # 7'
      write(io6,*)' ----------------------------------------------'
c --- Loop over cells
      do ic=1,npt2
         write(io6,*)
         write(io6,*)'     Source: ',cid(ic)
         write(io6,'(a13,6f7.2)')'Range  1-6 : ',(recirc(k,ic),k=1,6)
         write(io6,'(a13,6f7.2)')'Range  7-12: ',(recirc(k,ic),k=7,12)
         write(io6,'(a13,6f7.2)')'Range 13-18: ',(recirc(k,ic),k=13,18)
         write(io6,'(a13,6f7.2)')'Range 19-24: ',(recirc(k,ic),k=19,24)
         write(io6,'(a13,6f7.2)')'Range 25-30: ',(recirc(k,ic),k=25,30)
         write(io6,'(a13,6f7.2)')'Range 31-36: ',(recirc(k,ic),k=31,36)
      enddo

      write(io6,*)
      write(io6,*)' ----------------------------------------------'
      write(io6,*)' --- Building Downwash Data -- Sect # 8'
c      write(io6,*)' ----------------------------------------------'
	write(io6,'(a80)') sec8header
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,*) 'No. of sources with downwash: ',icount
	write(io6,*) 'No. of sources without downwash: ',npt2-icount
      return
      end
c----------------------------------------------------------------------
      subroutine rdisc(idathr,tmpk,rh,fvec,kdathr,leof)
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 031013                  RDISC
c                D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Read ISCMET.DAT meteorological data for one hour
c
c --- UPDATE:
c     5.1(020503) to 5.2(031013) --- DGS
c          (1)   Add Y2K check
c     5.0(990125) to 5.1(020503) --- DGS
c          (1)   Add wind flow direction to output for recirculation
c                feature (direction-dependent)
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                   FVEC - real    - flow direction (degrees)
c                 KDATHR - integer - Date of hour of data (YYJJJHH)
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.cte'
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 --- Enforce YYYY format for year
      call YR4(io6,iy,ierr)
      if(ierr.NE.0) then
         write(*,*)
         stop 'Halted in RDISC'
      endif

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(io6,iy,im,id,ijul)
         call INCR(io6,iy,ijul,ih,1)
         call GRDAY(io6,iy,ijul,im,id)
      else
         call JULDAY(io6,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 --- CTEMISS    Version: 5.3      Level: 061020                WRHDEM2
c                D. Strimaitis
c
c --- PURPOSE:  Write the header records for the PTEMARB point source
c               emissions file (arbitrarily-varying emissions)
c
c --- UPDATE:
c     5.0(990125) to 5.3(061020) --- DGS
c           -   Dataset 2.1 format
c           -   Drop option for Binary output file
c           -   Use MIDNITE to condition form of time reported
c
c --- INPUTS:
c            IUNIT - integer    - Fortran unit number of PTEMARB file
c     Common block /MAP/:
c           lutm,llcc,lps,lem,llaza,lttm,
c           iutmzn,feast,fnorth,
c           pmap,utmhem,datum,daten,clat0,clon0,clat1,clat2
c     Common Block /PT2/ variables:
c            NPT2, LASCII,
c            IBDAT2, IBHR2, IBSEC2
c            IEDAT2, IEHR2, IESEC2
c     Common Block /QA/ variables:
c            ver,level
c     Parameters:
c            IO6, MXSPEC
c
c
c --- OUTPUT:
c            none
c
c --- WRHDEM2 called by: SETUP
c --- WRHDEM2 calls:     MIDNITE
c----------------------------------------------------------------------
c
      include 'params.cte'
      include 'map.cte'
      include 'pt2.cte'
      include 'qa.cte'

      character*12 daten
      character*16 dset,dver
      character*64 dmod

      data dset/'PTEMARB.DAT     '/
      data dver/'2.1             '/
      data dmod/'Comments, times with seconds, time zone, coord info'/
      data daten/'02-21-2003  '/

      data ncom/1/
      data zero/0.0/

c --- Set species for FOG modeling
      nse2=2
      cslst2(1)='H2O'
      cslst2(2)='TXS'
      xmwem2(1)=1.0
      xmwem2(2)=1.0

c --- Extract year and Julian day from dates
      ibyr=ibdat2/1000
      ibjul=ibdat2-ibyr*1000
      ieyr=iedat2/1000
      iejul=iedat2-ieyr*1000

c --- Condition format of times written
      kbyr=ibyr
      kbjul=ibjul
      keyr=ieyr
      kejul=iejul
      kbhr=ibhr2
      kbsec=ibsec2
      kehr=iehr2
      kesec=iesec2
c --- Beginning time if midnight should be HR=00
      if(ibhr2.EQ.24) then
         kbhr=0
         call MIDNITE(io6,'TO 00h',ibyr,zero,zero,ibjul,
     &                             kbyr,zero,zero,kbjul)
      endif
c --- Ending time if midnight should be HR=24
      if(iehr2.EQ.0 .AND. iesec2.EQ.0) then
         kehr=24
         call MIDNITE(io6,'TO 24h',ieyr,zero,zero,iejul,
     &                             keyr,zero,zero,kejul)
      endif


c --- Write records in ASCII form
c -----------------------------------
c --- Record 1:  Dataset, Version, Modifier
      write(iunit,'(2a16,a64)') dset,dver,dmod
c --- Record 2:  Number of comment records
      write(iunit,'(i1)') ncom
c --- Record 3:  Comment (optional/repeatable)
      write(iunit,'(a59)') 'Created by CTEMISS Version: '//ver//
     &                     'Level: '//level
c --- Record 4:  Map projection
      write(iunit,'(a8)') pmap
c --- Record 5:  Map projection parameters
      if(LUTM) then
         write(iunit,'(i4,a4)') iutmzn,utmhem
      elseif(LLCC) then
         write(iunit,'(4a16)') clat0,clon0,clat1,clat2
      elseif(LPS) then
         write(iunit,'(3a16)') clat0,clon0,clat1
      elseif(LEM.or.LLAZA.or.LTTM) then
         write(iunit,'(2a16)') clat0,clon0
      endif
c --- Record 5a:  Map false Easting/Northing
      if(LLCC.or.LLAZA.or.LTTM) then
         write(iunit,*) feast,fnorth
      endif
c --- Record 6:  Map DATUM
      write(iunit,'(a8,a12)') datum,daten
c --- Record 7:  XYUNIT
      write(iunit,'(a4)') 'KM  '
c --- Record 8:  Time Zone
      write(iunit,'(a8)') atzone
c --- Record 9:  beginning & ending dates/times
      write(iunit,'(2(3i4,i6.4,2x))') kbyr,kbjul,kbhr,kbsec,
     &                                keyr,kejul,kehr,kesec
c --- Record 10:  No. sources, emitted species
      write(iunit,*) npt2,nse2

c --- Record #11:  Species list
      write(iunit,'(6a12)')(cslst2(n),n=1,nse2)

c --- Record #12:  Molecular weights of each emitted species
      write(iunit,'(6f12.3)')(xmwem2(n),n=1,nse2)

      return
      end
c----------------------------------------------------------------------
      subroutine WRTIEM2(iunit,lascii,nsrc2,cid,tiem2,bht2,bwd2)
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 061020                WRTIEM2
c                D. Strimaitis
c
c --- PURPOSE:  Write the time-invariant records to the PTEMARB point
c               source emissions file (arbitrarily-varying emissions)
c
c --- UPDATE:
c     5.0(990125) to 5.3(061020) --- DGS
c           -   Dataset 2.1 format (TIEM2 has 8 variables/source)
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(8,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(8,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,8)
            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,8)
            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,4f8.2)
102   format(6f9.3)

      return
      end
c----------------------------------------------------------------------
      subroutine WREM2(ibyr,ibjul,ibhr,ibsec,ieyr,iejul,iehr,iesec)
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 061020                  WREM2
c                D. Strimaitis
c
c --- PURPOSE:  Write the time-varying records of the PTEMARB point
c               source emissions file (arbitrarily-varying emissions)
c
c --- UPDATE:
c     5.0(990125) to 5.3(061020) --- DGS
c           -   Dataset 2.1 format: Begin/End times with seconds
c                                   Initial sigmas in data record
c
c --- INPUTS:
c     ibyr,ibjul,ibhr,ibsec - integer - Time at start of output period
c     ieyr,iejul,iehr,iesec - integer - Time at end of output period
c
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.cte'
      include 'pt2.cte'

      data zero/0.0/

      if(LASCII) then
      
c ---    Write records in ASCII form
c -----------------------------------
c ---    Write starting, ending date and time
         write(io16,101)ibyr,ibjul,ibhr,ibsec,ieyr,iejul,iehr,iesec
c ---    Loop over sources
         do i=1,npt2
            write(io16,102)cid(i),tstak2(i),exitw2(i),zero,zero,
     &                     (qstak2(k,i),k=1,2)
         enddo

      else
      
c ---    Write records in Binary form
c -----------------------------------
c ---    Write starting, ending date and time
         write(io16)ibyr,ibjul,ibhr,ibsec,ieyr,iejul,iehr,iesec
c ---    Loop over sources
         do i=1,npt2
            write(io16)cid(i),tstak2(i),exitw2(i),zero,zero,
     &                 (qstak2(k,i),k=1,2)
         enddo

      endif

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

      return
      end
c----------------------------------------------------------------------
      subroutine ctfluxa(iflow,tamb0,pamb0,rhamb0,month)
c----------------------------------------------------------------------
c --- CTEMISS    Version: 5.3      Level: 021031                CTFLUXA
c                D. Strimaitis,  Earth Tech
c
c --- PURPOSE:  ABATED Cooling Tower Emissions
c               Computes the water vapor emission rate (in addition to
c               ambient) and the exit temperature for a given heat
c               rejection rate, assuming that the effluent is at
c               saturation in the WET section, and additional heat is
c               added in the DRY section
c
c
c --- UPDATE:
c     5.1(020503) to 5.1(021031) --- DGS
c          (1)  Add eps2 for tolerance on dry section iterations
c          (2)  Reduce mandatory recirculation iterations
c          (3)  Add STOP for runaway iterations (iloop)
c
c     5.0(000203) to 5.1(020503) --- DGS
c          (1)  Add recirculation effect (estimated)

c     5.0(990223) to 5.0(000203) --- DGS
c          (1)  Exit velocity varies by EITHER month or temperature
c          (2)  Test for ZERO flux
c
c     5.0(990223) to 5.0(990325) --- DGS
c          (1)  Allow zero % ambient air drawn into DRY section to
c               override any non-zero % BTUs contributed by DRY section
c               so that tower becomes a pure wet tower if the air
c               intakes for the dry section are closed.
c
c --- INPUTS:
c     Common block /PT2/:
c           NPT2, TIEM2(7,mxpt2),
c           bturate(mxpt2),mbtu(mxpt2),btufac12(12,mxpt2),
c           wexit12(12,mxpt2),mwexit(mxpt2),
c           btudry12(12,mxpt2),airdry12(12,mxpt2),
c           mbtudry(mxpt2),tbin(12), recirc(36,mxpt2)
c     Parameters:
c           MXPT2, MXSPEC
c
c --- OUTPUT:
c     Common block /PT2/:
c           TSTAK2(mxpt2),EXITW2(mxpt2),QSTAK2(mxspec,mxpt2)
c
c --- CTFLUXA called by:  COMP
c --- CTFLUXA calls:      VSAT, HQFLUX, HQFLUXD
c----------------------------------------------------------------------
c
c --- Include parameter statements and commons
      include 'params.cte'
      include 'pt2.cte'

      logical labate,ldb

c --- Set constants
      data pi/3.1415927/
      data btu2cal/252.02/
      data hr2sec/3600./
      data eps/.00001/
      data eps2/.00005/
c --- CONST is MW(h2o) / Rgas = 18.016/.0831372
      data const/216.702/

c      ldb=.TRUE.
      ldb=.FALSE.

c --- Transfer ambient measures to working variables
      tamb=tamb0
      pamb=pamb0
      rhamb=rhamb0

c --- Compute background water vapor mixing ratio (g h2o/g air) which
c     is assumed constant for the surface layer; and report as
c     --- PBSAT   saturation vapor pressure (mb) at TAMB
c     --- PBACK0  background vapor pressure (mb) at TAMB
c     --- QBSAT   saturation vapor mass conc (g/m^3)
c     --- QBACK0  background vapor mass conc (g/m^3)
      call VSAT(tamb,pbsat,qbsat)
      frh=0.01*rhamb
      qback0=frh*qbsat
      pback0=frh*pbsat

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 --- Save variables based on true ambient conditions
      frha=frh
      qback0a=qback0
      pback0a=pback0
      faca=fac
      itcata=itcat

c --- Loop over cells (sources)
c -----------------------------
      do is=1,npt2

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

c ---    Set recirculation information for this cell
         rf=recirc(iflow,is)
         if(rf.GT.0.0) then
            irecirc=1
         else
            irecirc=0
         endif

c ---    Transfer stored ambient measures to working variables
         tamb=tamb0
         pamb=pamb0
         rhamb=rhamb0
         frh=frha
         qback0=qback0a
         pback0=pback0a
         fac=faca
         itcat=itcata
         itcatm1=itcat-1
         tlast=tamb

c ---    Initialize ambient water fluxes for the cell
         qamb0=0.0
         qambw=0.0
         qambd=0.0

c ---    Re-entry point for recirculation loop
c --------------------------------------------
4        if(irecirc.GT.1) then
c ---       Update inlet air properties as modified by recirculation
c ---       Approximate entrainment effect with simple averages
            tamb=(1.-rf)*tamb0+rf*texit
            qback0=qback0a+rf*qh2o/vf0
            call VSAT(tamb,pbsat,qbsat)
            frh=qback0/qbsat
            frh=AMIN1(frh,1.0)
            rhamb=100.*frh
            qback0=frh*qbsat
            pback0=frh*pbsat
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
         endif

c ---    Get exit velocity
         if(mwexit(is).EQ.1) then
            wexit=wexit12(month,is)
         elseif(mwexit(is).EQ.2) then
            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
         else
            stop 'CTFLUXA:  mwexit is not 1 or 2'
         endif

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

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

c ---    Get BTU/HR heat rejection rate
         if(mbtu(is).EQ.1) then
            btufac=btufac12(month,is)
         elseif(mbtu(is).EQ.2) then
            if(itcat.EQ.1) then               
               btufac=btufac12(1,is)
            elseif(itcat.EQ.13) then
               btufac=btufac12(12,is)
            else
               btufac=btufac12(itcatm1,is)+
     &                fac*(btufac12(itcat,is)-btufac12(itcatm1,is))
            endif
         else
            stop 'CTFLUXA:  mbtu is not 1 or 2'
         endif
         btu0=bturate(is)*btufac

c ---    Get FRACTION of BTU/HR heat rejection from DRY section
         if(mbtudry(is).EQ.1) then
            btudry=btudry12(month,is)
         elseif(mbtudry(is).EQ.2) then
            if(itcat.EQ.1) then               
               btudry=btudry12(1,is)
            elseif(itcat.EQ.13) then
               btudry=btudry12(12,is)
            else
               btudry=btudry12(itcatm1,is)+
     &                fac*(btudry12(itcat,is)-btudry12(itcatm1,is))
            endif
         else
            stop 'CTFLUXA:  mbtudry is not 1 or 2'
         endif
         fdbtu=0.01*AMAX1(0.0,btudry)

c ---    Get FRACTION of ambient air entering DRY section
         fdvol=0.01*AMAX1(0.0,airdry12(month,is))

c ---    Check for ZERO ambient airflow into DRY section and reset the
c ---    FRACTION of BTU/HR heat rejection from DRY section to ZERO
         if(fdvol.EQ.0.0) fdbtu=0.0

c ---    Check for zero flux and skip to end
         if(wexit.LE.0.0001 .OR. btu0.LE.0.0001) then
            tstak2(is)=tamb
            exitw2(is)=0.0
            qstak2(1,is)=0.0
            qstak2(2,is)=0.0
            goto 50
         endif

c ---    Convert heat rejection rate from BTU/Hr to cal/s
         heat0=btu0*btu2cal/hr2sec

c ---    Set heat rates for wet and dry sections
         dheat0=fdbtu*heat0
         wheat0=heat0-dheat0

c ---    Is there a dry section to this tower?
         labate=.TRUE.
         if(fdbtu.LE.0.0 .AND. fdvol.LE.0.0) labate=.FALSE.

c ---    Iteration loop for volume flux at top of wet section
c -----------------------------------------------------------
         iloop=0
c ---    First guess: use full flux times input fraction
         vfwet1=vf0*(1.-fdvol)
         vfwet2=vfwet1

5        iloop=iloop+1
         if(iloop.EQ.2) then
c ---       Second guess: use 90% full flux times input fraction
            vfwet2=0.9*vf0*(1.-fdvol)
         endif

c ---    Iterate to find temperature of effluent from WET section of CT
c ---    First guess: assume temperature rise of 10 K
         t1=tamb+10.
         call HQFLUX(t1,tamb,pamb,pback0,vfwet2,qh2o,qambw,psat,h1)
c ---    Second guess: assume temperature rise of 20 K
         t2=tamb+20.
         call HQFLUX(t2,tamb,pamb,pback0,vfwet2,qh2o,qambw,psat,h2)
c ---    Iterate: assume linear relationship between H and T
10       hdiff=wheat0-h2
         if(ABS(hdiff/wheat0) .LT. eps) goto 20
         t=t2+(t1-t2)*(hdiff)/(h1-h2)
         t1=t2
         h1=h2
         t2=t
         call HQFLUX(t2,tamb,pamb,pback0,vfwet2,qh2o,qambw,psat,h2)
         goto 10
      
20       qrate=qh2o
c ---    Set ambient vapor flux into DRY section as a fraction of
c ---    the total
         qambd=qambw*fdvol/(1.-fdvol)
         twet=t2
c ---    Compute temperature at top of CT
         texit=vf0*const/((qambw+qh2o)/psat+qambd/pback0)

c ---    Compute and add heat flux from in DRY section to obtain total
         call HQFLUXD(texit,twet,tamb,pamb,pback0,psat,
     &                qh2o,qambw,qambd,hd2)
         ht2=h2+hd2
         if(iloop.EQ.1) then
            ht1=ht2
            if(LABATE) then
               goto 5
            else
c ---          No iteration on wet section flux needed
               goto 40
            endif
         endif
c ---    Iterate: assume linear relationship between vol flux in wet
c ---             section and total computed flux
30       hdiff=heat0-ht2

c ---    Special DEBUG section
         if(ldb) write(io6,*)'heat0,ht2 = ',heat0,ht2
         if(iloop.GT.100) stop 'iloop GT 100'
c ---    Special DEBUG section

         if(ABS(hdiff/heat0) .LT. eps2) goto 40
         vfwet=vfwet2+(vfwet1-vfwet2)*(hdiff)/(ht1-ht2)
         vfwet1=vfwet2
         ht1=ht2
         vfwet2=vfwet
         goto 5

40       continue

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

c ---    Store ambient water fluxes from the first pass
         if(irecirc.LE.1) qamb0=qambw+qambd

c ---    DEBUG SECTION
         if(LDB) then
            write(io6,*)'CTFLUXA:    labate = ',labate
            write(io6,*)'       CELL Number = ',is
            write(io6,*)'  btu0,fdbtu,fdvol = ',btu0,fdbtu,fdvol
            write(io6,*)'      tamb0,rhamb0 = ',tamb0,rhamb0
            write(io6,*)'        rf,irecirc = ',rf,irecirc
            write(io6,*)'        tamb,rhamb = ',tamb,rhamb
            write(io6,*)'        qamb0,qh2o = ',qamb0,qh2o
            write(io6,*)'  texit,qamb0,qh2o = ',texit,qamb0,qh2o
            write(io6,*)'       qambw,qambd = ',qambw,qambd
         endif

c ---    Recirculation
c ---    Do the next pass if needed
         deltamb=ABS(tamb-tlast)
         tlast=tamb
         if(irecirc.LE.2 .AND. irecirc.NE.0) then
            irecirc=irecirc+1
            goto 4
         elseif(deltamb.GT.0.01 .AND. irecirc.NE.0) then
            irecirc=irecirc+1
            goto 4
         else
            delq=AMAX1(0.0,qambw+qambd-qamb0)
         endif

c ---    Pass results to output arrays
         tstak2(is)=texit
         exitw2(is)=wexit
c ---    (Add elevated ambient water flux from recirculation)
         qstak2(1,is)=qrate+delq
         qstak2(2,is)=dtemp*vf0

50       continue

      enddo
c --- End Loop over sources

      return
      end
c----------------------------------------------------------------------
      subroutine hqflux(t,tamb,pamb,pbkga,vf0,qh2o,qamb,psat,h)
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 990125                 HQFLUX
c                D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Compute flux of heat (latent and sensible) and the
c               corresponding mass flux of water vapor evaporated into
c               the plume to saturate it
c
c --- INPUTS:
c              t - real       - Temperature of plume (K)
c           tamb - real       - Temperature of ambient air (K)
c           pamb - real       - Pressure of ambient air (mb)
c          pbkga - real       - Ambient H2O vapor pressure at TAMB (mb)
c            vf0 - real       - Exit volume flux (m^3/s)
c
c --- OUTPUT:
c           qh20 - real       - Water vapor flux (g/s) added to plume
c           qamb - real       - Ambient water vapor flux (g/s)
c           psat - real       - Saturation vapor pressure at T (mb)
c              h - real       - Heat flux (cal/s) added to plume
c
c --- HQFLUX called by:   CTFLUX
c --- HQFLUX calls:       VSAT
c----------------------------------------------------------------------
      data h2omw/18.016/    !  g/mole
      data rgas/8.31372e-2/ !  mb m^3/mole/K
c         (cpdry/1.004/     !  j/g/K)
c         (cpwet/1.8765/    !  j/g/K)
      data cpdry/.23986/    !  cal/g/K
      data cpwet/.44837/    !  cal/g/K
      data heatl0/597.49/   !  cal/g (at 273.15K)
      data dheatldt/-.566/  !  cal/g/K

      const=h2omw/(rgas*t)
      camb=const*pbkga

      call VSAT(t,psat,csat)
      cplm=csat
      qh2o=(cplm-camb)*vf0
      qamb=camb*vf0

      rsat=.622*psat/(pamb-psat)
      rplm=rsat

      qheat=qh2o*(heatl0+dheatldt*(t-273.15))
      theat=(t-tamb)*(cpwet*camb+cpdry*cplm/rplm)*vf0
      h=qheat+theat

c !!! DEBUG
c      write(*,*)'HQFLUX: tamb,pamb,pbkga,vf0 = '
c      write(*,*) tamb,pamb,pbkga,vf0
c      write(*,*)'     t,psat,csat = ',t,psat,csat
c      write(*,*)'       rsat,rplm = ',rsat,rplm
c      write(*,*)'       cplm,camb = ',cplm,camb
c      write(*,*)'qh2o,qheat,theat = ',qh2o,qheat,theat
c !!! DEBUG

      return
      end
c----------------------------------------------------------------------
      subroutine hqfluxd(texit,twet,tamb,pamb,pback0,psat,
     &                   qh2o,qambw,qambd,h)
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 990125                HQFLUXD
c                D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Compute flux of heat (sensible) added by the dry section
c               of the Cooling Tower
c
c --- INPUTS:
c          texit - real       - Temperature of plume at exit (K)
c           twet - real       - Temperature of flow from WET section (K)
c           tamb - real       - Temperature of ambient air (K)
c           pamb - real       - Pressure of ambient air (mb)
c         pback0 - real       - Ambient H2O vapor pressure at TAMB (mb)
c           psat - real       - Saturation H2O vapor press at TWET (mb)
c           qh2o - real       - Water flux added in WET section (g/s)
c          qambw - real       - Ambient water flux into WET section (g/s)
c          qambd - real       - Ambient water flux into DRY section (g/s)
c
c --- OUTPUT:
c              h - real       - Heat flux (cal/s) added to plume in the
c                               DRY section
c
c --- HQFLUXD called by:   CTFLUXA
c --- HQFLUXD calls:       none
c----------------------------------------------------------------------
c --- dry air specific heat     (cpdry/1.004/     !  j/g/K)
c --- water vapor specific heat (cpwet/1.8765/    !  j/g/K)
      data cpdry/.23986/    !  cal/g/K
      data cpwet/.44837/    !  cal/g/K

c --- Compute mixing ratios (g h2o / g air) 
      rsat=.622*psat/(pamb-psat)
      ramb=.622*pback0/(pamb-pback0)

c --- Heat added to warm the flow from the WET section to exit temp
      theat1=(texit-twet)*(qh2o+qambw)*(cpwet+cpdry/rsat)

c --- Heat added to warm the ambient flow into the DRY section to the
c --- exit temp
      theat2=(texit-tamb)*qambd*(cpwet+cpdry/ramb)

      h=theat1+theat2

c !!! DEBUG
c      write(*,*)
c      write(*,*)'HQFLUXD: texit,twet,tamb = ',texit,twet,tamb
c      write(*,*)'        pamb,pback0,psat = ',pamb,pback0,psat
c      write(*,*)'        qh2o,qambw,qambd = ',qh2o,qambw,qambd
c      write(*,*)'               rsat,ramb = ',rsat,ramb
c      write(*,*)'           theat1,theat2 = ',theat1,theat2
c !!! DEBUG

      return
      end
c----------------------------------------------------------------------
      subroutine vsat(tdegk,psat,csat)
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 980918F                  VSAT
c                J. Scire, Earth Tech
c
c --- PURPOSE:  Compute saturation vapor pressure and corresponding
c               saturation mass concentration for given temperature
c
c --- INPUTS:
c          tdegk - real       - Temperature (K)
c
c --- OUTPUT:
c           psat - real       - Saturation water vapor pressure (mb)
c           csat - real       - Saturation water vapor mass
c                               concentration (g h2o/m**3 air)
c
c --- VSAT called by:   
c --- VSAT calls:       none
c----------------------------------------------------------------------
c
      data gascon/8.31372e-2/

c --- Compute saturation vapor pressure (Goff-Gratch eqns.)
      trat=273.16/tdegk
      if(tdegk.lt.273.16)then
c ---    temperature below 0. deg. C
         psat=exp(2.303*(-9.09685*(trat-1.0)+
     1   0.87682*(1.0-tdegk/273.16))-3.56654*(5.6100577-alog(tdegk))+
     2   1.81)
      else
c ---    temperature above 0. deg. C
         psat=exp(2.303*(10.79574*(1.0-trat)+
     1   1.50475e-4*(1.0-10.**(-8.2969*(tdegk/273.16-1.0)))+
     2   4.2873e-4*(10.**(4.76955*(1.0-trat))-1.0))-
     3   5.028*(alog(tdegk)-5.6100577)+1.81)
      endif

c --- Compute saturation water vapor mass concentration (g h2o/m**3 air)
c --- (gascon is the gas constant -- 8.31372e-2 m**3*mb/(mole*deg k)
      csat=18.016*psat/(gascon*tdegk)

      return
      end
c----------------------------------------------------------------------
      subroutine DIAGINI
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 990225                DIAGINI
c                D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Initialize diagnostic variables
c
c --- INPUTS:
c            none
c
c --- OUTPUT:
c     Common block /DIAG/:
c           NHOURS,QRATESUM,IRHDIST(100),IVIS(12),IVISHR(0:23,12),
c           I100(0:23),IFOGAMB
c                            
c --- OUTPUT:
c
c --- DIAGINI called by: COMP
c --- DIAGINI calls:     none
c----------------------------------------------------------------------
c
c --- Include parameter statements and commons
      include 'diag.cte'

      nhours=0
      qratesum=0.0
      ifogamb=0
      do i=1,100
         irhdist(i)=0
      enddo
      do i=0,23
         i100(i)=0
      enddo
      do i=1,12
         ivis(i)=0
         do j=0,23
            ivishr(j,i)=0
         enddo
      enddo

      return
      end
c----------------------------------------------------------------------
      subroutine DIAGUP(tamb,pamb,rhamb,ihr,imo,irh)
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 990225                 DIAGUP
c                D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Update diagnostic information
c
c --- INPUTS:
c            TAMB - real    - Ambient temperature (K)
c            PAMB - real    - Ambient pressure (mb)
c           RHAMB - real    - Ambient RH (%)
c             IHR - integer - Hour (0-23)
c             IMO - integer - Month (Jan = 1)
c
c     Common block /PT2/:
c           TSTAK2(mxpt2),QSTAK2(mxspec,mxpt2),EXITW2(mxpt2),
c           TIEM2(7,mxpt2) [(4,-) = Stack diameter (m)]
c                            
c     Parameters:
c           MXPT2, MXSPEC
c
c --- OUTPUT:
c             IRH - integer - Plume RH at exit (%)
c
c     Common block /DIAG/:
c           NHOURS,QRATESUM,IRHDIST(100),IVIS(12),IVISHR(0:23,12),
c           I100(0:23),IFOGAMB
c
c --- DIAGUP called by: COMP
c --- DIAGUP calls:     VSAT
c----------------------------------------------------------------------
c
c --- Include parameter statements and commons
      include 'params.cte'
      include 'diag.cte'
      include 'pt2.cte'

      logical lsat,ldb

      data pi/3.1415927/

      ldb=.FALSE.

c --- Update counter
      nhours=nhours+1

c --- Add current water vapor emission rate from FIRST source to sum
      qratesum=qratesum+qstak2(1,1)

c --- Compute background water vapor mixing ratio (g h2o/g air) which
c     is assumed constant for the surface layer; and report as
c     --- PBSAT   saturation vapor pressure (mb) at TAMB
c     --- PBACK0  background vapor pressure (mb) at TAMB
c     --- QBSAT   saturation vapor mass conc (g/m^3)
c     --- QBACK0  background vapor mass conc (g/m^3)
      call VSAT(tamb,pbsat,qbsat)
      frh=0.01*rhamb
      qback0=frh*qbsat
      pback0=frh*pbsat

c --- Compute relative humidity (%) of the plume as integer
      call VSAT(tstak2(1),p0sat,q0sat)
      qbypiu=qstak2(1,1)/(pi*exitw2(1))
      c0=qbypiu/(0.25*(tiem2(4,1)**2))
      rhplm=100.*(c0/q0sat+pback0/p0sat)
      irh=NINT(rhplm)

c --- Update RH distribution array (RH = 30 to 130 %)
      index=irh-29
      index=MAX(index,1)
      index=MIN(index,100)
      irhdist(index)=irhdist(index)+1

c --- Update monthly distribution of visible plume-hours
      lsat=.FALSE.
      camb=qback0
      cplm=q0sat*rhplm*.01
      tplm=tstak2(1)
      irhamb=NINT(rhamb)
      if(irhamb.GE.100) then
c ---    Do not include ambient fog hours
         lsat=.TRUE.
         i100(ihr)=i100(ihr)+1
         ifogamb=ifogamb+1
      elseif(rhplm.GE.100.5) then
         lsat=.TRUE.
         ivis(imo)=ivis(imo)+1
         ivishr(ihr,imo)=ivishr(ihr,imo)+1
      else
         delc=cplm-camb
         delt=tplm-tamb
c ---    Check for saturation in steps of 0.005 along mixing line
         step=0.0
10       step=step+0.005
         t=tamb+step*delt
         c=camb+step*delc
         call VSAT(t,psat,csat)
         if(csat.LE.c) then
            lsat=.TRUE.
         elseif(step.LT.0.98) then
            goto 10
         endif
         if(lsat) then
            ivis(imo)=ivis(imo)+1        
            ivishr(ihr,imo)=ivishr(ihr,imo)+1
         endif
      endif

c --- DEBUG OUTPUT
      if(LDB) then
         write(io6,'(i6,2(f10.3,f10.4),f8.2,5x,L1)')
     &         ihr,tamb,camb,tplm,cplm,rhplm,lsat
      endif

      return
      end
c----------------------------------------------------------------------
      subroutine DIAGOUT
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 990225                DIAGOUT
c                D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Write the diagnostic information to the list file
c
c --- INPUTS:
c     Common block /DIAG/:
c           NHOURS,QRATESUM,IRHDIST(100),IVIS(12),IVISHR(0:23,12),
c           I100(0:23),IFOGAMB
c                            
c --- OUTPUT:
c
c --- DIAGOUT called by: COMP
c --- DIAGOUT calls:     none
c----------------------------------------------------------------------
c
c --- Include parameter statements and commons
      include 'params.cte'
      include 'diag.cte'

      integer ihrtot(24)

c --- Set total water added by source
      qtotkg=qratesum*3.6

c --- Set average emission rate
      qaverage=qratesum/FLOAT(nhours)

c --- Output to file
      write(io6,*)
      write(io6,*)
      write(io6,*)'    DIAGNOSTICS for Source 1   -------------------'
      write(io6,*)
      write(io6,*)'    Number of hours             ',nhours
      write(io6,*)'    Total emission (kg)         ',qtotkg
      write(io6,*)'    Average emission rate (g/s) ',qaverage
      write(io6,*)
      write(io6,*)'    Plume RH(%) Bin      Number(hours)'
      write(io6,*)
      do i=100,1,-1
         irh=i+29
         write(io6,'(2(10x,i8))') irh,irhdist(i)       
      enddo
      write(io6,*)
      write(io6,*)
      write(io6,*)
     &'    Visible Plume Results for Hours With Ambient RH < 100%'
      write(io6,*)'    (number of hours excluded: ',ifogamb,' )'
      write(io6,*)
      write(io6,*)
      write(io6,*)'             Month        Visible Plumes'
      write(io6,*)
      ivistot=0
      do i=1,12
         write(io6,'(2(10x,i8))') i,ivis(i)       
         ivistot=ivistot+ivis(i)
      enddo
      write(io6,*)
      write(io6,'(28x,i8)') ivistot       
      write(io6,*)
      write(io6,*)'Visible Plumes (by hour)'
      write(io6,*)
      write(io6,'(2a)')' Mon h01 h02 h03 h04 h05 h06 h07 h08 h09 h10',
     &     ' h11 h12 h13 h14 h15 h16 h17 h18 h19 h20 h21 h22 h23 h24'
      write(io6,*)
      do j=1,24
         ihrtot(j)=0
      enddo
      do i=1,12
         write(io6,'(25i4)') i,(ivishr(j,i),j=1,23),ivishr(0,i)
         do j=0,23
            jhr=j
            if(j.EQ.0) jhr=24
            ihrtot(jhr)=ihrtot(jhr)+ivishr(j,i)
         enddo
      enddo
      write(io6,*)
      write(io6,'(4x,24i4)') ihrtot       
      write(io6,*)
      write(io6,*)'Ambient Fog (by hour)'
      write(io6,*)
         write(io6,'(4x,24i4)') (i100(j),j=1,23),i100(0)

      return
      end
c----------------------------------------------------------------------
      subroutine fin
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.3      Level: 061020                    FIN
c ---            J. Scire
c
c --- PURPOSE:  Run termination routine -- compute runtime,
c               write last day processed
c
c --- UPDATE:
c     5.2(031013) to 5.3(061020) --- DGS
c          (1)  Report start and end of last period with seconds
c
c     5.0(990125) to 5.2(031013) --- DGS
c          (1)  RDATE is c*10 with YYYY for year
c
c --- INPUTS:
c       Common block /QA/
c          rdate, rtime, rcpu, kyr, kjul, khr, ksec, isec
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.cte'
      include 'qa.cte'

      character*8 rtime2
      character*10 rdate2

      write(iomesg,*)'TERMINATION PHASE'

c --- Write last day/hour processed
      write(io6,5)kyr,kmonth,kday,kjul,khr,ksec
      write(io6,6)kyr,kmonth,kday,kjul,khr,isec
5     format(//2x,'LAST DAY/HOUR PROCESSED STARTS AT:'/5x,'Year: ',i4,
     1   2x,'Month: ',i2,3x,'Day: ',i2,3x,'Julian day: ',i3,3x,'Hour: ',
     2   i2,2x,'Second: ',i4)
6     format(//2x,'LAST DAY/HOUR PROCESSED ENDS AT  :'/5x,'Year: ',i4,
     1   2x,'Month: ',i2,3x,'Day: ',i2,3x,'Julian day: ',i3,3x,'Hour: ',
     2   i2,2x,'Second: ',i4)


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),10)iyr1
         call julday(io6,iyr1,imo1,iday1,ijul1)

         read(rdate2(1:2),10)imo2
         read(rdate2(4:5),10)iday2
         read(rdate2(7:10),10)iyr2
         call julday(io6,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 yr4(io,iyr,ierr)
c----------------------------------------------------------------------
c
c --- CALUTILS   Version: 2.1      Level: 991104                    YR4
c ---            D. Strimaitis,   Earth Tech
c
c --- PURPOSE: Checks/converts 2-digit year to 4-digit year
c
c --- INPUTS:
c                IO - integer    - Unit number for list file output
c               IYR - integer    - Year (YYYY or YY)
c
c     Common block /Y2K/:
c             IYYLO - integer    - Smallest 2-digit year for which
c                                  'old' century marker is used
c             ICCLO - integer    - 2-digit ('old') century
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----------------------------------------------------------------------
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(io,*)'ERROR in YR4 --- Year not recognized: ',iyr
         write(*,*)'ERROR in YR4 --- Year not recognized: ',iyr
      endif

      return
      end
c----------------------------------------------------------------------
      subroutine yr4c(iyr)
c----------------------------------------------------------------------
c
c --- CALUTILS   Version: 2.1      Level: 020828                   YR4C
c ---            D. Strimaitis,  Earth Tech
c
c --- PURPOSE: Checks/converts 2-digit year to 4-digit year (CURRENT)
c
c --- UPDATE
c --- V1.0-V1.1     020828  (DGS): Add check for YYYY on input
c
c --- INPUTS:
c               IYR - integer    - Year (YYYY or YY)
c
c --- OUTPUT:
c               IYR - integer    - Year (YYYY)
c
c --- YR4C called by:  host subroutines
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 --- CALUTILS   Version: 2.1      Level: 991104                  QAYR4
c ---            D. Strimaitis,   Earth Tech
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 - integer    - Smallest 2-digit year for which
c                                  'old' century marker is used
c             ICCLO - integer    - 2-digit ('old') century
c
c --- QAYR4 called by:  host subroutines
c --- QAYR4 calls:      none
c----------------------------------------------------------------------
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/50/

      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
         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 julday(io,iyr,imo,iday,ijuldy)
c----------------------------------------------------------------------
c
c --- CALUTILS   Version: 2.1      Level: 000602                 JULDAY
c ---            J. Scire, SRC
c
c --- PURPOSE:  Compute the Julian day number from the Gregorian
c               date (month, day)
c
c --- UPDATE
c ---               000602  (DGS): YYYY format for year
c
c --- INPUTS:
c            IO - integer      - Unit number for list file output
c           IYR - integer      - Year
c           IMO - integer      - Month
c          IDAY - integer      - Day
c
c --- OUTPUT:
c          IJUL - integer      - Julian day
c
c --- JULDAY called by:  host subroutines
c --- JULDAY calls:      none
c----------------------------------------------------------------------
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(io,*)
         write(io,*)'ERROR in SUBR. JULDAY'
         write(io,*)'Invalid date - IYR = ',iyr,' IMO = ',
     1    imo,' IDAY = ',iday
         write(*,*)
         stop 'Halted in JULDAY -- see list file.'
      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(io,iyr,ijul,imo,iday)
c----------------------------------------------------------------------
c
c --- CALUTILS   Version: 2.1      Level: 000602                  GRDAY
c                J. Scire, SRC
c
c --- PURPOSE:  Compute the Gregorian date (month, day) from the
c               Julian day
c
c --- UPDATE
c ---               000602  (DGS): YYYY format for year
c
c --- INPUTS:
c            IO - integer      - Unit number for list file output
c           IYR - integer      - Year
c          IJUL - integer      - Julian day
c
c --- OUTPUT:
c           IMO - integer      - Month
c          IDAY - integer      - Day
c
c --- GRDAY called by:  host subroutines
c --- GRDAY calls:      none
c----------------------------------------------------------------------
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(io,12)iyr,ijul
12    format(//2x,'ERROR in SUBR. GRDAY -- invalid Julian day '//2x,
     1 'iyr = ',i5,3x,'ijul = ',i5)
      write(*,*)
      stop 'Halted in GRDAY -- see list file.'
      end
c------------------------------------------------------------------------------
      subroutine deltt(j1yr,j1jul,j1hr,j2yr,j2jul,j2hr,jleng)
c------------------------------------------------------------------------------
c
c --- CALUTILS   Version: 2.1      Level: 941215                  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:  host subroutines
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 incr(io,iyr,ijul,ihr,nhrinc)
c----------------------------------------------------------------------
c
c --- CALUTILS   Version: 2.1      Level: 000602                   INCR
c                J. Scire, SRC
c
c --- PURPOSE:  Increment the time and date by "NHRINC" hours
c
c --- UPDATE
c ---               000602  (DGS): add message to "stop"
c ---               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       IO     - integer - Unit number for list file output
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 --- OUTPUT:
c       IYR    - integer - Updated year
c       IJUL   - integer - Updated Julian day
c       IHR    - integer - Updated hour (00-23)
c
c --- INCR called by: host subroutines
c --- INCR calls:     none
c----------------------------------------------------------------------
c
c --- Check nhrinc
      if(nhrinc.lt.-24) then
         write(io,*)'ERROR IN SUBR. INCR -- Invalid value of NHRINC ',
     1   '-- NHRINC = ',nhrinc
         write(*,*)
         stop 'Halted in INCR -- see list file.'
      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 midnite(io,ctrans,iyr,imo,iday,ijul,
     &                             kyr,kmo,kday,kjul)
c----------------------------------------------------------------------
c
c --- CALUTILS   Version: 2.54     Level: 041029                MIDNITE
c ---            D. Strimaitis, Earth Tech
c
c --- PURPOSE:  Converts date/time at midnight between day N, 0000
c               and day N-1, 2400.  Direction is determined by the
c               CTRANS instruction.
c
c --- INPUTS:
c            IO - integer      - Unit number for list file output
c        CTRANS - character    - Instruction 'TO 24h' or 'TO 00h'
c           IYR - integer      - Year
c           IMO - integer      - Month
c          IDAY - integer      - Day
c          IJUL - integer      - Julian day
c
c --- OUTPUT:
c           KYR - integer      - Year
c           KMO - integer      - Month
c          KDAY - integer      - Day
c          KJUL - integer      - Julian day
c
c --- MIDNITE called by:  host subroutines
c --- MIDNITE calls:      JULDAY, INCR, GRDAY
c----------------------------------------------------------------------
      character*6 ctrans

      ierr =0

c --- Get Julian day from month/day if needed
      if(ijul.LE.0) call JULDAY(io,iyr,imo,iday,ijul)

      kyr=iyr
      kmo=imo
      kday=iday
      kjul=ijul

      if(ctrans.EQ.'TO 24h') then
c ---    Convert from 0000 on ijul to 2400 on kjul
         ihr=0
         nhr=-1
         call INCR(io,kyr,kjul,ihr,nhr)
         call GRDAY(io,kyr,kjul,kmo,kday)
      elseif(ctrans.EQ.'TO 00h') then
c ---    Convert from 2400 on ijul to 0000 on kjul
         ihr=23
         nhr=1
         call INCR(io,kyr,kjul,ihr,nhr)
         call GRDAY(io,kyr,kjul,kmo,kday)
      else
         ierr=1
      endif

      if(ierr.eq.1)then
         write(io,*)
         write(io,*)'ERROR in SUBR. MIDNITE'
         write(io,*)'Invalid instruction: ',ctrans
         write(io,*)'           Expected: TO 24h'
         write(io,*)'              OR   : TO 00h'
         write(*,*)
         stop 'Halted in MIDNITE -- see list file.'
      endif

      return
      end
c----------------------------------------------------------------------
      subroutine datetm(rdate,rtime,rcpu)
c----------------------------------------------------------------------
c
c --- CALUTILS   Version: 2.1      Level: 010731                 DATETM
c ---            J. Scire, SRC
c
c --- PURPOSE:  Get system date and time from system clock
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, FIN
c --- DATETM calls:      DATE, TIME (Lahey F77 compiler utilities)
c                        DATE_AND_TIME (COMPAQ, LF95 compilers)
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!]

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

c --- COMPAQ DF90/95 and Lahey LF95 Compilers (begin)
c ---------------------------------------------------
c --- System date in CCYYMMDD
c --- System clock in HHMMSS.sss, where sss = thousandths of seconds
      call DATE_AND_TIME(sdate,stime)
c --- Pass to output formats (MM-DD-YYYY) and (HH:MM:SS)
      rdate='  -  -    '
      rdate(1:2)=sdate(5:6)
      rdate(4:5)=sdate(7:8)
      rdate(7:10)=sdate(1:4)
      rtime='  :  :  '
      rtime(1:2)=stime(1:2)
      rtime(4:5)=stime(3:4)
      rtime(7:8)=stime(5:6)
c --- COMPAQ DF90/95 and Lahey LF95 Compilers (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 --- CALPUFF    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 --- CALPUFF    Version: 5.3      Level: 991217                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 ---           Lahey version
c ---------------------------
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 'ctemiss.inp'
c
c --- COMLINE called by:  SETUP
c --- COMLINE calls:      GETCL - Lahey system routine
c                         NARGS, GETARG - Microsoft routines
c                         IARGC, GETARG - Sun routines
c                         IARGC, GETARG - HP routines
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/'ctemiss.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 --- If no command line arguments, use default
      if(ctext(1:1).eq.' ')ctext=cdeflt
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
      return
      end
