      FUNCTION STORE (X)
      REAL X
C
C***********************************************************
C
C                                    ROBERT RENKA
C                             NORTH TEXAS STATE UNIV.
C                                    (817) 565-2767
C
C   THIS FUNCTION FORCES ITS ARGMENT X TO BE STORED IN MAIN
C MEMORY.  THIS IS USEFUL FOR COMPUTING MACHINE DEPENDENT
C PARAMETERS (SUCH AS THE MACHINE PRECISION) WHERE IT IS
C NECESSARY TO AVOID COMPUTATION IN HIGH PRECISION REG-
C ISTERS.
C
C INPUT PARAMETER -
C
C      X - VALUE TO BE STORED.
C
C X IS NOT ALTERED BY THIS FUNCTION.
C
C OUTPUT PARAMETER -
C
C      STORE - VALUE OF X AFTER IT HAS BEEN STORED AND
C            (POSSIBLY) TRUNCATED OR ROUNDED TO THE
C            MACHINE WORD LENGTH.
C
C MODULES REQUIRED BY STORE - NONE
C
C***********************************************************
C
      COMMON/STCOM/Y
      Y = X
      STORE = Y
      RETURN
      END
      SUBROUTINE QSHEP3 (N,X,Y,Z,F,NQ,NW,NR, LCELL,LNEXT,
     .                   XYZMIN,XYZDEL,RMAX,RSQ,A,IER)
      INTEGER N, NQ, NW, NR, LCELL(NR,NR,NR), LNEXT(N), IER
      REAL    X(N), Y(N), Z(N), F(N), XYZMIN(3), XYZDEL(3),
     .            RMAX, RSQ(N), A(9,N)
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C                                        01/08/90
C
C   THIS SUBROUTINE COMPUTES A SET OF PARAMETERS A AND RSQ
C DEFINING A SMOOTH (ONCE CONTINUOUSLY DIFFERENTIABLE) TRI-
C VARIATE FUNCTION Q(X,Y,Z) WHICH INTERPOLATES DATA VALUES
C F AT SCATTERED NODES (X,Y,Z).  THE INTERPOLANT Q MAY BE
C EVALUATED AT AN ARBITRARY POINT BY FUNCTION QS3VAL, AND
C ITS FIRST DERIVATIVES ARE COMPUTED BY SUBROUTINE QS3GRD.
C   THE INTERPOLATION SCHEME IS A MODIFIED QUADRATIC SHEPARD
C METHOD --
C
C Q = (W(1)*Q(1)+W(2)*Q(2)+..+W(N)*Q(N))/(W(1)+W(2)+..+W(N))
C
C FOR TRIVARIATE FUNCTIONS W(K) AND Q(K).  THE NODAL FUNC-
C TIONS ARE GIVEN BY
C
C Q(K)(X,Y,Z) = A(1,K)*DX**2 + A(2,K)*DX*DY + A(3,K)*DY**2 +
C            A(4,K)*DX*DZ + A(5,K)*DY*DZ + A(6,K)*DZ**2 +
C            A(7,K)*DX + A(8,K)*DY + A(9,K)*DZ + F(K)
C
C WHERE DX = (X-X(K)), DY = (Y-Y(K)), AND DZ = (Z-Z(K)).
C THUS, Q(K) IS A QUADRATIC FUNCTION WHICH INTERPOLATES THE
C DATA VALUE AT NODE K.  ITS COEFFICIENTS A(,K) ARE OBTAINED
C BY A WEIGHTED LEAST SQUARES FIT TO THE CLOSEST NQ DATA
C POINTS WITH WEIGHTS SIMILAR TO W(K).      NOTE THAT THE RADIUS
C OF INFLUENCE FOR THE LEAST SQUARES FIT IS FIXED FOR EACH
C K, BUT VARIES WITH K.
C   THE WEIGHTS ARE TAKEN TO BE
C
C W(K)(X,Y,Z) = ( (R(K)-D(K))+ / R(K)*D(K) )**2
C
C WHERE (R(K)-D(K))+ = 0 IF R(K) .LE. D(K), AND D(K)(X,Y,Z)
C IS THE EUCLIDEAN DISTANCE BETWEEN (X,Y,Z) AND NODE K.  THE
C RADIUS OF INFLUENCE R(K) VARIES WITH K AND IS CHOSEN SO
C THAT NW NODES ARE WITHIN THE RADIUS.      NOTE THAT W(K) IS
C NOT DEFINED AT NODE (X(K),Y(K),Z(K)), BUT Q(X,Y,Z) HAS
C LIMIT F(K) AS (X,Y,Z) APPROACHES (X(K),Y(K),Z(K)).
C
C ON INPUT --
C
C      N = NUMBER OF NODES AND ASSOCIATED DATA VALUES.
C          N .GE. 10.
C
C      X,Y,Z = ARRAYS OF LENGTH N CONTAINING THE CARTESIAN
C            COORDINATES OF THE NODES.
C
C      F = ARRAY OF LENGTH N CONTAINING THE DATA VALUES
C          IN ONE-TO-ONE CORRESPONDENCE WITH THE NODES.
C
C      NQ = NUMBER OF DATA POINTS TO BE USED IN THE LEAST
C           SQUARES FIT FOR COEFFICIENTS DEFINING THE NODAL
C           FUNCTIONS Q(K).  A RECOMMENDED VALUE IS NQ =
C           17.  9 .LE. NQ .LE. MIN(40,N-1).
C
C      NW = NUMBER OF NODES WITHIN (AND DEFINING) THE RADII
C           OF INFLUENCE R(K) WHICH ENTER INTO THE WEIGHTS
C           W(K).  FOR N SUFFICIENTLY LARGE, A RECOMMENDED
C           VALUE IS NW = 32.      1 .LE. NW .LE. MIN(40,N-1).
C
C      NR = NUMBER OF ROWS, COLUMNS, AND PLANES IN THE CELL
C           GRID DEFINED IN SUBROUTINE STORE3.  A BOX CON-
C           TAINING THE NODES IS PARTITIONED INTO CELLS IN
C           ORDER TO INCREASE SEARCH EFFICIENCY.  NR =
C           (N/3)**(1/3) IS RECOMMENDED.  NR .GE. 1.
C
C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C      LCELL = ARRAY OF LENGTH .GE. NR**3.
C
C      LNEXT = ARRAY OF LENGTH .GE. N.
C
C      XYZMIN,XYZDEL = ARRAYS OF LENGTH .GE. 3.
C
C      RSQ = ARRAY OF LENGTH .GE. N.
C
C      A = ARRAY OF LENGTH .GE. 9N.
C
C ON OUTPUT --
C
C      LCELL = NR BY NR BY NR ARRAY OF NODAL INDICES ASSO-
C            CIATED WITH CELLS.  REFER TO STORE3.
C
C      LNEXT = ARRAY OF LENGTH N CONTAINING NEXT-NODE INDI-
C            CES.  REFER TO STORE3.
C
C      XYZMIN,XYZDEL = ARRAYS OF LENGTH 3 CONTAINING MINI-
C                  MUM NODAL COORDINATES AND CELL DIM-
C                  ENSIONS, RESPECTIVELY.      REFER TO
C                  STORE3.
C
C      RMAX = SQUARE ROOT OF THE LARGEST ELEMENT IN RSQ --
C             MAXIMUM RADIUS R(K).
C
C      RSQ = ARRAY CONTAINING THE SQUARES OF THE RADII R(K)
C            WHICH ENTER INTO THE WEIGHTS W(K).
C
C      A = 9 BY N ARRAY CONTAINING THE COEFFICIENTS FOR
C          QUADRATIC NODAL FUNCTION Q(K) IN COLUMN K.
C
C   NOTE THAT THE ABOVE OUTPUT PARAMETERS ARE NOT DEFINED
C UNLESS IER = 0.
C
C      IER = ERROR INDICATOR --
C            IER = 0 IF NO ERRORS WERE ENCOUNTERED.
C            IER = 1 IF N, NQ, NW, OR NR IS OUT OF RANGE.
C            IER = 2 IF DUPLICATE NODES WERE ENCOUNTERED.
C            IER = 3 IF ALL NODES ARE COPLANAR.
C
C MODULES REQUIRED BY QSHEP3 -- GETNP3, GIVENS, ROTATE,
C                          SETUP3, STORE3
C
C INTRINSIC FUNCTIONS CALLED BY QSHEP3 -- ABS, AMIN1, FLOAT,
C                                  MAX0, MIN0, SQRT
C
C***********************************************************
C
      INTEGER I, IB, IERR, IP1, IRM1, IROW, J, JP1, K, LMAX,
     .            LNP, NEQ, NN, NNQ, NNR, NNW, NP, NPTS(40),
     .            NQWMAX
      REAL    AV, AVSQ, B(10,10), C, DMIN, DTOL, FK, RQ, RS,
     .            RSMX, RSOLD, RTOL, RWS, S, SF, SUM, T, XK, YK,
     .            ZK, XYZDL(3), XYZMN(3)
C
      DATA    RTOL/1.E-5/, DTOL/.01/, SF/1./
C
C LOCAL PARAMETERS --
C
C AV =             ROOT-MEAN-SQUARE DISTANCE BETWEEN K AND THE
C             NODES IN THE LEAST SQUARES SYSTEM (UNLESS
C             ADDITIONAL NODES ARE INTRODUCED FOR STABIL-
C             ITY).      THE FIRST 6 COLUMNS OF THE MATRIX
C             ARE SCALED BY 1/AVSQ, THE LAST 3 BY 1/AV
C AVSQ =       AV*AV
C B =             TRANSPOSE OF THE AUGMENTED REGRESSION MATRIX
C C =             FIRST COMPONENT OF THE PLANE ROTATION USED TO
C             ZERO THE LOWER TRIANGLE OF B**T -- COMPUTED
C             BY SUBROUTINE GIVENS
C DMIN =       MINIMUM OF THE MAGNITUDES OF THE DIAGONAL
C             ELEMENTS OF THE REGRESSION MATRIX AFTER
C             ZEROS ARE INTRODUCED BELOW THE DIAGONAL
C DTOL =       TOLERANCE FOR DETECTING AN ILL-CONDITIONED
C             SYSTEM.  THE SYSTEM IS ACCEPTED WHEN DMIN
C             .GE. DTOL
C FK =             DATA VALUE AT NODE K -- F(K)
C I =             INDEX FOR A, B, NPTS, XYZMIN, XYZMN, XYZDEL,
C             AND XYZDL
C IB =             DO-LOOP INDEX FOR BACK SOLVE
C IERR =       ERROR FLAG FOR THE CALL TO STORE3
C IP1 =        I+1
C IRM1 =       IROW-1
C IROW =       ROW INDEX FOR B
C J =             INDEX FOR A AND B
C JP1 =        J+1
C K =             NODAL FUNCTION INDEX AND COLUMN INDEX FOR A
C LMAX =       MAXIMUM NUMBER OF NPTS ELEMENTS (MUST BE CON-
C             SISTENT WITH THE DIMENSION STATEMENT ABOVE)
C LNP =        CURRENT LENGTH OF NPTS
C NEQ =        NUMBER OF EQUATIONS IN THE LEAST SQUARES FIT
C NN,NNQ,NNR = LOCAL COPIES OF N, NQ, AND NR
C NNW =        LOCAL COPY OF NW
C NP =             NPTS ELEMENT
C NPTS =       ARRAY CONTAINING THE INDICES OF A SEQUENCE OF
C             NODES TO BE USED IN THE LEAST SQUARES FIT
C             OR TO COMPUTE RSQ.  THE NODES ARE ORDERED
C             BY DISTANCE FROM K AND THE LAST ELEMENT
C             (USUALLY INDEXED BY LNP) IS USED ONLY TO
C             DETERMINE RQ, OR RSQ(K) IF NW .GT. NQ
C NQWMAX =     MAX(NQ,NW)
C RQ =             RADIUS OF INFLUENCE WHICH ENTERS INTO THE
C             WEIGHTS FOR Q(K) (SEE SUBROUTINE SETUP3)
C RS =             SQUARED DISTANCE BETWEEN K AND NPTS(LNP) --
C             USED TO COMPUTE RQ AND RSQ(K)
C RSMX =       MAXIMUM RSQ ELEMENT ENCOUNTERED
C RSOLD =      SQUARED DISTANCE BETWEEN K AND NPTS(LNP-1) --
C             USED TO COMPUTE A RELATIVE CHANGE IN RS
C             BETWEEN SUCCEEDING NPTS ELEMENTS
C RTOL =       TOLERANCE FOR DETECTING A SUFFICIENTLY LARGE
C             RELATIVE CHANGE IN RS.  IF THE CHANGE IS
C             NOT GREATER THAN RTOL, THE NODES ARE
C             TREATED AS BEING THE SAME DISTANCE FROM K
C RWS =        CURRENT VALUE OF RSQ(K)
C S =             SECOND COMPONENT OF THE PLANE GIVENS ROTATION
C SF =             MARQUARDT STABILIZATION FACTOR USED TO DAMP
C             OUT THE FIRST 6 SOLUTION COMPONENTS (SECOND
C             PARTIALS OF THE QUADRATIC) WHEN THE SYSTEM
C             IS ILL-CONDITIONED.  AS SF INCREASES, THE
C             FITTING FUNCTION APPROACHES A LINEAR
C SUM =        SUM OF SQUARED EUCLIDEAN DISTANCES BETWEEN
C             NODE K AND THE NODES USED IN THE LEAST
C             SQUARES FIT (UNLESS ADDITIONAL NODES ARE
C             ADDED FOR STABILITY)
C T =             TEMPORARY VARIABLE FOR ACCUMULATING A SCALAR
C             PRODUCT IN THE BACK SOLVE
C XK,YK,ZK =   COORDINATES OF NODE K -- X(K), Y(K), Z(K)
C XYZDL =      LOCAL VARIABLES FOR XYZDEL
C XYZMN =      LOCAL VARIABLES FOR XYZMIN
C
      NN = N
      NNQ = NQ
      NNW = NW
      NNR = NR
      NQWMAX = MAX0(NNQ,NNW)
      LMAX = MIN0(40,NN-1)
      IF (9 .GT. NNQ  .OR.  1 .GT. NNW      .OR.  NQWMAX .GT.
     .        LMAX      .OR.  NNR .LT. 1) GO TO 20
C
C CREATE THE CELL DATA STRUCTURE, AND INITIALIZE RSMX.
C
      CALL STORE3 (NN,X,Y,Z,NNR, LCELL,LNEXT,XYZMN,XYZDL,
     .               IERR)
      IF (IERR .NE. 0) GO TO 22
      RSMX = 0.
C
C OUTER LOOP ON NODE K
C
      DO 18 K = 1,NN
      XK = X(K)
      YK = Y(K)
      ZK = Z(K)
      FK = F(K)
C
C MARK NODE K TO EXCLUDE IT FROM THE SEARCH FOR NEAREST
C   NEIGHBORS.
C
      LNEXT(K) = -LNEXT(K)
C
C INITIALIZE FOR LOOP ON NPTS.
C
      RS = 0.
      SUM = 0.
      RWS = 0.
      RQ = 0.
      LNP = 0
C
C COMPUTE NPTS, LNP, RWS, NEQ, RQ, AND AVSQ.
C
    1      SUM = SUM + RS
        IF (LNP .EQ. LMAX) GO TO 3
        LNP = LNP + 1
        RSOLD = RS
        CALL GETNP3 (XK,YK,ZK,X,Y,Z,NNR,LCELL,LNEXT,XYZMN,
     .                   XYZDL, NP,RS)
        IF (RS .EQ. 0.) GO TO 21
        NPTS(LNP) = NP
        IF ( (RS-RSOLD)/RS .LT.  RTOL ) GO TO 1
        IF (RWS .EQ. 0.  .AND.  LNP .GT. NNW) RWS = RS
        IF (RQ .NE. 0.  .OR.      LNP .LE. NNQ) GO TO 2
C
C   RQ = 0 (NOT YET COMPUTED) AND LNP .GT. NQ.      RQ =
C     SQRT(RS) IS SUFFICIENTLY LARGE TO (STRICTLY) INCLUDE
C     NQ NODES.  THE LEAST SQUARES FIT WILL INCLUDE NEQ =
C     LNP-1 EQUATIONS FOR 9 .LE. NQ .LE. NEQ .LT. LMAX
C     .LE. N-1.
C
        NEQ = LNP - 1
        RQ = SQRT(RS)
        AVSQ = SUM/FLOAT(NEQ)
C
C   BOTTOM OF LOOP -- TEST FOR TERMINATION.
C
    2        IF (LNP .GT. NQWMAX) GO TO 4
        GO TO 1
C
C ALL LMAX NODES ARE INCLUDED IN NPTS.      RWS AND/OR RQ**2 IS
C   (ARBITRARILY) TAKEN TO BE 10 PERCENT LARGER THAN THE
C   DISTANCE RS TO THE LAST NODE INCLUDED.
C
    3      IF (RWS .EQ. 0.) RWS = 1.1*RS
      IF (RQ .NE. 0.) GO TO 4
        NEQ = LMAX
        RQ = SQRT(1.1*RS)
        AVSQ = SUM/FLOAT(NEQ)
C
C STORE RSQ(K), UPDATE RSMX IF NECESSARY, AND COMPUTE AV.
C
    4      RSQ(K) = RWS
      IF (RWS .GT. RSMX) RSMX = RWS
      AV = SQRT(AVSQ)
C
C SET UP THE AUGMENTED REGRESSION MATRIX (TRANSPOSED) AS THE
C   COLUMNS OF B, AND ZERO OUT THE LOWER TRIANGLE (UPPER
C   TRIANGLE OF B) WITH GIVENS ROTATIONS -- QR DECOMPOSITION
C   WITH ORTHOGONAL MATRIX Q NOT STORED.
C
      I = 0
    5   I = I + 1
        NP = NPTS(I)
        IROW = MIN0(I,10)
        CALL SETUP3 (XK,YK,ZK,FK,X(NP),Y(NP),Z(NP),F(NP),
     .                   AV,AVSQ,RQ, B(1,IROW))
        IF (I .EQ. 1) GO TO 5
        IRM1 = IROW-1
        DO J = 1,IRM1
          JP1 = J + 1
          CALL GIVENS (B(J,J),B(J,IROW),C,S)
          CALL ROTATE (10-J,C,S,B(JP1,J),B(JP1,IROW))
        end do
        IF (I .LT. NEQ) GO TO 5
C
C TEST THE SYSTEM FOR ILL-CONDITIONING.
C
      DMIN = AMIN1( ABS(B(1,1)),ABS(B(2,2)),ABS(B(3,3)),
     .                  ABS(B(4,4)),ABS(B(5,5)),ABS(B(6,6)),
     .                  ABS(B(7,7)),ABS(B(8,8)),ABS(B(9,9)) )
      IF (DMIN*RQ .GE. DTOL) GO TO 13
      IF (NEQ .EQ. LMAX) GO TO 10
C
C INCREASE RQ AND ADD ANOTHER EQUATION TO THE SYSTEM TO
C   IMPROVE THE CONDITIONING.  THE NUMBER OF NPTS ELEMENTS
C   IS ALSO INCREASED IF NECESSARY.
C
    7      RSOLD = RS
      NEQ = NEQ + 1
      IF (NEQ .EQ. LMAX) GO TO 9
      IF (NEQ .EQ. LNP) GO TO 8
C
C   NEQ .LT. LNP
C
      NP = NPTS(NEQ+1)
      RS = (X(NP)-XK)**2 + (Y(NP)-YK)**2 + (Z(NP)-ZK)**2
      IF ( (RS-RSOLD)/RS .LT. RTOL ) GO TO 7
      RQ = SQRT(RS)
      GO TO 5
C
C   ADD AN ELEMENT TO NPTS.
C
    8      LNP = LNP + 1
      CALL GETNP3 (XK,YK,ZK,X,Y,Z,NNR,LCELL,LNEXT,XYZMN,
     .                 XYZDL, NP,RS)
      IF (NP .EQ. 0) GO TO 21
      NPTS(LNP) = NP
      IF ( (RS-RSOLD)/RS .LT. RTOL ) GO TO 7
      RQ = SQRT(RS)
      GO TO 5
C
    9      RQ = SQRT(1.1*RS)
      GO TO 5
C
C STABILIZE THE SYSTEM BY DAMPING SECOND PARTIALS -- ADD
C   MULTIPLES OF THE FIRST SIX UNIT VECTORS TO THE FIRST
C   SIX EQUATIONS.
C
   10 DO I = 1,6
        B(I,10) = SF
        IP1 = I + 1
        DO J = IP1,10
          B(J,10) = 0.
        end do
        DO J = I,9
          JP1 = J + 1
          CALL GIVENS (B(J,J),B(J,10),C,S)
          CALL ROTATE (10-J,C,S,B(JP1,J),B(JP1,10))
        end do
      end do
C
C TEST THE STABILIZED SYSTEM FOR ILL-CONDITIONING.
C
      DMIN = AMIN1( ABS(B(1,1)),ABS(B(2,2)),ABS(B(3,3)),
     .                  ABS(B(4,4)),ABS(B(5,5)),ABS(B(6,6)),
     .                  ABS(B(7,7)),ABS(B(8,8)),ABS(B(9,9)) )
      IF (DMIN*RQ .LT. DTOL) GO TO 22
C
C SOLVE THE 9 BY 9 TRIANGULAR SYSTEM FOR THE COEFFICIENTS
C
   13 DO IB = 1,9
        I = 10-IB
        T = 0.
        IF (I .EQ. 9) GO TO 15
        IP1 = I + 1
        DO J = IP1,9
          T = T + B(J,I)*A(J,K)
        end do
   15   continue
        A(I,K) = (B(10,I)-T)/B(I,I)
      end do
C
C SCALE THE COEFFICIENTS TO ADJUST FOR THE COLUMN SCALING.
C
      DO I = 1,6
        A(I,K) = A(I,K)/AVSQ
      end do
      A(7,K) = A(7,K)/AV
      A(8,K) = A(8,K)/AV
      A(9,K) = A(9,K)/AV
C
C UNMARK K AND THE ELEMENTS OF NPTS.
C
      LNEXT(K) = -LNEXT(K)
      DO I = 1,LNP
        NP = NPTS(I)
        LNEXT(NP) = -LNEXT(NP)
      end do

   18 CONTINUE
C
C NO ERRORS ENCOUNTERED.
C
      DO I = 1,3
        XYZMIN(I) = XYZMN(I)
        XYZDEL(I) = XYZDL(I)
      end do
      RMAX = SQRT(RSMX)
      IER = 0
      RETURN
C
C N, NQ, NW, OR NR IS OUT OF RANGE.
C
   20 IER = 1
      RETURN
C
C DUPLICATE NODES WERE ENCOUNTERED BY GETNP3.
C
   21 IER = 2
      RETURN
C
C NO UNIQUE SOLUTION DUE TO COLLINEAR NODES.
C
   22 DO I = 1,3
        XYZMIN(I) = XYZMN(I)
        XYZDEL(I) = XYZDL(I)
      end do
      IER = 3
      RETURN
      END
      FUNCTION QS3VAL (PX,PY,PZ,N,X,Y,Z,F,NR,LCELL,LNEXT,
     .                   XYZMIN,XYZDEL,RMAX,RSQ,A)
      INTEGER N, NR, LCELL(NR,NR,NR), LNEXT(N)
      REAL    PX, PY, PZ, X(N), Y(N), Z(N), F(N), XYZMIN(3),
     .            XYZDEL(3), RMAX, RSQ(N), A(9,N)
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C                                        10/28/87
C
C   THIS FUNCTION RETURNS THE VALUE Q(PX,PY,PZ) WHERE Q IS
C THE WEIGHTED SUM OF QUADRATIC NODAL FUNCTIONS DEFINED IN
C SUBROUTINE QSHEP3.  QS3GRD MAY BE CALLED TO COMPUTE A
C GRADIENT OF Q ALONG WITH THE VALUE, AND/OR TO TEST FOR
C ERRORS.
C
C ON INPUT --
C
C      PX,PY,PZ = CARTESIAN COORDINATES OF THE POINT P AT
C               WHICH Q IS TO BE EVALUATED.
C
C      N = NUMBER OF NODES AND DATA VALUES DEFINING Q.
C          N .GE. 10.
C
C      X,Y,Z,F = ARRAYS OF LENGTH N CONTAINING THE NODES
C              AND DATA VALUES INTERPOLATED BY Q.
C
C      NR = NUMBER OF ROWS, COLUMNS, AND PLANES IN THE CELL
C           GRID.  REFER TO STORE3.  NR .GE. 1.
C
C      LCELL = NR BY NR BY NR ARRAY OF NODAL INDICES ASSO-
C            CIATED WITH CELLS.  REFER TO STORE3.
C
C      LNEXT = ARRAY OF LENGTH N CONTAINING NEXT-NODE INDI-
C            CES.  REFER TO STORE3.
C
C      XYZMIN,XYZDEL = ARRAYS OF LENGTH 3 CONTAINING MINI-
C                  MUM NODAL COORDINATES AND CELL DIM-
C                  ENSIONS, RESPECTIVELY.      XYZDEL
C                  ELEMENTS MUST BE POSITIVE.  REFER TO
C                  STORE3.
C
C      RMAX = SQUARE ROOT OF THE LARGEST ELEMENT IN RSQ --
C             MAXIMUM RADIUS.
C
C      RSQ = ARRAY OF LENGTH N CONTAINING THE SQUARED RADII
C            WHICH ENTER INTO THE WEIGHTS DEFINING Q.
C
C      A = 9 BY N ARRAY CONTAINING THE COEFFICIENTS FOR THE
C          NODAL FUNCTIONS DEFINING Q.
C
C   INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION.      THE
C PARAMETERS OTHER THAN PX, PY AND PZ SHOULD BE INPUT UNAL-
C TERED FROM THEIR VALUES ON OUTPUT FROM QSHEP3.  THIS FUNC-
C TION SHOULD NOT BE CALLED IF A NONZERO ERROR FLAG WAS
C RETURNED BY QSHEP3.
C
C ON OUTPUT --
C
C      QS3VAL = FUNCTION VALUE Q(PX,PY,PZ) UNLESS N, NR,
C             XYZDEL, OR RMAX IS INVALID, IN WHICH CASE
C             NO VALUE IS RETURNED.
C
C MODULES REQUIRED BY QS3VAL -- NONE
C
C INTRINSIC FUNCTIONS CALLED BY QS3VAL -- IFIX, SQRT
C
C***********************************************************
C
      qs3val = 0.0

      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 .LT. 10  .OR.  NR .LT. 1  .OR.  DX .LE. 0.
     .        .OR.      DY .LE. 0.  .OR.  DZ .LE. 0.  .OR.
     .        RMAX .LT. 0.) RETURN
C
C SET IMIN, IMAX, JMIN, JMAX, KMIN, AND KMAX TO CELL INDICES
C   DEFINING THE RANGE OF THE SEARCH FOR NODES WHOSE RADII
C   INCLUDE P.      THE CELLS WHICH MUST BE SEARCHED ARE THOSE
C   INTERSECTED BY (OR CONTAINED IN) A SPHERE OF RADIUS RMAX
C   CENTERED AT P.
C
      IMIN = IFIX((XP-XMIN-RMAX)/DX) + 1
      IMAX = IFIX((XP-XMIN+RMAX)/DX) + 1
      IF (IMIN .LT. 1) IMIN = 1
      IF (IMAX .GT. NR) IMAX = NR
      JMIN = IFIX((YP-YMIN-RMAX)/DY) + 1
      JMAX = IFIX((YP-YMIN+RMAX)/DY) + 1
      IF (JMIN .LT. 1) JMIN = 1
      IF (JMAX .GT. NR) JMAX = NR
      KMIN = IFIX((ZP-ZMIN-RMAX)/DZ) + 1
      KMAX = IFIX((ZP-ZMIN+RMAX)/DZ) + 1
      IF (KMIN .LT. 1) KMIN = 1
      IF (KMAX .GT. NR) KMAX = NR
C
C THE FOLLOWING IS A TEST FOR NO CELLS WITHIN THE SPHERE
C   OF RADIUS RMAX.
C
      IF (IMIN .GT. IMAX  .OR.      JMIN .GT. JMAX      .OR.
     .        KMIN .GT. KMAX) GO TO 5
C
C ACCUMULATE WEIGHT VALUES IN SW AND WEIGHTED NODAL FUNCTION
C   VALUES IN SWQ.  THE WEIGHTS ARE W(L) = ((R-D)+/(R*D))**2
C   FOR R**2 = RSQ(L) AND D = DISTANCE BETWEEN P AND NODE L.
C
      SW = 0.
      SWQ = 0.
C
C OUTER LOOP ON CELLS (I,J,K).
C
      DO K = KMIN,KMAX
        DO J = JMIN,JMAX
          DO I = IMIN,IMAX
            L = LCELL(I,J,K)
            IF (L .EQ. 0) GO TO 3
C
C INNER LOOP ON NODES L.
C
    1       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 .GE. RS) GO TO 2
            IF (DS .EQ. 0.) GO TO 4
            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) )
C
C BOTTOM OF LOOP ON NODES IN CELL (I,J,K).
C
    2       LP = L
            L = LNEXT(LP)
            IF (L .NE. LP) GO TO 1
    3       CONTINUE
          end do
        end do
      end do
C
C SW = 0 IFF P IS NOT WITHIN THE RADIUS R(L) FOR ANY NODE L.
C
      IF (SW .EQ. 0.) GO TO 5
      QS3VAL = SWQ/SW
      RETURN
C
C (PX,PY,PZ) = (X(L),Y(L),Z(L))
C
    4 QS3VAL = F(L)
      RETURN
C
C ALL WEIGHTS ARE 0 AT P.
C
    5 QS3VAL = 0.
      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)
      INTEGER N, NR, LCELL(NR,NR,NR), LNEXT(N), IER
      REAL    PX, PY, PZ, X(N), Y(N), Z(N), F(N), XYZMIN(3),
     .            XYZDEL(3), RMAX, RSQ(N), A(9,N), Q, QX, QY, QZ
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C                                        10/28/87
C
C   THIS SUBROUTINE COMPUTES THE VALUE AND GRADIENT AT
C (PX,PY,PZ) OF THE INTERPOLATORY FUNCTION Q DEFINED IN SUB-
C ROUTINE QSHEP3.  Q(X,Y,Z) IS A WEIGHTED SUM OF QUADRATIC
C NODAL FUNCTIONS.
C
C ON INPUT --
C
C      PX,PY,PZ = CARTESIAN COORDINATES OF THE POINT P AT
C               WHICH Q AND ITS PARTIALS ARE TO BE EVAL-
C               UATED.
C
C      N = NUMBER OF NODES AND DATA VALUES DEFINING Q.
C          N .GE. 10.
C
C      X,Y,Z,F = ARRAYS OF LENGTH N CONTAINING THE NODES
C              AND DATA VALUES INTERPOLATED BY Q.
C
C      NR = NUMBER OF ROWS, COLUMNS AND PLANES IN THE CELL
C           GRID.  REFER TO STORE3.  NR .GE. 1.
C
C      LCELL = NR BY NR BY NR ARRAY OF NODAL INDICES ASSO-
C            CIATED WITH CELLS.  REFER TO STORE3.
C
C      LNEXT = ARRAY OF LENGTH N CONTAINING NEXT-NODE INDI-
C            CES.  REFER TO STORE3.
C
C      XYZMIN,XYZDEL = ARRAYS OF LENGTH 3 CONTAINING MINI-
C                  MUM NODAL COORDINATES AND CELL DIM-
C                  ENSIONS, RESPECTIVELY.      XYZDEL
C                  ELEMENTS MUST BE POSITIVE.  REFER TO
C                  STORE3.
C
C      RMAX = SQUARE ROOT OF THE LARGEST ELEMENT IN RSQ --
C             MAXIMUM RADIUS.
C
C      RSQ = ARRAY OF LENGTH N CONTAINING THE SQUARED RADII
C            WHICH ENTER INTO THE WEIGHTS DEFINING Q.
C
C      A = 9 BY N ARRAY CONTAINING THE COEFFICIENTS FOR THE
C          NODAL FUNCTIONS DEFINING Q.
C
C   INPUT PARAMETERS ARE NOT ALTERED BY THIS SUBROUTINE.
C THE PARAMETERS OTHER THAN PX, PY, AND PZ SHOULD BE INPUT
C UNALTERED FROM THEIR VALUES ON OUTPUT FROM QSHEP3.  THIS
C SUBROUTINE SHOULD NOT BE CALLED IF A NONZERO ERROR FLAG
C WAS RETURNED BY QSHEP3.
C
C ON OUTPUT --
C
C      Q = VALUE OF Q AT (PX,PY,PZ) UNLESS IER .EQ. 1, IN
C          WHICH CASE NO VALUES ARE RETURNED.
C
C      QX,QY,QZ = FIRST PARTIAL DERIVATIVES OF Q AT
C               (PX,PY,PZ) UNLESS IER .EQ. 1.
C
C      IER = ERROR INDICATOR
C            IER = 0 IF NO ERRORS WERE ENCOUNTERED.
C            IER = 1 IF N, NR, XYZDEL, OR RMAX IS INVALID.
C            IER = 2 IF NO ERRORS WERE ENCOUNTERED BUT
C                  (PX,PY,PZ) IS NOT WITHIN THE RADIUS
C                  R(K) FOR ANY NODE K (AND THUS Q = QX =
C                  QY = QZ = 0).
C
C MODULES REQUIRED BY QS3GRD -- NONE
C
C INTRINSIC FUNCTIONS CALLED BY QS3GRD -- IFIX, SQRT
C
C***********************************************************
C
      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 .LT. 10  .OR.  NR .LT. 1  .OR.  DX .LE. 0.
     .        .OR.      DY .LE. 0.  .OR.  DZ .LE. 0.  .OR.
     .        RMAX .LT. 0.) GO TO 5
C
C SET IMIN, IMAX, JMIN, JMAX, KMIN, AND KMAX TO CELL INDICES
C   DEFINING THE RANGE OF THE SEARCH FOR NODES WHOSE RADII
C   INCLUDE P.      THE CELLS WHICH MUST BE SEARCHED ARE THOSE
C   INTERSECTED BY (OR CONTAINED IN) A SPHERE OF RADIUS RMAX
C   CENTERED AT P.
C
      IMIN = IFIX((XP-XMIN-RMAX)/DX) + 1
      IMAX = IFIX((XP-XMIN+RMAX)/DX) + 1
      IF (IMIN .LT. 1) IMIN = 1
      IF (IMAX .GT. NR) IMAX = NR
      JMIN = IFIX((YP-YMIN-RMAX)/DY) + 1
      JMAX = IFIX((YP-YMIN+RMAX)/DY) + 1
      IF (JMIN .LT. 1) JMIN = 1
      IF (JMAX .GT. NR) JMAX = NR
      KMIN = IFIX((ZP-ZMIN-RMAX)/DZ) + 1
      KMAX = IFIX((ZP-ZMIN+RMAX)/DZ) + 1
      IF (KMIN .LT. 1) KMIN = 1
      IF (KMAX .GT. NR) KMAX = NR
C
C THE FOLLOWING IS A TEST FOR NO CELLS WITHIN THE SPHERE
C   OF RADIUS RMAX.
C
      IF (IMIN .GT. IMAX  .OR.      JMIN .GT. JMAX      .OR.
     .        KMIN .GT. KMAX) GO TO 6
C
C Q = SWQ/SW = SUM(W(L)*Q(L))/SUM(W(L)) WHERE THE SUM IS
C   FROM L = 1 TO N, Q(L) IS THE QUADRATIC NODAL FUNCTION,
C   AND W(L) = ((R-D)+/(R*D))**2 FOR RADIUS R(L) AND DIST-
C   ANCE D(L).      THUS
C
C       QX = (SWQX*SW - SWQ*SWX)/SW**2
C       QY = (SWQY*SW - SWQ*SWY)/SW**2
C       QZ = (SWQZ*SW - SWQ*SWZ)/SW**2
C
C   WHERE SWQX AND SWX ARE PARTIAL DERIVATIVES WITH RESPECT
C   TO X OF SWQ AND SW, RESPECTIVELY.  SWQY, SWY, SWQZ, AND
C   SWZ ARE DEFINED SIMILARLY.
C
      SW = 0.
      SWX = 0.
      SWY = 0.
      SWZ = 0.
      SWQ = 0.
      SWQX = 0.
      SWQY = 0.
      SWQZ = 0.
C
C OUTER LOOP ON CELLS (I,J,K).
C
      DO K = KMIN,KMAX
        DO J = JMIN,JMAX
          DO I = IMIN,IMAX
            L = LCELL(I,J,K)
            IF (L .EQ. 0) GO TO 3
C
C INNER LOOP ON NODES L.
C
    1       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 .GE. RS) GO TO 2
            IF (DS .EQ. 0.) GO TO 4
            RDS = RS*DS
            RD = SQRT(RDS)
            W = (RS+DS-RD-RD)/RDS
            T = 2.*(RD-RS)/(DS*RDS)
            WX = DELX*T
            WY = DELY*T
            WZ = DELZ*T
            QLX = 2.*A(1,L)*DELX + A(2,L)*DELY + A(4,L)*DELZ
            QLY = A(2,L)*DELX + 2.*A(3,L)*DELY + A(5,L)*DELZ
            QLZ = A(4,L)*DELX + A(5,L)*DELY + 2.*A(6,L)*DELZ
            QL = (QLX*DELX + QLY*DELY + QLZ*DELZ)/2. +
     .      A(7,L)*DELX + A(8,L)*DELY + A(9,L)*DELZ +
     .        F(L)
            QLX = QLX + A(7,L)
            QLY = QLY + A(8,L)
            QLZ = QLZ + A(9,L)
            SW = SW + W
            SWX = SWX + WX
            SWY = SWY + WY
            SWZ = SWZ + WZ
            SWQ = SWQ + W*QL
            SWQX = SWQX + WX*QL + W*QLX
            SWQY = SWQY + WY*QL + W*QLY
            SWQZ = SWQZ + WZ*QL + W*QLZ
C
C BOTTOM OF LOOP ON NODES IN CELL (I,J,K).
C
    2       LP = L
            L = LNEXT(LP)
            IF (L .NE. LP) GO TO 1
    3       CONTINUE
          end do
        end do
      end do
C
C SW = 0 IFF P IS NOT WITHIN THE RADIUS R(L) FOR ANY NODE L.
C
      IF (SW .EQ. 0.) GO TO 6
      Q = SWQ/SW
      SWS = SW*SW
      QX = (SWQX*SW - SWQ*SWX)/SWS
      QY = (SWQY*SW - SWQ*SWY)/SWS
      QZ = (SWQZ*SW - SWQ*SWZ)/SWS
      IER = 0
      RETURN
C
C (PX,PY,PZ) = (X(L),Y(L),Z(L))
C
    4 Q = F(L)
      QX = A(7,L)
      QY = A(8,L)
      QZ = A(9,L)
      IER = 0
      RETURN
C
C INVALID INPUT PARAMETER.
C
    5 IER = 1
      RETURN
C
C NO CELLS CONTAIN A POINT WITHIN RMAX OF P, OR
C   SW = 0 AND THUS DS .GE. RSQ(L) FOR ALL L.
C
    6 Q = 0.
      QX = 0.
      QY = 0.
      QZ = 0.
      IER = 2
      RETURN
      END
      SUBROUTINE GETNP3 (PX,PY,PZ,X,Y,Z,NR,LCELL,LNEXT,
     .                   XYZMIN,XYZDEL, NP,DSQ)
      INTEGER NR, LCELL(NR,NR,NR), LNEXT(*), NP
      REAL PX, PY, PZ, X(*), Y(*), Z(*), XYZMIN(3),
     .            XYZDEL(3), DSQ
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C
C   GIVEN A SET OF N NODES AND THE DATA STRUCTURE DEFINED IN
C SUBROUTINE STORE3, THIS SUBROUTINE USES THE CELL METHOD TO
C FIND THE CLOSEST UNMARKED NODE NP TO A SPECIFIED POINT P.
C NP IS THEN MARKED BY SETTING LNEXT(NP) TO -LNEXT(NP).  (A
C NODE IS MARKED IF AND ONLY IF THE CORRESPONDING LNEXT ELE-
C MENT IS NEGATIVE.  THE ABSOLUTE VALUES OF LNEXT ELEMENTS,
C HOWEVER, MUST BE PRESERVED.)      THUS, THE CLOSEST M NODES TO
C P MAY BE DETERMINED BY A SEQUENCE OF M CALLS TO THIS ROU-
C TINE.  NOTE THAT IF THE NEAREST NEIGHBOR TO NODE K IS TO
C BE DETERMINED (PX = X(K), PY = Y(K), AND PZ = Z(K)), THEN
C K SHOULD BE MARKED BEFORE THE CALL TO THIS ROUTINE.
C   THE SEARCH IS BEGUN IN THE CELL CONTAINING (OR CLOSEST
C TO) P AND PROCEEDS OUTWARD IN BOX-SHAPED LAYERS UNTIL ALL
C CELLS WHICH CONTAIN POINTS WITHIN DISTANCE R OF P HAVE
C BEEN SEARCHED, WHERE R IS THE DISTANCE FROM P TO THE FIRST
C UNMARKED NODE ENCOUNTERED (INFINITE IF NO UNMARKED NODES
C ARE PRESENT).
C
C ON INPUT --
C
C      PX,PY,PZ = CARTESIAN COORDINATES OF THE POINT P
C               WHOSE NEAREST UNMARKED NEIGHBOR IS TO BE
C               FOUND.
C
C      X,Y,Z = ARRAYS OF LENGTH N, FOR N .GE. 2, CONTAINING
C            THE CARTESIAN COORDINATES OF THE NODES.
C
C      NR = NUMBER OF ROWS, COLUMNS, AND PLANES IN THE CELL
C           GRID.  NR .GE. 1.
C
C      LCELL = NR BY NR BY NR ARRAY OF NODAL INDICES ASSO-
C            CIATED WITH CELLS.
C
C      LNEXT = ARRAY OF LENGTH N CONTAINING NEXT-NODE INDI-
C            CES (OR THEIR NEGATIVES).
C
C      XYZMIN,XYZDEL = ARRAYS OF LENGTH 3 CONTAINING MINI-
C                  MUM NODAL COORDINATES AND CELL DIM-
C                  ENSIONS, RESPECTIVELY.      XYZDEL
C                  ELEMENTS MUST BE POSITIVE.
C
C   INPUT PARAMETERS OTHER THAN LNEXT ARE NOT ALTERED BY
C THIS ROUTINE.  WITH THE EXCEPTION OF (PX,PY,PZ) AND THE
C SIGNS OF LNEXT ELEMENTS, THESE PARAMETERS SHOULD BE UNAL-
C TERED FROM THEIR VALUES ON OUTPUT FROM SUBROUTINE STORE3.
C
C ON OUTPUT --
C
C      NP = INDEX (FOR X, Y, AND Z) OF THE NEAREST UNMARKED
C           NODE TO P, OR 0 IF ALL NODES ARE MARKED OR NR
C           .LT. 1 OR AN ELEMENT OF XYZDEL IS NOT POSITIVE.
C           LNEXT(NP) .LT. 0 IF NP .NE. 0.
C
C      DSQ = SQUARED EUCLIDEAN DISTANCE BETWEEN P AND NODE
C            NP, OR 0 IF NP = 0.
C
C MODULES REQUIRED BY GETNP3 -- NONE
C
C INTRINSIC FUNCTIONS CALLED BY GETNP3 -- IABS, IFIX, SQRT
C
C***********************************************************
C
      LOGICAL FIRST
      XP = PX
      YP = PY
      ZP = PZ
      DX = XYZDEL(1)
      DY = XYZDEL(2)
      DZ = XYZDEL(3)
C
C TEST FOR INVALID INPUT PARAMETERS.
C
      IF (NR .LT. 1  .OR.  DX .LE. 0.  .OR.  DY .LE. 0.
     .        .OR.      DZ .LE. 0.) GO TO 10
C
C INITIALIZE PARAMETERS --
C
C   FIRST = TRUE IFF THE FIRST UNMARKED NODE HAS YET TO BE
C          ENCOUNTERED,
C   IMIN,...,KMAX = CELL INDICES DEFINING THE RANGE OF THE
C                SEARCH,
C   DELX,DELY,DELZ = PX-XYZMIN(1), PY-XYZMIN(2), AND
C                 PZ-XYZMIN(3),
C   I0,J0,K0 = CELL CONTAINING OR CLOSEST TO P,
C   I1,...,K2 = CELL INDICES OF THE LAYER WHOSE INTERSECTION
C            WITH THE RANGE DEFINED BY IMIN,...,KMAX IS
C            CURRENTLY BEING SEARCHED.
C
      FIRST = .TRUE.
      IMIN = 1
      IMAX = NR
      JMIN = 1
      JMAX = NR
      KMIN = 1
      KMAX = NR
      DELX = XP - XYZMIN(1)
      DELY = YP - XYZMIN(2)
      DELZ = ZP - XYZMIN(3)
      I0 = IFIX(DELX/DX) + 1
      IF (I0 .LT. 1) I0 = 1
      IF (I0 .GT. NR) I0 = NR
      J0 = IFIX(DELY/DY) + 1
      IF (J0 .LT. 1) J0 = 1
      IF (J0 .GT. NR) J0 = NR
      K0 = IFIX(DELZ/DZ) + 1
      IF (K0 .LT. 1) K0 = 1
      IF (K0 .GT. NR) K0 = NR
      I1 = I0
      I2 = I0
      J1 = J0
      J2 = J0
      K1 = K0
      K2 = K0
C
C OUTER LOOP ON LAYERS, INNER LOOP ON LAYER CELLS, EXCLUDING
C   THOSE OUTSIDE THE RANGE (IMIN,IMAX) X (JMIN,JMAX) X
C   (KMIN,KMAX).
C
    1 DO 7 K = K1,K2
      IF (K .GT. KMAX) GO TO 8
      IF (K .LT. KMIN) GO TO 7
      DO 6 J = J1,J2
        IF (J .GT. JMAX) GO TO 7
        IF (J .LT. JMIN) GO TO 6
        DO 5 I = I1,I2
          IF (I .GT. IMAX) GO TO 6
          IF (I .LT. IMIN) GO TO 5
          IF (K .NE. K1  .AND.  K .NE. K2  .AND.
     .            J .NE. J1  .AND.  J .NE. J2  .AND.
     .            I .NE. I1  .AND.  I .NE. I2) GO TO 5
C
C SEARCH CELL (I,J,K) FOR UNMARKED NODES L.
C
          L = LCELL(I,J,K)
          IF (L .EQ. 0) GO TO 5
C
C   LOOP ON NODES IN CELL (I,J,K).
C
    2          LN = LNEXT(L)
          IF (LN .LT. 0) GO TO 4
C
C   NODE L IS NOT MARKED.
C
          RSQ = (X(L)-XP)**2 + (Y(L)-YP)**2 + (Z(L)-ZP)**2
          IF (.NOT. FIRST) GO TO 3
C
C   NODE L IS THE FIRST UNMARKED NEIGHBOR OF P ENCOUNTERED.
C     INITIALIZE LMIN TO THE CURRENT CANDIDATE FOR NP, AND
C     RSMIN TO THE SQUARED DISTANCE FROM P TO LMIN.  IMIN,
C     IMAX, JMIN, JMAX, KMIN, AND KMAX ARE UPDATED TO DEFINE
C     THE SMALLEST RECTANGLE CONTAINING A SPHERE OF RADIUS
C     R = SQRT(RSMIN) CENTERED AT P, AND CONTAINED IN (1,NR)
C     X (1,NR) X (1,NR) (EXCEPT THAT, IF P IS OUTSIDE THE
C     BOX DEFINED BY THE NODES, IT IS POSSIBLE THAT IMIN
C     .GT. NR OR IMAX .LT. 1, ETC.).  FIRST IS RESET TO
C     FALSE.
C
          LMIN = L
          RSMIN = RSQ
          R = SQRT(RSMIN)
          IMIN = IFIX((DELX-R)/DX) + 1
          IF (IMIN .LT. 1) IMIN = 1
          IMAX = IFIX((DELX+R)/DX) + 1
          IF (IMAX .GT. NR) IMAX = NR
          JMIN = IFIX((DELY-R)/DY) + 1
          IF (JMIN .LT. 1) JMIN = 1
          JMAX = IFIX((DELY+R)/DY) + 1
          IF (JMAX .GT. NR) JMAX = NR
          KMIN = IFIX((DELZ-R)/DZ) + 1
          IF (KMIN .LT. 1) KMIN = 1
          KMAX = IFIX((DELZ+R)/DZ) + 1
          IF (KMAX .GT. NR) KMAX = NR
          FIRST = .FALSE.
          GO TO 4
C
C   TEST FOR NODE L CLOSER THAN LMIN TO P.
C
    3          IF (RSQ .GE. RSMIN) GO TO 4
C
C   UPDATE LMIN AND RSMIN.
C
          LMIN = L
          RSMIN = RSQ
C
C   TEST FOR TERMINATION OF LOOP ON NODES IN CELL (I,J,K).
C
    4          IF (IABS(LN) .EQ. L) GO TO 5
          L = IABS(LN)
          GO TO 2
    5          CONTINUE
    6        CONTINUE
    7      CONTINUE
C
C TEST FOR TERMINATION OF LOOP ON CELL LAYERS.
C
    8 IF (I1 .LE. IMIN      .AND.  I2 .GE. IMAX  .AND.
     .        J1 .LE. JMIN      .AND.  J2 .GE. JMAX  .AND.
     .        K1 .LE. KMIN      .AND.  K2 .GE. KMAX) GO TO 9
      I1 = I1 - 1
      I2 = I2 + 1
      J1 = J1 - 1
      J2 = J2 + 1
      K1 = K1 - 1
      K2 = K2 + 1
      GO TO 1
C
C UNLESS NO UNMARKED NODES WERE ENCOUNTERED, LMIN IS THE
C   CLOSEST UNMARKED NODE TO P.
C
    9 IF (FIRST) GO TO 10
      NP = LMIN
      DSQ = RSMIN
      LNEXT(LMIN) = -LNEXT(LMIN)
      RETURN
C
C ERROR -- NR OR XYZDEL IS INVALID OR ALL NODES ARE MARKED.
C
   10 NP = 0
      DSQ = 0.
      RETURN
      END
      SUBROUTINE GIVENS ( A,B, C,S)
      REAL A, B, C, S
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C
C   THIS ROUTINE CONSTRUCTS THE GIVENS PLANE ROTATION --
C     ( C  S)
C G = (     ) WHERE C*C + S*S = 1 -- WHICH ZEROS THE SECOND
C     (-S  C)
C ENTRY OF THE 2-VECTOR (A B)-TRANSPOSE.  A CALL TO GIVENS
C IS NORMALLY FOLLOWED BY A CALL TO ROTATE WHICH APPLIES
C THE TRANSFORMATION TO A 2 BY N MATRIX.  THIS ROUTINE WAS
C TAKEN FROM LINPACK.
C
C ON INPUT --
C
C      A,B = COMPONENTS OF THE 2-VECTOR TO BE ROTATED.
C
C ON OUTPUT --
C
C      A = VALUE OVERWRITTEN BY R = +/-SQRT(A*A + B*B)
C
C      B = VALUE OVERWRITTEN BY A VALUE Z WHICH ALLOWS C
C          AND S TO BE RECOVERED AS FOLLOWS --
C            C = SQRT(1-Z*Z), S=Z     IF ABS(Z) .LE. 1.
C            C = 1/Z, S = SQRT(1-C*C) IF ABS(Z) .GT. 1.
C
C      C = +/-(A/R)
C
C      S = +/-(B/R)
C
C MODULES REQUIRED BY GIVENS -- NONE
C
C INTRINSIC FUNCTIONS CALLED BY GIVENS - ABS, SQRT
C
C***********************************************************
C
      REAL AA, BB, R, U, V
C
C LOCAL PARAMETERS --
C
C AA,BB = LOCAL COPIES OF A AND B
C R =        C*A + S*B = +/-SQRT(A*A+B*B)
C U,V =   VARIABLES USED TO SCALE A AND B FOR COMPUTING R
C
      AA = A
      BB = B
      IF (ABS(AA) .LE. ABS(BB)) GO TO 1
C
C ABS(A) .GT. ABS(B)
C
      U = AA + AA
      V = BB/U
      R = SQRT(.25 + V*V) * U
      C = AA/R
      S = V * (C + C)
C
C NOTE THAT R HAS THE SIGN OF A, C .GT. 0, AND S HAS
C   SIGN(A)*SIGN(B).
C
      B = S
      A = R
      RETURN
C
C ABS(A) .LE. ABS(B)
C
    1 IF (BB .EQ. 0.) GO TO 2
      U = BB + BB
      V = AA/U
C
C STORE R IN A.
C
      A = SQRT(.25 + V*V) * U
      S = BB/A
      C = V * (S + S)
C
C NOTE THAT R HAS THE SIGN OF B, S .GT. 0, AND C HAS
C   SIGN(A)*SIGN(B).
C
      B = 1.
      IF (C .NE. 0.) B = 1./C
      RETURN
C
C A = B = 0.
C
    2 C = 1.
      S = 0.
      RETURN
      END
      SUBROUTINE ROTATE (N,C,S, X,Y )
      INTEGER N
      REAL    C, S, X(N), Y(N)
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C
C                                   ( C  S)
C   THIS ROUTINE APPLIES THE GIVENS ROTATION (         ) TO THE
C                                   (-S  C)
C            (X(1) ... X(N))
C 2 BY N MATRIX (            ).
C            (Y(1) ... Y(N))
C
C ON INPUT --
C
C      N = NUMBER OF COLUMNS TO BE ROTATED.
C
C      C,S = ELEMENTS OF THE GIVENS ROTATION.      THESE MAY BE
C            DETERMINED BY SUBROUTINE GIVENS.
C
C      X,Y = ARRAYS OF LENGTH .GE. N CONTAINING THE VECTORS
C            TO BE ROTATED.
C
C PARAMETERS N, C, AND S ARE NOT ALTERED BY THIS ROUTINE.
C
C ON OUTPUT --
C
C      X,Y = ROTATED VECTORS.
C
C MODULES REQUIRED BY ROTATE -- NONE
C
C***********************************************************
C
      INTEGER I
      REAL    XI, YI
C
C LOCAL PARAMETERS --
C
C I =        DO-LOOP INDEX
C XI,YI = X(I), Y(I)
C
      IF (N .LE. 0 .OR. (C .EQ. 1. .AND. S .EQ. 0.)) RETURN
      DO 1 I = 1,N
      XI = X(I)
      YI = Y(I)
      X(I) = C*XI + S*YI
      Y(I) = -S*XI + C*YI
    1      CONTINUE
      RETURN
      END
      SUBROUTINE SETUP3 (XK,YK,ZK,FK,XI,YI,ZI,FI,S1,S2,
     .                   R, ROW)
      REAL XK, YK, ZK, FK, XI, YI, ZI, FI, S1, S2,
     .         R, ROW(10)
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C
C   THIS ROUTINE SETS UP THE I-TH ROW OF AN AUGMENTED RE-
C GRESSION MATRIX FOR A WEIGHTED LEAST-SQUARES FIT OF A
C QUADRATIC FUNCTION Q(X,Y,Z) TO A SET OF DATA VALUES F,
C WHERE Q(XK,YK,ZK) = FK.  THE FIRST 6 COLUMNS (QUADRATIC
C TERMS) ARE SCALED BY 1/S2, AND COLUMNS 7, 8, AND 9 (LINEAR
C TERMS) ARE SCALED BY 1/S1.  THE WEIGHT IS (R-D)/(R*D) IF
C R .GT. D, AND 0 IF R .LE. D, WHERE D IS THE DISTANCE BE-
C TWEEN NODES I AND K.
C
C ON INPUT --
C
C      XK,YK,ZK,FK = COORDINATES AND DATA VALUE AT NODE K
C                  (INTERPOLATED BY Q).
C
C      XI,YI,ZI,FI = COORDINATES AND DATA VALUE AT NODE I.
C
C      S1,S2 = RECIPROCALS OF THE SCALE FACTORS.
C
C      R = RADIUS OF INFLUENCE ABOUT NODE K DEFINING THE
C          WEIGHT.
C
C      ROW = ARRAY OF LENGTH 10.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C ON OUTPUT --
C
C      ROW = ARRAY CONTAINING A ROW OF THE AUGMENTED
C            REGRESSION MATRIX.
C
C MODULES REQUIRED BY SETUP3 -- NONE
C
C INTRINSIC FUNCTION CALLED BY SETUP3 -- SQRT
C
C***********************************************************
C
      INTEGER I
      REAL    DX, DY, DZ, DXSQ, DYSQ, DZSQ, D, W, W1, W2
C
C LOCAL PARAMETERS -
C
C I =       DO-LOOP INDEX
C DX =       XI - XK
C DY =       YI - YK
C DZ =       ZI - ZK
C DXSQ = DX*DX
C DYSQ = DY*DY
C DZSQ = DZ*DZ
C D =       DISTANCE BETWEEN NODES K AND I
C W =       WEIGHT ASSOCIATED WITH THE ROW
C W1 =       W/S1
C W2 =       W/S2
C
      DX = XI - XK
      DY = YI - YK
      DZ = ZI - ZK
      DXSQ = DX*DX
      DYSQ = DY*DY
      DZSQ = DZ*DZ
      D = SQRT(DXSQ + DYSQ + DZSQ)
      IF (D .LE. 0.  .OR.  D .GE. R) GO TO 1
      W = (R-D)/R/D
      W1 = W/S1
      W2 = W/S2
      ROW(1) = DXSQ*W2
      ROW(2) = DX*DY*W2
      ROW(3) = DYSQ*W2
      ROW(4) = DX*DZ*W2
      ROW(5) = DY*DZ*W2
      ROW(6) = DZSQ*W2
      ROW(7) = DX*W1
      ROW(8) = DY*W1
      ROW(9) = DZ*W1
      ROW(10) = (FI - FK)*W
      RETURN
C
C NODES K AND I COINCIDE OR NODE I IS OUTSIDE OF THE RADIUS
C   OF INFLUENCE.  SET ROW TO THE ZERO VECTOR.
C
    1 DO I = 1,10
        ROW(I) = 0.
      end do

      RETURN
      END
      SUBROUTINE STORE3 (N,X,Y,Z,NR, LCELL,LNEXT,XYZMIN,
     .                   XYZDEL,IER)
      INTEGER N, NR, LCELL(NR,NR,NR), LNEXT(N), IER
      REAL    X(N), Y(N), Z(N), XYZMIN(3), XYZDEL(3)
C
C***********************************************************
C
C                                    ROBERT RENKA
C                              UNIV. OF NORTH TEXAS
C                                    (817) 565-2767
C
C   GIVEN A SET OF N ARBITRARILY DISTRIBUTED NODES IN THREE-
C SPACE, THIS SUBROUTINE CREATES A DATA STRUCTURE FOR A
C CELL-BASED METHOD OF SOLVING CLOSEST-POINT PROBLEMS.      THE
C SMALLEST BOX CONTAINING THE NODES IS PARTITIONED INTO AN
C NR BY NR BY NR UNIFORM GRID OF CELLS, AND NODES ARE AS-
C SOCIATED WITH CELLS.      IN PARTICULAR, THE DATA STRUCTURE
C STORES THE INDICES OF THE NODES CONTAINED IN EACH CELL.
C FOR A UNIFORM RANDOM DISTRIBUTION OF NODES, THE NEAREST
C NODE TO AN ARBITRARY POINT CAN BE DETERMINED IN CONSTANT
C EXPECTED TIME.
C
C ON INPUT --
C
C      N = NUMBER OF NODES.  N .GE. 2.
C
C      X,Y,Z = ARRAYS OF LENGTH N CONTAINING THE CARTESIAN
C            COORDINATES OF THE NODES.
C
C      NR = NUMBER OF ROWS, COLUMNS, AND PLANES IN THE
C           GRID.  THE CELL DENSITY (AVERAGE NUMBER OF
C           NODES PER CELL) IS D = N/(NR**3).      A RECOMMEND-
C           ED VALUE, BASED ON EMPIRICAL EVIDENCE, IS D = 3
C           -- NR = (N/3)**(1/3).  NR .GE. 1.
C
C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C      LCELL = ARRAY OF LENGTH .GE. NR**3.
C
C      LNEXT = ARRAY OF LENGTH .GE. N.
C
C      XYZMIN,XYZDEL = ARRAYS OF LENGTH .GE. 3.
C
C ON OUTPUT --
C
C      LCELL = NR BY NR BY NR CELL ARRAY SUCH THAT
C            LCELL(I,J,K) CONTAINS THE INDEX (FOR X, Y,
C            AND Z) OF THE FIRST NODE (NODE WITH SMALLEST
C            INDEX) IN CELL (I,J,K), OR LCELL(I,J,K) = 0
C            IF NO NODES ARE CONTAINED IN THE CELL.      THE
C            CORNER OF CELL (I,J,K) WHICH IS FARTHEST
C            FROM THE BOX CORNER DEFINED BY XYZMIN HAS
C            COORDINATES (XMIN+I*DX,YMIN+J*DY,ZMIN+K*DZ),
C            WHERE (XMIN,YMIN,ZMIN) ARE THE ELEMENTS OF
C            XYZMIN.  LCELL IS NOT DEFINED IF IER .NE. 0.
C
C      LNEXT = ARRAY OF NEXT-NODE INDICES SUCH THAT
C            LNEXT(L) CONTAINS THE INDEX OF THE NEXT NODE
C            IN THE CELL WHICH CONTAINS NODE L, OR
C            LNEXT(L) = L IF L IS THE LAST NODE IN THE
C            CELL FOR L = 1,...,N.  (THE NODES CONTAINED
C            IN A CELL ARE ORDERED BY THEIR INDICES.)
C            IF, FOR EXAMPLE, CELL (I,J,K) CONTAINS NODES
C            2, 3, AND 5 (AND NO OTHERS), THEN
C            LCELL(I,J,K) = 2, LNEXT(2) = 3, LNEXT(3) =
C            5, AND LNEXT(5) = 5.  LNEXT IS NOT DEFINED
C            IF IER .NE. 0.
C
C      XYZMIN = ARRAY OF LENGTH 3 CONTAINING THE MINIMUM
C             NODAL COORDINATES XMIN, YMIN, AND ZMIN (IN
C             THAT ORDER) UNLESS IER = 1.  THE OPPOSITE
C             CORNER OF THE BOX DEFINED BY THE NODES IS
C             (XMIN+NR*DX,YMIN+NR*DY,ZMIN+NR*DZ).
C
C      XYZDEL = ARRAY OF LENGTH 3 CONTAINING THE DIMENSIONS
C             OF THE CELLS UNLESS IER = 1.  XYZDEL(1) =
C             (XMAX-XMIN)/NR, XYZDEL(2) = (YMAX-YMIN)/NR,
C             AND XYZDEL(3) = (ZMAX-ZMIN)/NR, WHERE XMIN,
C             XMAX, YMIN, YMAX, ZMIN, AND ZMAX ARE THE
C             EXTREMA OF X, Y, AND Z.
C
C      IER = ERROR INDICATOR --
C            IER = 0 IF NO ERRORS WERE ENCOUNTERED.
C            IER = 1 IF N .LT. 2 OR NR .LT. 1.
C            IER = 2 IF A COMPONENT OF XYZDEL IS NOT
C                  POSITIVE.
C
C MODULES REQUIRED BY STORE3 -- NONE
C
C INTRINSIC FUNCTIONS CALLED BY STORE3 -- FLOAT, IFIX
C
C***********************************************************
C
      NN = N
      NNR = NR
      IF (NN .LT. 2  .OR.  NNR .LT. 1) GO TO 4
C
C COMPUTE THE DIMENSIONS OF THE RECTANGLE CONTAINING THE
C   NODES.
C
      XMN = X(1)
      XMX = XMN
      YMN = Y(1)
      YMX = YMN
      ZMN = Z(1)
      ZMX = ZMN
      DO L = 2,NN
        IF (X(L) .LT. XMN) XMN = X(L)
        IF (X(L) .GT. XMX) XMX = X(L)
        IF (Y(L) .LT. YMN) YMN = Y(L)
        IF (Y(L) .GT. YMX) YMX = Y(L)
        IF (Z(L) .LT. ZMN) ZMN = Z(L)
        IF (Z(L) .GT. ZMX) ZMX = Z(L)
      end do
      XYZMIN(1) = XMN
      XYZMIN(2) = YMN
      XYZMIN(3) = ZMN
C
C COMPUTE CELL DIMENSIONS AND TEST FOR ZERO AREA.
C
      DELX = (XMX-XMN)/FLOAT(NNR)
      DELY = (YMX-YMN)/FLOAT(NNR)
      DELZ = (ZMX-ZMN)/FLOAT(NNR)
      XYZDEL(1) = DELX
      XYZDEL(2) = DELY
      XYZDEL(3) = DELZ
      IF (DELX .EQ. 0. .OR.  DELY .EQ. 0.  .OR.
     .        DELZ .EQ. 0.) GO TO 5
C
C INITIALIZE LCELL.
C
      DO K = 1,NNR
        DO J = 1,NNR
          DO I = 1,NNR
            LCELL(I,J,K) = 0
          end do
        end do
      end do
C
C LOOP ON NODES, STORING INDICES IN LCELL AND LNEXT.
C
      NP1 = NN + 1
      DO LL = 1,NN
        LB = NP1 - LL
        I = IFIX((X(LB)-XMN)/DELX) + 1
        IF (I .GT. NNR) I = NNR
        J = IFIX((Y(LB)-YMN)/DELY) + 1
        IF (J .GT. NNR) J = NNR
        K = IFIX((Z(LB)-ZMN)/DELZ) + 1
        IF (K .GT. NNR) K = NNR
        L = LCELL(I,J,K)
        LNEXT(LB) = L
        IF (L .EQ. 0) LNEXT(LB) = LB
        LCELL(I,J,K) = LB
      end do
C
C NO ERRORS ENCOUNTERED
C
      IER = 0
      RETURN
C
C INVALID INPUT PARAMETER
C
    4 IER = 1
      RETURN
C
C NONPOSITIVE XYZDEL COMPONENT
C
    5 IER = 2
      RETURN
      END
