c----------------------------------------------------------------------
c --- POSTRM --  Postprocessor for receptor-mode FOG model output
c----------------------------------------------------------------------
c
c --- POSTRM     Version: 1.2      Level: 060214                   MAIN
c ---            J. Scire,  SRC
c
c --- Written by:  Joe Scire
c                  Sigma Research Corporation
c                  196 Baker Avenue
c                  Concord, MA  01742
c
c                  TEL. (978) 371-4200
c                  FAX  (978) 371-2468
c
c --- Copyright (c) 1998-2006 by TRC Environmental Corporation
c
c----------------------------------------------------------------------


c *** NOTE:  QA output for source contributions was turned off 9/4/2002

c
c --- Changes from 020228 to 060214:
c
c --- 060214   (DGS): Remove ijcode,lscale,iscode,lsum,nsum
c                     from NAMELIST
c              (DGS): Add fogmet, rhout,lstout, fogdat
c                     to NAMELIST, and change fogdat default from
c                     FOG.BIN to FOG.DAT
c
c --- Changes from 020202 to 020228:
c
c --- 020228   (DGS): Add change in RH and change in Temperature to
c                     output file(s) produced in LRHMODE
c              (DGS): Add Humidex output in LRHMODE
c
c --- Changes from 010924 to 020202:
c
c --- 020202   (DGS): Add output file that lists the weather codes for
c                     hours in which at least 1 receptor shows fog/ice
c                     (unit 9, 'foghrs.met').  Both LPRT and LMET
c                     namelist options must be true for this output.
c
c --- Changes from 001010 to 010924:
c
c --- 010905   (DGS): Revise date processing to allow CALPUFF output
c                     file to be processed completely (multiple years)
c              (DGS): Change version and level to character variables
c
c --- Changes from 000811 to 001010:
c
c --- 001010   (MNB): Include/update PROC144RM from 991203F:
c                     Add present weather procesing (CD144 metdat)
c                     (adapted from BdF changes in POSTRM)
c                     Add lmet & metdat to input file
c                     if lmet = TRUE: read file metdat (CD144)
c                     and print weather and cloud codes
c
c --- 001003   (DGS): XPND assumes that a single zero will not be coded
c                     as '-1.0', but CALPUFF output appears to pass
c                     a coded '-1.0'.  Remove this assumption.
c	       (DGS): Change receptor number format (label 1115) to
c                     write year as i4
c	       (DGS): Add error mesage if receptors exceed MXREC
c
c --- Changes from 870225 to 000811:
c
c --- 000811   (DGS): Use LCOMPRS logical in header record 2 from
c                     CALPUFF v5.4, Level 000602_2 in place of
c                     detecting compression from nfrec>1000
c              (DGS): Change year format to write i4 (YYYY)
c
c --- 990226F  (DGS): Add new output/processing features
c                      1 - Add LCALMPRT to control output of calm-hour
c                                records
c                      2 - Add LRHMODE to activate output of RH
c                                histogram for receptors
c                      3 - Add LDB for debug output
c
c --- 981116F  (DGS): Accept CALPUFF variation of FOG.BIN
c                      1 - IDAYS array zeroed (process all
c                                hours in file)
c                      2 - QSAT, QBACK provided for each receptor
c 
c----------------------------------------------------------------------
c --- File Unit Numbers
c      5   file='postrm2.inp'
c      6   file='postrm2.lst'
c      7   file='postrm2.rh'
c      9   file='foghrs.met'
c      10  file='fog.bin'
c      17  file= metdat
c----------------------------------------------------------------------
c
      parameter(mxrec=5000)
c
      real chi(mxrec),chiin(mxrec)

c --- Add arrays for qback and qsat at each receptor
      real qbackr(mxrec),qsatr(mxrec)

c --- Add array for counting RH frequency at each receptor (& ambient)
      integer numrh(20,0:mxrec)
      real rhavg(0:mxrec)

c --- Add arrays for counting delta-RH frequency and delta-T frequency
      integer numdrh(26,mxrec)
      integer numdt(21,mxrec)
      real percent(26,0:mxrec)
      real bindrh(27),bindt(22)

c --- Add arrays for counting HUMIDEX frequency
      integer numhdx(25,0:mxrec)
      real binhdx(26)

      integer ipcp(51),idays(366)
      integer nfog(mxrec),nice(mxrec)
      logical lpmode,lcompr,lend,lscale,lsum,lprt
      logical lfog,lice

      logical lcalpuff,lcalm,lcalmprt,lrhmode,ldb

c --- Define logicals for printing weather record for current hour
      logical lweather,lhour

      character*4 title(20)
      character*3 cfoglb
      character*1 alpar
      character*12 verson,level
c
      common/chiprt/nta,nsum,iscode(62),lsum,ijcode
      common/scale/ascale(62),bscale(62),iscale(62),lscale
      common/acode/iyr,iday,ihour,icode,istab,ws,iwd,idpbl
      common/bcode/icd,imet2,icd1,icd3,idayhr,tdegk,qback,qsat

C --- variables for CD144 data
      integer yy,mm,dd,hr,jd,jday
      logical lmet
      character*70 metdat, fogmet, rhout,lstout, fogdat
      character*8 pwth
      character*20 weather
      character*81 line
C
c      namelist/opts/ijcode,lscale,iscode,lprt,lsum,nsum,
c     &              lcalmprt,lrhmode,ldb,
c     &              lmet,metdat 
      namelist/opts/lprt,
     &              lcalmprt,lrhmode,ldb,
     &              lmet,metdat,fogmet,rhout,lstout,fogdat 
c
c --- initialize common variables in block data
      data nfog/mxrec*0/,nice/mxrec*0/
      data nbkfog/0/,nbkice/0/,ntothr/0/
      data npifog/0/,npiice/0/
      data nlines/0/,alpar/'('/
      data lend/.false./,lprt/.false./
      data lcalmprt/.false./,lrhmode/.false./,ldb/.false./
C      
      DATA lmet/.false./, metdat/'metdat.144'/
      data fogmet/'foghrs.met'/, rhout/'postrm.rh'/
      data lstout/'postrm.lst'/, fogdat/'fog.dat'/
      DATA yy/0/,mm/0/, dd/0/,hr/0/,jd/0/,jday/0/  

c --- Set lower boundaries for the delta-RH bins
      data bindrh/-5.,-4.,-3.,-2.,-1.,
     &             0., 1., 2., 3., 4.,
     &             5., 6., 7., 8., 9.,
     &            10.,11.,12.,13.,14.,
     &            15.,16.,17.,18.,19.,
     &            20.,2000./
c --- Set lower boundaries for the delta-temperature bins (F)
      data bindt/ 0.0,0.5,1.0,1.5,
     &            2.0,2.5,3.0,3.5,
     &            4.0,4.5,5.0,5.5,
     &            6.0,6.5,7.0,7.5,
     &            8.0,8.5,9.0,9.5,
     &           10.0,2000./
c --- Set lower boundaries for the HUMIDEX bins (F)
      data binhdx/ 81., 83., 85., 87., 89.,
     &             91., 93., 95., 97., 99.,
     &            101.,103.,105.,107.,109.,
     &            111.,113.,115.,117.,119.,
     &            121.,123.,125.,127.,129.,
     &            2000./

c --- Process control file
      open(5,file='postrm.inp',status='old')
      read(5,opts)
c
c --- open files
c
      open(6,file=lstout)
      open(10,file=fogdat,form='unformatted',status='old')
c
c --- compute runtime stats
      call runtime(6,dtsec)
c
c --- set version/level numbers
      verson='1.2'
      level='060214'
      write(6,1200)verson,level
1200  format(//11x,'POSTRM -- Receptor-mode Postprocessor',10x,
     1 'Version: ',a12,4x,'Level: ',a12/1x,13('**********')/)
c
c --- write inputs
      write(6,opts)
C
C --- open CD144 data:
      if (lmet) then
         write(6,*)
         write(6,*)'Opening', metdat
         write(6,*)'For Meteorological Data'
         open(17,file=metdat)
      endif
c
c --- Set up file for reporting weather codes (conditional)
      if(lprt .AND. lmet) then
         lweather=.true.
         open(9,file=fogmet,status='unknown')
         write(9,1200)verson,level
         write(9,opts)
         write(9,*)
      else
         lweather=.false.
      endif

c
c --- read scaling factors
      if(lscale)then
         if(lsum)then
            nscale=nsum
         else
            nscale=1
         endif
         do 15 i=1,nscale
         read(5,*)iscale(i),ascale(i),bscale(i)
15       continue
      endif
c
c --- check validity of nsum, iscode values
      if(lsum)then
         if(nsum.le.1)then
            write(6,22)nsum,lsum
22          format(//1x,'ERROR -- NSUM must be > 1 when LSUM is .true.'/
     1       2x,'NSUM = ',i3,3x,'LSUM = ',l1)
             stop
         endif
c
         nsumm1=nsum-1
         do 32 i=1,nsumm1
         if(iscode(i+1).le.iscode(i))then
            write(6,26)(n,iscode(n),n=1,nsum)
26          format(//1x,'ERROR  -- values of source codes (ISCODE) ',
     1      'must be entered in order of increasing numerical value'/
     2      1x,'For example, enter ISCODE = 1,2,101, not ISCODE = ',
     3      '101,2,1 '/5x,'I',6x,'ISCODE(I)'//62(3x,i2,9x,i3/))
            stop
         endif
32       continue
      endif
c
c --- check validity of iscale array values
      if(lscale)then
         if(lsum)then
            do 42 i=1,nscale
            if(iscale(i).eq.iscode(i))go to 42
            write(6,40)lscale,lsum,(n,iscode(n),iscale(n),n=1,nscale)
40          format(//1x,'ERROR -- When LSCALE and LSUM are true, ',
     1      'source code arrays (ISCALE & ISCODE) must be equal and '/
     2      1x,'source code values must be entered in order of ',
     3      ' increasing numerical value -- For example,'/
     4      1x,'enter iscale = 1,2,101, not iscale = 101,2,1'/
     5      1x,'LSCALE = ',l1,3x,'LSUM = ',l1/
     6      4x,'I',6x,'ISCODE(I)',5x,'ISCALE(I)'//
     7      62(3x,i2,8x,i4,10x,i4/))
            stop
42          continue
         else
            if(iscale(1).ne.ijcode)then
               write(6,44)lscale,lsum,ijcode,iscale(1)
44             format(//1x,'ERROR -- When LSCALE is true and LSUM is ',
     1         'false, iscale(1) must equal ijcode'/
     2         1x,'LSCALE = ',l1,3x,'LSUM = ',l1/
     3         1x,'IJCODE = ',i4,3x,'ISCALE(1) = ',i4)
               stop
            endif
         endif
      endif
c
c --- read/write header record # 1 of FOG output file
      read(10)title,nnrec,npts,ipcp,iyr,idays,lpmode

c --- Set logical to T if IDAYS array is filled with zeroes
      lcalpuff=.FALSE.
      ndays=0
      do i=1,366
         ndays=ndays+idays(i)
      enddo
      if(ndays.EQ.0) lcalpuff=.TRUE.

      write(6,16)title,nnrec,npts,iyr,lpmode,ipcp,idays
      if(lweather) write(9,16)title,nnrec,npts,iyr,lpmode,ipcp,idays
16    format(//1x,'HEADER RECORD DATA'//
     1         1x,'title:  ',20a4/
     1         1x,'nnrec:  ',i8/
     2         1x,'npts:   ',i8/
     3         1x,'iyr:    ',i8/
     4         1x,'lpmode: ',l8/
     5         1x,'ipcp:   ',51i2/
     6         1x,'idays:  ',10(10i1,1x)/
     7         1x,'        ',10(10i1,1x)/
     8         1x,'        ',10(10i1,1x)/
     8         1x,'        ',10(10i1,1x))
c
c --- check that output is receptor-mode output
      if(lpmode)then
         write(6,54)lpmode
54       format(//1x,'ERROR -- POSTRM can be used with RECEPTOR-MODE ',
     1   'output only'//1x,'lpmode: ',l1)
         stop
      endif

      if(lcalpuff) then
c ---    Read a second header record for compression logical
         read(10) lcompr
         nrec=nnrec
      else
c ---    decode compression-mode flag (not needed for CALPUFF file)
c ---       if nnrec > 1000, output arrays are compressed
c ---       if nnrec < 1000, output arrays are not compressed
         if(nnrec.gt.1000)then
            lcompr=.true.
            nrec=mod(nnrec,1000)
         else
            lcompr=.false.
            nrec=nnrec
         endif
      endif
c
c --- compute no. records each hour
      npp=0
      do 65 i=1,51
65    npp=npp+ipcp(i)
c --- nta is the number of records for each hour
      nta=1+npp
c
c --- compute no. days in this run
      ndys=0
      do 75 i=1,366
75    ndys=ndys+idays(i)
      write(6,77)ndys,idays
77    format(//1x,'Total number of days included in this run: ',i3//
     1 1x,'(0=not included; 1=included)'//
     2 3(/1x,10(10i1,3x)/),/1x,6(10i1,3x),6i1)
c
c --- REMOVE source contribution QA for CALPUFF applications
      if(.NOT.lcalpuff) then
      write(6,85)
85    format(///1x,'Source contributions from the following ',
     1 'point sources are available: '//1x,'(0=not available; ',
     2 '1=available)'//1x,'Point source number',5x,'Availability')
      do 95 i=1,npts
      write(6,90)i,ipcp(i)
90    format(9x,i2,19x,i1)
95    continue
      endif
c
      if(lsum)then
         write(6,105)nsum,(alpar,nxx,iscode(nxx),nxx=1,nsum)
105      format(///1x,'The concentration data used in this run of ',
     1   'POSTRM is the sum of the following ',i3,' sets'/2x,'of ',
     2   'source contribution data: '//1x,'ISCODE = ',
     3   9(a1,i1,')',i4,2x),6(/9x,10(a1,i2,')',i4,1x)))
         write(6,110)
110      format(/9x,'Where the concentration data associated with ',
     1   'each ISCODE is identified as follows: '//1x,9x,'ISCODE',
     2   18x,'Concentration data'//1x,9x,'101-150',2x,'Source ',
     3   'contribution - point source number "ISCODE - 100" (not ',
     4   'including background)'/
     5   10x,'999',6x,'Total concentration (including background))'/
     6   10x,'1000',5x,'Background')
      else
         write(6,115)ijcode
115      format(//'0',9x,'IJCODE = ',i4)
         write(6,117)
117      format(/1x,'IJCODE specifies the concentration ',
     1   'output data used in this run'//1x,'The concentration ',
     2   'data associated with IJCODE is as follows: '//
     3   10x,'IJCODE',18x,'Concentration data'//
     4   10x,'101-150',2x,'Source contribution - point source ',
     5   'number "IJCODE - 100" (not including background)'/
     6   10x,'999',6x,'Total concentration (including background))')
      endif
c
c --- write scaling factors
      if(lscale)then
         write(6,123)
123      format(///1x,'Concentration scaling option requested -- ',
     1    'concentrations will be scaled in the following way: '//20x,
     2    'CONC.(new) = ASCALE * CONC.(old) + BSCALE'//
     3    1x,'Where ASCALE and BSCALE are input by the user for each ',
     4    'source contribution group'///
     5    1x,'Source',10x,'Source code',6x,'ASCALE',7x,'BSCALE',
     5    /2x,'Group',12x,'(ISCALE)',19x,' (g/m**3)'/)
         do 140 i=1,nscale
         itmp=iscale(i)
         jtmp=itmp-100
         if(itmp.ge.101.and.itmp.le.150)write(6,132)jtmp,itmp,
     1    ascale(i),bscale(i)
132      format(1x,'point ',i2,10x,i4,5x,1pe12.5,3x,1pe12.5)
         if(itmp.eq.999)write(6,134)itmp,ascale(i),bscale(i)
134      format(1x,'all sources',7x,i4,5x,1pe12.5,3x,1pe12.5)
         if(itmp.eq.1000)write(6,136)itmp,ascale(i),bscale(i)
136      format(1x,'background',8x,i4,5x,1pe12.5,3x,1pe12.5)
140      continue
      endif
c
c --- write heading for tables
      if(lprt)then
         write(6,155)
155      format('1',35x,'Hours of plume-induced fog or ice'//1x,'Yr',
     1   2x,'Day',2x,'Hour',2x,'Receptor',2x,'Stab.',3x,'Wind',
     2   4x,'Vector',3x,'Temp.',3x,'Rel.',3x,'Background',2x,
     3   'Saturation',3x,'Total',4x,'Total-Sat.',3x,'Type',3x,'Weather'/
     4   19x,'No.',11x,'Speed',3x,'Wind Dir.',9x,'Hum.'/33x,'(m/s)',
     5   5x,'(deg)',2x,'(deg K)',2x,'(%)',5x,'(g/m**3)',4x,
     6   '(g/m**3)',3x,'(g/m**3)',3x,'(g/m**3)'/)
         if(lweather) write(9,156)
156      format('1',//5x,'Weather for hours of plume-induced fog or ice',
     1   //1x,'Yr',2x,'Day',2x,'Hour',1x,'Stab.',3x,'Wind',
     2   4x,'Vector',3x,'Temp.',3x,'Rel.',7x,'Weather'/
     3   22x,'Speed',3x,'Wind Dir.',9x,'Hum.'/22x,'(m/s)',
     4   5x,'(deg)',2x,'(deg K)',2x,'(%)',/)
      endif
c
c --- Initialize RH histogram processing (conditional)
c -----------------------------------------------------
      if(LRHMODE) then
         open(7,file=rhout)
         write(7,1200)verson,level
         write(7,16)title,nnrec,npts,iyr,lpmode,ipcp,idays
         write(7,*)
         do ir=0,mxrec
            rhavg(ir)=0.0
            do i=1,20
               numrh(i,ir)=0
            enddo
            do i=1,25
               numhdx(i,ir)=0
            enddo
         enddo
         do ir=1,mxrec
            do i=1,26
               numdrh(i,ir)=0
            enddo
            do i=1,21
               numdt(i,ir)=0
            enddo
         enddo
      endif

c --- Set counter for CALM hours
      ncalm=0
c
      do 1000 jday=1,366
c --- if(idays(jday).ne.1)go to 1000
      if(idays(jday).ne.1 .AND. .not.LCALPUFF)go to 1000
c
      do 500 jhr=1,24

c --- Set top of loop for processing all CALPUFF hours
c --- (this negates explicit loops over hours and days)
1     continue

c --- Initialize logical for writing hourly weather codes
      lhour=.false.

c --- Read in concentration (CHI is TOTAL, including background)
      if(LCALPUFF) then
         call getpuf(nrec,chi,chiin,qbackr,qsatr,lend,lcompr)
      else
         call getdat(nrec,chi,chiin,lend,lcompr)
      endif
      if(lend)go to 9999
c
      ntothr=ntothr+1
      lfog=.false.
      lice=.false.
      wd=iwd
      irh=NINT(100.*qback/qsat)

c --- Identify calm hours
      if(ws.LE.0.5) then
         lcalm=.TRUE.
      else
         lcalm=.FALSE.
      endif

c --- Debug Writes
      if(LDB) then
         write(6,*)'iyr,iday,ihour = ',iyr,iday,ihour
         write(6,*)'qback,qsat,irh = ',qback,qsat,irh
         do n=1,nrec
            write(6,*)'chi,qbackr,qsatr   = ',n,chi(n),qbackr(n),
     &                 qsatr(n)
         enddo
      endif
c
c --- Read CD144 data:--------------------------------------
c --- Current date is from FOG.DAT record (iyr,iday,ihour)
c --- (ihour is 1-24)
      if(lmet) then
         icur2 = iyr*100000+iday*100+ihour
         ibeg2 = yy*100000+jd*100+hr
 169     continue
         if(ibeg2.lt.icur2) then
            read(17,161,end=167) line
 161        format(A81)
            read(line,162) yy,mm,dd,hr,pwth
 162        format(5x,4i2,10x,a8)
c ---       Convert YY to YYYY year format
            if(yy.gt.48) then
               yy=yy+1900
            else 
               yy=yy+2000
            endif
            call JULDAY(yy,mm,dd,Jd)
            if(hr.eq.0) then
c ---          Roll back 1 hour to find previous day's date
               nhrinc=-1
               call INCR(yy,jd,hr,nhrinc)
c ---          Now add the hour back on to get hour 24
               hr = hr+1
            endif
            ibeg2 = yy*100000+jd*100+hr
            goto 169
         endif
 167     continue
         if(ibeg2.eq.icur2) then
c ---       Correct time found
c ---       Process weather code from cd144:
            call PROC144RM(pwth,weather)
         endif
      endif
C---------------------------------------------------------------------

c --- Update RH frequency array, using 5% bins (conditional), and sums
c --- for average RH
      if(LRHMODE) then
c ---    Set ambient
         k=(irh+4)/5
         k=MAX(k,1)
         numrh(k,0)=numrh(k,0)+1
         rhavg(0)=rhavg(0)+irh
c ---    Set each receptor
         do ir=1,nrec
            if(chi(ir).LE.0.0) then
c ---          No plume impact
               irhrec=irh
            elseif(LCALPUFF) then
               irhrec=NINT(100.*chi(ir)/qsatr(ir))
            else
               irhrec=NINT(100.*chi(ir)/qsat)
            endif
            irhrec=MIN(100,irhrec)
            irhrec=MAX(1,irhrec)
            k=(irhrec+4)/5
            numrh(k,ir)=numrh(k,ir)+1
            rhavg(ir)=rhavg(ir)+irhrec
         enddo
c ---    Also delta-RH and delta-T, and HUMIDEX
         call DELTA(nrec,chi,qsatr,numdrh,numdt,numhdx,
     &              bindrh,bindt,binhdx,LCALPUFF)
      endif
c
c --- compute background fog/icing hours
      if(irh.ge.100)then
         if(tdegk.lt.273.16)then
            nbkice=nbkice+1
         else
            nbkfog=nbkfog+1
         endif
         go to 290
      endif
c
c --- compute plume-induced fog/ice hours
      qsat0=qsat
      qback0=qback
      if(tdegk.lt.273.16)then
         do 180 i=1,nrec
         if(LCALPUFF) then
            qsat0=qsatr(i)
            qback0=qbackr(i)
         endif
         if(chi(i).gt.qsat0)then
            lhour=.true.
            lice=.true.
            nice(i)=nice(i)+1
            if(lprt .AND. (lcalmprt .OR. .NOT.lcalm))then
               cfoglb='ICE'
               cliq=chi(i)-qsat0
               nlines=nlines+1
               if(nlines.gt.50)then
                  nlines=1
                  write(6,155)
               endif
               write(6,175)iyr,iday,ihour,i,istab,ws,wd,tdegk,irh,
     1         qback0,qsat0,chi(i),cliq,cfoglb,pwth,weather
         
175            format(i4,1x,i3,3x,i2,4x,i4,7x,i1,4x,f5.1,5x,f5.1,3x,
     1         f5.1,3x,i3,5x,f7.3,5x,f7.3,4x,f7.3,4x,f7.3,5x,a3,
     &         3X,A8,1X,A16)
            endif
         endif
180      continue
      else
         do 280 i=1,nrec
         if(LCALPUFF) then
            qsat0=qsatr(i)
            qback0=qbackr(i)
         endif
         if(chi(i).gt.qsat0)then
            lhour=.true.
            lfog=.true.
            nfog(i)=nfog(i)+1
            if(lprt .AND. (lcalmprt .OR. .NOT.lcalm))then
               cfoglb='FOG'
               cliq=chi(i)-qsat0
               nlines=nlines+1
               if(nlines.gt.50)then
                  nlines=1
                  write(6,155)
               endif
               write(6,175)iyr,iday,ihour,i,istab,ws,wd,tdegk,irh,
     1         qback0,qsat0,chi(i),cliq,cfoglb,pwth,weather
            endif
         endif
280      continue
      endif

c --- Write weather information to foghour.met file if induced
c --- fog/ice is found this hour
      if(lweather .AND. lhour) then
         write(9,285) iyr,iday,ihour,istab,ws,wd,tdegk,irh,
     &                pwth,weather
285      format(i4,1x,i3,3x,i2,4x,i1,4x,f5.1,5x,f5.1,3x,
     1          f5.1,3x,i3,5x,3X,A8,1X,A16)
      endif

290   continue
c
c --- sum plume-induced fog/ice hours
      if(lfog)then
         npifog=npifog+1
         if(lcalm) ncalm=ncalm+1
      else if(lice)then
         npiice=npiice+1
         if(lcalm) ncalm=ncalm+1
      endif

      if(LCALPUFF) goto 1

500   continue
1000  continue
9999  continue

c
      if(lprt.and.npifog+npiice.eq.0)write(6,1002)
1002  format(25x,'No plume-induced fog/ice hours')
c
c --- write fog/ice table
      tothr=float(ntothr)
      xtmp=100./tothr
      nrec1=1
      nrec2=0
1100  continue
      nrec2=nrec2+50
      nrec2=min0(nrec2,nrec)
      write(6,1110)
1110  format('1','Receptor',4x,'Plume-induced',4x,'Plume-induced',4x,
     1 'Plume-induced'/17x,'fog',14x,'ice',14x,'total'/
     2 14x,3('(hrs)',3x,'(%)',6x)/)
c
      do 1120 i=nrec1,nrec2
      ntot=nfog(i)+nice(i)
      pc1=nfog(i)*xtmp
      pc2=nice(i)*xtmp
      pc3=ntot*xtmp
      write(6,1115)i,nfog(i),pc1,nice(i),pc2,ntot,pc3
1115  format(3x,i4,6x,3(i5,f7.2,5x))
1120  continue
      nrec1=nrec1+50
      if(nrec2.lt.nrec)go to 1100
c
c --- write summary statistics
      nbktot=nbkfog+nbkice
      npitot=npifog+npiice
      pc1=nbkfog*xtmp
      pc2=nbkice*xtmp
      pc3=nbktot*xtmp
      pc4=npifog*xtmp
      pc5=npiice*xtmp
      pc6=npitot*xtmp
      write(6,1130)ntothr,nbkfog,pc1,nbkice,pc2,nbktot,pc3,npifog,pc4,
     1 npiice,pc5,npitot,pc6,ncalm
1130  format('1','Total number of hours = ',i5//1x,
     1 '"Background fog" hours =',i5,3x,f6.2,' (%)'//1x,
     2 '"Background ice" hours =',i5,3x,f6.2,' (%)'//1x,
     3 'Total "background fog + ice" hours =',i5,3x,f6.2,' (%)'//1x,
     4 '"Plume-induced fog" hours =',i5,3x,f6.2,' (%)'//1x,
     5 '"Plume-induced ice" hours =',i5,3x,f6.2,' (%)'//1x,
     6 'Total "plume-induced fog + ice" hours =',i5,3x,f6.2,' (%)'//
     a 1x,'Total "calm plume-induced fog + ice" hours =',i5,//
     7 1x,'------------------------'//1x,'Notes:  ',
     8 '(1) "Background fog" hours --> background RH = 100%, ',
     9 'Temp > 32 deg. F'/9x,
     a '(2) "Background ice" hours --> background RH = 100%, ',
     b 'Temp <= 32 deg. F'/9x,
     c '(3) "Plume-induced fog" hours --> background RH < 100%, ',
     d 'Total Conc. > Saturation, Temp > 32 deg. F'/9x,
     e '(4) "Plume-induced ice" hours --> background RH < 100%, ',
     f 'Total Conc. > Saturation, Temp <= 32 deg. F'//)
c
c --- Write RH frequency tables with 20 receptors in each (conditional)
c --- using 5% binning
      if(LRHMODE) then
c ---    Report number of hours
         write(7,*)
         write(7,*)'Number of Hours Processed:',ntothr
         write(7,*)
c ---    Convert RH sums to averages
         do ir=0,nrec
            rhavg(ir)=rhavg(ir)/tothr
         enddo
c ---    Convert RH counts from hours to %
         do ir=0,nrec
            do i=1,20
               percent(i,ir)=100.*numrh(i,ir)/tothr
            enddo
         enddo
         ntab=nrec/20
         nrectab=ntab*20
         if(nrec.GT.nrectab) ntab=ntab+1
         do itab=1,ntab
            n1=(itab-1)*20+1
            n2=n1+19
            n2=MIN(n2,nrec)
            write(7,*)
            write(7,*)
            write(7,*)
            write(7,*)'Relative Humidity Histogram for Receptors',n1,n2
            write(7,*)'                     (%)'
            write(7,*)
            write(7,'(a,20i6)') ' RH Bin   Amb',(i,i=n1,n2)
            write(7,*)
            do i=20,1,-1
               l=5*i
               k=l-4
               write(7,'(i2,a,i3,21f6.1)') k,' -',l,percent(i,0),
     &                                   (percent(i,j),j=n1,n2)
            enddo
            write(7,*)
            write(7,'(a,21f6.2)') 'AVG RH%',rhavg(0),(rhavg(j),j=n1,n2)
         enddo

c ---    Set a large real for use when *'s are wanted in output field
         aster=10000.

c ---    Tables for the corresponding delta-RH
c ---    Convert delta-RH counts from hours to %
         do ir=1,nrec
            do i=1,26
               percent(i,ir)=100.*numdrh(i,ir)/tothr
            enddo
         enddo
         do itab=1,ntab
            n1=(itab-1)*20+1
            n2=n1+19
            n2=MIN(n2,nrec)
            write(7,*)
            write(7,*)
            write(7,*)
            write(7,*)'DELTA-RH(%) Histogram for Receptors',n1,n2
            write(7,*)'               (%)'
            write(7,*)
            write(7,'(a,20i6)') ' dRH(%) Bin ',(i,i=n1,n2)
            write(7,*)
            do i=26,2,-1
               write(7,'(f4.1,a4,f4.1,20f6.1)') bindrh(i),' to ',
     &                            bindrh(i+1),(percent(i,j),j=n1,n2)
            enddo
c ---       Lower limit of lowest bin range should be undefined
            write(7,'(f4.1,a4,f4.1,20f6.1)') aster,' to ',
     &                            bindrh(2),(percent(1,j),j=n1,n2)
            write(7,*)
         enddo

c ---    Tables for the corresponding delta-Temperature
c ---    Convert delta-T counts from hours to %
         do ir=1,nrec
            do i=1,21
               percent(i,ir)=100.*numdt(i,ir)/tothr
            enddo
         enddo
         do itab=1,ntab
            n1=(itab-1)*20+1
            n2=n1+19
            n2=MIN(n2,nrec)
            write(7,*)
            write(7,*)
            write(7,*)
            write(7,*)'DELTA-T(F) Histogram for Receptors',n1,n2
            write(7,*)'              (%)'
            write(7,*)
            write(7,'(a,20i6)') ' dT(F) Bin  ',(i,i=n1,n2)
            write(7,*)
            do i=21,2,-1
               write(7,'(f4.1,a4,f4.1,20f6.1)') bindt(i),' to ',
     &                            bindt(i+1),(percent(i,j),j=n1,n2)
            enddo
c ---       Lower limit of lowest bin range should be undefined
            write(7,'(f4.1,a4,f4.1,20f6.1)') aster,' to ',
     &                            bindt(2),(percent(1,j),j=n1,n2)
            write(7,*)
         enddo

c ---    Tables for the corresponding HUMIDEX
         do itab=1,ntab
            n1=(itab-1)*20+1
            n2=n1+19
            n2=MIN(n2,nrec)
            write(7,*)
            write(7,*)
            write(7,*)
            write(7,*)'HUMIDEX(F) Histogram for Receptors',n1,n2
            write(7,*)'            (hours)'
            write(7,*)
            write(7,'(a,20i6)') ' H Bin         Amb',(i,i=n1,n2)
            write(7,*)
            do i=25,2,-1
               nbin=NINT(binhdx(i))
               write(7,'(f4.0,a4,f4.0,21i6)') binhdx(i),' to ',
     &                            binhdx(i+1),numhdx(i,0),
     &                            (numhdx(i,j),j=n1,n2)
               if(nbin.EQ.103) write(7,*)'--------------------'
               if(nbin.EQ.113) write(7,*)'--------------------'
            enddo
c ---       Lower limit of lowest bin range should be undefined
            write(7,'(f4.0,a4,f4.0,21i6)') aster,' to ',
     &                            binhdx(2),numhdx(1,0),
     &                            (numhdx(1,j),j=n1,n2)
            write(7,*)
            write(7,*)'Less than 84F    --- No Discomfort'
            write(7,*)'84F  to 102F     --- Some Discomfort'
            write(7,*)'103F to 113F     --- Great Discomfort'
            write(7,*)'                     (Avoid Exertion)'
            write(7,*)'113F to 129F     --- Dangerous'
            write(7,*)'Above 129F       --- Heat Stroke Imminent'
            write(7,*)
         enddo

      endif
c
      call runtime(6,dtsec)
c
      stop
      end
      block data
c
      logical lsum,lscale
      common/chiprt/nta,nsum,iscode(62),lsum,ijcode
      common/scale/ascale(62),bscale(62),iscale(62),lscale
      data ijcode/999/,lsum/.false./,iscode/62*0/
      data nsum/0/,lscale/.false./
      end
      subroutine getdat(nrec,chi,chiin,lend,lcompr)
c
      parameter(mxrec=5000)
c
      real chia(mxrec)
      real chi(nrec),chiin(nrec)
      logical lsum,lend,lcompr,lscale
c
      common/chiprt/nta,nsum,iscode(62),lsum,ijcode
      common/acode/iyr,iday,ihour,icode,istab,ws,iwd,idpbl
      common/bcode/icd,imet2,icd1,icd3,idayhr,tdegk,qback,qsat
      common/scale/ascale(62),bscale(62),iscale(62),lscale
c
      if(.not.lsum)go to 50
      index=1
c --- nta is the number of records per hour (with lcompr=.false.)
c --- (2*nta records/hr with lcompr=.true.)
      do 10 i=1,nta
      if(lcompr)go to 101
c --- array compression option not used
      read(10,end=999)idayhr,icd,imet2,tdegk,qback,qsat,chi
      icd1=icd/1000
      icode=icd-1000*icd1
      if(icode.ne.iscode(index))go to 10
      go to 102
101   continue
c --- array compression option used
      read(10,end=999)ii
      call rd(idayhr,icd,imet2,tdegk,qback,qsat,ii,chiin)
      icd1=icd/1000
      icode=icd-1000*icd1
      if(icode.ne.iscode(index))go to 10
c --- expand the concentration array
      call xpnd(ii,nrec,chiin,chi)
102   continue
      if(lscale)call xscale(index,nrec,chi)
      index=index+1
      if(index.ne.2)go to 8
      do 7 j=1,nrec
7     chia(j)=chi(j)
      go to 10
8     do 9 j=1,nrec
9     chia(j)=chia(j)+chi(j)
10    continue
      xback=0.0
      if(index.eq.nsum+1)go to 15
c --- check if background should be added
      if(index.eq.nsum.and.iscode(index).eq.1000)then
         xback=qback
         if(lscale)xback=ascale(index)*xback+bscale(index)
         go to 15
      endif
      write(6,12)nsum,index,iscode
12    format(/////'0','ERROR -- invalid value of nsum or iscode'/
     1 /1x,'nsum = ',i5,5x,'index = ',i5,5x,
     2 6(/'0',10(i5,2x))/'0',2(i5,2x))
      stop
15    continue
c --- xback is zero except if iscode(nsum) is 1000
      do 17 j=1,nrec
17    chi(j)=chia(j)+xback
      go to 60
50    continue
      do 55 i=1,nta
      if(lcompr)go to 105
c     array compression option not used
      read(10,end=999)idayhr,icd,imet2,tdegk,qback,qsat,chi
      icd1=icd/1000
      icode=icd-1000*icd1
      if(icode.ne.ijcode)go to 55
      go to 106
105   continue
c     array compression option used
      read(10,end=999)ii
      call rd (idayhr,icd,imet2,tdegk,qback,qsat,ii,chiin)
      icd1=icd/1000
      icode=icd-1000*icd1
      if(icode.ne.ijcode)go to 55
c --- expand the concentration array
      call xpnd(ii,nrec,chiin,chi)
106   continue
      if(lscale)call xscale(1,nrec,chi)
      nta2=nta-i
      if(nta2.eq.0)go to 60
c
c --- skip remaining records in this hour
      if(lcompr)go to 153
c --- array compression not used
      do 54 j=1,nta2
      read(10)ni1,ni2,ni3,s1,s2,s3,chiin
54    continue
      go to 60
153   continue
c --- array compression used
      do 154 j=1,nta2
      read(10)ii
      call rd(ni1,ni2,ni3,s1,s2,s3,ii,chiin)
154   continue
      go to 60
55    continue
      write(6,56)ijcode
56    format(//1x,'ERROR -- error in the specification of ijcode'//
     1 1x,'ijcode = ',i5)
      stop
60    continue
c
c --- decode idayhr
c
      iday=idayhr/100
      ihour=idayhr-iday*100
c
c --- decode icd
c
      iws=icd1/10
      ws=float(iws)/10.
      istab=icd1-10*iws
      icd3=idayhr*10+istab
c
c --- decode imet2
c
      iwd=imet2/10000
      idpbl=imet2-iwd*10000
      return
999   lend=.true.
      return
      end
      subroutine rd(idayhr,icd,imet2,tdegk,qback,qsat,ii,chiin)
      real chiin(ii)
c
      read(10)idayhr,icd,imet2,tdegk,qback,qsat,chiin
      return
      end
      subroutine xpnd(ii,nrec,chiin,chi)
c
      real chi(nrec),chiin(ii)
c
c --- [DGS 001003]
c --- Some versions of the compression routine appear to replace a 
c --- single zero with the code '-1.0', so remove the (ii.eq.nrec) test
c ---     if(ii.eq.nrec)go to 200

      jj=0
      do 100 i=1,ii
      if(chiin(i).gt.0.0)go to 55
      nzero=-chiin(i)+0.0001
      do 40 j=1,nzero
      jj=jj+1
      chi(jj)=0.0
40    continue
      go to 100
55    continue
      jj=jj+1
      chi(jj)=chiin(i)
100   continue
      return
c --- 200  continue
c ---      do 210 i=1,ii
c ---      chi(i)=chiin(i)
c --- 210  continue
c ---      return
      end
      subroutine xscale(n,nrec,chi)
c
      real chi(nrec)
      logical lscale
      common/scale/ascale(62),bscale(62),iscale(62),lscale
c
      a=ascale(n)
      b=bscale(n)
      do 100 i=1,nrec
      chi(i)=a*chi(i)+b
100   continue
c
      return
      end
      SUBROUTINE RUNTIME(IUNIT,DTSEC)
C
      SAVE IJUL1,NUMCL
C
      CHARACTER*11 TSTR
      CHARACTER*8  DSTR
C
      DATA NUMCL/1/
C
C --- COMPUTES CUMULATIVE CLOCK TIME DIFFERENCE FROM FIRST CALL OF
C --- RUNTIME TO EACH SUBSEQUENT CALL OF RUNTIME
C
C --- DETERMINE JULIAN DAY
      CALL DATE(DSTR)
C --- NUMCT IS A FUNCTION WHICH CONVERTS A CHARACTER INTO AN INTEGER
      IMO=10*NUMCT(DSTR(1:1))+NUMCT(DSTR(2:2))
      IDAY=10*NUMCT(DSTR(4:4))+NUMCT(DSTR(5:5))
      IYR=10*NUMCT(DSTR(7:7))+NUMCT(DSTR(8:8))
      CALL JULDAY(IYR,IMO,IDAY,IJUL)
C
C --- DETERMINE CURRENT TIME
      CALL TIME(TSTR)
C
C --- WRITE DATE, TIME TO SCREEN AND "IUNIT"
      WRITE(*,131)NUMCL,DSTR,TSTR
      IF(IUNIT.NE.0)WRITE(IUNIT,131)NUMCL,DSTR,TSTR
131   FORMAT(1X,'RUNTIME CALL NO.: ',I8,3x,'DATE: ',A8,3X,'TIME: ',A11)
C
      IHR=10*NUMCT(TSTR(1:1))+NUMCT(TSTR(2:2))
      IMIN=10*NUMCT(TSTR(4:4))+NUMCT(TSTR(5:5))
      RSEC=10.*FLOAT(NUMCT(TSTR(7:7)))+FLOAT(NUMCT(TSTR(8:8)))+
     1 FLOAT(NUMCT(TSTR(10:10)))/10.+FLOAT(NUMCT(TSTR(11:11)))/100.
      RTSEC=3600.*IHR+60.*IMIN+RSEC
C
C --- IF FIRST CALL, JUST STORE DATA
      IF(NUMCL.EQ.1)THEN
         IJUL1=IJUL
         RTSEC1=RTSEC
         NUMCL=NUMCL+1
         RETURN
      ELSE
C ---    COMPUTE DELTA TIME (SECONDS) -- (ASSUMES SAME YEAR)
         IF(IJUL.EQ.IJUL1)THEN
            DTSEC=RTSEC-RTSEC1
         ELSE
            DDAY=IJUL-IJUL1
            DTSEC=(86400.-RTSEC1)+86400.*(DDAY-1.)+RTSEC
         ENDIF
      ENDIF
C
C --- WRITE ELAPSED TIME FROM FIRST CALL TO RTIME
      WRITE(*,133)DTSEC
      IF(IUNIT.NE.0)WRITE(IUNIT,133)DTSEC
133   FORMAT(1X,'DELTA TIME: ',F10.2,' (SEC)')
C
C --- UPDATE CALL COUNTER
      NUMCL=NUMCL+1
C
      RETURN
      END
      FUNCTION NUMCT(C1)
C
C --- IO6 = I/O UNIT FOR ERROR MESSAGE
      PARAMETER(IO6=6)
C
      CHARACTER*1 CNUM(10),C1
C
      DATA CNUM/'0','1','2','3','4','5','6','7','8','9'/
C
C --- CONVERT A CHARACTER VARIABLE INTO AN INTEGER
      DO 10 I=1,10
      IF(C1.EQ.CNUM(I))THEN
         NUMCT=I-1
         RETURN
      ENDIF
10    CONTINUE
C
      WRITE(IO6,12)C1
12    FORMAT(1X,'ERROR IN SUBR. NUMCT -- INVALID CHARACTER -- ',
     1 'C1 = ',A1)
      STOP1
      END
      SUBROUTINE JULDAY(IYR,IMO,IDAY,IJULDY)
C
      INTEGER KDAY(12)
      DATA KDAY/0,31,59,90,120,151,181,212,243,273,304,334/
C
C --- COMPUTE THE JULIAN DAY NUMBER FROM THE GREGORIAN DATE (MONTH, DAY)
C
      IJULDY=KDAY(IMO)+IDAY
      IF(IMO.LE.2)RETURN
      IF(MOD(IYR,4).EQ.0)IJULDY=IJULDY+1
C
      RETURN
      END

c----------------------------------------------------------------------
      subroutine getpuf(nrec,chi,chiin,qbackr,qsatr,lend,lcompr)
c----------------------------------------------------------------------
      parameter(mxrec=5000)
c
      real chi(nrec),qbackr(nrec),qsatr(nrec),chiin(nrec)
      logical lsum,lend,lcompr,lscale
c
      common/chiprt/nta,nsum,iscode(62),lsum,ijcode
      common/acode/iyr,iday,ihour,icode,istab,ws,iwd,idpbl
      common/bcode/icd,imet2,icd1,icd3,idayhr,tdegk,qback,qsat
      common/scale/ascale(62),bscale(62),iscale(62),lscale

c
c --- Read ambient data record
      read(10,end=999) nyrjulhr,icd,imet2,tdegk,frh,qsat

c --- Compute QBACK
      qback=frh*qsat

c --- Split year from day-hour
      iyr=nyrjulhr/100000
      idayhr=nyrjulhr-iyr*100000
c
c --- Decode idayhr
      iday=idayhr/100
      ihour=idayhr-iday*100
c
c --- Decode icd
      icd1=icd/1000
      icode=icd-1000*icd1
      iws=icd1/10
      ws=float(iws)/10.
      istab=icd1-10*iws
      icd3=idayhr*10+istab
c
c --- Decode imet2
      iwd=imet2/10000
      idpbl=imet2-iwd*10000

c --- Read data array records
      if(lcompr) then
         read(10)ii
         call rdpuf(ii,chiin)
         call xpnd(ii,nrec,chiin,chi)
         read(10)ii
         call rdpuf(ii,chiin)
         call xpnd(ii,nrec,chiin,qbackr)
         read(10)ii
         call rdpuf(ii,chiin)
         call xpnd(ii,nrec,chiin,qsatr)
      else
         call rdpuf(nrec,chi)
         call rdpuf(nrec,qbackr)
         call rdpuf(nrec,qsatr)
      endif

c --- Apply scaling to plume concentration
      if(lscale)call xscale(1,nrec,chi)

c --- Add background to plume concentration to obtain TOTAL
      do i=1,nrec
         chi(i)=chi(i)+qbackr(i)
      enddo

      return

999   lend=.true.
      return

      end
c -------------------------------------------------------------------
      subroutine rdpuf(ii,chiin)
c -------------------------------------------------------------------
      real chiin(ii)
      character*15 cname
      read(10)cname,chiin
      return
      end
C -------------------------------------------------------------------
      subroutine PROC144RM(pwth,weather)
c -----------------------------------------------------------------
C                                                         PROC144RM
c --- PURPOSE:  PROCESS THE WEATHER CODE FROM CD144 

c --- INPUTS:
c       PWTH    - CHARACTER - WEATHER CODE
c       WEATHER - CHARACTER - WEATHR CODE INTERPRETATION

c --- OUTPUT:
c       PWTH    - CHARACTER - WEATHER CODE
c       WEATHER - CHARACTER - WEATHR CODE INTERPRETATION

c --- PROC144RM called by: POSTRM2
C
c --- INCR calls:     none
C ----------------------------------------------------------------------

      character*8 pwth
      character*20 weather

      weather = '                    '

      if (pwth.ne.'        ') then
c not missing

         if (pwth(2:2).ne.'0'.or.pwth(3:3).ne.'0') then
c ---       Liquid Precipitation Present
            if (pwth(4:4).ne.'0'.or.pwth(5:5).ne.'0'.or.
     &          pwth(6:6).ne.'0') then
c ---           Frozen Precipitation Also Present
               if (pwth(7:7).ge.'1' .AND. pwth(7:7).le.'3') then
                  weather = 'Rain Snow Fog     '
               elseif (pwth(8:8).eq.'5') then
                  weather = 'Rain Snow BlSnow  '
               elseif (pwth(8:8).ge.'1' .AND. pwth(8:8).le.'4') then
                  weather = 'Rain Snow Haze    '
               else
                  weather = 'Rain Snow         '
               endif
            elseif (pwth(7:7).ge.'1' .AND. pwth(7:7).le.'3') then
               weather = 'Rain Fog          '
            elseif (pwth(8:8).eq.'5') then
               weather = 'Rain BlSnow       '
            elseif (pwth(8:8).ge.'1' .AND. pwth(8:8).le.'4') then
               weather = 'Rain Haze         '
            else
               weather = 'Rain               '
            endif
         elseif (pwth(4:4).ne.'0'.or.pwth(5:5).ne.'0'
     &           .or.pwth(6:6).ne.'0') then
c ---       Frozen Precipitation (No Liquid)
            if (pwth(7:7).ge.'1' .AND. pwth(7:7).le.'3') then
               weather = 'Snow Fog          '
            elseif (pwth(8:8).eq.'5') then
               weather = 'Snow BlSnow       '
            elseif (pwth(8:8).ge.'1' .AND. pwth(8:8).le.'4') then
               weather = 'Snow Haze         '
            else
               weather = 'Snow              '
            endif
         elseif (pwth(7:7).ge.'1' .AND. pwth(7:7).le.'3') then
c ---       No Precipitation
            if (pwth(8:8).ge.'1' .AND. pwth(8:8).le.'4') then
               weather = 'Fog Haze          '
            else
               weather = 'Fog               '
            endif
         elseif (pwth(8:8).eq.'5') then
c ---       No Precipitation
            weather = 'BlSnow            '
         elseif (pwth(8:8).ge.'1' .AND. pwth(8:8).le.'4') then
c ---       No Precipitation
            weather = 'Haze              '
         else
            weather = '                    '
         endif

      endif

      return
      end
c----------------------------------------------------------------------
      subroutine incr(iyr,ijul,ihr,nhrinc)
c----------------------------------------------------------------------
c
c --- CALPUFF    Version: 5.2      Level:  980304                  INCR
c                J. Scire, SRC
c
c --- PURPOSE:  Increment the time and date by "NHRINC" hours
c
c --- UPDATE
c --- V5.0-V5.0     980304  (DGS): Allow for a negative "increment" of
c                                  up to 24 hours
c                   980304  (DGS): Allow for arbitrarily large nhrinc
c
c --- INPUTS:
c       IYR    - integer - Current year
c       IJUL   - integer - Current Julian day
c       IHR    - integer - Current hour (00-23)
c       NHRINC - integer - Time increment (hours)
c
c               NOTE: "NHRINC" must >= -24
c                      Hour is between 00-23
c
c --- OUTPUT:
c       IYR    - integer - Updated year
c       IJUL   - integer - Updated Julian day
c       IHR    - integer - Updated hour (00-23)
c
c --- INCR called by: INCRS, MET1, MET2, MET3, MET4, RDSIG,
c                     RDMET, RDISC, RDPLM, RDMET4
c --- INCR calls:     none
c----------------------------------------------------------------------
c
c --- Check nhrinc
      if(nhrinc.lt.-24) then
         write(*,*)'ERROR IN SUBR. INCR -- Invalid value of NHRINC ',
     1   '-- NHRINC = ',nhrinc
         stop
      endif

c --- Save increment remaining (needed if nhrinc > 8760)
      nleft=nhrinc
c
c --- Process change in hour
      if(nhrinc.gt.0)then
c
10       ninc=MIN0(nleft,8760)
         nleft=nleft-ninc
c
c ---    Increment time
         ihr=ihr+ninc
         if(ihr.le.23)return
c
c ---    Increment day
         ijul=ijul+ihr/24
         ihr=mod(ihr,24)
c
c ---    ILEAP = 0 (non-leap year) or 1 (leap year)
         if(mod(iyr,4).eq.0)then
            ileap=1
         else
            ileap=0
         endif
c
         if(ijul.gt.365+ileap) then
c ---       Update year
            iyr=iyr+1
            ijul=ijul-(365+ileap)
         endif
c
c ---    Repeat if more hours need to be added
         if(nleft.GT.0) goto 10
c
      elseif(nhrinc.lt.0)then
c ---    Decrement time
         ihr=ihr+nhrinc
         if(ihr.lt.0)then
            ihr=ihr+24
            ijul=ijul-1
            if(ijul.lt.1)then
               iyr=iyr-1
               if(mod(iyr,4).eq.0)then
                  ijul=366
               else
                  ijul=365
               endif
            endif
         endif
      endif
c
      return
      end
c----------------------------------------------------------------------
      subroutine DELTA(nrec,chi,qsatr,numdrh,numdt,numhdx,
     &                 bindrh,bindt,binhdx,LCALPUFF)
c----------------------------------------------------------------------

      parameter(mxrec=5000)

      real chi(mxrec),qsatr(mxrec)
      real bindrh(27),bindt(22),binhdx(26)
      integer numdrh(26,mxrec)
      integer numdt(21,mxrec)
      integer numhdx(25,0:mxrec)

      logical LCALPUFF

      common/bcode/icd,imet2,icd1,icd3,idayhr,tdegk,qback,qsat

c --- Ambient RH(%)
      rhamb=100.*qback/qsat

c --- Ambient HUMIDEX (C)
c --- Get dew-point temperature to nearest 0.1 K
      tdewk=tdegk
      if(qback.LT.qsat) then
10       tdewk=tdewk-.1
         call VSAT(tdewk,psat,csat)
         if(csat.GE.qback) goto 10
      endif
      eterm=6.11*EXP(5417.7530*(1./273.16-1./tdewk))
      hamb=(tdegk-273.16)+(eterm-10.)/1.8
c --- Place HUMIDEX in its bin
c --- Convert from C to F
      hamb=hamb*1.8+32.
      ibin=1
      do ib=1,20
         if(hamb.GE.binhdx(ib)) ibin=ib
      enddo
      numhdx(ibin,0)=numhdx(ibin,0)+1

c --- Loop over receptors
      do ir=1,nrec
c ---    Compute change in relative humidity (%)
         if(chi(ir).LE.0.0) then
c ---       No plume impact
            drh=0.0
         elseif(LCALPUFF) then
            rh=100.*chi(ir)/qsatr(ir)
            drh=rh-rhamb
         else
            rh=100.*chi(ir)/qsat
            drh=rh-rhamb
         endif
c ---    Place delta RH in its bin
         ibin=1
         do ib=1,26
            if(drh.GE.bindrh(ib)) ibin=ib
         enddo
         numdrh(ibin,ir)=numdrh(ibin,ir)+1

c ---    Compute change in temperature (K)
         if(chi(ir).LE.0.0) then
c ---       No plume impact
            dtk=0.0
         elseif(LCALPUFF) then            
c ---       Find the plume temperature associated with qsatr
            tk=tdegk
            if(qsatr(ir).GE.qsat) then
c ---          Plume is warmer than ambient
20             tk=tk+.1
               call VSAT(tk,psat,csat)
               if(csat.LE.qsatr(ir)) goto 20
            else
c ---          Plume is cooler than ambient
30             tk=tk-.1
               call VSAT(tk,psat,csat)
               if(csat.GE.qsatr(ir)) goto 30
            endif
            dtk=tk-tdegk
         else
            dtk=0.0
         endif
c ---    Place delta Temperature in its bin
c ---    Convert from K to F
         dtf=dtk*1.8
         ibin=1
         do ib=1,21
            if(dtf.GE.bindt(ib)) ibin=ib
         enddo
         numdt(ibin,ir)=numdt(ibin,ir)+1

c ---    Compute HUMIDEX (C)
         if(chi(ir).LE.0.0) then
c ---       No plume impact
            hdex=hamb
         elseif(LCALPUFF) then            
c ---       Get dew-point temperature to nearest 0.1 K
            tdewk=tk
            if(chi(ir).LT.qsatr(ir)) then
40             tdewk=tdewk-.1
               call VSAT(tdewk,psat,csat)
               if(csat.GE.chi(ir)) goto 40
            endif
            eterm=6.11*EXP(5417.7530*(1./273.16-1./tdewk))
            hdex=(tk-273.16)+(eterm-10.)/1.8
c ---       Convert from C to F
            hdex=hdex*1.8+32.
         else
            hdex=hamb
         endif
c ---    Place HUMIDEX in its bin
         ibin=1
         do ib=1,20
            if(hdex.GE.binhdx(ib)) ibin=ib
         enddo
         numhdx(ibin,ir)=numhdx(ibin,ir)+1

      enddo

      return
      end
c----------------------------------------------------------------------
      subroutine vsat(tdegk,psat,csat)
c----------------------------------------------------------------------
c
c --- CTEMISS    Version: 5.0      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
