subroutine splint(x,y,y2,n,xx,yy,dydx)
c
c=======================================================================
c                                                                    ===
c  Given arrays X and Y of length N, which tabulate a function,      ===
c  Y=f(X), with the  X  in ascending order, and given the array      ===
c  Y2 which contains the second derivative of the interpolating      ===
c  function  at the  tabulated points X as computed  by routine      ===
c  SPLINE, and given a value  XX, this routine returns a cubic-      ===
c  spline interpolated value YY.  For a value of XX outside the      ===
c  range of X, constant extrapolation based on the nearest Y         ===
c  value is used.                                                    ===
c                                                                    ===
c  Reference :                                                       ===
c                                                                    ===
c  Press, W.H, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,   ===
c         1986: Numerical Recipes, the art of scientific computing.  ===
c         Cambridge University Press.                                ===
c                                                                    ===
c  Modified by H.G. Arango (1989) to output the first derivative     ===
c  DYDX at a given value XX.                                         ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <pconst.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer k,khi,klo,n
      FLOAT
     *      a,b,c,d,dydx,e,f,h,xx,yy
      FLOAT
     *      x(n),y(n),y2(n)
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c  Solve simple extrapolation cases first.
c-----------------------------------------------------------------------
c
      if (xx.le.x(1)) then
         yy   = y(1)
         dydx = c0
         return
      end if
c
      if (xx.ge.x(n)) then
         yy   = y(n)
         dydx = c0
         return
      end if
c
c-----------------------------------------------------------------------
c  Find the right place of XX in the table by means of bisection.
c-----------------------------------------------------------------------
c
      klo=1
      khi=n
  10  if((khi-klo).gt.1) then
        k=(khi+klo)/2
        if(x(k).gt.xx) then
          khi=k
        else
          klo=k
        endif
        goto 10
      endif
c
c-----------------------------------------------------------------------
c  KLO and KHI now bracket the input value XX.
c-----------------------------------------------------------------------
c
      h=x(khi)-x(klo)
      if(h.eq.c0) then
        print *, ' SPLINT: bad X input, they must be distinct.'
        call exitus('SPLINT')
      endif
c
c-----------------------------------------------------------------------
c  Evaluate cubic spline polynomial.
c-----------------------------------------------------------------------
c
      a=(x(khi)-xx)/h
      b=(xx-x(klo))/h
      c=(a*a*a-a)*(h*h)*r6
      d=(b*b*b-b)*(h*h)*r6
      e=(c3*(a*a)-c1)*h*r6
      f=(c3*(b*b)-c1)*h*r6
      yy=a*y(klo)+b*y(khi)+c*y2(klo)+d*y2(khi)
      dydx=(y(khi)-y(klo))/h-e*y2(klo)+f*y2(khi)
c
      return
      end