!!$---------------------------------------------------------------------
SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2)
!!$---------------------------------------------------------------------
!*$* NO CONSURRENTIZE
!!$---------------------------------------------------------------------
!!$  Given arrays x(1:n) and y(1:n) containing a tabulated function,
!!$  y(i)=f(x(i)), with x(1)<...<x(i)<...<x(n), and values yp1 and ypn
!!$  for the first derivative of the interpolating function at points 
!!$  1 and n respectively, this routine returns an array y2(1:n) which
!!$  contains the second derivatives of the interpolating function at
!!$  the tabulated points x(i).  If yp1 and/or ypn are equal to 1.0e+30
!!$  or larger the algorithm follws natural spline in which the second
!!$  derivative at the boundaries are set to zero.
!!$   [Adapted from Numerical Recipes in Fortran, S3.3, p109]
!!$---------------------------------------------------------------------
  USE PRCSN
!!$
  IMPLICIT NONE
  REAL(PRC), PARAMETER :: MAXVAL = 1.0E+30_PRC

  INTEGER, INTENT(IN) :: N
  REAL(PRC), INTENT(IN), DIMENSION(N) :: X,Y
  REAL(PRC), INTENT(IN) :: yp1,ypn
  REAL(PRC), INTENT(OUT), DIMENSION(N) :: y2

  INTEGER :: I,K
  REAL(PRC), DIMENSION(N) :: U
  REAL(PRC) :: SIG,P,QN,UN

!!$  do i=1,n
!!$     write(*,*) i,x(i),y(i)
!!$  end do
  IF (YP1 >= MAXVAL) THEN
     Y2(1)=0.0_PRC
     U(1)=0.0_PRC
  ELSE
     Y2(1)=-0.5_PRC
     U(1)=(3.0_PRC/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
  END IF
  DO I=2,N-1
     SIG   = (X(I)-X(I-1))/(X(I+1)-X(I-1))
     P     = SIG*Y2(I-1)+2.0_PRC
     Y2(I) = (SIG-1.0_PRC)/P
     U(I)  = (6.0_PRC*((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))-SIG*U(I-1))/P
  END DO
  IF (YPN >= MAXVAL) THEN
     QN=0.0_PRC
     UN=0.0_PRC
  ELSE
     QN=0.5_PRC
     UN=(3.0_PRC/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
  END IF
  Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.0_PRC)
  DO K=N-1,1,-1
     Y2(K)=Y2(K)*Y2(K+1)+U(K)
  END DO
  RETURN
END SUBROUTINE SPLINE
