      subroutine splint(xa,ya,y2a,n,x,y)
*
* given the input arrays xa & ya (y=f(x)) of length n, tabulating
* a function, and given the spline coefficients y2a (from spline)
* and given an evaluation point x, this routine returns a cubic-spline
* interpolated value y.
* Routine is from Numerical Recipies
*
      implicit double precision(a-h,o-z)

      dimension xa(n),ya(n),y2a(n)
*
* find place in table
*
      call hunt(xa,n,x,klo)
      khi = klo + 1
      h=xa(khi)-xa(klo)
      if (h.eq.0.)then
*
* xa points are not distinct
*
        stop
      endif
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
*
* evaluate spline
*
      y=a*ya(klo)+b*ya(khi)+
     1      ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.

      return
      end

************************************************************************

