C   
C     SPLINE INTERPOLATION PROGRAM FROM DOMINION ASTROPHYSICAL OBSERVATORY  
C     PUBLICATIONS, VOLUME XVI, NO.6
C   
      SUBROUTINE INTEP(XP,P,X,F,N,IER)  
      REAL LP1,LP2,L1,L2
      DIMENSION X(*),F(*)   
      IER=1 
      IO=1  
      IUP=0 
      IF(X(2).LT.X(1))IUP=1 
      N1=N-1
      IF((XP.GE.X(N).AND.IUP.EQ.0).OR.(XP.LE.X(N).AND.IUP.EQ.1))THEN
    5 P=F(N)
      GOTO 6
      ELSE IF((XP.LE.X(1).AND.IUP.EQ.0).OR. 
     * (XP.GE.X(1).AND.IUP.EQ.1))THEN   
      P=F(1)
    6 IER=2 
      RETURN
      END IF
      ENTRY EINTEP(XP,P,X,F,N,IER)  
      IER=1 
    8 DO 1 I=IO,N   
       IF(XP.LT.X(I).AND.IUP.EQ.0)GOTO 2
       IF(XP.GT.X(I).AND.IUP.EQ.1)GOTO 2
    1 CONTINUE  
      P=F(N)
      IER=2
      RETURN

    2 I=I-1 
      IF(I.EQ.IO-1)GOTO 4   
      IO=I+1
      LP1=1./(X(I)-X(I+1))  
      LP2=1./(X(I+1)-X(I))  
      IF(I.EQ.1)FP1=(F(2)-F(1))/(X(2)-X(1)) 
      IF(I.EQ.1)GOTO 3  
      FP1=(F(I+1)-F(I-1))/(X(I+1)-X(I-1))   
    3 IF(I.GE.N1)FP2=(F(N)-F(N-1))/(X(N)-X(N-1))
      IF(I.GE.N1)GOTO 4 
      FP2=(F(I+2)-F(I))/(X(I+2)-X(I))   
    4 XPI1=XP-X(I+1)
      XPI=XP-X(I)   
      L1=XPI1*LP1   
      L2=XPI*LP2
      P=F(I)*(1.-2.*LP1*XPI)*L1*L1+F(I+1)*(1.-2.*LP2*XPI1)  
     1*L2*L2+FP2*XPI1*L2*L2+FP1*XPI*L1*L1   
      RETURN
      END   
