c----------------------------------------------------------------------
c --- POSTPM --  Postprocessor for plume-mode FOG model output
c----------------------------------------------------------------------
c
c --- POSTPM     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 Exponent, Inc.
c
c----------------------------------------------------------------------
c --- Changes from  980514 to 060214:
c
c --- 060214   (DGS): Change name of input FOG data file produced by
c                     CALPUFF from FOG.BIN to FOG.DAT to be consistent
c                     with default name in CALPUFF, and add to control
c                     file
c              (DGS): Add list-file name to control file
c --- 010905   (DGS): Change output statement to write VPL0 rather than
c                     VPL (VPL is from the 'last trail'; VPL0 is the
c                     maximum over all receptor trails)
c              (DGS): Set plume height to zero when the length is less
c                     than the distance to the first receptor
c              (DGS): Fix in-plume RH reported (ambient water content
c                     had been added twice, resulting in RH=100%)
c              (DGS): Change version and level to character variables
c              (DGS): Increment year being processed when Julian day
c                     decreases
c --- 001020   (BdF): read present weather from HUSWO,
c                     replace lmet by iwx in namelist(1:CD144,2:HUSWO)
c --- 000811   (DGS): Add flow vector to output records
c              (DGS): Add Y2K changes to met data processing
c              (DGS): Restore printing logic for selected days/hours
c              (DGS): Add CALPUFF v5.4 (000602_2) changes
c                      1 - Accept LCOMPRS logical in header record 2
c                      2 - Each source produces a receptor trail, so
c                          plume length is the maximum among these
c
c --- 991203F  (DGS): Add present weather procesing (CD144 metdat)
c                     as coded by BdF 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 --- 980918F  (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
      parameter(mxrec=200)


c --- Add array for qsat at each receptor
c --- real chi(mxrec),xr(mxrec)
      real chi(mxrec),xr(mxrec),qsatr(mxrec),qbakr(mxrec),tempr(mxrec)
      real hr(mxrec)

c --- Variables for CD144 data
      logical lmet
      integer iwx
      character*70 metdat, lstout, fogdat
      integer yy144,mm144,dd144,hr144,jd144
      character*1 opaque
      character*8 pwth
      character*20 weather
      character*81 line
c --- Variables for HUSWO data
      character*114 linehus
c string describing metdata format
      character*5 smet(2)

      integer ipcp(51),idays(366)
      logical lpmode,lprt,lvpl,lprt1
      logical lcalpuff,lcompr
      character*4 title(20)
      character*3 cfoglb

      character*12 verson,level
c
      data lprt/.false./,lvpl/.true./
      data nday1/1/,nday2/366/,nhr1/1/,nhr2/24/
      data nbkfog/0/,nbkice/0/,nbktot/0/
      data nbksfog/0/,nbksice/0/
      data dum/0.0/
      data lmet/.false./,metdat/'metdat.144'/
      data fogdat/'fog.dat'/, lstout/'postpm.lst'/
      data iwx/0/
      data smet/'CD144','HUSWO'/
      data pwth/'        '/
      data weather/'                    '/
c     
      namelist/opts/lprt,lvpl,nday1,nday2,nhr1,nhr2,
     &              iwx,metdat,lstout,fogdat

c --- Set 'previous day' to zero
      jdayold=0

c --- Process control file
      open(5,file='postpm.inp',status='old')
      read(5,opts)
c
c --- open files
      open(6,file=lstout,status='unknown',carriagecontrol=
     &            'fortran')
      open(10,file=fogdat,status='old',form='unformatted')
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,'POSTPM -- Plume-mode Postprocessor',10x,
     1 'Version: ',a12,4x,'Level: ',a12/1x,13('**********')/)
c
c --- write control file inputs
      write(6,opts)
c
c --- Open met data if needed - iwx must be 1 or 2
      if (iwx.eq.1.or.iwx.eq.2) then
         lmet = .true.
      else
         lmet = .false.
         iwx = 0
      endif
      if (lmet) then
         write(6,*)
         write(6,*) 'Opening ',metdat
         write(6,*) 'For ',smet(iwx),' Meteorological Data'
         write(6,*) 'Weather entry in table below contains '
     &        ,smet(iwx),' code followed by decoded description'
         write(6,*) 'Missing skycover = 99, Missing Weather = Blank'
         write(6,*)
         open(17,file=metdat)         
      endif
c
c --- read header #1 of FOG output file
      read(10)title,nnrec,npts,ipcp,iyr,idays,lpmode

c --- Set CALPUFF 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.

c --- check that output is plume-mode output
      if(.not.lpmode)then
         write(6,12)title,nnrec,npts,iyr,lpmode,ipcp,idays
12       format(//1x,'ERROR -- POSTPM can be used with PLUME-MODE ',
     1   'output only'//
     2   1x,'title:  ',20a4/
     3   1x,'nnrec:  ',i8/
     4   1x,'npts:   ',i8/
     5   1x,'iyr:    ',i8/
     6   1x,'lpmode: ',l8/
     7   1x,'ipcp:   ',51i2/
     8   1x,'idays:  ',10(10i1,1x)/
     8   1x,'        ',10(10i1,1x)/
     8   1x,'        ',10(10i1,1x)/
     8   1x,'        ',10(10i1,1x))
         stop
      endif

c --- Process header record #2 differently if from CALPUFF
      if(LCALPUFF) then
         nrec=nnrec
         read(10)(lcompr,(xr(n),n=1,nrec))
         if(LCOMPR) then
            write(*,*)'FATAL Error --- Data Compression in FOG data'
            write(*,*)'file is not supported!'
            stop
         endif
      else
c ---    decode compression-mode flag
         nrec=mod(nnrec,1000)
         read(10)(xr(n),n=1,nrec)
      endif

c --- Obtain current century for this run
      icent=iyr/100

c --- write header data
      write(6,16)title,nnrec,npts,iyr,lpmode,ipcp,idays
16    format(//1x,'HEADER RECORD DATA'//
     1         1x,'title:  ',20a4/
     2         1x,'nnrec:  ',i8/
     3         1x,'npts:   ',i8/
     4         1x,'iyr:    ',i8/
     5         1x,'lpmode: ',l8/
     6         1x,'ipcp:   ',51i2/
     7         1x,'idays:  ',10(10i1,1x)/
     8         1x,'        ',10(10i1,1x)/
     9         1x,'        ',10(10i1,1x)/
     a         1x,'        ',10(10i1,1x))
c
      write(6,18)(xr(n),n=1,nrec)
18    format(//1x,'Receptor distances (m): ',10f10.1/
     1 100(25x,10f10.1/))
c
c --- receptors must be in order of increasing downwind distance
      do 72 i=2,nrec
      if(xr(i).ge.xr(i-1))go to 72
c --- write error message & terminate
      write(6,*)'ERROR -- plume-mode receptors must be in order of ',
     1 'increasing downwind distance -- i = ',i,' xr(i) = ',xr(i),
     2 ' xr(i-1) = ',xr(i-1)
      stop
72    continue
c
c --- loop over days
      do 1000 iday=1,366
c --- if(idays(iday).ne.1)go to 1000
      if(idays(iday).ne.1 .AND. .not.LCALPUFF)go to 1000
c
c --- loop over hours
      do 500 ihr=1,24

c --- Set top of loop for processing all CALPUFF hours
1     continue

c
c --- read hourly date/time/met. data
      if(LCALPUFF) then
         read(10,end=9999)jday,jhr,ws,fv,istab,tdegk,dpbl,qback,qsat,
     &                    numrec
      else
         read(10,end=9999)jday,jhr,ws,fv,istab,tdegk,dpbl,qback,qsat
         numrec=nrec
      endif
c --- Set flow vector index (1-36)
      ifv36=NINT(fv/10.)
      if(ifv36.EQ.0) ifv36=36

c --- If JDAY is less than the previous day processed, then year must
c --- be incremented
      if(jday.LT.jdayold) then
         iyr=iyr+1
      endif
      jdayold=jday

c --- Get additional met data from corresponding metdata file
      idatim1=iyr*100000+jday*100+jhr   
      if (LMET) then
         if (iwx.eq.1) then
c read CD144 file:
31       read(17,'(a81)',end=932) line
         read(line,'(5x,4i2,10x,a8,47x,a1)') yy144,mm144,dd144,
     &                                       hr144,pwth,opaque
         call JULDAY(yy144,mm144,dd144,jd144)
         if (hr144.eq.0) then
c ---       Roll back 1 hour to find previous day's date
            nhrinc=-1
            call INCR(yy144,jd144,hr144,nhrinc)
c ---       Now add the hour back on to get hour 24
            hr144 = hr144+1
         endif

c ---    Convert 2-digit year to 4-digit year
         if(yy144.LT.1000 .AND. iyr.GT.1000) then
c ---       Adopt nearest century
            if(ABS(((icent-1)*100+yy144)-iyr).LT.100) then
               yy144=(icent-1)*100+yy144
            elseif(ABS((icent*100+yy144)-iyr).LT.100) then
               yy144=icent*100+yy144
            elseif(ABS(((icent+1)*100+yy144)-iyr).LT.100) then
               yy144=(icent+1)*100+yy144
            else
               write(*,*)'Problem with dates in met file'
               write(*,*)'Current processing year:  ',iyr
               write(*,*)'          Met file year:  ',yy144
               stop
            endif
         endif
         idatim2=yy144*100000+jd144*100+hr144
         if(idatim2.LT.idatim1) goto 31
         if(idatim2.GT.idatim1) then
            write(*,*) 'Current Date-Time NOT Found in CD144 file!'
            write(*,*) 'Current Date-Time = ',idatim1
            write(*,*) '  CD144 Date-Time = ',idatim2
            write(6,*) 'Current Date-Time NOT Found in CD144 file!'
            write(6,*) 'Current Date-Time = ',idatim1
            write(6,*) '  CD144 Date-Time = ',idatim2
            stop
         endif
         elseif (iwx.eq.2) then
c read HUSWO file:
41       read(17,'(a114)',end=932) linehus
         read(linehus,'(6x,i4,3i2,13x,i2,46x,a8)') yy144,mm144,dd144,
     &        hr144,isky,pwth
         call JULDAY(yy144,mm144,dd144,jd144)
c HUSWO hour is 1-24, year is 4 digit already
         idatim2=yy144*100000+jd144*100+hr144
         if(idatim2.LT.idatim1) goto 41
         if(idatim2.GT.idatim1) then
            write(*,*) 'Current Date-Time NOT Found in HUSWO file!'
            write(*,*) 'Current Date-Time = ',idatim1
            write(*,*) '  HUSWO Date-Time = ',idatim2
            write(6,*) 'Current Date-Time NOT Found in HUSWO file!'
            write(6,*) 'Current Date-Time = ',idatim1
            write(6,*) '  HUSWO Date-Time = ',idatim2
            stop
         endif
         endif
      endif
c
c --- compute rh (%) and mixing ratio deficit
      irh=100.*qback/qsat+0.5
      qdef=qsat-qback
c
c --- save ambient temperature and RH for output
      tsfc=tdegk
      irhsfc=irh
c
c --- set print flag
      if(jday.ge.nday1.and.jday.le.nday2.and.
     1 jhr.ge.nhr1.and.jhr.le.nhr2)then
         lprt1=.true.
      else
         lprt1=.false.
      endif
c
c --- write hourly detailed plume data
      if(lprt.and.lprt1)then
         write(6,102)jday,jhr,ws,fv,istab,tdegk,dpbl,irh,qback,
     1   qsat,qdef
         write(6,104)
102      format(/1x,'Day=',i3,' Hr=',i2,/1x,'SURFACE Data:',' Ws=',
     1   f6.2,' Flow Dir=',f5.1,' Stab=',i1,' Temp=',f5.1,' Mix. Ht=',
     2   f5.0,' RH=',i3,' Back. Q=',f6.3,' Sat. Q=',f6.3,
     3   ' Q Def.=',f6.3)
104      format(' Recep.',2x,'Distance',
     1   2x,'---------- Sy -----------',
     1   4x,'---------- Sz -----------'/19x,'Bldg  Amb.   Buoy.  Total',
     2   4x,'Bldg  Amb.   Buoy.  Total',3x,'Plume',5x,'Centerline',4x,
     3   'Centerline'/19x,'Down  Turb.  Ind.',12x,'Down  Turb.  Ind.',
     4   12x,'Ht.',6x,'Plm. Conc.',4x,'Ambient  Sat.  Liquid  T-Plume')
      endif

c --- Process weather data
      if(LMET) then 
         if (iwx.eq.1) then
            call PROC144(pwth,opaque,weather,isky)
         elseif (iwx.eq.2) then
            call PROCHUS(pwth,weather)
c reset missing weather to spaces:
            do ii = 1,7,2
               if (pwth(ii:ii+1).eq.'99') then
                  pwth(ii:ii+1) = '  '
               endif
            enddo
         endif
      endif
c
c --- loop over receptors

c --- CALPUFF writes multiple receptor trails (one from each cell) so
c --- set up distinct outer loop point here
      jtrail=0
      vpl0=-1000
240   continue
      jtrail=jtrail+1

c --- loop over receptors
      do 250 irec=1,nrec
c
c --- read receptor data

c --- CALPUFF record contains qsat and QBACK for plume temperature
      if(LCALPUFF) then
         krec=(jtrail-1)*nrec+irec
         read(10)jrec,h,conc,qsat,qback,tdegk
      else
         krec=irec
         read(10)jrec,sy0,syturb,sybuoy,sigy,sz0,szturb,szbuoy,sigz,
     &           h,conc
      endif
c
      if(jrec.ne.krec)then
         write(6,*)'ERROR -- receptor numbers do not match -- '
         write(6,*)'         Expected: ',krec
         write(6,*)'            Found: ',jrec
         stop
      endif
c
      qsatr(irec)=qsat
      qbakr(irec)=qback
      tempr(irec)=tdegk
      hr(irec)=h
      qtot=qback+conc
      chi(irec)=qtot
      cliq=amax1(qtot-qsat,0.0)
c
c --- label plume-induced fog/ice hours
      if(qtot.gt.qsat)then
         if(tdegk.lt.273.16)then
            cfoglb='ICE'
         else
            cfoglb='FOG'
         endif
      else
         cfoglb='   '
      endif
c
c --- write detailed plume information
      
      if(lprt.and.lprt1) then
         if(LCALPUFF) then
            write(6,106)krec,xr(irec),dum,dum,dum,dum,dum,
     1                  dum,dum,dum,h,conc,qback,qsat,cliq,
     2                  tdegk,cfoglb
         else
            write(6,106)irec,xr(irec),sy0,syturb,sybuoy,sigy,sz0,
     1                  szturb,szbuoy,sigz,h,conc,qback,qsat,cliq,
     2                  tdegk,cfoglb
         endif
      endif
106   format(i4,5x,f6.0,2x,4(f6.1,1x),1x,4(f6.1,1x),2x,f6.1,3x,
     1 1pe10.3,5x,0p3f7.3,f7.1,3x,a3)
250   continue
c
c --- compute and report visible plume length
      if(.not.lvpl) go to 400
c
c --- Find the receptor nearest the source where conc<qsat
      isavp1=nrec+1
      if(ws.GT.0.1) then
         do i=nrec,1,-1
            if(chi(i).LT.qsatr(i)) isavp1=i
         enddo
      else
         isavp1=0
      endif
      isav=isavp1-1
c
c --- Find the receptor farthest from the source where conc>qsat
      jsav=0
      if(ws.GT.0.1) then
         do j=1,nrec
            if(chi(j).GE.qsatr(j)) jsav=j
         enddo
      else
         jsav=-1
      endif
      jsavp1=jsav+1
c
c --- Use the far edge only when NO fog is found at first distance
      if(isav.NE.jsav) then
         if(isav.EQ.0) then
            isav=jsav
            isavp1=jsavp1
         endif
      endif
c
      if(isav.lt.0)then
c ---    visible plume length not defined in CALM
         vpl=-999.0
      elseif(isav.eq.0)then
c ---    visible plume length < xr(1)
         vpl=0.0
      elseif(isav.eq.nrec)then
c ---    visible plume length > xr(nrec)
         vpl=9999.
      else
c ---    visible plume length between xr(1) and xr(nrec) -- interpolate
         xess=chi(isav)-qsatr(isav)
         xessp1=chi(isavp1)-qsatr(isavp1)
         vpl=xr(isav)+(xr(isavp1)-xr(isav))*xess/(xess-xessp1)
      endif
      if(vpl.GT.vpl0) then
c ---    Update values reported to those just after plume dissipates
         i=MAX0(1,isavp1)
         i=MIN0(i,nrec)
         vqdef=qsatr(i)-qbakr(i)
         vqsat=qsatr(i)
         vqback=qbakr(i)
         vtdegk=tempr(i)
         if(vpl.GT.0.0) then
            vh=hr(i)
         else
            vh=0.0
         endif
         kk=(jtrail-1)*nrec+i

c ---    Compute in-plume relative humidity just after it dissipates
         irh=NINT(100.*chi(i)/vqsat)
         irh=MIN0(100,irh)
         vpl0=vpl
      endif
c
      if(krec.EQ.numrec) then
c ---    write visible plume length data
         if(lvpl.and.lprt1)write(6,306)jday,jhr,irh,vtdegk,vqback,vqsat,
     1                                 vqdef,vh,vpl0,tsfc,irhsfc,ifv36,
     2                                 isky,pwth,weather,kk
      endif
 306  format(1x,'Day= ',i3,'   Hr= ',i2,'   RH= ',i3,'   Temp= ',f5.1,
     1 '  Q_Backgrd.= ',f6.3,'   Q_Sat.=',f6.3,'   Q_Deficit=',f6.3,
     2 '  Plume_ht. = ',f6.1,'  Visible_Length= ',f6.0,
     3 '  Sfc_T= ',f5.1,'  Sfc_RH= ',i3,'  FV36= ',i2,
     4 '  Cloud= ',i2,'  Weather= ',a8,1x,a16,'  Rec= ',i3)

400   continue

c --- Return to receptor loop if another trail of receptors is in file
      if(krec.LT.numrec) goto 240

      nbktot=nbktot+1

      
c --- Keep track of no. background fog/ice hours at plume centerline
      if(irh.ge.100)then
         if(tdegk.lt.273.16)then
            nbkice=nbkice+1
         else
            nbkfog=nbkfog+1
         endif
      endif

c --- Keep track of no. background fog/ice hours at surface
      if(irhsfc.ge.100)then
         if(tsfc.lt.273.16)then
            nbksice=nbksice+1
         else
            nbksfog=nbksfog+1
         endif
      endif

      if(LCALPUFF) goto 1

500   continue
1000  continue
c
9999  continue
c
c --- write summary on background fog/ice hours at plume ht
      nbksum=nbkfog+nbkice
      write(6,1002)nbktot,nbkfog,nbkice,nbksum
1002  format(//1x,'Plume-Height:',/
     1 1x,'Total no. hours = ',i5/
     2 1x,'"Background fog" hours = ',i5/
     3 1x,'"Background ice" hours = ',i5/
     4 1x,'Total "background fog + ice" hours = ',i5)
c
c --- write summary on background fog/ice hours at surface
      nbksums=nbksfog+nbksice
      write(6,1003)nbktot,nbksfog,nbksice,nbksums
1003  format(//1x,'Surface:',/
     1 1x,'Total no. hours = ',i5/
     1 1x,'"Background fog" hours = ',i5/
     2 1x,'"Background ice" hours = ',i5/
     3 1x,'Total "background fog + ice" hours = ',i5//
     4 1x,'------------------------'//1x,'Notes:  ',
     5 '(1) "Background fog" hours --> background RH = 100%, ',
     6 'Temp > 32 deg. F'/9x,
     5 '(2) "Background ice" hours --> background RH = 100%, ',
     6 'Temp <= 32 deg. F'/)
c
      call runtime(6,dtsec)
c
      stop
c
932   write(*,*)'EOF in reading met data'
      write(*,*)'Searching for period:  ',idatim1
      write(*,*)'     Last met period:  ',idatim2
      stop

      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 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 proc144(pwth,opaque,weather,isky)
c -----------------------------------------------------------------
      character*1 opaque
      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

      isky=0
      if(opaque.EQ.'-') then
         isky=10
      elseif (opaque.eq.' ') then
         isky = 99
      else
         read(opaque,'(i1)') isky
      endif

      return
      end
c -----------------------------------------------------------------
      subroutine prochus(pwth,weather)
c
c     convert HUSWO weather code to text description for list file
c     Benjamin de Foy, 20th October 2000
c -----------------------------------------------------------------
      character*8 pwth
      character*20 weather
      integer ipwth(4)

      weather = '                    '

      read(pwth,'(4i2)') ipwth(1),ipwth(2),ipwth(3),ipwth(4)

c keep track of which events are reported to avoid duplication:
      istorm = 0
      irain = 0
      isnow = 0
      idust = 0
      isand = 0
      ifog = 0
      iblsnow = 0
      ihaze = 0
      iice = 0

      np = 1
      do i = 1,4
         iw = ipwth(i)
         if (iw.ne.99.and.iw.ge.10) then
            if (iw.lt.20) then
c thunderstorm, squall, spout, tornado = storm
               if (istorm.eq.0) then
                  weather(np:np+6) = 'Storm '
                  np = np+6
                  istorm = 1
               endif
            elseif (iw.lt.30) then
c rain, rain showers, freezing rain
               if (irain.eq.0) then
                  weather(np:np+5) = 'Rain '
                  np = np+5
                  irain = 1
               endif
            elseif (iw.lt.40) then
c Squall/drizzle = rain
               if (irain.eq.0) then
                  weather(np:np+5) = 'Rain '
                  np = np+5
                  irain = 1
               endif
            elseif (iw.lt.50) then
c snow, snow pellets, ice crystals = snow
               if (isnow.eq.0) then
                  weather(np:np+5) = 'Snow '
                  np = np+5
                  isnow = 1
               endif
            elseif (iw.lt.60) then
c snow showers, snow squalls, snow grains = snow
               if (isnow.eq.0) then
                  weather(np:np+5) = 'Snow '
                  np = np+5
                  isnow = 1
               endif
            elseif (iw.lt.70) then
c sleet, sleet showers, hail = snow
               if (isnow.eq.0) then
                  weather(np:np+5) = 'Snow '
                  np = np+5
                  isnow = 1
               endif
            elseif (iw.eq.73) then
c dust
               if (idust.eq.0) then
                  weather(np:np+5) = 'Dust '
                  np = np+5
                  idust = 1
               endif
            elseif (iw.eq.74) then
c sand
               if (isand.eq.0) then
                  weather(np:np+5) = 'Sand '
                  np = np+5
                  isand = 1
               endif
            elseif (iw.lt.80) then
c fog
               if (ifog.eq.0) then
                  weather(np:np+4) = 'Fog '
                  np = np+4
                  ifog = 1
               endif
            elseif (iw.eq.84) then
c blowing snow
               if (iblsnow.eq.0) then
                  weather(np:np+7) = 'BlSnow '
                  np = np+7
                  iblsnow = 1
               endif
            elseif (iw.lt.90) then
c smoke, haze, dust, spray, volcanic ash = haze
               if (ihaze.eq.0) then
                  weather(np:np+5) = 'Haze '
                  np = np+5
                  ihaze = 1
               endif
            elseif (iw.lt.99) then
c ice pellets = Ice
               if (iice.eq.0) then
                  weather(np:np+4) = 'Ice '
                  np = np+4
                  iice = 1
               endif
            endif
         endif
      enddo


      return
      end







