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

C     + + + PURPOSE + + +
C
C     DEPTH H IS THE SOLUTION TO THE PARTIAL DIFFERENTIAL EQUATION:
C
C     H1+ALPHA*M*H**(M-1)*H2=S ,T AND X NON-NEGATIVE    (1)
C
C     WHERE:
C
C     H1    =  PARTIAL OF H WITH RESPECT TO T,
C     H2    =  PARTIAL OF H WITH RESPECT TO X,
C     ALPHA =  CHEZY DEPTH-DISCHARGE COEFFICIENT,
C     M     =  CHEZY DEPTH-DISCHARGE EXPONENT,
C     S     =  STEP FUNCTION WITH NS INITIAL NON-ZERO STEPS THAT IS
C              ZERO AFTER THE LAST STEP, AND
C     H IS ZERO WHEN EITHER T OR X IS ZERO.
C
C     LET PHI(T) BE PSI(T,0), WHERE PSI(T,U) IS DEFINED TO BE THE 
C     INTEGRAL FROM U TO T WRT L OF THE M-1TH POWER OF THE INTEGRAL 
C     OF S FROM TIME U TO L  
C
C     PSI(T,U) STRICTLY DECREASES FROM PHI(T) TO 0 FOR FIXED T AND
C     FOR U BETWEEN 0 AND AM1N(T,TSTAR) (S>0 BEFORE TSTAR,S=0 AFTER) 
C 
C     LET PSIINV(T,-) BE THE INVERSE OF PSI(T,-) - THEN THE GENERAL 
C     SOLUTION TO (1) IS GIVEN BY:
C
C     FOR X >= X1, H(T,X)=INTEGRAL OF S FROM 0 TO T             (2)
C     FOR X <= X1, H(T,X)=INTEGRAL OF S FROM PSIINV(T,X2) TO T  (3)
C
C     WHERE X1 = ALPHA*M*PHI(T) AND X2 = X/(ALPHA*M).
C
C     SUBROUTINE PHI(-) COMPUTES PHI, SUBROUTINE PSIINV(T,-) COMPUTES 
C     THE INVERSE OF PSI(T,-), AND SUBROUTINE SINT(T) COMPUTES THE 
C     INTEGRAL OF S FROM 0 TO TIME
C
C     CALLED FROM HDRIVE
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 + + +

      INTEGER II, NS
      REAL A1, A2, X, S(MXTIME), T(MXTIME), SI(MXTIME+1),
     1     TSTAR, M, HDPTHO, A, MRND
      DOUBLE PRECISION T2

C
C     + + + ARGUMENT DEFINITIONS + + +
C
C     A1      -
C     A2      -
C     T2      -
C     X       -
C     S       -
C     T       -
C     SI      -
C     TSTAR   -
C     M       -
C     HDPTHO  -
C
C     + + + LOCAL VARIABLES + + +
C  
      REAL TERM1, TERM2, OPHI, OSINT, XX
      DOUBLE PRECISION PSI, DPSI, OPSII

C     + + + END SPECIFICATIONS + + +
C
      XX = X / A1
C
      CALL PHI(T2, TSTAR, SI, NS, II, M, A2, S, T, OPHI)
      IF (XX.GE.OPHI) THEN
         CALL SINT(T2, T, TSTAR, II, S, SI, NS, OSINT)
         HDPTHO = OSINT
         RETURN
      END IF
C
C     REPLACE ALL FUNCTIONS WITH SUBROUTINE CALLS
C     HDEPTH = SINT(T2) - SINT(PSIINV(T2,XX))
C
      CALL SINT(T2, T, TSTAR, II, S, SI, NS, OSINT)
      TERM1 = OSINT
      CALL PSIINV(T2, XX, TSTAR, T, S, SI, OSINT, NS, A2,
     1            II, M, PSI, DPSI, OPSII, A, MRND)
      CALL SINT(OPSII, T, TSTAR, II, S, SI, NS, OSINT)
      TERM2 = OSINT

      HDPTHO = TERM1 - TERM2
C
      RETURN
      END