!+ Subroutine to set coefficients for a spline fit.
!
! Method:
!	The second derivatives of the cubic spline are found by solving
!	a matrix equation using natural boundary conditions.
!
! 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_FIT(N_DATA, X, Y, Y2, U)
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/error.pcf'
!
!     DUMMY ARGUMENTS
      INTEGER	!, INTENT(IN)
     &     N_DATA
!             NUMBER OF DATA POINTS
      REAL*8	!, INTENT(IN)
     &     X(N_DATA)
!             ABSICISSAE OF DATA
     &   , Y(N_DATA)
!             ORDINATES OF DATA
      REAL*8	!, INTENT(OUT)
     &     Y2(N_DATA)
!             CALCULATED SECOND DERIVATIVES
     &   , U(N_DATA)
!             TEMPORARY WORKING VARIABLE
!
!     LOCAL ARGUMENTS.
      INTEGER
     &     I
!             LOOP VARIABLE
      REAL*8
     &     SIGMA
!             TEMPORARY VARIABLE
     &   , RHO
!             "
!
!
!     NO FIT IS POSSIBLE WITH ONE TERM.
      IF (N_DATA.EQ.1) RETURN
!
!     THIS ROUTINE IS ESSENTIALLY TRIDIAGONAL SOLUTION OF A MATRIX
!     EQUATION. THE STARTING POINT IS THE NATURAL BOUNDARY CONDITION.
      Y2(1)=0.0D+00
      U(1)=0.0D+00
      DO I=2, N_DATA-1
         SIGMA=(X(I)-X(I-1))/(X(I+1)-X(I-1))
         RHO=SIGMA*Y2(I-1)+2.0D+00
         Y2(I)=(SIGMA-1.0D+00)/RHO
         U(I)=(6.0D+00*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))
     &      /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIGMA*U(I-1))/RHO
      ENDDO
      Y2(N_DATA)=0.0D+00
!
!     BACKSUBSTITUTION.
      DO I=N_DATA-1, 1, -1
         Y2(I)=Y2(I)*Y2(I+1)+U(I)
      ENDDO
!
!
      RETURN
      END
