      PROGRAM AggregateModelResults
!
!     This program is a DOS program to read daily model flows and storages for
!     all subbasins about a lake and write daily aggregate model flows and
!     storages.  Compile and link as follows (Lahey Fortran 95 used here):
!
!     lf95 ErrorProcessDOS.f95             -c
!     lf95 GLERLutil.f95                   -c
!     lf95 AggregateModelResults.f95       -c
!     lf95 AggregateModelResults.obj GLERLutil.obj ErrorProcessDOS -out AggregateModelResults.EXE
!
      USE GLERLUTIL
      IMPLICIT NONE
      CHARACTER :: BnS*3
      PARAMETER (BnS = 'SUP')
      INTEGER :: NoSubbasins, MaxNoDays
      PARAMETER (NoSubbasins = 22)
      PARAMETER (MaxNoDays = 36525)
      INTEGER :: iD, iM, iY, I, J, K, M, MinM, MaxM
      REAL :: R(MaxNoDays), Area(NoSubbasins), TotArea, U(MaxNoDays),   &
     & L(MaxNoDays), G(MaxNoDays), S(MaxNoDays), P(MaxNoDays), Runoff,  &
     & Uszm, Lszm, Ground, Surface, Snow
      CHARACTER :: N*2, CR*1

      CR = CHAR(13)
      DO M = 1, MaxNoDays
       R(M) = 0.
       U(M) = 0.
       L(M) = 0.
       G(M) = 0.
       S(M) = 0.
       P(M) = 0.
      ENDDO
      DO K = 1, NoSubbasins
       Area(K) = -9999.
      ENDDO
      DO K = 1, NoSubbasins
       N = CHAR(K / 10 + 48) // CHAR(K - (K / 10) * 10 + 48)
       OPEN(UNIT=1,FILE='Param'//N//'.'//BnS,STATUS='OLD',              &
     &  CARRIAGECONTROL='LIST',IOSTAT=I)
       IF (I .EQ. 0) THEN
        READ(1, 1, IOSTAT=I)
        READ(1, 1, IOSTAT=I) Area(K)
        IF (I .NE. 0) STOP 'Could not read area !'
        J = J + 1
       ENDIF
       CLOSE(1, IOSTAT=I)
      ENDDO

      K = 0
      I = 1
      DO WHILE (I .NE. 0 .AND. K .LT. NoSubbasins)
       CLOSE(1, IOSTAT=I)
       K = K + 1
       N = CHAR(K / 10 + 48) // CHAR(K - (K / 10) * 10 + 48)
       OPEN(UNIT=1,FILE='OUTPUT'//N//'.'//BnS,STATUS='OLD',             &
     &  CARRIAGECONTROL='LIST', IOSTAT=I)
      ENDDO
      IF (I .NE. 0) STOP 'No valid OUTPUTxx files present'
      READ(1, 1)
      READ(1, 1)
      READ(1, 1)
      READ(1, 2, IOSTAT=I) iD, iM, iY
      CALL DateSequence(iD, iM, iY, MinM)
      DO WHILE (I .EQ. 0)
       CALL DateSequence(iD, iM, iY, MaxM)
       READ(1, 2, IOSTAT=I) iD, iM, iY
      ENDDO
      CLOSE(1)

      TotArea = 0.
      DO K = 1, NoSubbasins
       N = CHAR(K / 10 + 48) // CHAR(K - (K / 10) * 10 + 48)
       WRITE(6, 4) CR, 'Reading ' // 'OUTPUT' // N // '.' // BnS
       OPEN(UNIT=1,FILE='OUTPUT'//N//'.'//BnS,STATUS='OLD',             &
     &  CARRIAGECONTROL='LIST',IOSTAT=I)
       IF (I .EQ. 0) THEN
        IF (Area(K) .LT. 0.) STOP 'Missing area (parameter file)'
        TotArea = TotArea + Area(K)
        READ(1, 1)
        READ(1, 1)
        READ(1, 1)
        READ(1, 2, IOSTAT=I) iD, iM, iY, Runoff, Uszm, Lszm, Ground,     &
     &   Surface, Snow
        CALL DateSequence(iD, iM, iY, J)
        IF (J .NE. MinM) STOP 'Mismatch in start date'
        DO WHILE (I .EQ. 0)
         CALL DateSequence(iD, iM, iY, J)
         M = J - MinM + 1
         IF (M .LE. MaxNoDays) THEN
          R(M) = R(M) + Runoff * Area(K)
          U(M) = U(M) + Uszm * Area(K)
          L(M) = L(M) + Lszm * Area(K)
          G(M) = G(M) + Ground * Area(K)
          S(M) = S(M) + Surface * Area(K)
          P(M) = P(M) + Snow * Area(K)
         ELSE
          STOP 'Insufficient array dimension'
         ENDIF
         READ(1, 2, IOSTAT=I) iD, iM, iY, Runoff, Uszm, Lszm, Ground,    &
     &    Surface, Snow
        ENDDO
        IF (J .NE. MaxM) STOP 'Mismatch in end date'
       ENDIF
       CLOSE(1)
      ENDDO

      WRITE(6, 4) CR, 'Writing ' // 'OUTPUTDM.' // BnS
      OPEN(UNIT=1,FILE='OUTPUTDM'//'.'//BnS,STATUS='OLD',               &
     & CARRIAGECONTROL='LIST',IOSTAT=I)
      CLOSE(1, STATUS='DELETE',IOSTAT=I)
      OPEN(UNIT=1,FILE='OUTPUTDM'//'.'//BnS,STATUS='NEW',               &
     & CARRIAGECONTROL='LIST')
      WRITE(1, 3) BnS
      DO M = MinM, MaxM
       CALL SequenceDate(iD, iM, iY, M)
       Runoff = R(M - MinM + 1) / TotArea
       Uszm = U(M - MinM + 1) / TotArea
       Lszm = L(M - MinM + 1) / TotArea
       Ground = G(M - MinM + 1) / TotArea
       Surface = S(M - MinM + 1) / TotArea
       Snow = P(M - MinM + 1) / TotArea
       WRITE(1, 2) iD, iM, iY, Runoff, Uszm, Lszm, Ground, Surface, Snow
      ENDDO
      CLOSE(1)
      WRITE(6, 4) CR, '                    '
    1 FORMAT(e13.6e2)
    2 FORMAT(2I2, I4, 6e11.4e2)
    3 FORMAT('This is the distributed model output lumped together ',   &
     & 'for ', A3, '.', /, 'Output from LBRM.  Runoff and storages ',   &
     & 'are in cm over the basin.', /, 'DDMMYYYY     RUNOFF       USZM',&
     & '       LSZM     GROUND    SURFACE       SNOW')
    4 FORMAT(A1, A, $)
      END PROGRAM AggregateModelResults
