sw = 0.0D+00
  swx = 0.0D+00
  swy = 0.0D+00
  swq = 0.0D+00
  swqx = 0.0D+00
  swqy = 0.0D+00
!
!  Outer loop on cells (I,J).
!
  do j = jmin, jmax

    do i = imin, imax

      k = lcell(i,j)
!
!  Inner loop on nodes K.
!
      if ( k /= 0 ) then

        do

          delx = xp - x(k)
          dely = yp - y(k)
          ds = delx * delx + dely * dely
          rs = rsq(k)

          if ( ds == 0.0D+00 ) then
            q = f(k)
            qx = a(4,k)
            qy = a(5,k)
            ier = 0
            return
          end if

          if ( ds < rs ) then

            rds = rs * ds
            rd = sqrt ( rds )
            w = ( rs + ds - rd - rd ) / rds
            t = 2.0D+00 * ( rd - rs ) / ( ds * rds )
            wx = delx * t
            wy = dely * t
            qkx = 2.0D+00 * a(1,k) * delx + a(2,k) * dely
            qky = a(2,k) * delx + 2.0D+00 * a(3,k) * dely
            qk = ( qkx * delx + qky * dely ) / 2.0D+00
            qkx = qkx + a(4,k)
            qky = qky + a(5,k)
            qk = qk + a(4,k) * delx + a(5,k) * dely + f(k)
            sw = sw + w
            swx = swx + wx
            swy = swy + wy
            swq = swq + w * qk
            swqx = swqx + wx * qk + w * qkx
            swqy = swqy + wy * qk + w * qky

          end if

          kp = k
          k = lnext(kp)

          if ( k == kp ) then
            exit
          end if

        end do

      end if

    end do

  end do
!
!  SW = 0 if and only if P is not within the radius R(K) for any node K.
!
  if ( sw /= 0.0D+00 ) then

    q = swq / sw
    sws = sw * sw
    qx = ( swqx * sw - swq * swx ) / sws
    qy = ( swqy * sw - swq * swy ) / sws
    ier = 0

  else

    q = 0.0D+00
    qx = 0.0D+00
    qy = 0.0D+00
    ier = 2

  end if

  return
end
subroutine qshep2 ( n, x, y, f, nq, nw, nr, lcell, lnext, xmin, &
  ymin, dx, dy, rmax, rsq, a, ier )

!***********************************************************************
!
!! QSHEP2 computes an interpolant to scattered data in the plane.
!
!  Discussion:
!
!    QSHEP2 computes a set of parameters A and RSQ defining a smooth, 
!    once continuously differentiable, bi-variate function Q(X,Y) which 
!    interpolates given data values F at scattered nodes (X,Y).  
!
!    The interpolant function Q(X,Y) may be evaluated at an arbitrary point 
!    by passing the parameters A and RSQ to the function QS2VAL.  The
!    first derivatives dQdX(X,Y) and dQdY(X,Y) may be evaluated by 
!    subroutine QS2GRD.
!
!    The interpolation scheme is a modified quadratic Shepard method:
!
!      Q = ( W(1) * Q(1) + W(2) * Q(2) + .. + W(N) * Q(N) ) 
!        / ( W(1)        + W(2)        + .. + W(N) )
!
!    for bivariate functions W(K) and Q(K).  The nodal functions are given by
!
!      Q(K)(X,Y) = 
!          F(K)
!        + A(1,K) * ( X - X(K) )**2 
!        + A(2,K) * ( X - X(K) ) * ( Y - Y(K) )
!        + A(3,K) * ( Y - Y(K) )**2 
!        + A(4,K) * ( X - X(K) )
!        + A(5,K) * ( Y - Y(K) ).
!
!    Thus, Q(K) is a quadratic function which interpolates the
!    data value at node K.  Its coefficients A(*,K) are obtained
!    by a weighted least squares fit to the closest NQ data
!    points with weights similar to W(K).  Note that the radius
!    of influence for the least squares fit is fixed for each
!    K, but varies with K.
!
!    The weights are taken to be
!
!      W(K)(X,Y) = ( (R(K)-D(K))+ / R(K) * D(K) )**2
!
!    where (R(K)-D(K))+ = 0 if R(K) <= D(K) and D(K)(X,Y) is
!    the euclidean distance between (X,Y) and (X(K),Y(K)).  The
!    radius of influence R(K) varies with K and is chosen so
!    that NW nodes are within the radius.  Note that W(K) is
!    not defined at node (X(K),Y(K)), but Q(X,Y) has limit F(K)
!    as (X,Y) approaches (X(K),Y(K)).
!
!  Author:
!
!    Robert Renka,
!    University of North Texas
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 660: QSHEP2D, Quadratic Shepard method for bivariate
!    interpolation of scattered data,
!    ACM Transactions on Mathematical Software,
!    Volume 14, 1988, pages 149-150.
!
!  Parameters:
!
!    Input, integer N, the number of nodes (X,Y) at which data values
!    are given.  N must be at least 6.
!
!    Input, real ( kind = 8 ) X(N), Y(N), the coordinates of the nodes at which
!    data has been supplied.
!
!    Input, real ( kind = 8 ) F(N), the data values.
!
!    Input, integer NQ, the number of data points to be used in the least
!    squares fit for coefficients defining the nodal functions Q(K).  
!    A highly recommended value is NQ = 13.  
!    NQ must be at least 5, and no greater than the minimum of 40 and N-1.
!
!    Input, integer NW, the number of nodes within (and defining) the radii
!    of influence R(K) which enter into the weights W(K).  For N 
!    sufficiently large, a recommended value is NW = 19.   NW must be
!    at least 1, and no greater than the minimum of 40 and N-1.
!
!    Input, integer NR, the number of rows and columns in the cell grid 
!    defined in subroutine STORE2.  A rectangle containing the nodes 
!    is partitioned into cells in order to increase search efficiency.  
!    NR = SQRT(N/3) is recommended.  NR must be at least 1.
!
!    Output, integer LCELL(NR,NR), array of nodal indices associated
!    with cells.
!
!    Output, integer LNEXT(N), contains next-node indices ( or their 
!    negatives ).
!
!    Output, real ( kind = 8 ) XMIN, YMIN, DX, DY, the minimum nodal X, Y
!    coordinates, and the X, Y dimensions of a cell.
!
!    Output, real ( kind = 8 ) RMAX, the square root of the largest element
!    in RSQ, the maximum radius of influence.
!
!    Output, real ( kind = 8 ) RSQ(N), the squared radii which enter into 
!    the weights defining the interpolant Q.
!
!    Output, real ( kind = 8 ) A(5,N), the coefficients for the nodal functions 
!    defining the interpolant Q.
!
!    Output, integer IER, error indicator.
!    0, if no errors were encountered.
!    1, if N, NQ, NW, or NR is out of range.
!    2, if duplicate nodes were encountered.
!    3, if all nodes are collinear.
!
!  Local parameters:
!
! av =        root-mean-square distance between k and the
!             nodes in the least squares system (unless
!             additional nodes are introduced for stabil-
!             ity).      the first 3 columns of the matrix
!             are scaled by 1/avsq, the last 2 by 1/av
! avsq =      av*av
! b =         transpose of the augmented regression matrix
! c =         first component of the plane rotation used to
!             zero the lower triangle of b**t -- computed
!             by subroutine givens
! ddx,ddy =   local variables for dx and dy
! dmin =      minimum of the magnitudes of the diagonal
!             elements of the regression matrix after
!             zeros are introduced below the diagonal
! DTOL =      tolerance for detecting an ill-conditioned
!             system.  the system is accepted when DTOL <= DMIN.
! fk =        data value at node k -- f(k)
! i =         index for a, b, and npts
! ib =        do-loop index for back solve
! ierr =      error flag for the call to store2
! irow =      row index for b
! j =         index for a and b
! jp1 =       j+1
! k =         nodal function index and column index for a
! lmax =      maximum number of npts elements (must be con-
!             sistent with the dimension statement above)
! lnp =       current length of npts
! neq =       number of equations in the least squares fit
! nn,nnr =    local copies of n and nr
! np =        npts element
! npts =      array containing the indices of a sequence of
!             nodes to be used in the least squares fit
!             or to compute rsq.  the nodes are ordered
!             by distance from k and the last element
!             (usually indexed by lnp) is used only to
!             determine rq, or rsq(k) if NQ < NW.
! nqwmax =    max(nq,nw)
! rq =        radius of influence which enters into the
!             weights for q(k) (see subroutine setup2)
! rs =        squared distance between k and npts(lnp) --
!             used to compute rq and rsq(k)
! rsmx =      maximum rsq element encountered
! rsold =     squared distance between k and npts(lnp-1) --
!             used to compute a relative change in rs
!             between succeeding npts elements
! RTOL =      tolerance for detecting a sufficiently large
!             relative change in rs.  if the change is
!             not greater than RTOL, the nodes are
!             treated as being the same distance from k
! rws =       current value of rsq(k)
! s =         second component of the plane givens rotation
! SF =        marquardt stabilization factor used to damp
!             out the first 3 solution components (second
!             partials of the quadratic) when the system
!             is ill-conditioned.  as SF increases, the
!             fitting function approaches a linear
! sum2 =      sum of squared euclidean distances between
!             node k and the nodes used in the least
!             squares fit (unless additional nodes are
!             added for stability)
! t =         temporary variable for accumulating a scalar
!             product in the back solve
! xk,yk =     coordinates of node k -- x(k), y(k)
! xmn,ymn =   local variables for xmin and ymin
!
  implicit none

  integer n
  integer nr

  real ( kind = 8 ) a(5,n)
  real ( kind = 8 ) av
  real ( kind = 8 ) avsq
  real ( kind = 8 ) b(6,6)
  real ( kind = 8 ) c
  real ( kind = 8 ) ddx
  real ( kind = 8 ) ddy
  real ( kind = 8 ) dmin
  real ( kind = 8 ), parameter :: dtol = 0.01D+00
  real ( kind = 8 ) dx
  real ( kind = 8 ) dy
  real ( kind = 8 ) f(n)
  real ( kind = 8 ) fk
  integer i
  integer ier
  integer ierr
  integer irow
  integer j
  integer jp1
  integer k
  integer lcell(nr,nr)
  integer lmax
  integer lnext(n)
  integer lnp
  integer neq
  integer nn
  integer nnr
  integer np
  integer npts(40)
  integer nq
  integer nqwmax
  integer nw
  real ( kind = 8 ) rmax
  real ( kind = 8 ) rq
  real ( kind = 8 ) rs
  real ( kind = 8 ) rsmx
  real ( kind = 8 ) rsold
  real ( kind = 8 ) rsq(n)
  real ( kind = 8 ), parameter :: rtol = 1.0D-05
  real ( kind = 8 ) rws
  real ( kind = 8 ) s
  real ( kind = 8 ), parameter :: SF = 1.0D+00
  real ( kind = 8 ) sum2
  real ( kind = 8 ) t
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) xk
  real ( kind = 8 ) xmin
  real ( kind = 8 ) xmn
  real ( kind = 8 ) y(n)
  real ( kind = 8 ) yk
  real ( kind = 8 ) ymin
  real ( kind = 8 ) ymn

  nn = n
  nnr = nr
  nqwmax = max ( nq, nw )
  lmax = min ( 40, n-1 )

  if ( nq < 5 ) then
    ier = 1
    return
  else if ( nw < 1 ) then
    ier = 1
    return
  else if ( lmax < nqwmax ) then
    ier = 1
    return
  else if ( nr < 1 ) then
    ier = 1
    return
  end if
!
!  Create the cell data structure, and initialize RSMX.
!
  call store2 ( nn, x, y, nnr, lcell, lnext, xmn, ymn, ddx, ddy, ierr )

  if ( ierr /= 0 ) then
    xmin = xmn
    ymin = ymn
    dx = ddx
    dy = ddy
    ier = 3
    return
  end if

  rsmx = 0.0D+00
!
!  Outer loop on node K.
!
  do k = 1, nn

    xk = x(k)
    yk = y(k)
    fk = f(k)
!
!  Mark node K to exclude it from the search for nearest neighbors.
!
    lnext(k) = -lnext(k)
!
!  Initialize for loop on NPTS.
!
    rs = 0.0D+00
    sum2 = 0.0D+00
    rws = 0.0D+00
    rq = 0.0D+00
    lnp = 0
!
!  Compute NPTS, LNP, RWS, NEQ, RQ, and AVSQ.
!
1   continue

    sum2 = sum2 + rs

    if ( lnp == lmax ) then
      go to 3
    end if

    lnp = lnp + 1
    rsold = rs

    call getnp2 ( xk, yk, x, y, nnr, lcell, lnext, xmn, ymn, ddx, ddy, np, rs )

    if ( rs == 0.0D+00 ) then
      ier = 2
      return
    end if

    npts(lnp) = np

    if ( ( rs - rsold ) / rs < RTOL ) then
      go to 1
    end if

    if ( rws == 0.0D+00 .and. nw < lnp ) then
      rws = rs
    end if
!
!  RQ = 0 (not yet computed) and NQ < lnp.     
!
!  RQ = sqrt(RS) is sufficiently large to (strictly) include NQ nodes.  
!
!  The least squares fit will include NEQ = LNP - 1 equations for 
!  5 <= NQ <= NEQ < LMAX <= N-1.
!
    if ( rq == 0.0D+00 .and. nq < lnp ) then
      neq = lnp - 1
      rq = sqrt ( rs )
      avsq = sum2 / real ( neq, kind = 8 )
    end if

    if ( nqwmax < lnp ) then
      go to 4
    else
      go to 1
    end if
!
!  All LMAX nodes are included in NPTS.   RWS and/or RQ**2 is
!  (arbitrarily) taken to be 10 percent larger than the
!  distance RS to the last node included.
!
3   continue

    if ( rws == 0.0D+00 ) then
      rws = 1.1D+00 * rs
    end if

    if ( rq == 0.0D+00 ) then
      neq = lmax
      rq = sqrt ( 1.1D+00 * rs )
      avsq = sum2 / real ( neq, kind = 8 )
    end if

4   continue
!
!  Store RSQ(K), update RSMX if necessary, and compute AV.
!
    rsq(k) = rws
    rsmx = max ( rsmx, rws )
    av = sqrt ( avsq )
!
!  Set up the augmented regression matrix (transposed) as the
!  columns of B, and zero out the lower triangle (upper
!  triangle of B) with Givens rotations -- QR decomposition
!  with orthogonal matrix Q not stored.
!
    i = 0

5   continue

    i = i + 1
    np = npts(i)
    irow = min ( i, 6 )

    call setup2 ( xk, yk, fk, x(np), y(np), f(np), av, avsq, rq, b(1,irow) )

    if ( i == 1 ) then
      go to 5
    end if

    do j = 1, irow-1
      jp1 = j + 1
      call givens ( b(j,j), b(j,irow), c, s )
      call rotate ( 6-j, c, s, b(jp1,j), b(jp1,irow) )
    end do

    if ( i < neq ) then
      go to 5
    end if
!
!  Test the system for ill-conditioning.
!
    dmin =  min ( abs ( b(1,1) ), abs ( b(2,2) ), abs ( b(3,3) ), &
      abs ( b(4,4) ), abs ( b(5,5) ) )

    if ( DTOL <= dmin * rq ) then
      go to 13
    end if

    if ( neq == lmax ) then
      go to 10
    end if
!
!  Increase RQ and add another equation to the system to improve conditioning.  
!  The number of NPTS elements is also increased if necessary.
!
7   continue

    rsold = rs
    neq = neq + 1

    if ( neq == lmax ) then
      go to 9
    end if
!
!   NEQ < LNP.
!
    if ( neq /= lnp ) then
      np = npts(neq+1)
      rs = ( x(np) - xk )**2 + ( y(np) - yk )**2
      if ( ( rs - rsold ) / rs < rtol ) then
        go to 7
      end if
      rq = sqrt(rs)
      go to 5
    end if
!
!  Add an element to NPTS.
!
    lnp = lnp + 1
    call getnp2 ( xk, yk, x, y, nnr, lcell, lnext, xmn, ymn, ddx, ddy, np, rs )

    if ( np == 0 ) then
      ier = 2
      return
    end if

    npts(lnp) = np

    if ( ( rs - rsold ) / rs < rtol ) then
      go to 7
    end if

    rq = sqrt ( rs )
    go to 5

9   continue

    rq = sqrt ( 1.1D+00 * rs )
    go to 5
!
!  Stabilize the system by damping second partials.  Add multiples of the 
!  first three unit vectors to the first three equations.
!
10  continue

    do i = 1, 3

      b(i,6) = sf
      b(i+1:6,6) = 0.0D+00

      do j = i, 5
        jp1 = j + 1
        call givens ( b(j,j), b(j,6), c, s )
        call rotate ( 6-j, c, s, b(jp1,j), b(jp1,6) )
      end do

    end do
!
!  Test the stabilized system for ill-conditioning.
!
    dmin = min ( abs ( b(1,1) ), abs ( b(2,2) ), abs ( b(3,3) ), &
      abs ( b(4,4) ), abs ( b(5,5) ) )

    if ( dmin * rq < dtol ) then
      xmin = xmn
      ymin = ymn
      dx = ddx
      dy = ddy
      ier = 3
      return
    end if
!
!  Solve the 5 by 5 triangular system for the coefficients.
!
13  continue