C      ALGORITHM 661, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 14, NO. 2, P.151.
C
C			   QS3TEST
C
C   THIS PROGRAM TESTS THE SCATTERED DATA INTERPOLATION
C PACKAGE QSHEP3D BY PRINTING THE MAXIMUM ERRORS ASSOCIATED
C WITH INTERPOLATED VALUES AND GRADIENTS ON A 5 BY 5 BY 5
C UNIFORM GRID IN THE UNIT CUBE.  THE DATA SET CONSISTS
C OF 64 NODES WITH DATA VALUES TAKEN FROM A QUADRATIC FUNC-
C TION FOR WHICH THE METHOD IS EXACT.  THE RATIO OF MAXIMUM
C INTERPOLATION ERROR RELATIVE TO THE MACHINE PRECISION IS
C ALSO PRINTED.  THIS SHOULD BE O(1).  THE INTERPOLATED
C VALUES FROM QS3VAL AND QS3GRD ARE COMPARED FOR AGREEMENT.
C
      INTEGER LCELL(3,3,3), LNEXT(64)
      REAL X(64), Y(64), Z(64), F(64), RSQ(64), A(9,64),
     .	   XYZMIN(3), XYZDEL(3)
      REAL P(5)
C
C QSHEP3 PARAMETERS AND LOGICAL UNIT FOR OUTPUT
C
      DATA N/64/,  NQ/17/,  NW/32/,  NR/3/,  LOUT/6/
C
C QUADRATIC TEST FUNCTION AND PARTIAL DERIVATIVES
C
      FQ(XX,YY,ZZ) = ((XX + 2.*YY + 3.*ZZ)/6.)**2
      FX(XX,YY,ZZ) = (XX + 2.*YY + 3.*ZZ)/18.
      FY(XX,YY,ZZ) = (XX + 2.*YY + 3.*ZZ)/9.
      FZ(XX,YY,ZZ) = (XX + 2.*YY + 3.*ZZ)/6.
C
C GENERATE A 4 BY 4 BY 4 GRID OF NODES IN THE UNIT CUBE.
C
      L = 0
      DO 1 K = 1,4
	ZL = FLOAT(K-1)/3.
	DO 1 J = 1,4
	  YL = FLOAT(J-1)/3.
	  DO 1 I = 1,4
	    L = L + 1
	    X(L) = FLOAT(I-1)/3.
	    Y(L) = YL
    1	    Z(L) = ZL
C
C COMPUTE THE DATA VALUES.
C
      DO 2 L = 1,N
    2	F(L) = FQ(X(L),Y(L),Z(L))
C
C COMPUTE PARAMETERS DEFINING THE INTERPOLANT Q.
C
      CALL QSHEP3 (N,X,Y,Z,F,NQ,NW,NR, LCELL,LNEXT,XYZMIN,
     .		   XYZDEL,RMAX,RSQ,A,IER)
      IF (IER .NE. 0) GO TO 6
C
C GENERATE A 5 BY 5 BY 5 UNIFORM GRID OF INTERPOLATION
C   POINTS (P(I),P(J),P(K)) IN THE UNIT CUBE.  THE EIGHT
C   CORNERS COINCIDE WITH NODES.
C
      DO 3 I = 1,5
    3	P(I) = FLOAT(I-1)/4.
C
C COMPUTE THE MACHINE PRECISION EPS.
C
      EPS = 1.
    4 EPS = EPS/2.
      EP1 = EPS + 1.
      IF (STORE(EP1) .GT. 1.) GO TO 4
      EPS = EPS*2.
C
C COMPUTE INTERPOLATION ERRORS AND TEST FOR AGREEMENT IN THE
C   Q VALUES RETURNED BY QS3VAL AND QS3GRD.
C
      EQ = 0.
      EQX = 0.
      EQY = 0.
      EQZ = 0.
      DO 5 K = 1,5
	PZ = P(K)
	DO 5 J = 1,5
	  PY = P(J)
	  DO 5 I = 1,5
	    PX = P(I)
	    Q1 = QS3VAL (PX,PY,PZ,N,X,Y,Z,F,NR,LCELL,LNEXT,
     .			 XYZMIN,XYZDEL,RMAX,RSQ,A)
	    CALL QS3GRD (PX,PY,PZ,N,X,Y,Z,F,NR,LCELL,LNEXT,
     .			 XYZMIN,XYZDEL,RMAX,RSQ,A, Q,QX,QY,
     .			 QZ,IER)
	    IF (IER .NE. 0) GO TO 7
	    IF (ABS(Q1-Q) .GT. 3.*ABS(Q)*EPS) GO TO 8
	    EQ = AMAX1(EQ,ABS(FQ(PX,PY,PZ)-Q))
	    EQX = AMAX1(EQX,ABS(FX(PX,PY,PZ)-QX))
	    EQY = AMAX1(EQY,ABS(FY(PX,PY,PZ)-QY))
    5	    EQZ = AMAX1(EQZ,ABS(FZ(PX,PY,PZ)-QZ))
C
C PRINT ERRORS AND THE RATIO EQ/EPS.
C
      RQ = EQ/EPS
      WRITE (LOUT,100)
      WRITE (LOUT,110) EQ, RQ
      WRITE (LOUT,120) EQX
      WRITE (LOUT,130) EQY
      WRITE (LOUT,140) EQZ
      STOP
  100 FORMAT (///1H ,31HMAXIMUM ABSOLUTE ERRORS IN THE ,
     .	      25HINTERPOLANT Q AND PARTIAL/
     .	      1H ,32HDERIVATIVES (QX,QY,QZ) RELATIVE ,
     .	      24HTO MACHINE PRECISION EPS//
     .	      1H ,10X,8HFUNCTION,3X,9HMAX ERROR,3X,
     .	      13HMAX ERROR/EPS/)
  110 FORMAT (1H ,13X,1HQ,7X,E9.3,7X,F4.2)
  120 FORMAT (1H ,13X,2HQX,6X,E9.3)
  130 FORMAT (1H ,13X,2HQY,6X,E9.3)
  140 FORMAT (1H ,13X,2HQZ,6X,E9.3)
C
C ERROR IN QSHEP3
C
    6 WRITE (LOUT,200) IER
      STOP
  200 FORMAT (///1H ,28H*** ERROR IN QSHEP3 -- IER =,I2,
     .	      4H ***)
C
C ERROR IN QS3GRD
C
    7 WRITE (LOUT,210) IER
      STOP
  210 FORMAT (///1H ,28H*** ERROR IN QS3GRD -- IER =,I2,
     .	      4H ***)
C
C VALUES RETURNED BY QS3VAL AND QS3GRD DIFFER BY A RELATIVE
C   AMOUNT GREATER THAN 3*EPS.
C
    8 WRITE (LOUT,220) Q1, Q
      STOP
  220 FORMAT (///1H ,33H*** ERROR -- INTERPOLATED VALUES ,
     .	      37HQ1 (QS3VAL) AND Q2 (QS3GRD) DIFFER --//
     .	      1H ,5X,5HQ1 = ,E21.14,5X,5HQ2 = ,E21.14)
      END
      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 6 J = 1,IRM1
	    JP1 = J + 1
	    CALL GIVENS (B(J,J),B(J,IROW),C,S)
    6	    CALL ROTATE (10-J,C,S,B(JP1,J),B(JP1,IROW))
	  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 12 I = 1,6
	  B(I,10) = SF
	  IP1 = I + 1
	  DO 11 J = IP1,10
   11	    B(J,10) = 0.
	  DO 12 J = I,9
	    JP1 = J + 1
	    CALL GIVENS (B(J,J),B(J,10),C,S)
   12	    CALL ROTATE (10-J,C,S,B(JP1,J),B(JP1,10))
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 15 IB = 1,9
	  I = 10-IB
	  T = 0.
	  IF (I .EQ. 9) GO TO 15
	  IP1 = I + 1
	  DO 14 J = IP1,9
   14	    T = T + B(J,I)*A(J,K)
   15	  A(I,K) = (B(10,I)-T)/B(I,I)
C
C SCALE THE COEFFICIENTS TO ADJUST FOR THE COLUMN SCALING.
C
	DO 16 I = 1,6
   16	  A(I,K) = A(I,K)/AVSQ
	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 17 I = 1,LNP
	  NP = NPTS(I)
   17	  LNEXT(NP) = -LNEXT(NP)
   18	CONTINUE
C
C NO ERRORS ENCOUNTERED.
C
      DO 19 I = 1,3
	XYZMIN(I) = XYZMN(I)
   19	XYZDEL(I) = XYZDL(I)
      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 23 I = 1,3
	XYZMIN(I) = XYZMN(I)
   23	XYZDEL(I) = XYZDL(I)
      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
      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 3 K = KMIN,KMAX
	DO 3 J = JMIN,JMAX
	  DO 3 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
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 3 K = KMIN,KMAX
	DO 3 J = JMIN,JMAX
	  DO 3 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
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(1), NP
      REAL    PX, PY, PZ, X(1), Y(1), Z(1), 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 2 I = 1,10
    2	ROW(I) = 0.
      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 1 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)
    1	IF (Z(L) .GT. ZMX) ZMX = Z(L)
      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 2 K = 1,NNR
	DO 2 J = 1,NNR
	  DO 2 I = 1,NNR
    2	    LCELL(I,J,K) = 0
C
C LOOP ON NODES, STORING INDICES IN LCELL AND LNEXT.
C
      NP1 = NN + 1
      DO 3 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
    3	LCELL(I,J,K) = LB
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