c   Polynomial interpolation routines from Numerical Recipes
c   For use with dfrmint.f, frmint.f

c   Ken Van Riper, kvr@lanl.gov

      subroutine polynomial_int (xa, ya, n, no, x, y, dy)

c                          Taken from Numerical Recipes
      integer nmax, n, no
      parameter (nmax=7)

      real*8 xa(n), ya(n), c(nmax), d(nmax)
      real*8 x, y, dy, dif, dift, w, den, ho, hp
      integer ns, i, m

      ns = 1
      dif = abs(x - xa(1))
      do i = 1,no
          dift = abs(x - xa(i))
          if (dift .lt. dif) then
              ns = i
              dif = dift
          endif
          c(i) = ya(i)
          d(i) = ya(i)
      enddo
      y = ya(ns)
      ns = ns - 1
      do m = 1,no-1
          do i = 1,no-m
              ho = xa(i) - x
              hp = xa(i+m) - x
              w = c(i+1) - d(i)
              den = ho - hp
              if (den .eq. 0) then
                  write (6,'(a)') 'Trouble in polynomial_int'
                  stop
              endif
              den = w/den
              d(i) = hp*den
              c(i) = ho*den
          enddo
          if (2*ns .lt. no-m) then
              dy = c(ns+1)
          else
              dy = d(ns)
              ns = ns -1
          endif
          y = y + dy
      enddo
      return
      end

      subroutine TwoDpolyint (etas, betas, ya, m, n, mo, no,
     $                         eta, beta, y, dy)

c                          Taken from Numerical Recipes
c  M for eta direction, N for beta direction

      integer nmax, mmax, m, n, mo, no
      parameter (nmax=7, mmax=7)
      real*8 etas(m), betas(n), ya(m,n), yntmp(nmax), ymtmp(mmax),
     $       eta, beta, dy, y
      integer j, k

      do j = 1,mo
          do k = 1,no
              yntmp(k) = ya(j,k)
          enddo
          call polynomial_int (betas, yntmp, n, no, beta, ymtmp(j), dy)
      enddo
      call polynomial_int (etas, ymtmp, m, mo, eta, y, dy)
      return
      end