!!$---------------------------------------------------------------------
SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y)
!!$---------------------------------------------------------------------
!!$  Given the arrays xa(1:n) and ya(1:n) containing a tabulated function,
!!$  ya(i)=f(xa(i)), with x(1)<...<x(i)<...<x(n), and y2a(1:n), which is 
!!$  the output from spline, and a value x, this routine returns a cubic-
!!$  spline interpolated value y.
!!$   [Adapted from Numerical Recipes in Fortran, S3.3, p110]
!!$---------------------------------------------------------------------
  USE PRCSN
!!$
  IMPLICIT NONE

  INTEGER, INTENT(IN) :: N
  REAL(PRC),INTENT(IN), DIMENSION(N) :: XA,YA,Y2A
  REAL(PRC), INTENT(IN) :: X
  REAL(PRC), INTENT(OUT) :: Y

  INTEGER :: k,khi,klo
  REAL(PRC) :: a,b,h

  KLO = 1
  KHI = N

  DO
     IF (KHI-KLO <= 1) EXIT
     K = (KHI+KLO)/2
     IF(XA(K) > X)THEN
        KHI = K
     ELSE
        KLO = K
     END IF
  END DO

  H=XA(KHI)-XA(KLO)
  IF (H == 0.0_PRC) THEN
     WRITE(*,'(" Bad XA input in splint.f90.  Quiting... ")')
     STOP
  END IF
  A = (XA(KHI)-X)/H
  B = (X-XA(KLO))/H
  Y = A*YA(KLO)+B*YA(KHI) &
       + ((A**3.0_PRC - A)*Y2A(KLO)+(B**3.0_PRC - B)*Y2A(KHI)) &
       * (H**2.0_PRC)/6.0_PRC
  RETURN
END SUBROUTINE SPLINT
