* ==================================================================== *
      PROGRAMM SUMSNOT

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  *
* This program calculates analyzes daily climate data of SnoTel-origin *
* and calculates monthly averages (temperature) and sums (precip) per  *
* water year (starting in autumn). There is an option wether to check  *
* for errors or not, before generating monthly values.                 *
*                                                                      *
*      written  (08/15/97):   by: Dave W. Roberts                      *
*                                 Utah State University                *
*                                 Logan, UT 84322-5215                 *
*                                 dvrbts@nr.usu.edu                    *
*                                                                      *
*     modified and much                                                *
*     extended  (11/21/97)                                             *  
*     finalized (04/03/99):   by: Niklaus E. Zimmermann                *
*                                 Utah State University                *
*                                 Logan, UT 84322-5215                 *
*                                 nez@nr.usu.edu                       *
*     Version: 1.2                                                     *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  *
c
c* common stats
c
        real*4 prec(12)
        real*4 avepcp(12)
        real*4 absmax(12)
        real*4 avemax(12)
        real*4 aveave(12)
        real*4 avemin(12)
        real*4 absmin(12)
        integer*2 numday(12)
c
      common / stats / prec,avepcp,absmax,avemax,aveave,avemin,
     +                 numday,absmin
c
c* common dates
c
        integer*2 yr
        integer*2 oldyr
        integer*2 mth
        integer*2 oldmth
        integer*2 day
c
        common / dates / yr,oldyr,mth,oldmth,day
c
c* common station
c
        character*12 stat
        character*6 stanum
        integer*2 type
c
        common / station / stat,stanum,type
c
c* common datval
c
        real*4 tmax
        real*4 tave
        real*4 tmin
        real*4 pcp
        integer*2 ok(0:9)  ! check temper-data
        integer*2 okpcp2   ! check precip-data ....
        integer*2 okpcp1   ! ... and hold it to check for errors at end-of-month
c
        common / datval / tmax,tave,tmin,pcp,ok
c
c* common queue
c
        integer*2 qmax   ! tests for monotonic de/increase in tmax
        integer*2 qave   ! tests for monotonic de/increase in tave
        integer*2 qmin   ! tests for monotonic de/increase in tmin
c
        common / queue / qmax,qave,qmin
c
c* common arglst
c
        integer*2    count
        character*40 argmnt(2)
c
      COMMON / arglst / argmnt,count
c
c* common frost
c
        integer*2    sfrost   ! Date of last spring frost
        integer*2    ffrost   ! date of first fall frost
c
      COMMON / frost / sfrost,ffrost
c
c* local
c
        character*80 args    ! corresponds to "line" in GETARGS
        character*80 line    ! used to read the data-lines
        character*80 pref    ! local variable cont. first argument [prefix of output-file]
        real*4 ttmax         ! used to hold previous value
        real*4 ttave         ! used to hold previous value
        real*4 ttmin         ! used to hold previous value
        integer*2 correct    ! 1=correction, 2=no correction
        integer*2 proc       ! 1=errorfree or correction=2
                             ! 2=T-errors and correction=1
* ------------------- analysis of arguments -------------------------- *
*
** This is UNIX ***
*      call getarg(1,pref)
** This is DOS  ***
      call getcl(args)
      call getargs(args)
      pref = argmnt(1)
      read(argmnt(2),1) correct
    1 format(I3)
      if (correct .ne. 2) correct=1
**
c
      do 10 i=80,1,-1
         if (pref(i:i) .ne. ' ') then
            open (unit=1,file=pref(1:i)//'.txt',status='old')
            open (unit=2,file=pref(1:i)//'.dat',status='unknown')
            open (unit=3,file=pref(1:i)//'.err',status='unknown')
            write (3,12) 
            open (unit=4,file=pref(1:i)//'.cor',status='unknown')
            goto 11
         endif
   10 continue
   11 continue
   12 format('        Error-Codes    Queue     #day',/,
     +       'StatID  123456789 10   Mx Av Mn   TDM   Original ',
     +       'Data Line',/,'------  ------------  ---------  ',
     +       '----   ------------------->>>')

* -------------- Read dataline and check for errors ------------------ *
*
   20 read(1,200,end=99) line              ! test for new station
      if (line(1:4) .eq. 'Stat')  then
          call clear
          call setup(line)                 ! check for format of data
          goto 22
      endif
      goto 20
*
   22 read(1,200,end=99) line              ! test for end of year
        okpcp1 = okpcp2
        if (line(3:6) .eq. ' 930') then
          call endmth(ok,okpcp1)           ! summarize September
          call endyr                       ! summarize water year
          goto 20                          ! .. and start reading data again
      endif
*
      if (line(3:6) .eq. '    ') then      ! test for empty line
          goto 20
      endif
*
      if (type .eq. 1) then
         read(line,220) yr,mth,day,tmax,tave,tmin,pcp
         call chkprec(line,pcp,okpcp2,ok)
         call chktemp(line,tmax,tave,tmin,ok)
*
      else if (type .eq. 2) then
         read(line,220) yr,mth,day,tmax,tmin,tave,pcp
         call chkprec(line,pcp,okpcp2,ok)
         call chktemp(line,tmax,tave,tmin,ok)
*
      else
         read(line,220) yr,mth,day,pcp
         if (line(21:26).eq.'      ' .or. pcp.eq.-99.9) then 
            okpcp2 = 1
         elseif (pcp .lt. 0.0) then
            okpcp2 = 2
         else
            okpcp2 = 0
         endif
*
         tmin  = -99.9
         tave  = -99.9
         tmax  = -99.9
         ok(2) = 2
         ok(0) = 1
*
      endif
*
* -------------- If no errors, analyze ------------------------------- *
*
      proc = 1
*
      if (ok(0).ge.1 .and. CORRECT.eq.1) then
         write(3,240) stanum,ok(1),ok(2),ok(3),ok(4),ok(5),ok(6),ok(7),
     +                ok(8),ok(9),okpcp2,qmax,qave,qmin,
     +                numday(mth),line(1:50)
         proc = 2
*
      elseif (okpcp2.ne.0 .and. ok(0).eq.0) then
         write(3,240) stanum,ok(1),ok(2),ok(3),ok(4),ok(5),ok(6),ok(7),
     +                ok(8),ok(9),okpcp2,qmax,qave,qmin,numday(mth),
     +                line(1:50)
      endif
*
      if (okpcp1.ne.0 .and. oldmth.ne.mth) then
         write(4,250) mth,day,yr,stanum,stat
      endif
*
      if (proc.eq.1) then
*
         if (oldmth .ne. mth) then
            call endmth(ok,okpcp1)
         endif
*
         numday(mth) = numday(mth) + 1
         oldmth = mth
         oldyr = yr
*
         avemin(mth) = avemin(mth) + tmin
         avemax(mth) = avemax(mth) + tmax
         aveave(mth) = aveave(mth) + tave
         avepcp(mth) = pcp
*
         if (numday(mth).eq.1) then
            absmax(mth)=tmax
            absmin(mth)=tmin
         endif
         if (tmax .gt. absmax(mth)) absmax(mth)=tmax
         if (tmin .lt. absmin(mth)) absmin(mth)=tmin
*
         if (tmin.lt.0.0 .and. tmin.ne.-99.9 .and. tmin.ne.-51.3 .and.
     +       tmin.ne.-49.9) then
             call frost(mth,day)
         endif
*
         if (tmax.ne.0.0 .and. tmax.ge.ttmax) qmax = qmax+1
         if (tave.ne.0.0 .and. tave.ge.ttave) qave = qave+1
         if (tmin.ne.0.0 .and. tmin.ge.ttmin) qmin = qmin+1
         if (tmax.ne.0.0 .and. tmax.lt.ttmax) qmax = qmax-1
         if (tave.ne.0.0 .and. tave.lt.ttave) qave = qave-1
         if (tmin.ne.0.0 .and. tmin.lt.ttmin) qmin = qmin-1
         ttmax=tmax
         ttave=tave
         ttmin=tmin
*
      elseif (proc.eq.2) then
*
         if (oldmth .ne. mth) then
            call endmth(ok,okpcp1)
         endif
*
         oldmth = mth
         oldyr = yr
*
         avepcp(mth) = pcp
*
      endif
*
      goto 22
*
   99 if (line(1:10) .eq. '----------') then
         goto 100
      elseif (line(1:8) .eq. '------- ') then
         goto 100
      elseif (line(1:4) .eq. 'Stat') then
         goto 100
      elseif (line(1:1) .eq. ' ') then
         goto 100
      elseif (line(1:1) .eq. '*') then
         goto 100
      else
         call endyr
      endif
  100 continue
*
  200 format (a80)
  220 format (i2,i2,i2,t21,f6.1,t29,f6.1,t37,f6.1,t45,f6.1)
  240 format(a6,2x,9(I1),2x,I1,2x,i3,i3,i3,4x,i2,3x,a50)
  250 format('Check month before ',I2,'/',I2,'/19',I2,
     +       ' for errors at station: ',a6,1x,a12)
*
      end

* ==================================================================== *
        subroutine clear
c
c* common stats
c
        real*4 prec(12)
        real*4 avepcp(12)
        real*4 absmax(12)
        real*4 avemax(12)
        real*4 aveave(12)
        real*4 avemin(12)
        real*4 absmin(12)
        integer*2 numday(12)
c
      common / stats / prec,avepcp,absmax,avemax,aveave,avemin,
     +                 numday,absmin
c
c* common queue
c
        integer*2 qmax   ! tests for monotonic de/increase in tmax
        integer*2 qave   ! tests for monotonic de/increase in tave
        integer*2 qmin   ! tests for monotonic de/increase in tmin
c
        common / queue / qmax,qave,qmin
*
*  common frost
*
        integer*2    sfrost
        integer*2    ffrost
*
      COMMON / frost / sfrost,ffrost
c
* -------------------------------------------------------------------- *
c
        do 10 i=1,12
           avemax(i) = 0.0
           avemin(i) = 0.0
           aveave(i) = 0.0
           avepcp(i) = 0.0
           prec(i) = 0.0
           numday(i) = 0
           absmax(i) = 0.0
           absmin(i) = 0.0
   10   continue
        qmax = 0
        qave = 0
        qmin = 0
        sfrost = 1
        ffrost = 365
c
        return
c
        end

* ==================================================================== *
        subroutine chktemp(line,tmax,tave,tmin,ok)
c
c* common queue
c
        integer*2 qmax   ! tests for monotonic de/increase in tmax
        integer*2 qave   ! tests for monotonic de/increase in tave
        integer*2 qmin   ! tests for monotonic de/increase in tmin
c
        common / queue / qmax,qave,qmin
c
        real*4 tmax,tave,tmin
        integer*2 ok(0:9)
        character*80 line    ! used to read the data-lines
c
* -------------------------------------------------------------------- *
c
        if (line(21:26) .eq. '      ' .or.
     +      line(29:34) .eq. '      ' .or.
     +      line(37:42) .eq. '      ') then
           ok(1) = 1
           ok(0) = ok(0)+1
        else if (tmax.eq.0.0 .and. tmin.eq.0.0 .and. tave.eq.0.0) then
           ok(2) = 2
           ok(0) = ok(0)+1
        else if (tmax .eq. -99.9 .or. tmin .eq. -99.9 .or.
     +           tave .eq. -99.9) then
           ok(3) = 3
           ok(0) = ok(0)+1
        else if (tmax.le.tave .or. tave.le.tmin .or.tmax.le.tmin) then
           ok(4) = 4
           ok(0) = ok(0)+1
        endif
*
        if (tmax.eq.-51.3.or.tave.eq.-51.3.or.tmin.eq.-51.3) then
          ok(5) = 5
           ok(0) = ok(0)+1
        endif
*
        if (tmax.eq.-49.9.or.tave.eq.-49.9.or.tmin.eq.-49.9) then
           ok(6) = 6
           ok(0) = ok(0)+1
        endif
*
        if (tmax.eq.95.8 .or. tave.eq.95.8 .or.tmin.eq.95.8) then
           ok(7) = 7
           ok(0) = ok(0)+1
        endif
*
        if (tmax.eq.95.7 .or. tave.eq.95.7. or.tmin.eq.95.7) then
           ok(8) = 8
           ok(0) = ok(0)+1
        endif
*
        if (abs(qmax).ge.17 .or. abs(qave).ge.17. .or.
     +           abs(qmin).ge.17) then
           ok(9) = 9
           ok(0) = ok(0)+1
        endif
*
c
        return
c
        end

* ==================================================================== *
        subroutine chkprec(line,pcp,okpcp2,ok)
c
        character*80 line
        integer*2 okpcp2
        integer*2 ok(0:9)
        real*4 pcp
c
* -------------------------------------------------------------------- *
c
        if (line(45:50) .eq. '      ') then
           okpcp2 = 1
        elseif (pcp .eq. -99.9) then
           okpcp2 = 1
        elseif (pcp .lt. 0.0) then
           okpcp2 = 2
        else
           okpcp2 = 0
           do 99 i=0,9
              ok(i) = 0
   99      continue
        endif
c
        return
c
        end

* ==================================================================== *
        subroutine endmth(ok,okpcp1)
c
c* common stats
c
        real*4 prec(12)
        real*4 avepcp(12)
        real*4 absmax(12)
        real*4 avemax(12)
        real*4 aveave(12)
        real*4 avemin(12)
        real*4 absmin(12)
        integer*2 numday(12)
c
      common / stats / prec,avepcp,absmax,avemax,aveave,avemin,
     +                 numday,absmin
c
c* common queue
c
        integer*2 qmax   ! tests for monotonic de/increase in tmax
        integer*2 qave   ! tests for monotonic de/increase in tave
        integer*2 qmin   ! tests for monotonic de/increase in tmin
c
        common / queue / qmax,qave,qmin
c
c
c* common dates
c
        integer*2 yr
        integer*2 oldyr
        integer*2 mth
        integer*2 oldmth
        integer*2 day
c
        common / dates / yr,oldyr,mth,oldmth,day
c
c* common station
c
        character*12 stat
        character*6 stanum
        integer*2 type
c
        common / station / stat,stanum,type
c
        integer*2 ok(0:9)
        integer*2 okpcp1
c
* -------------------------------------------------------------------- *
c
        if (ok(9).eq.9) then
           avemax(oldmth) = -99.9
           avemin(oldmth) = -99.9
           aveave(oldmth) = -99.9
           absmax(oldmth) = -99.9
           absmin(oldmth) = -99.9
           if (okpcp1 .ne. 0) then
              avepcp(oldmth) = -99.9
           endif
        elseif  (numday(oldmth) .lt. 15) then
           avemin(oldmth) = -99.9
           avemax(oldmth) = -99.9
           aveave(oldmth) = -99.9
           absmax(oldmth) = -99.9
           absmin(oldmth) = -99.9
           if (okpcp1 .ne. 0) then
              avepcp(oldmth) = -99.9
           endif
        else
           avemax(oldmth) = avemax(oldmth) / numday(oldmth)
           avemin(oldmth) = avemin(oldmth) / numday(oldmth)
           aveave(oldmth) = aveave(oldmth) / numday(oldmth)
           if (okpcp1 .ne. 0) then
              avepcp(oldmth) = -99.9
           endif
        endif
c
        if (okpcp1.ne.0 .and. oldmth.eq.9) then
            write(4,99) oldmth,yr,stanum,stat
        endif
c
        qmax = 0
        qave = 0
        qmin = 0
c
   99 format('Check the month of ',I2,'/30/19',I2,
     +       ' for errors at station: ',a6,1x,a12) 
c
        return
c
        end

* ==================================================================== *
        subroutine endyr
c
c* common stats
c
        real*4 prec(12)
        real*4 avepcp(12)
        real*4 absmax(12)
        real*4 avemax(12)
        real*4 aveave(12)
        real*4 avemin(12)
        real*4 absmin(12)
        integer*2 numday(12)
c
      common / stats / prec,avepcp,absmax,avemax,aveave,avemin,
     +                 numday,absmin
c
c* common dates
c
        integer*2 yr
        integer*2 oldyr
        integer*2 mth
        integer*2 oldmth
        integer*2 day
c
        common / dates / yr,oldyr,mth,oldmth,day
c
c* common station
c
        character*12 stat
        character*6 stanum
        integer*2 type
c
        common / station / stat,stanum,type
*
*  common frost
*
        integer*2    sfrost
        integer*2    ffrost
*
      COMMON / frost / sfrost,ffrost
c
* -------------------------------------------------------------------- *
c
        do 50 i=9,2,-1
           if (avepcp(i-1).ne.-99.9 .and. avepcp(i).ne.-99.9) then
              prec(i) = avepcp(i) - avepcp(i-1)
           else
              prec(i) = -99.9
           endif
   50   continue
*
        if (avepcp(1).ne.-99.9 .and. avepcp(12).ne.-99.9) then
           prec(1) = avepcp(1) - avepcp(12)
        else
           prec(1) = -99.9
        endif
*
        if (avepcp(12).ne.-99.9 .and. avepcp(11).ne.-99.9) then
           prec(12) = avepcp(12) - avepcp(11)
        else
           prec(12) = -99.9
        endif
*
        if (avepcp(11).ne.-99.9 .and. avepcp(10).ne.-99.9) then
           prec(11) = avepcp(11) - avepcp(10)
        else
           prec(11) = -99.9
        endif
*
        if (avepcp(10).ne.-99.9) then
           prec(10) = avepcp(10)
        else
           prec(10) = -99.9
        endif
c
        write(2,100) stanum,stat,oldyr,(absmax(j),j=1,12)
        write(2,110) stanum,stat,oldyr,(avemax(j),j=1,12)
        write(2,120) stanum,stat,oldyr,(aveave(j),j=1,12)
        write(2,130) stanum,stat,oldyr,(avemin(j),j=1,12)
        write(2,140) stanum,stat,oldyr,(absmin(j),j=1,12)
        write(2,150) stanum,stat,oldyr,(prec(j),j=1,12)
        write(2,160) stanum,stat,oldyr,(numday(j),j=1,12)
        write(2,170) stanum,stat,oldyr,sfrost,ffrost
c
        return
c
  100 format (a6,2x,a12,1x,'19',i2,1x,'AMX',12(f8.1))
  110 format (a6,2x,a12,1x,'19',i2,1x,'MAX',12(f8.1))
  120 format (a6,2x,a12,1x,'19',i2,1x,'AVE',12(f8.1))
  130 format (a6,2x,a12,1x,'19',i2,1x,'MIN',12(f8.1))
  140 format (a6,2x,a12,1x,'19',i2,1x,'AMN',12(f8.1))
  150 format (a6,2x,a12,1x,'19',i2,1x,'PRE',12(f8.1))
  160 format (a6,2x,a12,1x,'19',i2,1x,'DAY',12(i8))
  170 format (a6,2x,a12,1x,'19',i2,1x,'SFF',2(i8))

        end

* ==================================================================== *
        subroutine setup(line)
c
c* common dates
c
        integer*2 yr
        integer*2 oldyr
        integer*2 mth
        integer*2 oldmth
        integer*2 day
c
        common / dates / yr,oldyr,mth,oldmth,day
c
c* common station
c
        character*12 stat
        character*6 stanum
        integer*2 type
c
        common / station / stat,stanum,type
c
c* passed
c
        character*80 line
c
* -------------------------------------------------------------------- *
c
        oldmth = 10
        stat = line(19:31)
        stanum = line(11:16)
*
        do 10 i=1,4               ! jump the next 4 unimportant lines
           read(1,100) line
   10   continue
*
        read(1,100) line          ! start reading the first line of this station and determine the data format
        if (line(23:26) .eq. 'tmax') then
           if (line(31:34) .eq. 'tavg') then
              type = 1
           else
              type = 2
           endif
        else
              type = 3
        endif
c
        return
c
  100 format (a80)
c
        end

* ==================================================================== *
      SUBROUTINE FROST(mth,day)  

*   SUBROUTINE FROST determines for each station the number of days    *
*   between the last spring frost and the first frost in fall. It      *
*   consideres the 15th of July as the dividing line between spring    *
*   and fall, and calculates the number of days from that date.        *

*
*  common frost
*
        integer*2    sfrost
        integer*2    ffrost
*
      COMMON / frost / sfrost,ffrost
*
        integer*2    mth
        integer*2    day
c
* -------------------------------------------------------------------- *
c
      if (mth .eq. 1) then 
         sfrost = day
      elseif (mth .eq. 2) then
         sfrost = day + 31
      elseif (mth .eq. 3) then
         sfrost = day + 59
      elseif (mth .eq. 4) then
         sfrost = day + 90
      elseif (mth .eq. 5) then
         sfrost = day + 120
      elseif (mth .eq. 6) then
         sfrost = day + 151
      elseif (mth.eq.7 .and. day.le.15) then
         sfrost = day + 181
      elseif (mth.eq.7 .and. day.gt.15 .and. (day+181).lt.ffrost) then
         ffrost = day + 181
      elseif (mth.eq.8 .and. (day+212).lt.ffrost) then
         ffrost = day + 212
      elseif (mth.eq.9 .and. (day+243).lt.ffrost) then
         ffrost = day + 243
      elseif (mth.eq.10 .and. (day+273).lt.ffrost) then
         ffrost = day + 273
      elseif (mth.eq.11 .and. (day+304).lt.ffrost) then
         ffrost = day + 304
      elseif (mth.eq.12 .and. (day+334).lt.ffrost) then
         ffrost = day + 334
      endif
*
      return
* 
      end

* ==================================================================== *
      SUBROUTINE GETARGS(line)

*   SUBROUTINE GETCL takes a line of text and breaks it down into      *
*   arguments by checking for blanks between words.  The arguments     *
*   are ordered and stored in array ARGMNT.  The number of charac-     *
*   ters in each argument is stored in array SIZE, and the number      *
*   of arguments passed is stored in COUNT.                            *


      CHARACTER*40 argmnt(2)    ! list of command line arguments
      INTEGER*2 count           ! number of arguments passed
*
      COMMON / arglst / argmnt,count
*
***DOS and VMS only
*
       CHARACTER*80  line
       INTEGER*2     left,right
c
* -------------------------------------------------------------------- *
c
      left = 1
      right = 0
      count = 0
*
      DO 10 i=1,2
      argmnt(i) = ' '
   10 CONTINUE
*
      DO 11 I=1,32
      IF (line(i:i) .NE. ' ') THEN
        RIGHT = i
      ELSE IF (i .NE. 1 .AND. line(i-1:i-1) .NE. ' ') THEN
        count = count + 1
        argmnt(count) = line(left:right)
        left = I+1
      ELSE
        left = I+1
      ENDIF
   11 CONTINUE
*
      RETURN
*
      END

* ==================================================================== *