SUBROUTINE PSIINV(TIME, X, TSTAR, T, S, SI, OSINT, NS, A2,
     1                  II, M, PSI, DPSI, OPSII, A, MRND)

C
C     + + + PURPOSE + + +
C
C     SUBROUTINE PSIINV IS USED TO COMPUTE BOTH PSI AND ITS DERIVATIVE
C     DPSI
C
C     THIS SUBROUTINES RETURNS THE INVERSE OF THE PSI FUNCTION
C
C     SINCE PSI IS CONVEX FOR FIXED TIME, A SIMPLE APPLICATION OF 
C     NEWTON'S METHOD WILL SUFFICE TO INVERT PSI(TIME,U)=X 
C
C     IF TIME IS CHANGED BY A SMALL AMOUNT BETWEEN CALLS, THE OLD 
C     PSIINV VALUE PROVIDES A GOOD STARTING POINT FOR THE NEXT VALUE
C
C     CALLED FROM HDEPTH
C     AUTHOR(S): D. FLANAGAN, J. ASCOUGH
C     VERSION: THIS MODULE TAKEN FROM ASCOUGH STANDALONE IRS CODE
C     DATE CODED:  3-24-2005
C     CODED BY: D. FLANAGAN
C
C     + + + PARAMETER DECLARATIONS + + +
C
      INTEGER MXTIME
      PARAMETER (MXTIME = 1000)
C
C     + + + ARGUMENT DECLARATIONS + + +
C
      INTEGER NS, II
      REAL X, TSTAR, T(MXTIME), S(MXTIME), SI(MXTIME+1),
     1     OSINT, A2, M, A, MRND
      DOUBLE PRECISION PSI, DPSI, OPSII, TIME    

C     + + + ARGUMENT DEFINITIONS + + +
C
C     NS      -
C     II      -
C     TIME    -
C     X       -
C     TSTAR   -
C     T       -
C     S       -
C     SI      -
C     OSINT    -
C     A2      -
C     M       -
C     PSI     -
C     DPSI    -
C     OPSII   - OUTPUT VALUE FROM PSIINV (FORMERLY FUNCTION ARGUMENT)
C
C     + + + LOCAL VARIABLES + + +
C
      REAL TO, TECO, RNUMB
      DOUBLE PRECISION U, UM, UB, UT

C     + + + END SPECIFICATIONS + + +
C
      U = 0.
      TO = 0.
      TECO = MIN(TIME,TSTAR)
      UB = TO
      UT = TECO
C
C     U INITIALLY DOES NOT SEEM TO HAVE A VALUE JCAII 10-13-04
C     INITIALIZED THROUGH SALFORD FORTRAN
C
      IF (U.LE.TO.OR.TECO.LE.U) U = (TO+TECO) / 2.D0
C
   10 CALL PSIS(TIME, U, TSTAR, T, S, SI, OSINT, NS, A2, II,
     1          M, PSI, DPSI)
C
      IF (ABS(PSI-X).GE..005*X) THEN
C          
         IF (DPSI.EQ.0.D0) THEN
C         
            CALL RANDM(X, A, MRND, RNUMB)
            U = UB + (UT-UB) * RNUMB / 2.D0
         ELSE
            UM = U - (PSI-X) / DPSI
C
            IF (UM.LT.U) THEN
               UT = U
            ELSE
               UB = U
            END IF
C
            IF ((UT-UB).LT.0.00001D0) GO TO 20
C
            IF (UM.LE.UB) THEN
               U = (UB+U) / 2.D0
            ELSE IF (UM.LT.UT) THEN
               U = UM
            ELSE
               U = (UT+U) / 2.D0
            END IF
C
         END IF
C        
         GO TO 10
C
      END IF
C     
   20 OPSII = U
C     
      RETURN
      END