SUBROUTINE SINTRP(X,Y,XOUT,YOUT,YPOUT,NEQN,KOLD,PHI,IVC,IV,KGI,GI,
& ALPHA,G,W,XOLD,P)
C
C***BEGIN PROLOGUE SINTRP
C***DATE WRITTEN 740101 (YYMMDD)
C***REVISION DATE 840201 (YYMMDD)
C***CATEGORY NO. D2A2
C***KEYWORDS INITIAL VALUE ORDINARY DIFFERENTIAL EQUATIONS,
C VARIABLE ORDER ADAMS METHODS, SMOOTH INTERPOLANT FOR
C DEABM IN THE DEPAC PACKAGE
C***AUTHOR SHAMPINE, L.F., SNLA
C GORDON, M.K.
C MODIFIED BY H.A. WATTS
C***PURPOSE APPROXIMATES THE SOLUTION AT XOUT BY EVALUATING THE
C POLYNOMIAL COMPUTED IN STEPS AT XOUT. MUST BE USED IN
C CONJUNCTION WITH STEPS.
C***DESCRIPTION
C
C WRITTEN BY L. F. SHAMPINE AND M. K. GORDON
C
C ABSTRACT
C
C
C THE METHODS IN SUBROUTINE STEPS APPROXIMATE THE SOLUTION NEAR X
C BY A POLYNOMIAL. SUBROUTINE SINTRP APPROXIMATES THE SOLUTION AT
C XOUT BY EVALUATING THE POLYNOMIAL THERE. INFORMATION DEFINING THIS
C POLYNOMIAL IS PASSED FROM STEPS SO SINTRP CANNOT BE USED ALONE.
C
C THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT,
C COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS, THE INITIAL
C VALUE PROBLEM BY L. F. SHAMPINE AND M. K. GORDON.
C FURTHER DETAILS ON USE OF THIS CODE ARE AVAILABLE IN *SOLVING
C ORDINARY DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND INTRP*,
C BY L. F. SHAMPINE AND M. K. GORDON, SLA-73-1060.
C
C INPUT TO SINTRP --
C
C THE USER PROVIDES STORAGE IN THE CALLING PROGRAM FOR THE ARRAYS IN
C THE CALL LIST
C DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),P(NEQN),
C ALPHA(12),G(13),W(12),GI(11),IV(10)
C AND DEFINES
C XOUT -- POINT AT WHICH SOLUTION IS DESIRED.
C THE REMAINING PARAMETERS ARE DEFINED IN STEPS AND PASSED TO
C SINTRP FROM THAT SUBROUTINE.
C
C OUTPUT FROM SINTRP --
C
C YOUT(*) -- SOLUTION AT XOUT
C YPOUT(*) -- DERIVATIVE OF SOLUTION AT XOUT
C
C THE REMAINING PARAMETERS ARE RETURNED UNALTERED FROM THEIR INPUT
C VALUES. INTEGRATION WITH STEPS MAY BE CONTINUED.
C
C***REFERENCES SHAMPINE L.F., GORDON M.K., *SOLVING ORDINARY
C DIFFERENTIAL EQUATIONS WITH ODE, STEP, AND INTRP*,
C SLA-73-1060, SANDIA LABORATORIES, 1973.
C WATTS H.A., SHAMPINE L.F., *A SMOOTHER INTERPOLANT FOR
C DE/STEP,INTRP : II*, SAND84-0293, SANDIA LABORATORIES,
C 1984.
C***ROUTINES CALLED (NONE)
C***END PROLOGUE SINTRP
C
USE REAL_PRECISION
REAL (KIND=R8):: ALP,ALPHA,C,G,GAMMA,GDI,GDIF,GI,GTEMP,
& H,HI,HMU,P,PHI,RMU,SIGMA,TEMP1,TEMP2,TEMP3,W,WTEMP,
& X,XI,XIM1,XIQ,XOLD,XOUT,Y,YOUT,YPOUT
INTEGER I,IQ,IV,IVC,IW,J,JQ,KGI,KOLD,KP1,KP2,L,M,NEQN
C
DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),P(NEQN)
DIMENSION GTEMP(13),C(13),WTEMP(13),G(13),W(12),ALPHA(12),
& GI(11),IV(10)
C
C***FIRST EXECUTABLE STATEMENT
KP1 = KOLD + 1
KP2 = KOLD + 2
C
HI = XOUT - XOLD
H = X - XOLD
XI = HI/H
XIM1 = XI - 1.
C
C INITIALIZE WTEMP(*) FOR COMPUTING GTEMP(*)
C
XIQ = XI
DO IQ = 1,KP1
XIQ = XI*XIQ
TEMP1 = IQ*(IQ+1)
WTEMP(IQ) = XIQ/TEMP1
END DO
C
C COMPUTE THE DOUBLE INTEGRAL TERM GDI
C
IF (KOLD .LE. KGI) THEN
GDI = GI(KOLD)
ELSE
IF (IVC .GT. 0) THEN
IW = IV(IVC)
GDI = W(IW)
M = KOLD - IW + 3
ELSE
GDI = 1.0/TEMP1
M = 2
END IF
IF (M .LE. KOLD) THEN
DO I = M,KOLD
GDI = W(KP2-I) - ALPHA(I)*GDI
END DO
END IF
END IF
C
C COMPUTE GTEMP(*) AND C(*)
C
GTEMP(1) = XI
GTEMP(2) = 0.5*XI*XI
C(1) = 1.0
C(2) = XI
IF (KOLD .GE. 2) THEN
DO I = 2,KOLD
ALP = ALPHA(I)
GAMMA = 1.0 + XIM1*ALP
L = KP2 - I
DO JQ = 1,L
WTEMP(JQ) = GAMMA*WTEMP(JQ) - ALP*WTEMP(JQ+1)
END DO
GTEMP(I+1) = WTEMP(1)
C(I+1) = GAMMA*C(I)
END DO
END IF
C
C DEFINE INTERPOLATION PARAMETERS
C
SIGMA = (WTEMP(2) - XIM1*WTEMP(1))/GDI
RMU = XIM1*C(KP1)/GDI
HMU = RMU/H
C
C INTERPOLATE FOR THE SOLUTION -- YOUT
C AND FOR THE DERIVATIVE OF THE SOLUTION -- YPOUT
C
YOUT = 0.0
YPOUT = 0.0
DO J = 1,KOLD
I = KP2 - J
GDIF = G(I) - G(I-1)
TEMP2 = (GTEMP(I) - GTEMP(I-1)) - SIGMA*GDIF
TEMP3 = (C(I) - C(I-1)) + RMU*GDIF
YOUT = YOUT + TEMP2*PHI(:,I)
YPOUT = YPOUT + TEMP3*PHI(:,I)
END DO
YOUT = ((1.0 - SIGMA)*P + SIGMA*Y) +
& H*(YOUT + (GTEMP(1) - SIGMA*G(1))*PHI(:,1))
YPOUT = HMU*(P - Y) + (YPOUT + (C(1) + RMU*G(1))*PHI(:,1))
C
RETURN
END SUBROUTINE SINTRP