i = 10-ib
      t = 0.0D+00
      do j = i+1, 9
        t = t + b(j,i)*a(j,k)
      end do
      a(i,k) = (b(10,i)-t)/b(i,i)
    end do
!
!  Scale the coefficients to adjust for the column scaling.
!
    a(1:6,k) = a(1:6,k) / avsq
    a(7,k) = a(7,k) / av
    a(8,k) = a(8,k) / av
    a(9,k) = a(9,k) / av
!
!  Unmark K and the elements of NPTS.
!
    lnext(k) = -lnext(k)
    do i = 1, lnp
      np = npts(i)
      lnext(np) = -lnext(np)
    end do

  end do
!
!  No errors encountered.
!
  xyzmin(1:3) = xyzmn(1:3)
  xyzdel(1:3) = xyzdl(1:3)
  rmax = sqrt ( rsmx )
  ier = 0
  return
!
!  N, NQ, NW, or NR is out of range.
!
20 continue

  ier = 1
  return
!
!  Duplicate nodes were encountered by GETNP3.
!
21 ier = 2
  return
end
function qs3val ( px, py, pz, n, x, y, z, f, nr, lcell, lnext, xyzmin, &
  xyzdel, rmax, rsq, a )

!***********************************************************************
!
!! QS3VAL evaluates the interpolant function Q(X,Y,Z) created by QSHEP3.
!
!  Discussion:
!
!    This function returns the value Q(PX,PY,PZ) where Q is
!    the weighted sum of quadratic nodal functions defined in
!    subroutine QSHEP3.  QS3GRD may be called to compute a
!    gradient of Q along with the value, or to test for errors.
!
!    This function should not be called if a nonzero error flag was
!    returned by QSHEP3.
!
!  Author:
!
!    Robert Renka,
!    University of North Texas,
!    (817) 565-2767.
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 661: QSHEP3D, Quadratic Shepard method for trivariate
!    interpolation of scattered data,
!    ACM Transactions on Mathematical Software,
!    Volume 14, 1988, pages 151-152.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) PX, PY, PZ, the point P at which Q is 
!    to be evaluated.
!
!    Input, integer N, the number of nodes and data values defining Q.
!    10 <= N.
!
!    Input, real ( kind = 8 ) X(N), Y(N), Z(N), F(N), the node coordinates
!    and data values interpolated by Q.
!
!    Input, integer NR, the number of rows, columns and planes in the cell
!    grid.  Refer to STORE3.  1 <= NR.
!
!    Input, integer LCELL(NR,NR,NR), nodal indices associated with cells.  
!    Refer to STORE3.
!
!    Input, integer LNEXT(N), the next-node indices.  Refer to STORE3.
!
!    Input, real ( kind = 8 ) XYZMIN(3), XYZDEL(3), the minimum nodal
!    coordinates and cell dimensions, respectively.  XYZDEL elements 
!    must be positive.  Refer to STORE3.
!
!    Input, real ( kind = 8 ) RMAX, the square root of the largest 
!    element in RSQ, the maximum radius.
!
!    Input, real ( kind = 8 ) RSQ(N), the squared radii which enter 
!    into the weights defining Q.
!
!    Input, real ( kind = 8 ) A(9,N), the coefficients for the nodal 
!    functions defining Q.
!
!    Output, real ( kind = 8 ) QS3VAL, the function value Q(PX,PY,PZ)
!    unless N, NR, XYZDEL, or RMAX is invalid, in which case the 
!    value 0 is returned.
!
  implicit none

  integer n
  integer nr

  real ( kind = 8 ) a(9,n)
  real ( kind = 8 ) delx
  real ( kind = 8 ) dely
  real ( kind = 8 ) delz
  real ( kind = 8 ) dxsq
  real ( kind = 8 ) dysq
  real ( kind = 8 ) dzsq
  real ( kind = 8 ) ds
  real ( kind = 8 ) dx
  real ( kind = 8 ) dy
  real ( kind = 8 ) dz
  real ( kind = 8 ) f(n)
  integer i
  integer imax
  integer imin
  integer j
  integer jmax
  integer jmin
  integer k
  integer kmax
  integer kmin
  integer l
  integer lcell(nr,nr,nr)
  integer lnext(n)
  integer lp
  real ( kind = 8 ) px
  real ( kind = 8 ) py
  real ( kind = 8 ) pz
  real ( kind = 8 ) qs3val
  real ( kind = 8 ) rd
  real ( kind = 8 ) rds
  real ( kind = 8 ) rmax
  real ( kind = 8 ) rs
  real ( kind = 8 ) rsq(n)
  real ( kind = 8 ) sw
  real ( kind = 8 ) swq
  real ( kind = 8 ) w
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) xp
  real ( kind = 8 ) xyzdel(3)
  real ( kind = 8 ) xyzmin(3)
  real ( kind = 8 ) y(n)
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin
  real ( kind = 8 ) yp
  real ( kind = 8 ) z(n)
  real ( kind = 8 ) zmax
  real ( kind = 8 ) zmin
  real ( kind = 8 ) zp

  xp = px
  yp = py
  zp = pz
  xmin = xyzmin(1)
  ymin = xyzmin(2)
  zmin = xyzmin(3)
  dx = xyzdel(1)
  dy = xyzdel(2)
  dz = xyzdel(3)

  if ( n < 10 ) then
    qs3val = 0.0D+00
    return
  end if

  if ( nr < 1  .or.  dx <= 0.0 &
         .or.  dy <= 0.0  .or.  dz <= 0.0  .or. &
         rmax < 0.0 ) then
    qs3val = 0.0D+00
    return
  end if
!
!  Set IMIN, imax, jmin, jmax, kmin, and kmax to cell indices
!  defining the range of the search for nodes whose radii
!  include P.  The cells which must be searched are those
!  intersected by (or contained in) a sphere of radius rmax
!  centered at P.
!
  imin = int((xp-xmin-rmax)/dx) + 1
  imin = max ( imin, 1 )

  imax = int((xp-xmin+rmax)/dx) + 1
  imax = min ( imax, nr )

  jmin = int((yp-ymin-rmax)/dy) + 1
  jmin = max ( jmin, 1 )

  jmax = int((yp-ymin+rmax)/dy) + 1
  jmax = min ( jmax, nr )

  kmin = int((zp-zmin-rmax)/dz) + 1
  kmin = max ( kmin, 1 )

  kmax = int((zp-zmin+rmax)/dz) + 1
  kmax = min ( kmax, nr )
!
!  Test for no cells within the sphere of radius RMAX.
!
  if ( imax < imax .or. &
       jmax < jmin .or. &
       kmax < kmin ) then
    qs3val = 0.0D+00
    return
  end if
!
!  Accumulate weight values in SW and weighted nodal function
!  values in SWQ.  The weights are w(l) = ((r-d)+/(r*d))**2
!  for r**2 = rsq(l) and d = distance between P and node L.
!
  sw = 0.0D+00
  swq = 0.0D+00
!
!  Outer loop on cells (i,j,k).
!
  do k = kmin, kmax
    do j = jmin, jmax
      do i = imin, imax

        l = lcell(i,j,k)

        if ( l == 0 ) then
          cycle
        end if
!
!  Inner loop on nodes L.
!
        do

          delx = xp - x(l)
          dely = yp - y(l)
          delz = zp - z(l)

          dxsq = delx * delx
          dysq = dely * dely
          dzsq = delz * delz

          ds = dxsq + dysq + dzsq
          rs = rsq(l)

          if ( ds < rs ) then

            if ( ds == 0.0D+00 ) then
              qs3val = f(l)
              return
            end if

            rds = rs * ds
            rd = sqrt ( rds )
            w = ( rs + ds - rd - rd ) / rds
            sw = sw + w

            swq = swq + w *( a(1,l) * dxsq + a(2,l)*delx*dely + &
              a(3,l) * dysq + a(4,l)*delx*delz + &
              a(5,l) * dely*delz + a(6,l)*dzsq + &
              a(7,l) * delx + a(8,l)*dely + &
              a(9,l) * delz + f(l) )

          end if

          lp = l
          l = lnext(lp)

          if ( l == lp ) then
            exit
          end if

        end do

      end do
    end do
  end do
!
!  SW = 0 iff P is not within the radius R(L) for any node L.
!
  if ( sw == 0.0D+00 ) then
    qs3val = 0.0D+00
  else
    qs3val = swq / sw
  end if

  return
end
subroutine qs3grd ( px, py, pz, n, x, y, z, f, nr, lcell, lnext, xyzmin, &
  xyzdel, rmax, rsq, a, q, qx, qy, qz, ier )

!***********************************************************************
!
!! QS3GRD computes the value and gradient of the interpolant function.
!
!  Discussion:
!
!    This subroutine computes the value and gradient at (PX,PY,PZ) of 
!    the interpolatory function Q defined in subroutine QSHEP3.  
!
!    Q(X,Y,Z) is a weighted sum of quadratic nodal functions.
!
!  Author:
!
!    Robert Renka,
!    University of North Texas,
!    (817) 565-2767.
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 661: QSHEP3D, Quadratic Shepard method for trivariate
!    interpolation of scattered data,
!    ACM Transactions on Mathematical Software,
!    Volume 14, 1988, pages 151-152.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) PX, PY, PZ, the point P at which Q and
!    its partials are to be evaluated.
!
!    Input, integer N, the number of nodes and data values defining Q.
!    10 <= N.
!
!    Input, real ( kind = 8 ) X(N), Y(N), Z(N), F(N), the node coordinates and
!    data values interpolated by Q.
!
!    Input, integer NR, the number of rows, columns and planes in the cell
!    grid.  Refer to STORE3.  1 <= NR.
!
!    Input, integer LCELL(NR,NR,NR), nodal indices associated with cells.  
!    Refer to STORE3.
!
!    Input, integer LNEXT(N), the next-node indices.  Refer to STORE3.
!
!    Input, real ( kind = 8 ) XYZMIN(3), XYZDEL(3), the minimum nodal
!    coordinates and cell dimensions, respectively.  XYZDEL elements must
!    be positive.  Refer to STORE3.
!
!    Input, real ( kind = 8 ) RMAX, the square root of the largest element
!    in RSQ, the maximum radius.
!
!    Input, real ( kind = 8 ) RSQ(N), the squared radii which enter into
!    the weights defining Q.
!
!    Input, real ( kind = 8 ) A(9,N), the coefficients for the nodal 
!    functions defining Q.
!
!    Output, real ( kind = 8 ) Q, the value of Q at (PX,PY,PZ) unless
!    IER == 1, in which case no values are returned.
!
!    Output, real ( kind = 8 ) QX, QY, QZ, the first partial derivatives of Q at
!    (PX,PY,PZ) unless IER == 1.
!
!    Output, integer IER, error indicator
!    0, if no errors were encountered.
!    1, if N, NR, XYZDEL, or RMAX is invalid.
!    2, if no errors were encountered but (PX.PY.PZ) is not within the 
!       radius R(K) for any node K (and thus Q = QX = QY = QZ = 0).
!
  implicit none

  integer n
  integer nr

  real ( kind = 8 ) a(9,n)
  real ( kind = 8 ) delx
  real ( kind = 8 ) dely
  real ( kind = 8 ) delz
  real ( kind = 8 ) ds
  real ( kind = 8 ) dx
  real ( kind = 8 ) dxsq
  real ( kind = 8 ) dy
  real ( kind = 8 ) dysq
  real ( kind = 8 ) dz
  real ( kind = 8 ) dzsq
  real ( kind = 8 ) f(n)
  integer i
  integer ier
  integer imax
  integer imin
  integer j
  integer jmax
  integer jmin
  integer k
  integer kmax
  integer kmin
  integer l
  integer lcell(nr,nr,nr)
  integer lnext(n)
  integer lp
  real ( kind = 8 ) px
  real ( kind = 8 ) py
  real ( kind = 8 ) pz
  real ( kind = 8 ) q
  real ( kind = 8 ) ql
  real ( kind = 8 ) qlx
  real ( kind = 8 ) qly
  real ( kind = 8 ) qlz
  real ( kind = 8 ) qx
  real ( kind = 8 ) qy
  real ( kind = 8 ) qz
  real ( kind = 8 ) rd
  real ( kind = 8 ) rds
  real ( kind = 8 ) rmax
  real ( kind = 8 ) rs
  real ( kind = 8 ) rsq(n)
  real ( kind = 8 ) sw
  real ( kind = 8 ) swq
  real ( kind = 8 ) swqx
  real ( kind = 8 ) swqy
  real ( kind = 8 ) swqz
  real ( kind = 8 ) sws
  real ( kind = 8 ) swx
  real ( kind = 8 ) swy
  real ( kind = 8 ) swz
  real ( kind = 8 ) t
  real ( kind = 8 ) w
  real ( kind = 8 ) wx
  real ( kind = 8 ) wy
  real ( kind = 8 ) wz
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) xmax
  real ( kind = 8 ) xmin
  real ( kind = 8 ) xp
  real ( kind = 8 ) xyzdel(3)
  real ( kind = 8 ) xyzmin(3)
  real ( kind = 8 ) y(n)
  real ( kind = 8 ) ymax
  real ( kind = 8 ) ymin
  real ( kind = 8 ) yp
  real ( kind = 8 ) z(n)
  real ( kind = 8 ) zmax
  real ( kind = 8 ) zmin
  real ( kind = 8 ) zp

  xp = px
  yp = py
  zp = pz
  xmin = xyzmin(1)
  ymin = xyzmin(2)
  zmin = xyzmin(3)
  dx = xyzdel(1)
  dy = xyzdel(2)
  dz = xyzdel(3)

  if ( n < 10  .or.  nr < 1  .or.  dx <= 0. &
         .or.  dy <= 0.0D+00  .or.  dz <= 0.0D+00  .or. &
         rmax < 0.0D+00 ) then
    ier = 1
    return
  end if
!
!  Set IMIN, IMAX, jmin, jmax, kmin, and kmax to cell indices
!  defining the range of the search for nodes whose radii
!  include P.  The cells which must be searched are those
!  intersected by or contained in a sphere of radius RMAX
!  centered at P.
!
  imin = int((xp-xmin-rmax)/dx) + 1
  imin = max ( imin, 1 )

  imax = int((xp-xmin+rmax)/dx) + 1
  imax = min ( imax, nr )

  jmin = int((yp-ymin-rmax)/dy) + 1
  jmin = max ( jmin, 1 )

  jmax = int((yp-ymin+rmax)/dy) + 1
  jmax = min ( jmax, nr )

  kmin = int((zp-zmin-rmax)/dz) + 1
  kmin = max ( kmin, 1 )

  kmax = int((zp-zmin+rmax)/dz) + 1
  kmax = min ( kmax, nr )
!
!  Test for no cells within the sphere of radius RMAX.
!
  if ( imax < imin .or. &
       jmax < jmin .or. &
       kmax < kmin ) then
    q = 0.0D+00
    qx = 0.0D+00
    qy = 0.0D+00
    qz = 0.0D+00
    ier = 2
    return
  end if
!
!  Q = swq/sw = sum(w(l)*q(l))/sum(w(l)) where the sum is
!  from l = 1 to N, q(l) is the quadratic nodal function,