subroutine getnp2 ( px, py, x, y, nr, lcell, lnext, xmin, ymin, &
  dx, dy, np, dsq )

!***********************************************************************
!
!! GETNP2 seeks the closest unmarked node to a point.
!
!  Discussion:
!
!    GETNP2 uses the cell method to find the closest unmarked node NP
!    to a specified point P, given a set of N nodes and the data structure 
!    defined by STORE2.
!
!    NP is then marked by negating LNEXT(NP).  Thus, the closest M nodes to
!    P may be determined by a sequence of M calls to this routine.  
!
!    If the point P is itself actually a node K, and you want to find the
!    nearest point to P that is not node K, then you must be sure to mark
!    node K before calling.
!
!    The search is begun in the cell containing or closest to P and proceeds 
!    outward in rectangular layers until all cells which contain points 
!    within distance R of P have been searched.  R is the distance from P to 
!    the first unmarked node encountered, or infinite if no unmarked nodes
!    are present.
!
!    Input parameters other than LNEXT are not altered by this routine.  
!    With the exception of ( PX, PY ) and the signs of LNEXT elements, 
!    these parameters should be unaltered from their values on output 
!    from subroutine STORE2.
!
!  Modified:
!
!    10 July 1999
!
!  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, real ( kind = 8 ) PX, PY, the (X,Y) coordinates of the point P
!    whose nearest unmarked neighbor is to be found.
!
!    Input, real ( kind = 8 ) X(N), Y(N), the coordinates of the nodes at which
!    data has been supplied.
!
!    Input, integer NR, the number of rows and columns in the cell grid.
!    NR must be at least 1.
!
!    Input, integer LCELL(NR,NR), array of nodal indices associated
!    with cells.
!
!    Input/output, integer LNEXT(N), contains next-node indices ( or their 
!    negatives ).  On return, if the output value of NP is nonzero, then
!    LNEXT(NP) will be negative.
!
!    Input, real ( kind = 8 ) XMIN, YMIN, DX, DY, the minimum nodal X, Y
!    coordinates, and the X, Y dimensions of a cell.  DX and DY must be
!    positive.
!
!    Output, integer NP, the index into the vectors X and Y of the nearest
!    unmarked node to the point P.  NP will be 0 if all nodes are marked 
!    or if the values of NR, DX, DY are illegal.  LNEXT(NP) will be less
!    than 0 if NP is nonzero (this marks node NP as being used now).
!
!    Output, real ( kind = 8 ) DSQ, if NP is nonzero, then DSQ is the 
!    squared distance between P and node NP.
!
!  Local Parameters:
!
!    first = true iff the first unmarked node has yet to be encountered,
!
!    imin, imax, jmin, jmax = cell indices defining the range of the search,
!
!    delx, dely = px-xmin and py-ymin,
!
!    i0, j0 = cell containing or closest to P,
!
!    i1, i2, j1, j2 = cell indices of the layer whose intersection with 
!    the range defined by imin,...,jmax is currently being searched.
!
  implicit none

  integer nr

  real ( kind = 8 ) delx
  real ( kind = 8 ) dely
  real ( kind = 8 ) dsq
  real ( kind = 8 ) dx
  real ( kind = 8 ) dy
  logical first
  integer i
  integer i0
  integer i1
  integer i2
  integer imax
  integer imin
  integer j
  integer j0
  integer j1
  integer j2
  integer jmax
  integer jmin
  integer l
  integer lcell(nr,nr)
  integer lmin
  integer ln
  integer lnext(*)
  integer np
  real ( kind = 8 ) px
  real ( kind = 8 ) py
  real ( kind = 8 ) r
  real ( kind = 8 ) rsmin
  real ( kind = 8 ) rsq
  real ( kind = 8 ) x(*)
  real ( kind = 8 ) xmin
  real ( kind = 8 ) xp
  real ( kind = 8 ) y(*)
  real ( kind = 8 ) ymin
  real ( kind = 8 ) yp

  xp = px
  yp = py
!
!  Test for invalid input parameters.
!
  if ( nr < 1 .or. dx <= 0.0D+00 .or. dy <= 0.0D+00 ) then
    np = 0
    dsq = 0.0D+00
  end if
!
!  Initialize parameters:
!
  first = .true.
  imin = 1
  imax = nr
  jmin = 1
  jmax = nr
  delx = xp - xmin
  dely = yp - ymin

  i0 = int ( delx / dx ) + 1
  i0 = max ( i0, 1 )
  i0 = min ( i0, nr )

  j0 = int ( dely / dy ) + 1
  j0 = max ( j0, 1 )
  j0 = min ( j0, nr )

  i1 = i0
  i2 = i0
  j1 = j0
  j2 = j0
!
!  Outer loop on layers, inner loop on layer cells, excluding
!  those outside the range (imin,imax) x (jmin,jmax).
!
1 continue

  do j = j1, j2

    if ( jmax < j ) then
      exit
    end if

    if ( j < jmin ) then
      cycle
    end if

    do i = i1, i2

      if ( imax < i ) then
        exit
      end if

      if ( i < imin ) then
        cycle
      end if

      if ( j /= j1 .and. j /= j2 .and. i /= i1 .and. i /= i2 ) then
        cycle
      end if
!
!  Search cell (i,j) for unmarked nodes l.
!
      l = lcell(i,j)

      if ( 0 < l ) then
!
!  Loop on nodes in cell (i,j).
!
2       continue

        ln = lnext(l)
!
!  Node L is the first unmarked neighbor of P encountered.
!
!  Initialize LMIN to the current candidate for np, and
!  rsmin to the squared distance from p to lmin.  imin,
!  imax, jmin, and jmax are updated to define the smal-
!  lest rectangle containing a circle of radius r =
!  sqrt(rsmin) centered at p, and contained in (1,nr) x
!  (1,nr) (except that, if p is outside the rectangle
!  defined by the nodes, it is possible that imin .gt.
!  nr, imax < 1, nr < jmin, or jmax < 1).
!
        if ( 0 <= ln ) then

          rsq = ( x(l) - xp )**2 + ( y(l) - yp )**2

          if ( first ) then

            lmin = l
            rsmin = rsq
            r = sqrt ( rsmin )

            imin = int ( ( delx - r ) / dx ) + 1
            imin = max ( imin, 1 )

            imax = int ( ( delx + r ) / dx ) + 1
            imax = min ( imax, nr )

            jmin = int ( ( dely - r ) / dy ) + 1
            jmin = max ( jmin, 1 )

            jmax = int ( ( dely + r ) / dy ) + 1
            jmax = min ( jmax, nr )

            first = .false.

          else

            if ( rsq < rsmin ) then
              lmin = l
              rsmin = rsq
            end if
 
          end if

        end if

        if ( abs ( ln ) /= l ) then
          l = abs ( ln )
          go to 2
        end if

      end if

    end do

  end do
!
!  Test for termination of loop on cell layers.
!
  if ( imin < i1 .or. i2 < imax .or. jmin < j1 .or. j2 < jmax ) then

    i1 = i1 - 1
    i2 = i2 + 1
    j1 = j1 - 1
    j2 = j2 + 1
    go to 1

  end if

  if ( first ) then
    np = 0
    dsq = 0.0D+00
  else
    np = lmin
    dsq = rsmin
    lnext(lmin) = -lnext(lmin)
  end if

  return
end
subroutine givens ( a, b, c, s )

!***********************************************************************
!
!! GIVENS constructs a Givens plane rotation.
!
!  Discussion:
!
!    The transformation has the form of a 2 by 2 matrix G(C,S):
!
!      (   C  S )
!      ( - S  C )
!
!    where C*C + S*S = 1, which zeroes the second entry of the
!    the column vector ( A, B ) when C and S are properly chosen.
!    A call to GIVENS is normally followed by a call to ROTATE
!    which computes the product of G(C,S) with a 2 by N matrix.
!
!  Modified:
!
!    10 July 1999
!
!  Parameters:
!
!    Input/output, real ( kind = 8 ) A, B.
!
!    On input, A and B define the 2-vector whose second entry (B) is
!    to be annihilated by a Givens rotation.
!
!    On output, A has been overwritten by a value
!      R = +/- SQRT ( A*A + B*B )
!    and B has been overwritten by a value Z which allows C
!    and S to be recovered as:
!
!      if | Z | <= 1, then
!        C = SQRT ( 1 - Z*Z ), 
!        S = Z
!      else if 1 < | Z | then
!        C = 1 / Z, 
!        S = SQRT ( 1 - C*C ).
!
!    Output, real ( kind = 8 ) C, S, the components of the Givens
!    transformation, which may be computed by:
!
!      C = +/- A / SQRT ( A*A + B*B )
!      S = +/- B / SQRT ( A*A + B*B )
!
!  Local parameters:
!
!  r =        c*a + s*b = +/-sqrt(a*a+b*b)
!  u,v =   variables used to scale a and b for computing r
!
!  abs(b) < abs(a)
!
!  Note that r has the sign of a, 0 < c, and s has
!  sign(a)*sign(b).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) c
  real ( kind = 8 ) r
  real ( kind = 8 ) s
  real ( kind = 8 ) u
  real ( kind = 8 ) v

  if ( abs ( b ) < abs ( a ) ) then

    u = 2.0D+00 * a
    v = b / u
    r = sqrt ( 0.25D+00 + v * v ) * u
    c = a / r
    s = 2.0D+00 * v * c
    b = s
    a = r
!
!  abs(a) <= abs(b)
!
!  Store r in a.
!  Note that r has the sign of b, 0 < s, and c has sign(a)*sign(b).
!
  else if ( b /= 0.0D+00 ) then

    u = 2.0D+00 * b
    v = a / u
    a = sqrt ( 0.25D+00 + v * v ) * u
    s = b / a
    c = 2.0D+00 * v * s

    if ( c /= 0.0D+00 ) then
      b = 1.0D+00 / c
    else
      b = 1.0D+00
    end if
!
!  a = b = 0.
!
  else

    c = 1.0D+00
    s = 0.0D+00

  end if

  return
end
subroutine qs2grd ( px, py, n, x, y, f, nr, lcell, lnext, xmin, &
  ymin, dx, dy, rmax, rsq, a, q, qx, qy, ier )

!***********************************************************************
!
!! QS2GRD evaluates the interpolant and its first spatial derivatives.
!
!  Discussion:
!
!    QS2GRD computes the value and the gradient at the point (PX,PY) 
!    of the interpolatory function Q, defined by QSHEP2 for a given set
!    of scattered data.  Q(X,Y) is a weighted sum of quadratic
!    nodal functions.
!
!    Input parameters are not altered by this subroutine.  The parameters 
!    other than PX and PY should be input unaltered from their values 
!    on output from QSHEP2.  This subroutine should not be called if a 
!    nonzero error flag was returned by QSHEP2.
!
!  Modified:
!
!    10 July 1999
!
!  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, real ( kind = 8 ) PX, PY, the coordinates of the point at which the
!    interpolant and its derivatives are to be evaluated.
!
!    Input, integer N, the number of nodes and data values which
!    are to be interpolated.  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 at the nodes.
!
!    Input, integer NR, the number of rows and columns in the cell 
!    grid.  Refer to subroutine STORE2 for details.  NR must be at least 1.
!
!    Input, integer LCELL(NR,NR), array of nodal indices associated
!    with cells.
!
!    Input, integer LNEXT(N), contains next-node indices.
!
!    Input, real ( kind = 8 ) XMIN, YMIN, DX, DY, the minimum nodal X, Y
!    coordinates, and the X, Y dimensions of a cell.  Computed by QSHEP2.
!
!    Input, real ( kind = 8 ) RMAX, the square root of the largest element
!    in RSQ, the maximum radius of influence.  Computed by QSHEP2.
!
!    Input, real ( kind = 8 ) RSQ(N), the squared radii which enter into the
!    weights defining the interpolant Q.  Computed by QSHEP2.
!
!    Input, real ( kind = 8 ) A(5,N), the coefficients for the nodal functions 
!    defining the interpolant Q.  Computed by QSHEP2.
!
!    Output, real ( kind = 8 ) Q, QX, QY, the value of the interpolant, and
!    its derivatives with respect to X and Y, at (PX,PY).
!
!    Output, integer IER, error indicator.
!    0, if no errors were encountered.
!    1, if N, NR, DX, DY or RMAX is invalid.
!    2, if no errors were encountered but (PX,PY) is not within the 
!       radius R(K) for any node K and thus Q = QX = QY = 0.
!
  implicit none

  integer n
  integer nr

  real ( kind = 8 ) a(5,n)
  real ( kind = 8 ) delx
  real ( kind = 8 ) dely
  real ( kind = 8 ) ds
  real ( kind = 8 ) dx
  real ( kind = 8 ) dy
  real ( kind = 8 ) f(n)
  integer i
  integer ier
  integer imax
  integer imin
  integer j
  integer jmax
  integer jmin
  integer k
  integer kp
  integer lcell(nr,nr)
  integer lnext(n)
  real ( kind = 8 ) px
  real ( kind = 8 ) py
  real ( kind = 8 ) q
  real ( kind = 8 ) qk
  real ( kind = 8 ) qkx
  real ( kind = 8 ) qky
  real ( kind = 8 ) qx
  real ( kind = 8 ) qy
  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 ) sws
  real ( kind = 8 ) swx
  real ( kind = 8 ) swy
  real ( kind = 8 ) t
  real ( kind = 8 ) w
  real ( kind = 8 ) wx
  real ( kind = 8 ) wy
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) xmin
  real ( kind = 8 ) xp
  real ( kind = 8 ) y(n)
  real ( kind = 8 ) ymin
  real ( kind = 8 ) yp

  xp = px
  yp = py

  if ( n < 6 ) then
    ier = 1
    return
  else if ( nr < 1 ) then
    ier = 1
    return
  else if ( dx <= 0.0D+00 ) then
    ier = 1
    return
  else if ( dy <= 0.0D+00 ) then
    ier = 1
    return
  else if ( rmax < 0.0D+00 ) then
    ier = 1
    return
  end if
!
!  Set imin, imax, jmin, and jmax to cell indices defining
!  the range of the search for nodes whose radii include P.
!  The cells which must be searched are those inter-
!  sected by (or contained in) a circle 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 )
!
!  Test for no cells within the circle of radius RMAX.
!
  if ( imax < imin .or. jmax < jmin ) then
    q = 0.0D+00
    qx = 0.0D+00
    qy = 0.0D+00
    ier = 2
    return
  end if
!
!  Q = swq/sw = sum(w(k)*q(k))/sum(w(k)) where the sum is
!  from k = 1 to n, q(k) is the quadratic nodal function,
!  and w(k) = ((r-d)+/(r*d))**2 for radius r(k) and distance d(k).  Thus
!
!    qx = (swqx*sw - swq*swx)/sw**2  and
!    qy = (swqy*sw - swq*swy)/sw**2
!
!  where swqx and swx are partial derivatives with respect
!  to x of swq and sw, respectively.  swqy and swy are 
!  defined similarly.
!