C#######################################################################
C
C  Routines to compute cubic-splines for an array of data and to
C  evaluate the spline.
C  Note: Instead of defining a parameter NMAX, the array 'u'
C        appears in the variable list of the subroutine 'splinenm'.
C
C#######################################################################


C***********************************************************************
C Given arrays x(1:n) and y(1:n) containing a tabulated function,
C i.e., y_i=f(x_i), with x_1<x_2<...<x_n, and given values yp1 and
C ypn for the first derivative of the interpolating function at
C points 1 and n, respectively, this routine returns an array
C y2(1:n) of length n which contains the second derivatives of the
C interpolating function at the tabulated points x_i. If yp1 and/or ypn
C are equal to 1 times 10^30 or larger, the routine is signaled to
C set the corresponding boundary condition for a natural spline, with
C zero second derivative on that boundary.x
C***********************************************************************
      
      subroutine splinenm(x,y,n,yp1,ypn,y2,u)

      integer n
      real*8  yp1,ypn,x(n),y(n),y2(n),u(n)
      integer i,k
      real*8  p,qn,sig,un

      if (yp1.gt.0.99d30) then
         y2(1)=0d0
         u(1)=0d0
      else
         y2(1)=-0.5d0
         u(1)=(3d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif

      do i=2,n-1
         sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
         p=sig*y2(i-1)+2d0
         y2(i)=(sig-1d0)/p
         u(i)=(6d0*((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
      enddo

      if (ypn.gt.0.99d30) then
         qn=0d0
         un=0d0
      else
         qn=0.5d0
         un=(3d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif

      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1d0)

      do k=n-1,1,-1
         y2(k)=y2(k)*y2(k+1)+u(k)
      enddo

      return
      end


C***********************************************************************
C Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a 
C function (with the xa_i in order), and given the array y2a(1:n),
C which is the output from the subroutine 'splinenm' above, and given
C a value of x, this routine returns a cubic-spline interpolated
C value of y.
C***********************************************************************

      subroutine splintnm(xa,ya,y2a,n,x,y)

      integer n
      real*8 x,y,xa(n),y2a(n),ya(n)
      integer k,khi,klo
      real*8 a,b,h

      klo=1
      khi=n

 1    if (khi-klo.gt.1) then
         k=(khi+klo)/2
         if (xa(k).gt.x) then
            khi=k
         else
            klo=k
         endif
         goto 1
      endif

      h=xa(khi)-xa(klo)
      if (h.eq.0d0) then
         write(6,'(/a)') 'splintnm: Bad xa input'
         write(6,'(a,e15.6)') 'xa(khi) = xa(klo) =', xa(khi)
         stop
      endif
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+
     +     ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6d0

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     splin2nm: two-dimensional interpolation with cubic splines
c
c     routine corresponds to splintnm in the 1d case, i.e. it 
c     estimates the function value y for a given point (x1,x2) after
c     an array y2a of coefficients has been calculated with another routine.
c
c     input:
c        x1a   x values for the first dimension
c        x2a   x values for the second dimension
c        ya    2d array with tabulated function
c        y2a   2d array with spline coefficients
c        m     number of points in first dimension
c        n     number of points in second dimension
c        y1p1  first derivative at first point in first dimension
c        y1pm  first derivative at mth point in first dimension
c        x1,x2 point for which function value is to be estimated
c     output:
c        y     estimated function value
c     temporary memory:
c        y2tmp vector of length nt
c        ytmp  dito
c        yytmp dito
c        nt    length of temporary memory; nt>=max(m,n)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      
      subroutine splin2nm(x1a,x2a,ya,y2a,m,n,y1p1,y1pm,
     &     x1,x2,y,y2tmp,ytmp,yytmp,nt)

      implicit none

      integer m,n,nt
      real*8 x1,x2,y,x1a(m),x2a(n),y2a(m,n),ya(m,n)
      real*8 y1p1,y1pm,y2tmp(nt),ytmp(nt),yytmp(nt)
      integer j,k

      if(nt.lt.m.or.nt.lt.n) then
         write(6,'(/a)') 'splin2nm: temporary arrays too short'
         stop
      endif
      do j=1,m
         do k=1,n
            ytmp(k)=ya(j,k)
            y2tmp(k)=y2a(j,k)
         enddo
         call splintnm(x2a,ytmp,y2tmp,n,x2,yytmp(j))
      enddo
      call splinenm(x1a,yytmp,m,y1p1,y1pm,y2tmp,ytmp)
      call splintnm(x1a,yytmp,y2tmp,m,x1,y)
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c   spline2: calculates the auxiliary y2a array for splin2nm.
c
c   input: x2a(n)  - x-array for second dimension
c          ya(m,n) - 2D y-array
c          m,n - dimensions
c   output:y2a(m,n) - second derivative array
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
      subroutine spline2(x2a,ya,m,n,y2a,ytmp,y2tmp,u,nt)
 
      implicit none
 
      integer  m,n,nt,j,k
      real*8   x2a(n),ya(m,n),y2a(m,n),ytmp(nt),y2tmp(nt),u(nt)
 
 
      do j=1,m
         do k=1,n
            ytmp(k)=ya(j,k)
         enddo
         call splinenm(x2a,ytmp,n,1.d30,1.d30,y2tmp,u)
         do k=1,n
            y2a(j,k)=y2tmp(k)
         enddo
      enddo
 
      return
      end
