 
      SUBROUTINE DB3INK(X,NX,Y,NY,Z,NZ,FCN,LDF1,LDF2,KX,KY,KZ,TX,TY,TZ,
     *  BCOEF,WORK,IFLAG)
C***BEGIN PROLOGUE  DB3INK
C***DATE WRITTEN   25 MAY 1982
C***REVISION DATE  25 MAY 1982
C***CATEGORY NO.  E1A
C***KEYWORDS  INTERPOLATION, THREE-DIMENSIONS, GRIDDED DATA, SPLINES,
C             PIECEWISE POLYNOMIALS
C***AUTHOR  BOISVERT, RONALD, NBS
C             SCIENTIFIC COMPUTING DIVISION
C             NATIONAL BUREAU OF STANDARDS
C             WASHINGTON, DC 20234
C***PURPOSE  DOUBLE PRECISION VERSION OF DB3INK
C            DB3INK DETERMINES A PIECEWISE POLYNOMIAL FUNCTION THAT
C            INTERPOLATES THREE-DIMENSIONAL GRIDDED DATA. USERS SPECIFY
C            THE POLYNOMIAL ORDER (DEGREE+1) OF THE INTERPOLANT AND
C            (OPTIONALLY) THE KNOT SEQUENCE.
C***DESCRIPTION
C
C   DB3INK determines the parameters of a  function  that  interpolates
C   the three-dimensional gridded data (X(i),Y(j),Z(k),FCN(i,j,k))  for
C   i=1,..,NX, j=1,..,NY, and k=1,..,NZ. The interpolating function and
C   its derivatives may  subsequently  be  evaluated  by  the  function
C   DB3VAL.
C
C   The interpolating  function  is  a  piecewise  polynomial  function
C   represented as a tensor product of one-dimensional  B-splines.  The
C   form of this function is
C
C                      NX   NY   NZ
C        S(x,y,z)  =  SUM  SUM  SUM  a   U (x) V (y) W (z)
C                     i=1  j=1  k=1   ij  i     j     k
C
C   where the functions U(i), V(j), and  W(k)  are  one-dimensional  B-
C   spline basis functions. The coefficients a(i,j) are chosen so that
C
C   S(X(i),Y(j),Z(k)) = FCN(i,j,k)  for i=1,..,NX, j=1,..,NY, k=1,..,NZ
C
C   Note that for fixed values of y  and  z  S(x,y,z)  is  a  piecewise
C   polynomial function of x alone, for fixed values of x and z  S(x,y,
C   z) is a piecewise polynomial function of y  alone,  and  for  fixed
C   values of x and y S(x,y,z)  is  a  function  of  z  alone.  In  one
C   dimension a piecewise polynomial may be created by  partitioning  a
C   given interval into subintervals and defining a distinct polynomial
C   piece on each one. The points where adjacent subintervals meet  are
C   called knots. Each of the functions U(i), V(j), and W(k) above is a
C   piecewise polynomial.
C
C   Users of DB3INK choose  the  order  (degree+1)  of  the  polynomial
C   pieces used to define the piecewise polynomial in each of the x, y,
C   and z directions (KX, KY, and KZ). Users also may define their  own
C   knot sequence in x, y, and z separately (TX, TY, and TZ). If IFLAG=
C   0, however, DB3INK will choose sequences of knots that result in  a
C   piecewise  polynomial  interpolant  with  KX-2  continuous  partial
C   derivatives in x, KY-2 continuous partial derivatives in y, and KZ-
C   2 continuous partial derivatives in z. (KX  knots  are  taken  near
C   each endpoint in x, not-a-knot end conditions  are  used,  and  the
C   remaining knots are placed at data points  if  KX  is  even  or  at
C   midpoints between data points if KX is odd. The y and z  directions
C   are treated similarly.)
C
C   After a call to DB3INK, all information  necessary  to  define  the
C   interpolating function are contained in the parameters NX, NY,  NZ,
C   KX, KY, KZ, TX, TY, TZ, and BCOEF. These quantities should  not  be
C   altered until after the last call of the evaluation routine DB3VAL.
C
C
C   I N P U T
C   ---------
C
C   X       Double precision 1D array (size NX)
C           Array of x abcissae. Must be strictly increasing.
C
C   NX      Integer scalar (.GE. 3)
C           Number of x abcissae.
C
C   Y       Double precision 1D array (size NY)
C           Array of y abcissae. Must be strictly increasing.
C
C   NY      Integer scalar (.GE. 3)
C           Number of y abcissae.
C
C   Z       Double precision 1D array (size NZ)
C           Array of z abcissae. Must be strictly increasing.
C
C   NZ      Integer scalar (.GE. 3)
C           Number of z abcissae.
C
C   FCN     Double precision 3D array (size LDF1 by LDF2 by NY)
C           Array of function values to interpolate. FCN(I,J,K) should
C           contain the function value at the point (X(I),Y(J),Z(K))
C
C   LDF1    Integer scalar (.GE. NX)
C           The actual first dimension of FCN used in the
C           calling program.
C
C   LDF2    Integer scalar (.GE. NY)
C           The actual second dimension of FCN used in the calling
C           program.
C
C   KX      Integer scalar (.GE. 2, .LT. NX)
C           The order of spline pieces in x.
C           (Order = polynomial degree + 1)
C
C   KY      Integer scalar (.GE. 2, .LT. NY)
C           The order of spline pieces in y.
C           (Order = polynomial degree + 1)
C
C   KZ      Integer scalar (.GE. 2, .LT. NZ)
C           The order of spline pieces in z.
C           (Order = polynomial degree + 1)
C
C
C   I N P U T   O R   O U T P U T
C   -----------------------------
C
C   TX      Double precision 1D array (size NX+KX)
C           The knots in the x direction for the spline interpolant.
C           If IFLAG=0 these are chosen by DB3INK.
C           If IFLAG=1 these are specified by the user.
C                      (Must be non-decreasing.)
C
C   TY      Double precision 1D array (size NY+KY)
C           The knots in the y direction for the spline interpolant.
C           If IFLAG=0 these are chosen by DB3INK.
C           If IFLAG=1 these are specified by the user.
C                      (Must be non-decreasing.)
C
C   TZ      Double precision 1D array (size NZ+KZ)
C           The knots in the z direction for the spline interpolant.
C           If IFLAG=0 these are chosen by DB3INK.
C           If IFLAG=1 these are specified by the user.
C                      (Must be non-decreasing.)
C
C
C   O U T P U T
C   -----------
C
C   BCOEF   Double precision 3D array (size NX by NY by NZ)
C           Array of coefficients of the B-spline interpolant.
C           This may be the same array as FCN.
C
C
C   M I S C E L L A N E O U S
C   -------------------------
C
C   WORK    Double precision 1D array (size NX*NY*NZ + max( 2*KX*(NX+1),
C                             2*KY*(NY+1), 2*KZ*(NZ+1) )
C           Array of working storage.
C
C   IFLAG   Integer scalar.
C           On input:  0 == knot sequence chosen by B2INK
C                      1 == knot sequence chosen by user.
C           On output: 1 == successful execution
C                      2 == IFLAG out of range
C                      3 == NX out of range
C                      4 == KX out of range
C                      5 == X not strictly increasing
C                      6 == TX not non-decreasing
C                      7 == NY out of range
C                      8 == KY out of range
C                      9 == Y not strictly increasing
C                     10 == TY not non-decreasing
C                     11 == NZ out of range
C                     12 == KZ out of range
C                     13 == Z not strictly increasing
C                     14 == TY not non-decreasing
C
C***REFERENCES  CARL DE BOOR, A PRACTICAL GUIDE TO SPLINES,
C                 SPRINGER-VERLAG, NEW YORK, 1978.
C               CARL DE BOOR, EFFICIENT COMPUTER MANIPULATION OF TENSOR
C                 PRODUCTS, ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C                 VOL. 5 (1979), PP. 173-182.
C***ROUTINES CALLED  DBTPCF,DBKNOT
C***END PROLOGUE  DB3INK
 
 
