!$Author: fredfox $
!$Date: 2006-02-25 00:21:16 $
!$Revision: 1.2 $
!$Source: /weru/cvs/wepp/wepp.watbal/hdrive.for,v $

      SUBROUTINE HDRIVE(NQT, DURPQ, QTP, TPEE,
     1                  U, ISEED, II, DT, NS, NT, NQ, QTOT, Q, TQ1,
     1                  RECUM, T, S, SI, SLEN, ALPHA, M, DUREXR, A1, A2,
     1                  TSTAR, PEAKRO, DURRUN)

C     + + + PURPOSE + + +
C
C     HDRIVE COMPUTES KINEMATIC FLOW DEPTHS ON A PLANE AT SELECTED
C     DISTANCES DOWN THE PLANE AND AT SELECTED TIME SPACING  
C
C     THE LATERAL INFLOW (RAINFALL EXCESS) IS PRESENTED AS A POSITIVE 
C     STEP FUNCTION UP TO A GIVEN TIME AND IS ZERO THEREAFTER (IF 
C     STEP VALUE IS GIVEN AS ZERO, IT IS SET TO 1.E-8)
C
C     SUBROUTINE HDEPTH COMPUTES DEPTH OF FLOW ON THE PLANE
C
C     CALLED FROM MAIN
C     AUTHOR(S): D. FLANAGAN, J. ASCOUGH
C     VERSION: THIS MODULE TAKEN FROM ASCOUGH STANDALONE IRS CODE
C     DATE CODED:  3-28-2005
C     CODED BY: D. FLANAGAN
C
C     + + + PARAMETER DECLARATIONS + + +
C      
      INTEGER MXTIME
      PARAMETER (MXTIME = 1000)

C     + + + ARGUMENT DECLARATIONS + + +
C
      INTEGER II, NQ, NS, NT, NQT
      REAL DURPQ, QTP, TPEE, ISEED, DT, QTOT(MXTIME),
     1     Q(MXTIME+1), TQ1(MXTIME), RECUM(MXTIME), T(MXTIME), 
     1     S(MXTIME), SI(MXTIME+1), SLEN, ALPHA, M, DUREXR, A1, A2, 
     1     TSTAR, PEAKRO, DURRUN
      DOUBLE PRECISION U

C     + + + ARGUMENT DEFINITIONS + + +
C
C     II      -
C     NQ      -
C     NS      -
C     NT      -
C     NQT     -
C     DURPQ   -
C     QTP     -
C     TPEE    -
C     Z       -
C     ISEED   -
C     DT      -
C     QTOT    -
C     Q       -
C     TQ1     -
C     RECUM   -
C     T       -
C     S       -
C     SI      -
C     SLEN     -
C     ALPHA   -
C     M       -
C     DUREXR  -
C     A1      -
C     A2      -
C     TSTAR   -
C     PEAKRO  -
C     DURRUN  -
C     U       -
C
C     + + + LOCAL VARIABLES + + +
C
      INTEGER BEGRUN, I, NQI, IQT
      REAL I1, LQ, QTMAX, BEGTIM, D, QMAX, QMAX10, T1
      REAL X0, X, A, MRND, HDPTHO
      DOUBLE PRECISION T2, TQNEW
      
C     + + + END SPECIFICATIONS + + +
C
      X0 = 1948.
      CALL BGNRND(X0, X, A, MRND)
C
      U = 0.D0
      ISEED = 3734923
      BEGRUN = 0
      II = 1
      DT = 1.
      NT = NS + 10
      NQ = NT
      I1 = 0.D0
      LQ = 0.
C     
      DO I = 1, 1000
         QTOT(I) = 0.
         Q(I) = 0.
         TQ1(I) = 0.
      END DO
C     
      QTMAX = 0.97 * RECUM(NS+1)
C     
      I = 0
C
   20 I = I + 1
C
      IF (I.LE.(NS+1)) THEN
         T2 = T(I)
      ELSE
         T2 = T2 + DT
      END IF
C     

      CALL HDEPTH(T2, SLEN, A1, A2, TSTAR, T, S, SI, NS, II,
     1            M, HDPTHO, A, MRND)

      D = HDPTHO
      TQ1(I) = T2
C
      IF (BEGRUN.EQ.0.AND.D.NE.0.) THEN
         BEGRUN = 1
         BEGTIM = T2
      END IF
C
      Q(I) = ALPHA * D ** M
C     
      IF (I.GT.1) THEN
         I1 = I1 + (Q(I)+LQ) * (T2-T1) / 2.0
         QTOT(I+1) = I1
      END IF
C     
      LQ = Q(I)
      T1 = T2
      IF (I.EQ.1) QMAX = Q(I)
C     
      IF (Q(I).GT.QMAX) THEN
         NQI = I
         QMAX = Q(I)
         QTP = T2
         IQT = I
         QMAX10 = .1 * QMAX
      END IF
C     
      IF (I.LT.999) THEN
         IF ((QTOT(I+1)*(1000.0/SLEN)).LT.QTMAX) THEN
            GO TO 20
         ELSE
            IF (Q(I).GT.QMAX10) GO TO 20
         END IF
      END IF
C     
      NQT = I + 1
      TQ1(NQT) = T2 + DT
      TQNEW = TQ1(NQT)
      CALL HDEPTH(TQNEW, SLEN, A1, A2, TSTAR, T, S, SI, NS, II,
     1            M, HDPTHO, A, MRND)
      D = HDPTHO
      Q(NQT) = ALPHA * D ** M
C     
      PEAKRO = Q(NQI) * (60000.0/SLEN)
C     
      DO I = IQT, NQT
         IF (Q(I+1).LT.Q(I)) THEN
            TPEE = TQ1(I)
            DURPQ = TPEE - QTP
            GO TO 40
         END IF
      END DO 
C     
   40 IF ((TQ1(NQT)/60.0).GT.DUREXR) THEN
         DURRUN = (TQ1(NQT)-BEGTIM) / 60.0
      ELSE
         DURRUN = DUREXR
      END IF
C     
      RETURN
      END