c2--567--1---------2---------3---------4---------5---------6---------712
      Subroutine spline(x0, y0, n0, x, y, n, lam1, lamn0)
c
c     -- Spline/Spline interpolation/N. Nath/April '01.  (Adapted from 
c        a previous version, this routine has the following features: 
c        (1) it incorporates different types of end conditions; (2) the 
c        argument array {x0} may be an increasing or a decreasing 
c        sequence; (3) {x} need not be ordered; (4) for x lying outside
c        {x0}, the cubic in the nearest interval is used; this feature 
c        should be used with extreme caution.)
c
c     -- Function:  Performs a cubic spline interpolation for function 
c        values (y) at an array of arguments (x) of dimension n, with 
c        different end parameters lam1 & lamn0.  x0 & y0 (of dimension 
c        n0) are the input argument & associated function array.  lam1 
c        [alternately lamn0] may take the following values:  0.0, -1.0,
c        -2.0, 1.0, which correspond respectively to y20(1) = (0.0, 0.5,
c        1.0)*y20(2), & y10(1) = 0.0 [y20(n0) = (0.0, 0.5, 1.0)*
c        y20(n0-1), & y10(n0) = 0.0]; y20 and y10 are respectively the
c        second and first derivative at x0.
c
c     -- Usage Notes:  (1) Choice of lam1 and lamn0 depends on the 
c        knowledge of the function.  If nothing is known, use one of the
c        first two values.  Note also that the second & third values 
c        correspond respectively to a single cubic in the domain: 
c        (2.0*x0(1) - x0(2), x0(2)), with y20 at the left boundary = 
c        0.0, & a quadratic in the domain (x0(1), x0(2)), and similarly
c        for the other end.  (2) x0 must be monotonically increasing or 
c        decreasing.  (3) n0 must be > 2; if it exceeds 256, change the
c        value of n0max.
c
            implicit none
c
c     -- Input & Output Variables:
c
            real        x0(n0), y0(n0),
                              !base argument & function array (input),
     1                  x(n), y(n),
                              !argument array at which function values
                              !  are desired (input), & interpolated
                              !  function array (output)
     2                  lam1, lamn0
                              !end parameters at x0(1) & x0(n0) (input)

            integer     n0, n !dimensions of (x0, y0) & (x, y) (input)
c
c     -- Local Variables:
c
            integer     n0max
            parameter   (n0max = 256)
c
            real        y20(n0max) !second derivative array

            real        h(n0max), a(n0max), b(n0max), d(n0max),
     1                  phi, lam, dy0, dy00, c, c1, f

            integer     n01, l, u, m, i, j, k, i1

            logical     reverse                
c
c     -- Subprogram Called:  order
c
c     -- Ref.:  Ahlberg et al, Theory Of Splines & Their Applications,
c        Academic Press, 1967; Hamming, Numerical Methods, Mcgraw-Hill,
c        1973.
c
      phi(lam) = lam*(lam + 1.0)*(lam + 2.0)
c
c           Order x0 in an increasing sequence; change x accordingly
c                 !note:  this ordering, viz., x0 -> -x0 while y0 -> y0, 
c                 !  only changes the signs of the odd derivatives; thus
c                 !  y20 remains unchanged
c
      reverse =  (x0(1) .gt. x0(n0))
      if (reverse) call order(x0, n0, x, n)
c
c           Determine y20 by solving tri-diagonal linear equations
c
c           o  Define h, a, & d (quantities independent of y0)
c
      n01 = n0 - 1
      do i = 1,n01
         h(i) = x0(i+1) - x0(i)
      end do
      a(1) = -lam1/2.0
      do i = 2,n01
         i1 = i - 1
         d(i) = h(i1)*a(i1) + 2.0*(h(i1) + h(i))
         a(i) = -h(i)/d(i)
      end do
      d(n0) = h(n01)*(lamn0*a(n01) + 2.0)
c
c           o  Define b, & complete computations
c
      dy0 = (y0(2) - y0(1))/h(1)
      b(1) = (phi(lam1)*dy0)/h(1)/2.0
      do i = 2,n01
         dy00 = dy0
         dy0 = (y0(i+1) - y0(i))/h(i)
         b(i) = (6.0*(dy0 - dy00) - h(i-1)*b(i-1))/d(i)
      end do
      b(n0) = (-phi(lamn0)*dy0 - lamn0*h(n01)*b(n01))/d(n0)
      y20(n0) = b(n0)
      k = n0
      do i = 2,n0
         k = k - 1
         y20(k) = a(k)*y20(k+1) + b(k)
      end do
c
c           Interpolate using binary domain search
c
      do j = 1,n
         l = 1
         u = n0
         m = (l + u)/2
         do while (m .gt. l)
            if (x(j) .lt. x0(m)) then
                  u = m
               else
                  l = m
            end if
            m = (l + u)/2
         end do
         c = (x0(l+1) - x(j))/h(l)
         c1 = 1.0 - c
         f = (h(l)**2)/6.0
         y(j) = c*(y0(l) - f*y20(l)*(1.0 - c**2))
     1            + c1*(y0(l+1) - f*y20(l+1)*(1.0 - c1**2))
      end do
c
c           Reorder x0 and x
c
      if (reverse) call order(x0, n0, x, n)
c
      return
      end

c     ---
      Subroutine order(x0, n0, x, n)
c
c     -- Order/Adjunct to Spline/N. Nath/April '01
c
c     -- Input & Output Variables:
c
            real        x0(n0), x(n)

            integer     n0, n
c
c     -- Local Variables:
c
            integer     i, j     
c
c           Change signs of x0 and x
c
      do i = 1,n0
         x0(i) = -x0(i)
      end do
      do j = 1,n
         x(j) = -x(j)
      end do
c
      return
      end