SUBROUTINE PSIS(TIME, UU, TSTAR, T, S, SI, OSINT, NS, A2, II,
     1                M, PSI, DPSI)

C     + + + PURPOSE + + +
C
C     THIS ROUTINE COMPUTES FUNCTION PSI AND ITS DERIVATIVE DPSI  
C
C     ESSENTIALLY, PSI AND DPSI CAN BE CALCULATED FROM THE INTEGRAL OF 
C     S WRT TIME (THE LATTER IS CALCULATED BY SINT)
C
C     CALLED FROM PSIINV
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
C     + + + ARGUMENT DECLARATIONS + + +
C
      INTEGER NS, II, IU, IL, K
      REAL TSTAR, T(MXTIME), S(MXTIME), SI(MXTIME+1),
     1     A2, M, OSINT
      DOUBLE PRECISION UU, PSI, DPSI, TIME    

C     + + + ARGUMENT DEFINITIONS + + +
C
C     NS     -
C     II     -
C     TIME   -
C     TSTAR  -
C     T      -
C     S      -
C     SI     -
C     A2     -
C     M      -
C     OSINT  -
C     UU     -
C     PSI    -
C     DPSI   -
C
C     + + + LOCAL VARIABLES + + +
C
      DOUBLE PRECISION S1, S2, S1TOA2, S2TOA2, A, B, XU
C
C     + + + END SPECIFICATIONS + + +

C
      IF (TIME.GE.TSTAR) THEN
         B = SI(NS+1)
         IU = NS + 1
      ELSE
         CALL SINT(TIME, T, TSTAR, II, S, SI, NS, OSINT)
         B = OSINT
         IU = II + 1
      END IF
C
      XU = UU
      CALL SINT(XU, T, TSTAR, II, S, SI, NS, OSINT)
      A = OSINT
      IL = II + 1
      K = IL
      PSI = 0.D0
      DPSI = 0.D0
      S1 = 0.D0
      S1TOA2 = 0.D0
C
   10 IF (K.NE.IU) THEN
         S2 = MAX(SI(K)-A,0.D0)
         S2TOA2 = S2 ** A2
         PSI = PSI + (S2*S2TOA2-S1*S1TOA2) / S(K-1)
         DPSI = DPSI + (S2TOA2-S1TOA2) / S(K-1)
         K = K + 1
         S1 = S2
         S1TOA2 = S2TOA2
         GO TO 10
      END IF
C
      S2 = MAX(B-A,0.D0)
      S2TOA2 = S2 ** A2
      
      PSI = (PSI+(S2*S2TOA2-S1*S1TOA2)/S(K-1)) / M
      DPSI = -S(IL-1) * (DPSI+(S2TOA2-S1TOA2)/S(K-1))
      
      IF (TIME.LE.TSTAR.OR.S2.EQ.0.) RETURN
C   
      PSI = PSI + S2TOA2 * (TIME-TSTAR)
      DPSI = DPSI - S(IL-1) * A2 * S2TOA2 / S2 * (TIME-TSTAR)
C
      RETURN
      END