c----------------------------------------------------------------------
c --- SUMPOST --  Summary processor for plume-mode FOG model output
c----------------------------------------------------------------------
c
c --- SUMPOST    Version: 2.1      Level: 060215                   MAIN
c ---            J. Scire, D. Strimaitis    Earth Tech, Inc.
c
c --- Written by:  Joe Scire, David Strimaitis
c                  Earth Tech. Inc
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 --- PURPOSE    Summarize statistics from a POSTPM list file,
c                with special treatment for ambient fogging
c                (Version of SUMSTAT modified by D. Strimaitis)
c
c
c --- Update 02/15/06: Add command-line input of POSTPM list file name
c
c --- Update 11/28/01: Begin tabulations at first plume-length provided 
c                      in the control file (height is set to ZERO
c                      if computed length is less than or equal to
c                      the first distance)
c --- Update 11/28/01: Add banner to list file for version/level
c --- Update 9/7/01:   REMOVE merged plume height calculation
c --- Update 9/7/01:   Drop output of length & ht by flow direction
c --- Update 9/7/01:   Add report for daylight-only hours
c                      using sunrise/sunset calculation
c                      Requires new variables:
c                      IDAYLITE Flag for daylight module (0:off, 1:on)
c                      ITZONE	Time zone (e.g. 5 for EST)
c                      XLAT	N. Latitude (decimal degrees)
c                      XLON	W. Longitude (decimal degrees)
c                      ITWILITE Number of hours to add for twilight
c --- Update 10/20/00: read weather descriptor rather than code, BdF
c                      remove only rain, fog and snow
c --- Update 8/14/00:  Break out results by flow direction, reported 
c                      in 30-degree blocks (simulation period),and
c                      test only the ambient RH listed in POSTPM
c                      record when screening for ambient fog.
c --- Update 3/09/00:  Break out results by season as well as for the
c ---                  entire simulation period
c
c----------------------------------------------------------------------

c --- Set parameters
      parameter(mxhtbin=20,mxlgbin=20)
      parameter(mxhtp1=mxhtbin+1,mxlgp1=mxlgbin+1)
      parameter(mxhtbins=5*mxhtp1,mxlgbins=5*mxlgp1)
      parameter(io5=5,io6=6,io10=10)
c --- Add flow vector bins in 30-degree increments
      parameter(mxfvp1=13)

c --- Number of passes through data with different selection criteria
      parameter(mxpass=7)
c
      character*234 crec
      character*70 procfile
      character*30 seasons(5),flowdir(mxfvp1)
      character*3 fvdir(mxfvp1)
      character*16 cwx,fog,rain,snow
      character*12 verson,level
      logical strfind

      real thrshht(mxhtbin),thrshlg(mxlgbin)
      integer nhtbin(mxhtp1,5),nlgbin(mxlgp1,5)
      integer nhthr(24,mxhtp1,5),nlghr(24,mxlgp1,5)
      integer nhrs(5)
c --- Note:  elements 1-4 are winter/spring/summer/fall
c ---        element 5 is full period of simulation

c --- Array for joint occurence of heights, lengths by flow direction
      integer nhtlgfv(mxhtp1,mxlgp1,mxfvp1)
      integer nlgfv(mxlgp1,mxfvp1)
      integer nfv(mxfvp1)

c --- Arrays for daylight processing
      integer irise(366),iset(366),kk(6)
      logical lday
c
c --- Key flow-vector tabulation output to logical
      logical lflowvec
c
      data nhtbin/mxhtbins*0/,nlgbin/mxlgbins*0/
      data nhrs/5*0/
      data irhfog/98/
      data ipass/0/
      data seasons/'WINTER Season ----------------',
     &             'SPRING Season ----------------',
     &             'SUMMER Season ----------------',
     &             'FALL   Season ----------------',
     &             'SIMULATION PERIOD ------------'/
      data flowdir/' 15 - 45 deg  ----------------',
     &             ' 45 - 75 deg  ----------------',
     &             ' 75 - 105 deg (East) ---------',
     &             '105 - 135 deg  ---------------',
     &             '135 - 165 deg  ---------------',
     &             '165 - 195 deg (South) --------',
     &             '195 - 225 deg  ---------------',
     &             '225 - 255 deg  ---------------',
     &             '255 - 285 deg (West) ---------',
     &             '285 - 315 deg  ---------------',
     &             '315 - 345 deg  ---------------',
     &             '345 - 15 deg (North) ---------',
     &             'All Directions ---------------'/
      data fvdir/'030','060','090','120','150','180',
     &           '210','240','270','300','330','360','All'/
c weather events to be removed: (see strfind below)
      data Fog /'Fog             '/
      data Rain/'Rain            '/
      data Snow/'Snow            '/
c
c --- Daylight option
      data idaylite/0/,itzone/0/,xlat/0.0/,xlon/0.0/,itwilite/0/
c
c --- Set flow-vector tabulation output logical
      data lflowvec/.FALSE./
c --- Set default filename for POSTPM list file
      data procfile/'postpm.lst'/

c --- Extract name of POSTPM output list file from command line
      call COMLINE(procfile)

c
c --- Open files
      open(io5,file='sumpost.inp',status='old')
      open(io6,file='sumpost.lst',status='unknown')
      open(io10,file=procfile,status='old')
c
c --- set version/level numbers
      verson='2.1'
      level='060215'
      write(6,1200)verson,level
1200  format(//11x,'SUMPOST -- Plume-mode Summary',10x,
     1 'Version: ',a12,4x,'Level: ',a12/1x,13('**********')/)

c --- Report file being processed
      write(6,*)
      write(6,'(a18,a70)')'Processing File:  ',procfile
      write(6,*)

c
c --- Read control file
      read(io5,*)nht
      if(nht.gt.mxhtbin)then
         write(io6,*)'Error -- too many height bins specified ',
     1   'NHT = ',nht,' MXHTBIN = ',mxhtbin
         write(*,*)'Error -- too many height bins specified ',
     1   'NHT = ',nht,' MXHTBIN = ',mxhtbin
         stop
      endif
      read(io5,*)(thrshht(n),n=1,nht)
c
      read(io5,*)nlg
      if(nlg.gt.mxlgbin)then
         write(io6,*)'Error -- too many plume length bins specified ',
     1   'NLG = ',nlg,' MXLGBIN = ',mxlgbin
         write(*,*)'Error -- too many plume length bins specified ',
     1   'NLG = ',nlg,' MXLGBIN = ',mxlgbin
         stop
      endif
      read(io5,*)(thrshlg(n),n=1,nlg)

c --- Section for defining MERGED PLUME parameters
c --- Set Defaults
      nsrc=1
      dsep=0.
      hsrc=0.
c --- Disable merged plume ht option here!
c --- read(io5,*,end=2) nsrc,dsep,hsrc
2     cbrtn=(FLOAT(nsrc))**(1./3.)
      if(nsrc.GT.1) then
         write(io6,*)
         write(io6,*)
         write(io6,*)'-------------------------------------------'
         write(io6,*)'  Plume Rise Enhanced for Merging Plumes'
         write(io6,*)'-------------------------------------------'
         write(io6,*)'  Number of Plumes      : ',nsrc
         write(io6,*)'  Separation of Sources : ',dsep
         write(io6,*)'  Height of Sources     : ',hsrc
         write(io6,*)'-------------------------------------------'
         write(io6,*)
      endif

c --- Add section for daylight processing information
      read(io5,*,end=3) idaylite
      if(idaylite.EQ.1) then
         read(io5,*) itzone
         read(io5,*) xlat
         read(io5,*) xlon
         read(io5,*) itwilite
         write(io6,*)
         write(io6,*)
         write(io6,*)'-------------------------------------------'
         write(io6,*)'  Daylight Module Inputs'
         write(io6,*)'-------------------------------------------'
         write(io6,*)'  Time Zone             : ',itzone
         write(io6,*)'  N. Latitude (deg)     : ',xlat
         write(io6,*)'  W. Longitude (deg)    : ',xlon
         write(io6,*)'  Twilight (hours)      : ',itwilite
         write(io6,*)'-------------------------------------------'
         write(io6,*)
c ---    Set up sunrise/sunset arrays if requested, using mid-pt of hr
c ---    as criterion for calling hour day or night
         do kday=1,366
            call SUN(kday,xlat,xlon,itzone,riselst,setlst)
            irise(kday)=NINT(riselst)+1-itwilite
            iset(kday)=NINT(setlst)+itwilite
         enddo
c ---    Write last result to list file as check
         write(io6,*)'Example sunrise/sunset ---'
         write(io6,*)'Day 366:  rise(LST),irise = ',riselst,irise(366)
         write(io6,*)'            set(LST),iset = ',setlst,iset(366)
c ---    Report hours to list file
         write(io6,*)
         write(io6,1205)
         do i=1,61
            kk(1)=i
            do j=2,6
               kk(j)=kk(j-1)+61
            enddo
            write(io6,1206) (kk(j),irise(kk(j)),iset(kk(j)),j=1,6)
         enddo
      endif
1205  format(1x,'Computed sunrise/sunset hours for site location',//,
     &       1x,6('Jday Rise Set '),/)
1206  format(1x,6(i4,1x,i4,1x,i3,1x))
3     continue

c
c --- Define intermediate parameters
      nhtp1=nht+1
      nlgp1=nlg+1

c --- Initialize arrays
      do k=1,5
         nhrs(k)=0
         do i=1,mxhtp1
            nhtbin(i,k)=0
            do j=1,24
               nhthr(j,i,k)=0
            enddo
         enddo
         do i=1,mxlgp1
            nlgbin(i,k)=0
            do j=1,24
               nlghr(j,i,k)=0
            enddo
         enddo
      enddo
      do k=1,mxfvp1
         do j=1,mxlgp1
            do i=1,mxhtp1
               nhtlgfv(i,j,k)=0
            enddo
         enddo
      enddo
      do k=1,mxfvp1
         nfv(k)=0
      enddo

c --- Get run identifier from POST-file
4     continue
      read(io10,'(a234)',end=999)crec
      if(crec(1:6).NE.'title:') goto 4
      write(io6,'(a234)')crec

c --- Set top of loop for making additional passes through data
5     ipass=ipass+1
      nread=0

c --- Set constraints and write output identifier for this pass
      if(ipass.EQ.1) then
         irh0=100
         vpl0=-1000.
         isky0=10
         iprec0=99999
         ifog0=1
         write(io6,*)
         write(io6,*)
         write(io6,*)'----------------------------'
         write(io6,*)'  Tabulation for all hours'
         write(io6,*)'----------------------------'
         write(io6,*)
      elseif(ipass.EQ.2) then
         irh0=irhfog
         vpl0=-1000.
         isky0=10
         iprec0=99999
         ifog0=1
         write(io6,*)
         write(io6,*)
         write(io6,*)'----------------------------------------------'
         write(io6,'(''  Ambient FOG hours removed (RH > '',i2,''%)'')')
     &             irh0
         write(io6,*)'----------------------------------------------'
         write(io6,*)
      elseif(ipass.EQ.3) then
         irh0=100
         vpl0=-1.
         isky0=10
         iprec0=99999
         ifog0=1
         write(io6,*)
         write(io6,*)
         write(io6,*)'----------------------'
         write(io6,*)'  CALM hours removed'
         write(io6,*)'----------------------'
         write(io6,*)
      elseif(ipass.EQ.4) then
         irh0=irhfog
         vpl0=-1.
         isky0=10
         iprec0=99999
         ifog0=1
         write(io6,*)
         write(io6,*)
         write(io6,*)'--------------------------------------'
         write(io6,*)'  CALM and Ambient FOG hours removed'
         write(io6,*)'--------------------------------------'
         write(io6,*)
      elseif(ipass.EQ.5) then
         irh0=irhfog
         vpl0=-1.
         isky0=9
         iprec0=99999
         ifog0=1
         write(io6,*)
         write(io6,*)
         write(io6,*)'----------------------------------------------'
         write(io6,*)'  CALM/Ambient FOG/Overcast(10) hours removed'
         write(io6,*)'----------------------------------------------'
         write(io6,*)
      elseif(ipass.EQ.6) then
         irh0=irhfog
         vpl0=-1.
         isky0=10
         iprec0=0
         ifog0=4
         write(io6,*)
         write(io6,*)
         write(io6,*)'------------------------------------------'
         write(io6,*)'  CALM/Ambient FOG/Weather hours removed'
         write(io6,*)'------------------------------------------'
         write(io6,*)
      elseif(ipass.EQ.7) then
         irh0=irhfog
         vpl0=-1.
         isky0=9
         iprec0=0
         ifog0=4
         write(io6,*)
         write(io6,*)
         write(io6,*)'-------------------------------------------------'
         write(io6,*)'  CALM/Ambient FOG/Overcast/Weather hours removed'
         write(io6,*)'-------------------------------------------------'
         write(io6,*)
      else
         stop 'Invalid IPASS'
      endif      
c
c --- Skip to a data record in POSTPM list file
10    continue
      read(io10,'(a234)',end=999)crec
c
      if(crec(1:4).ne.'Day=')go to 10
      nread=nread+1
c
c --- Read the day (1-366)
      read(crec(6:8),'(i3)') iday

c --- Read the hour of the day (1-24)
      read(crec(16:17),'(i2)') ihr

c --- Read the plume height (m)
      read(crec(111:116),14)h
14    format(f6.1)

c --- Read the visible plume length (m)
      read(crec(135:140),16)vpl
16    format(f6.3)

c --- Read ambient Relative Humidity (%) at plume centerline
      read(crec(25:27),'(i3)') irh

c --- Read ambient Relative Humidity (%) at surface
      read(crec(165:167),'(i3)') irhss

c --- Read Flow Vector index (1-36)
      read(crec(176:177),'(i2)') ifv36

c --- Read the opaque sky cover code
      read(crec(187:188),'(i2)') isky

c --- Read the weather description (NB: 16chars (was 20 in postpm2))
      read(crec(209:226),'(a16)') cwx

c iprec and ifog relate to CD144 only. read weather descriptor instead
cc --- Read the precipitation codes
c      read(crec(201:205),'(i5)') iprec
c
cc --- Read the fogging code
c      read(crec(206:206),'(i1)') ifog
c
c --- Apply constraints here
      if(irhss.GT.irh0) goto 10
      if(vpl.LT.vpl0) goto 10

      if (ipass.eq.5.or.ipass.eq.7) then
c throw out overcast:
         if (isky.eq.10) goto 10
      endif
      if (ipass.eq.6.or.ipass.eq.7) then
c remove fog, rain and snow events.
c to add to the list, make new character*16 variables with the event to remove
c the 3rd parameter when calling strfind is the no. characters in that event
         if (strfind(cwx,Fog,3)) goto 10
         if (strfind(cwx,Rain,4)) goto 10
         if (strfind(cwx,Snow,4)) goto 10
c         if (cwx.ne.'                ') then
c            write(*,*) 'Weather event discarded: ',cwx
c         endif
      endif
c use ipass if loops above instead
c      if(iprec.GT.iprec0) goto 10
c      if(ifog.GT.0 .AND. ifog.LT.ifog0) goto 10
c      if(isky.gt.isky0) goto 10

c --- Compute revised plume height if merged plume data are provided
      if(nsrc.GT.1 .AND. dsep.GT.0.) then
         rise0=h-hsrc
         ss=6.*(FLOAT(nsrc-1)*dsep/(cbrtn*rise0))**1.5
         h=hsrc+rise0*((FLOAT(nsrc)+ss)/(1.+ss))**(1./3.)
      endif

c --- Identify the 30-degree flow direcion bin for this record
      iflo=NINT(FLOAT(ifv36)/3.)
      if(iflo.EQ.0) iflo=12

c --- Identify the season for this record
      if(iday.GE.336) then
         iseas=1
      elseif(iday.GE.245) then
         iseas=4
      elseif(iday.GE.153) then
         iseas=3
      elseif(iday.GE.61) then
         iseas=2
      else
         iseas=1
      endif

      lday=.FALSE.
      if(idaylite.EQ.1) then
c ---    Identify day/night for this record
         if(ihr.GE.irise(iday) .AND. ihr.LE.iset(iday)) lday=.TRUE.
      endif

c --- Screen out (set to zero) any plume height entry associated with
c --- plume lengths that are <= the first plume length tabulated
      if(vpl.LE.thrshlg(1)) h=0.0

c --- Put this record into its appropriate flow direction bin
      nfv(iflo)=nfv(iflo)+1
      nfv(mxfvp1)=nfv(mxfvp1)+1

c --- Put this record into its appropriate plume height bins
      do i=1,nht
         if(h.le.thrshht(i))then
c ---       Add screen for daylight hours
            if(LDAY .OR. idaylite.EQ.0) then
               nhtbin(i,5)=nhtbin(i,5)+1
               nhthr(ihr,i,5)=nhthr(ihr,i,5)+1
               nhtbin(i,iseas)=nhtbin(i,iseas)+1
               nhthr(ihr,i,iseas)=nhthr(ihr,i,iseas)+1
            endif
c ---       Add joint tabulation with length and direcion
            do j=1,nlg
               if(vpl.le.thrshlg(j))then
                  nhtlgfv(i,j,iflo)=nhtlgfv(i,j,iflo)+1
                  nhtlgfv(i,j,mxfvp1)=nhtlgfv(i,j,mxfvp1)+1
                  nlgfv(j,iflo)=nlgfv(j,iflo)+1
                  nlgfv(j,mxfvp1)=nlgfv(j,mxfvp1)+1
                  go to 96
               endif
            enddo
            nhtlgfv(i,nlgp1,iflo)=nhtlgfv(i,nlgp1,iflo)+1
            nhtlgfv(i,nlgp1,mxfvp1)=nhtlgfv(i,nlgp1,mxfvp1)+1
            nlgfv(nlgp1,iflo)=nlgfv(nlgp1,iflo)+1
            nlgfv(nlgp1,mxfvp1)=nlgfv(nlgp1,mxfvp1)+1
            go to 96
         endif
      enddo
c --- Add screen for daylight hours
      if(LDAY .OR. idaylite.EQ.0) then
         nhtbin(nhtp1,5)=nhtbin(nhtp1,5)+1
         nhthr(ihr,nlgp1,5)=nhthr(ihr,nhtp1,5)+1
         nhtbin(nhtp1,iseas)=nhtbin(nhtp1,iseas)+1
         nhthr(ihr,nlgp1,iseas)=nhthr(ihr,nhtp1,iseas)+1
      endif
c --- Add joint tabulation with length and direcion
      do j=1,nlg
         if(vpl.le.thrshlg(j))then
            nhtlgfv(nhtp1,j,iflo)=nhtlgfv(nhtp1,j,iflo)+1
            nhtlgfv(nhtp1,j,mxfvp1)=nhtlgfv(nhtp1,j,mxfvp1)+1
            nlgfv(j,iflo)=nlgfv(j,iflo)+1
            nlgfv(j,mxfvp1)=nlgfv(j,mxfvp1)+1
            go to 96
         endif
      enddo
      nhtlgfv(nhtp1,nlgp1,iflo)=nhtlgfv(nhtp1,nlgp1,iflo)+1
      nhtlgfv(nhtp1,nlgp1,mxfvp1)=nhtlgfv(nhtp1,nlgp1,mxfvp1)+1
      nlgfv(nlgp1,iflo)=nlgfv(nlgp1,iflo)+1
      nlgfv(nlgp1,mxfvp1)=nlgfv(nlgp1,mxfvp1)+1
96    continue
c
c --- Put this record into its appropriate plume length bin
      do i=1,nlg
         if(vpl.le.thrshlg(i))then
c ---       Add screen for daylight hours
            if(LDAY .OR. idaylite.EQ.0) then
               nlgbin(i,5)=nlgbin(i,5)+1
               nlghr(ihr,i,5)=nlghr(ihr,i,5)+1
               nlgbin(i,iseas)=nlgbin(i,iseas)+1
               nlghr(ihr,i,iseas)=nlghr(ihr,i,iseas)+1
            endif
            go to 196
         endif
      enddo
c --- Add tabulation for daylight hours
      if(LDAY .OR. idaylite.EQ.0) then
         nlgbin(nlgp1,5)=nlgbin(nlgp1,5)+1
         nlghr(ihr,nlgp1,5)=nlghr(ihr,nlgp1,5)+1
         nlgbin(nlgp1,iseas)=nlgbin(nlgp1,iseas)+1
         nlghr(ihr,nlgp1,iseas)=nlghr(ihr,nlgp1,iseas)+1
      endif
196   continue
c
c --- Repeat and read next record
      go to 10
999   continue
c

c --- Set base for fraction
      fract=1./FLOAT(nread)

c --- Write summary tables
      do k=1,5
         nhrs(k)=0
         do i=1,nhtp1
            nhrs(k)=nhrs(k)+nhtbin(i,k)
         enddo
      enddo
c --- QA Check
      nhr1=nhrs(5)
      nhr2=0
      do i=1,nlgp1
         nhr2=nhr2+nlgbin(i,5)
      enddo
c
c --- nhr1 and nhr2 should be the same
      if(nhr1.ne.nhr2)then
         write(io6,*)'Error - NHR1 is not equal to NHR2 -- NHR1 = ',
     1   nhr1,' NHR2 = ',nhr2
         write(*,*)'Error - NHR1 is not equal to NHR2 -- NHR1 = ',
     1   nhr1,' NHR2 = ',nhr2
         stop
      endif

c --- Report hours read and included
      write(io6,*)
      write(io6,*)'Total Number of Hours: ',nread
      write(io6,*)'      Hours in Tables: ',nhrs(5)
      write(io6,*)'              ( Winter ',nhrs(1),')'
      write(io6,*)'              ( Spring ',nhrs(2),')'
      write(io6,*)'              ( Summer ',nhrs(3),')'
      write(io6,*)'              ( Fall   ',nhrs(4),')'


c --- Loop over seasons (#5 is full simulation period)
      do k=1,5
c
c --- Summary table for plume heights
      write(io6,*)
      if(idaylite.EQ.1) then
         write(io6,*) seasons(k),'  (DAYLIGHT)'
      else
         write(io6,*) seasons(k)
      endif
      write(io6,'(''PLUME HEIGHT Summary Table'',t68,
     &            ''Diurnal Distribution (hours 01-24)'')')
      write(io6,'(32x,''Hours        %     Cum. %'',8x,24(2x,i2.2))')
     &           (n,n=1,24)
      cum=0.0
      do i=1,nht
         pcnt=100.*float(nhtbin(i,k))*fract
         cum=cum+pcnt
         write(io6,201)thrshht(i),nhtbin(i,k),pcnt,cum,
     &                 (nhthr(j,i,k),j=1,24)
      enddo
201   format(1x,' Plume height <= ',f7.1,' (m)   ',i5,2f10.1,8x,24i4)
      pcnt=100.*float(nhtbin(nhtp1,k))*fract
      cum=cum+pcnt
      write(io6,202)thrshht(nht),nhtbin(nhtp1,k),pcnt,cum,
     &              (nhthr(j,nhtp1,k),j=1,24)
202   format(1x,' Plume height  > ',f7.1,' (m)   ',i5,2f10.1,8x,24i4)
c
c --- Summary table for visible plume lengths
      write(io6,*)
      if(idaylite.EQ.1) then
         write(io6,*) seasons(k),'  (DAYLIGHT)'
      else
         write(io6,*) seasons(k)
      endif
      write(io6,'(''VISIBLE PLUME LENGTH Summary Table'',t68,
     &            ''Diurnal Distribution (hours 01-24)'')')
      write(io6,'(40x,''Hours        %     Cum. %'',24(2x,i2.2))')
     &           (n,n=1,24)
      cum=0.0
      do i=1,nlg
         pcnt=100.*float(nlgbin(i,k))*fract
         cum=cum+pcnt
         write(io6,301)thrshlg(i),nlgbin(i,k),pcnt,cum,
     &                 (nlghr(j,i,k),j=1,24)
      enddo
301   format(1x,' Visible plume length <= ',f7.1,' (m)   ',i5,2f10.1,
     &          24i4)
      pcnt=100.*float(nlgbin(nlgp1,k))*fract
      cum=cum+pcnt
      write(io6,302)thrshlg(nlg),nlgbin(nlgp1,k),pcnt,cum,
     &              (nlghr(j,nlgp1,k),j=1,24)
302   format(1x,' Visible plume length  > ',f7.1,' (m)   ',i5,2f10.1,
     &       24i4)

c --- End loop over seasons
      enddo

c
c --- Summary tables by flow direction
c-------------------------------------

c --- Skip this output section for most applications!
      if(LFLOWVEC) then     

c ---    Tabulate lengths by flow direction (all heights)
         write(io6,*)
         write(io6,*)
         write(io6,*) seasons(5)
         write(io6,'(''VISIBLE PLUME LENGTH Summary Table'',t68,
     &            ''Flow Direction (30-deg sectors centered @)'')')
         write(io6,'(40x,''Hours        %     Cum. %'',12(2x,a3))')
     &              (fvdir(n),n=1,12)
         cum=0.0
         do i=1,nlg
            pcnt=100.*float(nlgfv(i,mxfvp1))*fract
            cum=cum+pcnt
            write(io6,401)thrshlg(i),nlgfv(i,13),pcnt,cum,
     &                    (nlgfv(i,k),k=1,12)
         enddo
401   format(1x,' Visible plume length <= ',f7.1,' (m)   ',i5,2f10.1,
     &       12i5)

         pcnt=100.*float(nlgfv(nlgp1,mxfvp1))*fract
         cum=cum+pcnt
         write(io6,402)thrshlg(nlg),nlgfv(nlgp1,13),pcnt,cum,
     &              (nlgfv(nlgp1,k),k=1,12)
402   format(1x,' Visible plume length  > ',f7.1,' (m)   ',i5,2f10.1,
     &       12i5)

         pcnt=100.*float(nfv(mxfvp1))*fract
         cum=pcnt
         write(io6,403) nfv(mxfvp1),pcnt,cum,
     &              (nfv(k),k=1,12)
403   format(1x,' ALL Visible plume lengths ',12x,i5,2f10.1,
     &       12i5)

c ---    Tabulate height by flow direction for each length interval
c ---    Loop over visible lengths
         do j=1,nlgp1
            write(io6,*)
            write(io6,*)
            write(io6,*) seasons(5)
            if(j.EQ.1) then
               write(io6,'(''Visible PLUME LENGTH (m) <= '',f7.1)')
     &                 thrshlg(1)
            elseif(j.EQ.nlgp1) then
               write(io6,'(''Visible PLUME LENGTH (m) >  '',f7.1)')
     &                 thrshlg(j-1)
            else
             write(io6,
     &            '(f7.1,'' < Visible PLUME LENGTH (m) <='',f7.1)')
     &             thrshlg(j-1),thrshlg(j)
            endif
            write(io6,*)
            write(io6,'(''PLUME HEIGHT Summary Table'',t68,
     &                ''Flow Direction (30-deg sectors centered @)'')')
            write(io6,'(32x,''Hours        %     Cum. %'',
     &                  8x,12(2x,a3))') (fvdir(n),n=1,12)
            cum=0.0
            do i=1,nht
               pcnt=100.*float(nhtlgfv(i,j,mxfvp1))*fract
               cum=cum+pcnt
               write(io6,411)thrshht(i),nhtlgfv(i,j,13),pcnt,cum,
     &                       (nhtlgfv(i,j,k),k=1,12)
            enddo
411      format(1x,' Plume height <= ',f7.1,' (m)   ',i5,2f10.1,
     &          8x,12i5)

            pcnt=100.*float(nhtlgfv(nhtp1,j,mxfvp1))*fract
            cum=cum+pcnt
            write(io6,412)thrshht(nht),nhtlgfv(nhtp1,j,13),pcnt,cum,
     &                 (nhtlgfv(nhtp1,j,k),k=1,12)
         enddo
412   format(1x,' Plume height  > ',f7.1,' (m)   ',i5,2f10.1,
     &       8x,12i5)

      endif

c --- Check for next pass
      if(ipass.LT.mxpass) then
c ---    Clear arrays
         do k=1,5
            do i=1,mxhtp1
               nhtbin(i,k)=0
               do j=1,24
                  nhthr(j,i,k)=0
               enddo
            enddo
            do i=1,mxlgp1
               nlgbin(i,k)=0
               do j=1,24
                  nlghr(j,i,k)=0
               enddo
            enddo
         enddo
         do k=1,mxfvp1
            do j=1,mxlgp1
               nlgfv(j,k)=0
               do i=1,mxhtp1
                  nhtlgfv(i,j,k)=0
               enddo
            enddo
         enddo
         do k=1,mxfvp1
            nfv(k)=0
         enddo
c ---    Restore file
         REWIND(io10)
c ---    Start over again
         goto 5
      endif

c
      stop
      end

c -----------------------------------------------------------------
      logical function strfind(str,find,l2)
c
c returns true if the first l2 characters of find are present
c as a separate word in str.
c Benjamin de Foy, 20th October 2000
c -----------------------------------------------------------------
c string to look in:
      character*16 str
c length of string to look in:     
      integer l1
c string to be found:
      character*16 find
c length of string to be find:
      integer l2

      strfind = .false.
      l1 = 16
      j = 1
      newword = 1
      do i = 1,l1-l2
         if (str(i:i+l2).eq.find(j:j+l2)) then
c check that we have a single word, not part of a word:
            if (newword.eq.1.and.
     &           (str(i+l2+1:i+l2+1).eq.' '.or.(i+l2).eq.l1)) then
               strfind = .true.
c               write(*,*) 'Found ',find,' in ',str
               goto 900
            endif
         endif
c are we beginning a new word?
         if (str(i:i).eq.' ') then
            newword = 1
         else
            newword = 0
         endif
      enddo
 900  continue
      return
      end
c ----------------------------------------------------------------------
      subroutine sun(julday,xlat,xlon,itzone,riselst,setlst)
c ----------------------------------------------------------------------
c --- Adopted From PCRAMMET SUBROUTINE SUN(JULDAY) ---
C***********************************************************************
C*    SUN Module of PCRAMMET Meteorological Pre-processor
C*
C*    PURPOSE: This routine calculates the sunrise and sunset times
C*             of the day from the date, latitude, longitude, and time
C*             zone
C*
C*    PROGRAMMER: Jayant Hardikar
C*                PES Inc.
C*                Adapted From MPRM
C*
C*    DATE:       February 23, 1995
C*
C*    INPUTS:     Latitude, longitude, Julian day
C*
C*    OUTPUTS:    Sunrise, sunset, solar elevation for each hour of day
C*
C*    CALLED FROM: MAIN program
C***********************************************************************
C*
C*    Variable Declarations

C*    ZONE 05 = EASTERN
C*    ZONE 06 = CENTRAL
C*    ZONE 07 = MOUNTAIN
C*    ZONE 08 = PACIFIC

      data const/57.29578/                                                                                                             MET40800
C*
      DUM    = XLON/15.0 - ITZONE
      TEMPZ  = 15.0*ITZONE - XLON
      SINLAT = SIN(XLAT/CONST )
      COSLAT = COS(XLAT/CONST )

C---- Determine the angular(radians) fraction of a year for this date;
C     the constant 0.0172028=360./365.242*57.29578
      DAYNO  = (JULDAY-1.0)*0.0172028
      TDAYNO = 2.0*DAYNO
      SIND   = SIN(DAYNO)
      COSD   = COS(DAYNO)
      SINTD  = SIN(TDAYNO)
      COSTD  = COS(TDAYNO)
C
C---- Account for ellipticity of earth's orbit.

      SIGMA = 279.9348 + (DAYNO*CONST) + 1.914827*SIND -
     &        0.079525*COSD + 0.019938*SINTD - 0.00162*COSTD

C---- Find the sine of the solar declination;
C     the constant 0.39785=sin(0.4091720193=23.44383/57.29578)

      DSIN = 0.39785*SIN(SIGMA/CONST)
      DCOS = SQRT(1.0-DSIN*DSIN)

C---- Determine time(hrs) of meridian passage

      AMM  = 12.0 + 0.12357*SIND - 0.004289*COSD + 0.153809*SINTD +
     &                                             0.060783*COSTD
      HCOS = (-SINLAT*DSIN)/(COSLAT*DCOS)

C---- Determine solar hour angle of sunrise-sunset.
      H2 = (ATAN2(SQRT(1.0-HCOS*HCOS),HCOS)/15.0)*CONST

C---- Time of sunrise(riselst) and time of sunset(setlst) are expressed
C     in local standard time since the zone correction has already
C     been made; otherwise they would be in Greenwich mean time.
      riselst = AMM - H2 + DUM
      setlst  = AMM + H2 + DUM

      RETURN
      END
c----------------------------------------------------------------------
      subroutine comline(ctext)
c----------------------------------------------------------------------
c
c --- CALUTILS   Version: 2.2      Level: 960521                COMLINE
c                J. Scire, SRC
c
c --- PURPOSE:  Call the compiler-specific system routine that will
c               pass back the command line argument after the text
c               that executed the program
c
c               This routine contains calls for several different
c               compilers, but only one should be active at any one
c               time.
c
c --- INPUTS:
c
c          CTEXT - character*70  - Default command line argument #1
c
c --- OUTPUT:
c
c          CTEXT - character*70  - Command line argument #1
c                                  If command line argument is
c                                  missing, CTEXT is not changed
c
c --- COMLINE called by:  SETUP
c --- COMLINE calls:      GETCL - Lahey system routine
c                         NARGS, GETARG - Microsoft routines
c                         NARGS, GETARG - Compaq DF routines
c                         IARGC, GETARG - Sun routines
c                         IARGC, GETARG - HP routines
c----------------------------------------------------------------------
c
c ----------------------------------------
c --- COMPAQ DF compiler library directive
c *** USE DFLIB
c ----------------------------------------
c
      character*70 ctext,cdeflt
c
c --- Microsoft variables
c *** integer*2 iarg,istat
c
c --- HP declaration
c *** external getarg
c *** external iargc
c
c --- The following is for any system without a command line routine
c --- and is also used as a default
      cdeflt=ctext
c
c ------------------
c --- Lahey compiler
c ------------------
      call getcl(ctext)
c
c ----------------------
c --- COMPAQ DF compiler
c ----------------------
c ***  numargs=NARGS()
c ***  if(numargs.ge.1)then
c ***     call GETARG(1,ctext)
c ***  endif
c
c ----------------------
c --- Microsoft compiler
c ----------------------
c *** numargs=NARGS()
c *** if(numargs.ge.1)then
c ***    iarg=1
c ***    call GETARG(iarg,ctext,istat)
c *** endif
c
c ----------------
c --- Sun compiler
c ----------------
c *** numargs=IARGC()
c *** if(numargs.ge.1)then
c ***    call GETARG(1,ctext)
c *** endif
c
c -----------------------------------------------------
c --- HP compiler -- NOTE: needs +U77 switch on compile
c -----------------------------------------------------
c *** numargs=IARGC()
c *** if(numargs.ge.1)then
c ***    call GETARG(1,ctext)
c *** endif
c
c
c --- If no command line arguments, use default
      if(ctext(1:1).eq.' ')ctext=cdeflt

      return
      end
