!+ Subroutine to evaluate a cubic spline.
!
! Method:
!	Straightforward.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date			Comment
!	1.0		12-04-95		First Version under RCS
!						(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77  with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      SUBROUTINE SPLINE_EVALUATE(IERR, N_DATA, X, Y, Y2
     &   , X_POINT, Y_POINT)
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/error.pcf'
!
!     DUMMY ARGUMENTS.
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER	!, INTENT(IN)
     &     N_DATA
!             NUMBER OF DATA POINTS
      REAL*8	!, INTENT(IN)
     &     X(N_DATA)
!             ABSCISSAE OF DATA
     &   , Y(N_DATA)
!             ORDINATES OF DATA
     &   , Y2(N_DATA)
!             SECOND DERIVATIVE
     &   , X_POINT
!             POINT FOR EVALUATION
      REAL*8	!, INTENT(OUT)
     &     Y_POINT
!             VALUE OF SPLINE
!
!     LOCAL ARGUMENTS.
      INTEGER
     &     I_LOW
!             LOWER LIMIT OF BRACKETING PAIR
     &   , I_HIGH
!             UPPER LIMIT OF BRACKETING PAIR
     &   , I
!             INDEX OF NEW POINT
      REAL*8
     &     DELTA
!             WIDTH OF INTERVAL
     &   , A
!             WEIGHT OF LOWER POINT
     &   , B
!             WEIGHT OF UPPER POINT
     &   , C
!             TEMPORARY STORE
!
!
      IF (N_DATA.EQ.1) THEN
!        RETURN THE ONLY AVAILABLE VALUE.
         Y_POINT=Y(1)
      ELSE IF ( (X_POINT.GT.X(N_DATA)).OR.(X_POINT.LT.X(1)) ) THEN
!        OUTSIDE THE RANGE: RETURN AN ERROR FLAG FOR CORRECTIVE ACTION.
         IERR=I_ERR_RANGE
         RETURN
      ELSE
!        FIND THE BRACKETING POINTS BY A BINARY CHOP.
         I_LOW=1
         I_HIGH=N_DATA
1        IF ((I_HIGH-I_LOW).GT.1) THEN
            I=(I_LOW+I_HIGH)/2
            IF (X_POINT.GT.X(I)) THEN
               I_LOW=I
            ELSE
               I_HIGH=I
            ENDIF
            GOTO 1
         ENDIF
         DELTA=X(I_HIGH)-X(I_LOW)
         A=(X(I_HIGH)-X_POINT)/DELTA
         B=(X_POINT-X(I_LOW))/DELTA
         C=(DELTA**2/6.0D+00)
         Y_POINT=A*(Y(I_LOW)+(A*A-1)*Y2(I_LOW)*C)
     &      +B*(Y(I_HIGH)+(B*B-1.0D+00)*Y2(I_HIGH)*C)
      ENDIF
!
!
      RETURN
      END
