      PROGRAM ETRC_AGGREGATE_WEEKLY
C     *****************************
C     Aggregates weekly time-series from daily NLDAS-derived ASCE ETrs, for a
C     specific date. Output is a weekly total. Includes February 29s as a
C     normal date.

C     Author: Mike Hobbins
C     Date: January 30, 2018

      IMPLICIT  NONE
      INTEGER   mlen(12),nrows,ncols,numcrit,bnd_s,bnd_n,bnd_w,bnd_e,
     +          bnd_lo,bnd_hi,yearin,monin,dayin,row,col,colin,rowin,
     +          day,mon,year,dayct,j,lpdy,count
      REAL      Lat,Lon,Elev
      CHARACTER path*99,temppath*99,timefile*99,infile*99,outfile*99,
     +          header*19
      DATA      mlen/31,28,31,30,31,30,31,31,30,31,30,31/

C     Define number of rows and columns, variables, and years of analysis.
      PARAMETER (nrows    = 224,
     +           ncols    = 464,
     +           numcrit  = 7,
     +           bnd_s    = 25,
     +           bnd_n    = 53,
     +           bnd_w    = -125,
     +           bnd_e    = -67,
     +           bnd_lo   = -90,
     +           bnd_hi   = 4500,
     +           path     = '/Volumes/Drobo1/ETrc_data/',
     +           temppath =
     +        '/Users/mhobbins/scripts/EDDI/Automated_delivery/tempry/')
      INTEGER    Mask(nrows,ncols)
      REAL       in(7,nrows,ncols),out(ncols)


C     Read specifications file for date
      WRITE(timefile,'(A55,A13)') temppath,'Timestamp.tmp'
      OPEN(UNIT=11,FILE=timefile,STATUS='OLD',ERR=99)
      READ(11,'(I4,2(X,I2))') yearin,monin,dayin
      CLOSE(11)
      PRINT*,'ETrc_aggregate_weekly:   ',
     +       'Year',yearin,'; Month',monin,'; Day',dayin

C     Read in geographic data to set mask for analysis
      OPEN(UNIT=12,FILE='/Users/mhobbins/scripts/gtopomean15k.asc',
     +     STATUS='OLD',ERR=99)
      DO row = nrows, 1, -1
        DO col = 1, ncols
C         - read input data: Col; Row; Lat [deg], latitude; Lon [deg],
C         longitude; Elev [m], elevation
          READ(12,*) colin,rowin,Lat,Lon,Elev
          IF((rowin.NE.nrows-row+1).OR.(colin.NE.col)) THEN
            PRINT*,'Static grid input error: col ',colin,'; row ',rowin
            GOTO 99
          ENDIF
C         Generate geographic analysis mask
          Mask(row,col) = 0
          IF((Lat.GT.bnd_s).AND.(Lat.LT.bnd_n).AND.
     +       (Lon.GT.bnd_w).AND.(Lon.LT.bnd_e).AND.
     +       (Elev.GT.bnd_lo).AND.(Elev.LT.bnd_hi))
     +       Mask(row,col) = 1
        ENDDO                                    ! end col loop
      ENDDO                                      ! end row loop
      CLOSE(12)

      day  = dayin
      mon  = monin
      year = yearin

C     Open the weekly output file
      WRITE(outfile,'(A26,A16,I4,2I2.2,A4)')
     +      path,'weekly_agg/ETrs_',year,mon,day,'.asc'
      OPEN(UNIT=13,FILE=outfile,STATUS='UNKNOWN')

C     Set start date for week
      DO dayct = 1,7
C       Read day's input data
        WRITE(infile,'(A26,A11,I4,2I2.2,A4)')
     +        path,'daily/ETrs_',year,mon,day,'.asc'
        OPEN(UNIT=13+dayct,FILE=infile,STATUS='OLD',ERR=99)
        DO j = 1, 6
          READ(13+dayct,'(A19)') header
          IF(dayct.EQ.1) WRITE(13,'(A19)') header
        ENDDO
C       Fill input array for all rows, columns, and days
        DO row = 1, nrows
          READ(13+dayct,*) (in(dayct,row,col),col=1,ncols)
        ENDDO
        CLOSE(13+dayct)
C       - go back one day
        day = day - 1
        IF(day.LT.1) THEN
          mon = mon - 1
          IF(mon.LT.1) THEN
            mon  = 12
            year = year - 1
          ENDIF
          lpdy = 0                               ! adjust month-length for Feb 29s
          IF((MOD(year,4).EQ.0).AND.(mon.EQ.2)) lpdy = 1
          day = mlen(mon) + lpdy
        ENDIF
      ENDDO

C     Process the data arrays for each pixel in each day
      DO 40 row = 1, nrows
        DO 30 col = 1, ncols
          out(col) = 0.
          count    = 0
          DO 20 dayct = 1,7
            IF(Mask(row,col).EQ.1) THEN
              IF(in(dayct,row,col).NE.-9999.) THEN
                count = count + 1
                out(col) = out(col) + in(dayct,row,col)
              ENDIF
            ELSE
              out(col) = -9999.
            ENDIF
20        ENDDO                                  ! end dayct loop
          IF(count.GE.numcrit) THEN
            out(col) = out(col) * 7. / count
          ELSE
            out(col) = -9999.
          ENDIF
30      ENDDO                                    ! end column loop

C       Write the pixel's data to the output file
        WRITE(13,*) (out(col),col=1,ncols)
40    ENDDO                                      ! end row loop

C     Close output file
      CLOSE(13)

99    END