!+ Subroutine to evaluate the integral of a cubic spline.
!
! Method:
!	The range of integration is checked against the range of the
!	spline. If the range of integration lies within the range of 
!	the spline the integral is evaluated noting that between each
!	pair of knots the spline is a cubic which can be integrated
!	directly.
!
! 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 INTEGRATE_SPLINE(IERR, X1, X2
     &   , N_DATA, X, Y, Y2, INTEGRAL)
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/error.pcf'
      INCLUDE 'incl/def_std_io.icf'
!
!     DUMMY ARGUMENTS.
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER	!, INTENT(IN)
     &     N_DATA
!             NUMBER OF DATA POINTS
      REAL*8	!, INTENT(IN)
     &     X1
!             LOWER LIMIT OF INTEGRATION
     &   , X2
!             UPPER LIMIT OF INTEGRATION
     &   , X(N_DATA)
!             SPLINE KNOTTING POINTS
     &   , Y(N_DATA)
!             VALUES OF SPLINE
     &   , Y2(N_DATA)
!             SECOND DERIVATIVE
      REAL*8	!, INTENT(OUT)
     &     INTEGRAL
!             EVALUATED INTEGRAL
!
!     LOCAL VARIABLES.
      INTEGER
     &     I_LOW
!             LOWER KNOT POINT
     &   , I_HIGH
!             UPPER KNOT POINT
     &   , I_MID
!             BISECTING POINT
      REAL*8
     &     XL
!             LOWER LIMIT OF INTEGRAL
     &   , XH
!             UPPER LIMIT OF INTEGRAL
     &   , X_BAR
!             MID-POINT OF RANGE OF INTEGRATION
     &   , DELTA
!             WIDTH OF INTERVAL
     &   , DELTA_1
!             WIDTH OF INTEGRATION
     &   , IA
!             INTEGRAL OF COEFFICIENT
     &   , IB
!             INTEGRAL OF COEFFICIENT
     &   , I3A
!             INTEGRAL OF COEFFICIENT
     &   , I3B
!             INTEGRAL OF COEFFICIENT
!
!
      INTEGRAL=0.0D+00
!
!     AN INTEGRAL WITH ONE DATA POINT CANNOT BE PERFORMED.
      IF (N_DATA.EQ.1) THEN
         WRITE(IU_ERR, '(/A)') 
     &      '*** ERROR: SPLINE HAS ONLY ONE KNOT POINT.'
         IERR=I_ABORT_CALCULATION
         RETURN
      ENDIF
!
!     CHECK THAT THE INTEGRAND LIES WITHIN THE DOMAIN OF THE SPLINE
!     AND FIND THE FIRST INTERVAL.
      IF ( (X1.GT.X(N_DATA)).OR.(X2.LT.X(1)) ) THEN
!        DOMAIN OF INTEGRAND LIES OUTSIDE RANGE OF SPLINE.
         IERR=I_ERR_RANGE
         RETURN
      ELSE IF (X1.LE.X(1)) THEN
         XL=X(1)
         XH=MIN(X2, X(2))
         I_LOW=1
      ELSE
!        FIND THE INTERVAL IN WHICH X1 LIES.
         I_LOW=1
         I_HIGH=N_DATA
1           I_MID=(I_LOW+I_HIGH)/2
            IF (X1.GT.X(I_MID)) THEN
               I_LOW=I_MID
            ELSE
               I_HIGH=I_MID
            ENDIF
            IF ( (I_HIGH-I_LOW).GT.1) GOTO 1
!
         XL=X1
         XH=MIN(X(I_HIGH), X2)
      ENDIF
      I_HIGH=I_LOW+1
!
2     DELTA=(X(I_HIGH)-X(I_LOW))
      DELTA_1=XH-XL
      X_BAR=0.5D+00*(XL+XH)
      IA=(DELTA_1/DELTA)*(X(I_HIGH)-X_BAR)
      IB=(DELTA_1/DELTA)*(X_BAR-X(I_LOW))
      I3A=0.25D+00*((X(I_HIGH)-XL)**4-(X(I_HIGH)-XH)**4)/DELTA**3
      I3B=0.25D+00*((X(I_LOW)-XH)**4-(X(I_LOW)-XL)**4)/DELTA**3
      INTEGRAL=INTEGRAL+Y(I_LOW)*IA+Y(I_HIGH)*IB
     &   +Y2(I_LOW)*DELTA**2*(I3A-IA)/6.0D+00
     &   +Y2(I_HIGH)*DELTA**2*(I3B-IB)/6.0D+00
!     DETERMINE NEXT INTERVAL.
      IF (XH.LT.X2) THEN
         I_LOW=I_LOW+1
         I_HIGH=I_LOW+1
         IF (I_HIGH.GT.N_DATA) THEN
!           RANGE OF SPLINE IS EXHAUSTED, SO STOP ADDING TERMS.
            RETURN
         ENDIF
         XL=X(I_LOW)
         XH=MIN(X2, X(I_HIGH))
         GOTO 2
      ENDIF
!
!
      RETURN
      END
