C 
C Subroutines of the software package 'FITPACK' by A.K. Cline
C used to specify the model for the complete ray tracing algorithm.
C
C This file consists of the following parts:
C     (0) Auxiliary subroutine
C             SNHCSH
C             SNHCSH
C         common to all the following parts.
C     (1) The subroutines preparing the parameters necessary to compute
C         an interpolatory function:
C             CURVN1 (Hermite representation of 1-D function),
C             CURVB1 (B-spline representation of 1-D function),
C             SURFB1 (B-spline representation of 2-D function),
C             VAL3B1 (B-spline representation of 3-D function),
C             VGEN (auxiliary subroutine),
C             TERMS (auxiliary subroutine),
C             TRIDEC (auxiliary subroutine),
C             TRISOL (auxiliary subroutine).
C             CURVN1
C             CURVB1
C             SURFB1
C             VAL3B1
C             VGEN  
C             TERMS 
C             TRIDEC
C             TRISOL
C         Subroutines CURVN1 and CURVB1 are alternatives.
C     (2) The subroutines evaluating the value, first and second partial
C         derivatives of the interpolatory function at a given point:
C             CURV2D (Hermite representation of 1-D function),
C             CURVBD (B-spline representation of 1-D function),
C             SURFBD (B-spline representation of 2-D function),
C             VAL3BD (B-spline representation of 3-D function),
C             DSPLNZ (auxiliary subroutine),
C             INTRVL (auxiliary external function).
C             CURV2D
C             CURVBD
C             SURFBD
C             VAL3BD
C             DSPLNZ
C             INTRVL
C         Subroutines CURV2D and CURVBD are alternatives.
C
C Taken from:
C     FITPACK - A Software Package for Curve and Surface Fitting
C     Employing Splines under Tension
C     by Alan Kaylor Cline, Department of Computer Sciences,
C     The University of Texas at Austin, August 31, 1981.
C Note 1:
C     To conform with the FORTRAN77 standard, dummy array dimensions (1)
C     have been changed to (*), and subroutine TRISOL has been revised.
C Note 2:
C     Subroutines CURVB1 and CURVBD do not belong to the original
C     version of FITPACK.
C Note 3
C     The lines denoted by '*V' in the first two columns of file
C     'fit.for' calculate the model variations with respect to the model
C     parameters.
C     File 'fitv.for', intended for the model inversion, is created
C     from 'fit.for' by replacing each '*V' in the first two columns
C     by spaces using program 'clean.for'.  Subroutines VAR4 and VAR5
C     of file 'var.for' may then be called to handle the variations.
C Note 4:
C     To get the original versions of the subroutines SURFBD and VAL3BD,
C     the statement with 'CALL VAR2' must be removed from each of them.
C     The statements have been added by L.Klimes for the sake of inverse
C     modelling to the subroutines CURVBD, SURFBD, and VAL3BD.
C     The three appearances of the statements 'CALL VAR2' are denoted by
C     '*V' in the first 2 columns.  The three lines should be removed or
C     modified before compilation.
C
C=======================================================================
C
C Part 0:
C
C=======================================================================
C
C     
C
      SUBROUTINE SNHCSH (SINHM,COSHM,X,ISW)
C
      INTEGER ISW
      REAL SINHM,COSHM,X
C
C                            From FITPACK -- August 31, 1981
C                       Coded by A. K. Cline and R. J. Renka
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This subroutine returns approximations to
C       SINHM(X) = SINH(X)-X
C       COSHM(X) = COSH(X)-1
C and
C       COSHMM(X) = COSH(X)-1-X*X/2
C with relative error less than 6.16e-6
C
C On input--
C
C   X contains the value of the independent variable.
C
C   ISW indicates the function desired
C           = -1 if only SINHM is desired,
C           =  0 if both SINHM and coshm are desired,
C           =  1 if only COSHM is desired,
C           =  2 if only COSHMM is desired,
C           =  3 if both SINHM and COSHMM are desired.
C
C On output--
C
C   SINHM contains the value of SINHM(X) if ISW .LE. 0 or
C   ISW .EQ. 3 (SINHM is unaltered if ISW .EQ.1 or ISW .EQ.
C   2).
C
C   COSHM contains the value of COSHM(X) if ISW .EQ. 0 or
C   ISW .EQ. 1 and contains the value of COSHMM(X) if ISW
C   .GE. 2 (COSHM is unaltered if ISW .EQ. -1).
C
C And
C
C   X and ISW are unaltered.
C
C-----------------------------------------------------------
C
      DATA SP2/5.04850926418006E-04/,
     *     SP1/3.62841692246321E-02/,
     *     SQ1/-1.37157937097122E-02/
      DATA CP2/1.31625490355985E-03/,
     *     CP1/6.57866547762733E-02/,
     *     CQ1/-1.75465241841312E-02/
      DATA ZP2/1.40048186158693E-04/,
     *     ZP1/1.67309141907440E-02/,
     *     ZQ2/9.82154460147143E-05/,
     *     ZQ1/-1.66024148976133E-02/
      XX = X
      AX = ABS(XX)
      XS = XX*XX
      IF ((AX .GE. 2.20) .OR. (AX .GE. 1.17 .AND.
     *     ISW .NE. 2)) EXPX = EXP(AX)
C
C SINHM approximation
C
      IF (ISW .EQ. 1 .OR. ISW .EQ. 2) GO TO 2
      IF (AX .GE. 1.17) GO TO 1
      SINHM = (((SP2*XS+SP1)*XS+1.)*XS*XX)/((SQ1*XS+1.)*6.)
      GO TO 2
    1 SINHM = (EXPX-1./EXPX)/2.-AX
      IF (XX .LT. 0.) SINHM = -SINHM
C
C COSHM approximation
C
    2 IF (ISW .NE. 0 .AND. ISW .NE. 1) GO TO 4
      IF (AX .GE. 1.17) GO TO 3
      COSHM = (((CP2*XS+CP1)*XS+1.)*XS)/((CQ1*XS+1.)*2.)
      GO TO 4
    3 COSHM = (EXPX+1./EXPX)/2.-1.
C
C COSHMM approximation
C
    4 IF (ISW .LE. 1) RETURN
      IF (AX .GE. 2.20) GO TO 5
      COSHM = (((ZP2*XS+ZP1)*XS+1.)*XS*XS)/(((ZQ2*XS+ZQ1)*XS
     *        +1.)*24.)
      RETURN
    5 COSHM = (EXPX+1./EXPX)/2.-1.-XS/2.
      RETURN
      END
C
C=======================================================================
C
C Part 1:
C
C=======================================================================
C
C     
C
      SUBROUTINE CURVN1 (N,X,Y,YP,TEMP,SIGMA,IERR)
C
      INTEGER N,IERR
      REAL X(N),Y(N),YP(N),TEMP(N),SIGMA
C
C                            From FITPACK -- August 31, 1981
C                    Coded by a. K. Cline and s. E. Galinsky
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This subroutine determines the parameters necessary to
C compute a natural interpolatory spline under tension
C through a sequence of functional values. For actual
C computation of points on the curve it is necessary to call
C the function CURV2.
C
C On input--
C
C   N is the number of values to be interpolated (N.GE.2).
C
C   X is an array of the N increasing abscissae of the
C   functional values.
C
C   Y is an array of the N ordinates of the values, (i. e.
C   Y(K) is the functional value corresponding to X(K) ).
C
C   YP is an array of length at least N.
C
C   TEMP is an array of length at least N which is used for
C   scratch storage.
C
C And
C
C   SIGMA contains the tension factor. This value indicates
C   the curviness desired. If ABS(SIGMA) is nearly zero
C   (e.g. .001) the resulting curve is approximately a
C   cubic spline. If ABS(SIGMA) is large (e.g. 50.) the
C   resulting curve is nearly a polygonal line. If SIGMA
C   equals zero a cubic spline results.  A standard value
C   for SIGMA is approximately 1. In absolute value.
C
C On output--
C
C   YP contains the values of the second derivative of the
C   curve at the given nodes.
C
C   IERR contains an error flag,
C        = 0 for normal return,
C        = 1 if N is less than 2,
C        = 2 if X-values are not strictly increasing.
C
C And
C
C   N, X, Y, and SIGMA are unaltered.
C
C This subroutine references package modules SNHCSH.
C
C-----------------------------------------------------------
C
      NM1 = N-1
      NP1 = N+1
      IERR = 0
      IF (N .LE. 1) GO TO 4
      IF (X(N) .LE. X(1)) GO TO 5
C
C Denormalize tension factor
C
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
C
C Set up right hand side and tridiagonal system for YP and
C perform forward elimination
C
      DELX1 = X(2)-X(1)
      IF (DELX1 .LE. 0.) GO TO 5
      DX1 = (Y(2)-Y(1))/DELX1
      CALL TERMS (DIAG1,SDIAG1,SIGMAP,DELX1)
      SDIAG1 = 0.
      YP(1) = 0.
      TEMP(1) = 0.
      IF (N .EQ. 2) GO TO 2
      DO 1 I = 2,NM1
        DELX2 = X(I+1)-X(I)
        IF (DELX2 .LE. 0.) GO TO 5
        DX2 = (Y(I+1)-Y(I))/DELX2
        CALL TERMS (DIAG2,SDIAG2,SIGMAP,DELX2)
        DIAG = DIAG1+DIAG2-SDIAG1*TEMP(I-1)
        YP(I) = (DX2-DX1-SDIAG1*YP(I-1))/DIAG
        TEMP(I) = SDIAG2/DIAG
        DX1 = DX2
        DIAG1 = DIAG2
    1   SDIAG1 = SDIAG2
   2  YP(N) = 0.
      TEMP(N-1) = 0.
C
C Perform back substitution
C
      DO 3 I = 2,N
        IBAK = NP1-I
    3   YP(IBAK) = YP(IBAK)-TEMP(IBAK)*YP(IBAK+1)
      RETURN
C
C Too few points
C
    4 IERR = 1
      RETURN
C
C X-values not strictly increasing
C
    5 IERR = 2
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE CURVB1 (NX,X,W,C,VX,TEMP,SIGMA,IERR)
C
      INTEGER NX,IERR
      REAL X(NX),W(NX),C(NX),VX(5,NX),TEMP(*),SIGMA
C
C                                      Complement to FITPACK
C                                       by Alan Kaylor Cline
C                                   Coded -- October 9, 1986
C                                            by Ludek Klimes
C                                      Inst. Geol. Geotechn.
C                               Czechosl. Acad. Sci., Prague
C
C This subroutine determines the parameters necessary to
C compute an interpolatory function on a one dimensional
C grid. The function determined can be
C represented by splines under tension.  For actual
C mapping of points it is necessary to call the subroutine
C CURVBD, which also returns first and second derivatives.
C
C On input--
C
C   NX is the number of grid points.
C   (NX should be at least 2)
C
C   X is array of the NX coordinates of the grid points.
C   These should be strictly increasing.
C
C   W is an array of the NX functional values at the
C   the grid points, i. e. W(I,J) contains the functional
C   value at X(I) for I = 1,...,NX .
C
C   C is an array of at least NX locations. This
C   parameter may coincide with W in which case W is
C   destroyed on output.
C
C   VX is the array of at least 5 * NX  locations.
C
C   TEMP is an array of at least 3 * NX locations
C   which is used for scratch storage.
C
C   SIGMA contains the tension factor. This value indicate
C   the curviness desired. If ABS(SIGMA) is nearly zero
C   (e. g. .001) the resulting surface is approximately the
C   tensor product of cubic splines. If ABS(SIGMA) is large
C   (e. g. 50.) the resulting surface is approximately
C   bi-linear. If SIGMA equals zero tensor products of cubic
C   splines result. A standard value for SIGMA is
C   approximately 1. In absolute value.
C
C On output--
C
C   C contains the coefficients of a representation of the
C   interpolated function in a B-spline form.
C
C   VX contains B-spline under tension basis data.
C
C   IERR contains an error flag.
C        = 0 for normal return,
C        = 1 if NX is less than 2,
C        = 2 if the X-array is not strictly
C            increasing.
C
C And
C
C   None of the input parameters are altered (except W if
C   this parameter and C are identical in the calling
C   sequence).
C
C This subroutine references package modules VGEN, TERMS,
C SNHCSH, TRIDEC, and TRISOL.
C
C-----------------------------------------------------------
C
C Copy W into C
C
      DO 1 I = 1,NX
    1   C(I) = W(I)
C
C Generate basis functions along X-grid
C set up tridiagonal system and solve
C
      CALL VGEN (NX,X,SIGMA,VX,IERR)
      IF (IERR .NE. 0) RETURN
      DO 2 I = 2,NX
    2   TEMP(I) = VX(5,I-1)
      NXPI = NX
      DO 3 I = 1,NX
        NXPI = NXPI+1
    3   TEMP(NXPI) = 1.
      DO 4 I = 2,NX
        NXPI = NXPI+1
    4   TEMP(NXPI) = VX(4,I)
      CALL TRIDEC (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),
     *             TEMP(1),TEMP(NX+1),IERR)
      CALL TRISOL (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),C,NX,
     *             1,1)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SURFB1 (NX,NY,X,Y,W,NW1,C,VX,VY,TEMP,SIGMA,
     *                   IERR)
C
      INTEGER NX,NY,NW1,IERR
      REAL X(NX),Y(NY),W(NW1,NY),C(NX,NY),VX(5,NX),VY(5,NY),
     *     TEMP(*),SIGMA
C
C                            From FITPACK -- August 31, 1981
C                                 Coded by Alan Kaylor Cline
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This subroutine determines the parameters necessary to
C compute an interpolatory function on a two dimensional
C rectangular grid. The function determined can be
C represented as a tensor product of splines under tension
C for actual mapping of points it is necessary to call the
C subroutine SURFBD, which also returns first and second
C partial derivatives.
C
C On input--
C
C   NX and NY are the number of grid lines in the X- and Y
C   directions, respectively, of the rectangular grid. (NX
C   and NY should be at least 2.)
C
C   X and Y are arrays of the NX and NY coordinates of the
C   grid lines in X- and Y-directions, respectively. These
C   should be strictly increasing.
C
C   W is an array of the NX * NY functional values at the
C   the grid points, i. e. W(I,J) contains the functional
C   value at (X(I),Y(J)) for I = 1,...,NX, and J = 1,...,NY.
C
C   NW1 is the first dimension of the array W used in the
C   calling program (NW1 .GE. NX).
C
C   C is an array of at least NX * NY locations. This
C   parameter may coincide with W in which case W is
C   destroyed on output.
C
C   VX and VY are arrays of at least 5 * NX and 5 * NY
C   locations, respectively.
C
C   Temp is an array of at least 3 * MAX(NX, NY) locations
C   which is used for scratch storage.
C
C And
C
C   SIGMA contains the tension factor. This value indicate
C   the curviness desired. If ABS(SIGMA) is nearly zero
C   (e. G. .001) the resulting surface is approximately the
C   tensor product of cubic splines. If ABS(SIGMA) is large
C   (e. G. 50.) the resulting surface is approximately
C   bi-linear. If SIGMA equals zero tensor products of cubic
C   splines result. A standard value for SIGMA is
C   approximately 1. In absolute value.
C
C On output--
C
C   C contains the coefficients of a representation of the
C   interpolated function in a B-spline tensor production
C   form.
C
C   VX and VY contain B-spline under tension basis data.
C
C   IERR contains an error flag.
C        = 0 for normal return,
C        = 1 if NX or NY is less than 2,
C        = 2 if the X- or Y-arrays are not strictly
C            increasing.
C
C And
C
C   None of the input parameters are altered (except W if
C   this parameter and C are identical in the calling
C   sequence).
C
C This subroutine references package modules VGEN, TERMS,
C SNHCSH, TRIDEC, and TRISOL.
C
C--------------------------------------------------------- -
C
C Copy W into C
C
      DO 1 J = 1,NY
        DO 1 I = 1,NX
    1     C(I,J) = W(I,J)
C
C Generate basis functions along X-grid
C set up tridiagonal system and solve
C
      CALL VGEN (NX,X,SIGMA,VX,IERR)
      IF (IERR .NE. 0) RETURN
      DO 2 I = 2,NX
    2   TEMP(I) = VX(5,I-1)
      NXPI = NX
      DO 3 I = 1,NX
        NXPI = NXPI+1
    3   TEMP(NXPI) = 1.
      DO 4 I = 2,NX
        NXPI = NXPI+1
    4   TEMP(NXPI) = VX(4,I)
      CALL TRIDEC (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),
     *             TEMP(1),TEMP(NX+1),IERR)
      CALL TRISOL (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),C,NX,
     *            NY,1)
C
C Generate basis functions along Y-grid
C set up tridiagonal system and solve
C
      CALL VGEN (NY,Y,SIGMA,VY,IERR)
      IF (IERR .NE. 0) RETURN
      DO 5 J = 2,NY
    5   TEMP(J) = VY(5,J-1)
      NYPJ = NY
      DO 6 J = 1,NY
        NYPJ = NYPJ+1
    6   TEMP(NYPJ) = 1.
      DO 7 J = 2,NY
        NYPJ = NYPJ+1
    7   TEMP(NYPJ) = VY(4,J)
      CALL TRIDEC (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1),
     *             TEMP(1),TEMP(NY+1),IERR)
      CALL TRISOL (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1),C,1,
     *            NX,NX)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE VAL3B1 (NX,NY,NZ,X,Y,Z,W,NW1,NW2,C,VX,VY,
     *                   VZ,TEMP,SIGMA,IERR)
C
      INTEGER NX,NY,NZ,NW1,NW2,IERR
      REAL X(NX),Y(NY),Z(NZ),W(NW1,NW2,NZ),C(NX,NY,NZ),
     *     VX(5,NX),VY(5,NY),VZ(5,NZ),TEMP(*),SIGMA
C
C                            From FITPACK -- August 31, 1981
C                                 Coded by Alan Kaylor Cline
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This subroutine determines the parameters necessary to
C compute an interpolatory function on a three dimensional
C rectangular grid. The function determined can be
C represented as a tensor product of splines under tension.
C For actual mapping of points it is necessary to call the
C subroutine VAL3BD, which also returns first and second
C partial derivatives.
C
C On input--
C
C   NX, NY, and NZ are the number of grid lines in the X-,
C   Y-, and Z-directions, respectively, of the rectangular
C   grid. (NX, NY, and NZ should be at least 2.)
C
C   X, Y, and Z are arrays of the NX, NY, and NZ coordinates
C   of the grid lines in the X-, Y-, and Z-directions,
C   respectively. These should be strictly increasing.
C
C   W is an array of the NX * NY * NZ functional values at
C   the grid points, i. e. W(I,J,K) contains the functional
C   value at (X(I),Y(J),Z(K)) for I = 1,...,NX,
C   J = 1,...,NY, and K = 1,...,NZ.
C
C   NW1 and NW2 are the first two dimensions of the array W
C   used in the calling program (NW1 .GE. NX AND NW2 .GE.
C   NY).
C
C   C is an array of at least NX * NY * NZ locations. This
C   parameter may coincide with W in which case W is
C   destroyed on output.
C
C   VX, VY, and VZ are arrays of at least 5 * NX, 5 * NY,
C   and 5 * NZ locations, respectively.
C
C   Temp is an array of at least 3 * MAX(NX, NY, NZ)
C   locations which is used for scratch storage.
C
C And
C
C   SIGMA contains the tension factor. This value indicates
C   the curviness desired. If ABS(SIGMA) is nearly zero
C   (e. g. .001) the resulting surface is approximately the
C   tensor product of cubic splines. If ABS(SIGMA) is large
C   (e. g. 50.) the resulting surface is approximately
C   tri-linear. If SIGMA equals zero tensor products of
C   cubic splines result. A standard value for SIGMA is
C   approximately 1. In absolute value.
C
C On output--
C
C   C contains the coefficients of a representation of the
C   interpolated function in a B-spline tensor production
C   form.
C
C   VX, VY, and VZ contain B-spline under tension basis
C   data.
C
C   IERR contains an error flag.
C        = 0 for normal return,
C        = 1 if NX, NY, or NZ is less than 2,
C        = 2 if the X-, Y-, or Z-arrays are not strictly
C            increasing.
C
C And
C
C   None of the input parameters are altered (except W if
C   this parameter and C are identical in the calling
C   sequence).
C
C This subroutine references package modules VGEN, TERMS,
C SNHCSH, TRIDEC, and TRISOL.
C
C-----------------------------------------------------------
C
C Copy W into C
C
      DO 1 K = 1,NZ
        DO 1 J = 1,NY
          DO 1 I = 1,NX
    1       C(I,J,K) = W(I,J,K)
C
C Generate basis functions along X-grid
C set up tridiagonal system and solve
C
      CALL VGEN (NX,X,SIGMA,VX,IERR)
      IF (IERR .NE. 0) RETURN
      DO 2 I = 2,NX
    2   TEMP(I) = VX(5,I-1)
      NXPI = NX
      DO 3 I = 1,NX
        NXPI = NXPI+1
    3   TEMP(NXPI) = 1.
      DO 4 I = 2,NX
        NXPI = NXPI+1
    4   TEMP(NXPI) = VX(4,I)
      CALL TRIDEC (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),
     *             TEMP(1),TEMP(NX+1),IERR)
      CALL TRISOL (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),C,NX,
     *            NY*NZ,1)
C
C Generate basis functions along Y-grid
C set up tridiagonal system and solve
C
      CALL VGEN (NY,Y,SIGMA,VY,IERR)
      IF (IERR .NE. 0) RETURN
      DO 5 J = 2,NY
    5   TEMP(J) = VY(5,J-1)
      NYPJ = NY
      DO 6 J = 1,NY
        NYPJ = NYPJ+1
    6   TEMP(NYPJ) = 1.
      DO 7 J = 2,NY
        NYPJ = NYPJ+1
    7   TEMP(NYPJ) = VY(4,J)
      CALL TRIDEC (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1),
     *             TEMP(1),TEMP(NY+1),IERR)
      DO 8 K = 1,NZ
    8   CALL TRISOL (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1),C(1,1,K),
     *               1,NX,NX)
C
C Generate basis functions along Z-grid
C set up tridiagonal system and solve
C
      CALL VGEN (NZ,Z,SIGMA,VZ,IERR)
      IF (IERR .NE. 0) RETURN
      DO 9 K = 2,NZ
    9   TEMP(K) = VZ(5,K-1)
      NZPK = NZ
      DO 10 K = 1,NZ
        NZPK = NZPK+1
   10   TEMP(NZPK) = 1.
      DO 11 K = 2,NZ
        NZPK = NZPK+1
   11   TEMP(NZPK) = VZ(4,K)
      CALL TRIDEC (NZ,TEMP(1),TEMP(NZ+1),TEMP(2*NZ+1),
     *             TEMP(1),TEMP(NZ+1),IERR)
      CALL TRISOL (NZ,TEMP(1),TEMP(NZ+1),TEMP(2*NZ+1),C,1,
     *            NX*NY,NX*NY)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE VGEN (N,X,SIGMA,V,IERR)
C
      INTEGER N,IERR
      REAL X(N),SIGMA,V(5,N)
C
C                            From FITPACK -- August 31, 1981
C                                 Coded by Alan Kaylor Cline
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This subroutine generates an array of coefficients used by
C other subroutines for the determination of a B-spline
C under tension basis.
C
C On input--
C
C   N is the number of knots defining the basis (N .GE. 2).
C
C   X is the array of the N increasing knots. Any linear
C   combination of the resulting basis will have third
C   derivative discontinuities only at the interior knots,
C   (i. e. X(2),...,X(N-1) ).
C
C   SIGMA contains the tension factor. This value indicates
C   the curviness desired. If ABS(SIGMA) is nearly zero
C   (e. g. .001) the basis functions are approximately cubic
C   splines. If ABS(SIGMA) is large (e. g. 50.) the basis
C   functions are nearly piecewise linear. If SIGMA equals
C   zero a cubic spline basis results. A standard value for
C   SIGMA is approximately 1. In absolute value.
C
C And
C
C   V is an array of at least 5*N locations.
C
C On output--
C
C   V contains certain coefficients to be used by other
C   subprograms for the determination of the B-spline under
C   tension basis. Considered as a 5 by N array, for I = 1,
C   ... , N, B-spline basis function I is specified by--
C     V(1,I) = second derivative at X(I-1), for I .NE. 1,
C     V(2,I) = second derivative at X(I),   for all I,
C     V(3,I) = second derivative at X(I+1), for I .NE. N,
C     V(4,I) = function value at X(I-1),    for I .NE. 1,
C     V(5,I) = function value at X(I+1),    for I .NE. N,
C   and the properties that it has--
C     1. Function value 1 at X(I),
C     2. Function value and second derivative = 0 at
C        X(1), ... , X(I-2), and X(I+2), ... , X(N).
C   In V(5,N) and V(3,N) are contained function value and
C   second derivative of basis function zero at X(1),
C   respectively. In V(4,1) and V(1,1) are contained
C   function value and second derivative of basis function
C   N+1 at X(N), respectively. Function value and second
C   derivative of these two basis functions are zero at all
C   other knots. Only basis function zero has non-zero
C   second derivative value at X(1) and only basis
C   function N+1 has non-zero second derivative at X(N).
C
C   IERR contains an error flag,
C        = 0 for normal return,
C        = 1 if N is less than 2,
C        = 2 if X-values are not strictly increasing.
C
C And
C
C   N, X, and SIGMA are unaltered.
C
C This subroutine references package modules TERMS and
C SNHCSH.
C
C-----------------------------------------------------------
C
      NM1 = N-1
      IERR = 0
      IF (N .LE. 1) GO TO 3
      IF (X(N) .LE. X(1)) GO TO 4
C
C Denormalize tension factor
C
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
C
C Generate coefficients for left end basis functions
C
      D3 = X(2)-X(1)
      IF (D3 .LE. 0.) GO TO 4
      CALL TERMS (DIAG3,SDIAG3,SIGMAP,D3)
      D4 = D3
      IF (N .GE. 3) D4 = X(3)-X(2)
      IF (D4 .LE. 0.) GO TO 4
      CALL TERMS (DIAG4,SDIAG4,SIGMAP,D4)
      A22 = DIAG3+SDIAG3
      A23 = DIAG3+DIAG4+SDIAG3+SDIAG4
      V(2,1) = 0.
      V(3,1) = 1./(D3*(DIAG3+DIAG4)+(D3+D4)*SDIAG4)
      V(5,1) = SDIAG4*D4*V(3,1)
      IF (N .EQ. 2) GO TO 2
      A22 = 2.*A22
      D1 = D3
      D2 = D3
      D3 = D4
      DIAG1 = DIAG3
      DIAG2 = DIAG3
      DIAG3 = DIAG4
      SDIAG1 = SDIAG3
      SDIAG2 = SDIAG3
      SDIAG3 = SDIAG4
C
C Generate coefficients for interior basis functions
C
      DO 1 I = 2,NM1
        IF (I .NE. NM1) D4 = X(I+2)-X(I+1)
        IF (D4 .LE. 0.) GO TO 4
        IF (D4 .NE. D3) CALL TERMS (DIAG4,SDIAG4,SIGMAP,D4)
        A11 = DIAG1+DIAG2+SDIAG1*(1.+D1/D2)
        A12 = SDIAG2/A11
        B1 = 1./(D2*A11)
        A33 = DIAG3+DIAG4+SDIAG4*(1.+D4/D3)
        A32 = SDIAG3/A33
        B3 = 1./(D3*A33)
        A21 = A22
        A22 = A23
        A23 = DIAG3+DIAG4+SDIAG3+SDIAG4
        V(2,I) = -(A21*B1+A23*B3)/(A22-A21*A12-A23*A32)
        V(1,I) = B1-A12*V(2,I)
        V(3,I) = B3-A32*V(2,I)
        V(4,I) = SDIAG1*D1*V(1,I)
        V(5,I) = SDIAG4*D4*V(3,I)
C
C Save constants for next iteration
C
        D1 = D2
        D2 = D3
        D3 = D4
        DIAG1 = DIAG2
        DIAG2 = DIAG3
        DIAG3 = DIAG4
        SDIAG1 = SDIAG2
        SDIAG2 = SDIAG3
    1   SDIAG3 = SDIAG4
C
C Generate coefficients for right end basis functions
C
      V(2,N) = 0.
      V(1,N) = 1./(D2*(DIAG1+DIAG2)+(D2+D1)*SDIAG1)
      V(4,N) = SDIAG1*D1*V(1,N)
      V(3,N) = V(1,3)
      V(5,N) = V(4,3)
      V(1,1) = V(3,N-2)
      V(4,1) = V(5,N-2)
C
C Adjust basis for natural end conditions
C
      V(4,2) = V(4,2)-V(1,2)*V(5,N)/V(3,N)
      V(1,2) = 0.
      V(5,NM1) = V(5,NM1)-V(3,NM1)*V(4,1)/V(1,1)
      V(3,NM1) = 0.
      RETURN
C
C N equal to 2
C
    2 V(4,1) = V(5,1)
      V(1,1) = V(3,1)
      V(3,1) = 0.
      V(5,1) = 0.
      V(1,2) = 0.
      V(2,2) = 0.
      V(3,2) = V(1,1)
      V(4,2) = 0.
      V(5,2) = V(4,1)
      RETURN
C
C Too few knots
C
    3 IERR = 1
      RETURN
C
C X-values not strictly increasing
C
    4 IERR = 2
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE TERMS (DIAG,SDIAG,SIGMA,DEL)
C
      REAL DIAG,SDIAG,SIGMA,DEL
C
C                            From FITPACK -- August 31, 1981
C                       Coded by A. K. Cline and R. J. Renka
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This subroutine computes the diagonal and superdiagonal
C terms of the tridiagonal linear system associated with
C spline under tension interpolation.
C
C On input--
C
C   SIGMA contains the tension factor.
C
C And
C
C   DEL contains the step size.
C
C On output--
C
C               (SIGMA*DEL*COSH(SIGMA*DEL) - SINH(SIGMA*DEL)
C   DIAG = DEL*--------------------------------------------.
C                     (SIGMA*DEL)**2 * SINH(SIGMA*DEL)
C
C                   SINH(SIGMA*DEL) - SIGMA*DEL
C   SDIAG = DEL*----------------------------------.
C                (SIGMA*DEL)**2 * SINH(SIGMA*DEL)
C
C And
C
C   SIGMA and DEL are unaltered.
C
C This subroutine references package module SNHCSH.
C
C-----------------------------------------------------------
C
      IF (SIGMA .NE. 0.) GO TO 1
      DIAG = DEL/3.
      SDIAG = DEL/6.
      RETURN
    1 SIGDEL = SIGMA*DEL
      CALL SNHCSH (SINHM,COSHM,SIGDEL,0)
      DENOM = DEL/((SINHM+SIGDEL)*SIGDEL*SIGDEL)
      DIAG = DENOM*(SIGDEL*COSHM-SINHM)
      SDIAG = DENOM*SINHM
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE TRIDEC (N,SUBDI,DIAGI,SUPD,SUBDO,DIAGO,
     *                   IERR)
C
      INTEGER N,IERR
      REAL SUBDI(N),DIAGI(N),SUPD(N),SUBDO(N),DIAGO(N)
C
C                            From FITPACK -- August 31, 1981
C                                 Coded by Alan Kaylor Cline
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This subroutine factorizes a tridiagonal matrix in order
C to solve systems of linear equations. The factorization
C employs gaussian elimination without any interchanging of
C columns or rows. The subroutine trisol may be called to
C actually solve the system once the factorization has been
C performed.
C
C On input--
C
C   N contains the order of the matrix (N .GE. 1).
C
C   SUBDI is an array containing the subdiagonal elements of
C   the matrix in positions 2, ... , N.
C
C   DIAGI is an array containing the diagonal elements of
C   the matrix.
C
C   SUPD is an array containing the superdiagonal elements
C   of the matrix in positions 1, ... , N-1.
C
C And
C
C   SUBDO and DIAGO are arrays of length N. (The storage
C   for these may coincide with that for SUBDI and DIAGI,
C   respectively, in which case the original contents of
C   SUBDI and DIAGI will be destroyed.)
C
C On output--
C
C   SUBDO and DIAGO contain the subdiagonal and diagonal of
C   the factorization matrix.
C
C   IERR contains an error flag,
C        = 0 for normal return,
C        = 1 if N is less than 1,
C        = 2 if the system is singular.
C
C And
C
C   N, SUBDI, DIAGI, and SUPD are unaltered (unless storage
C   for SUBDI or DIAGI coincided with that for SUBDO
C   or DIAGO, respectively).
C
C-----------------------------------------------------------
C
      IF (N .LE. 0) GO TO 3
      IERR = 2
      DIAGO(1) = DIAGI(1)
      IF (N .EQ. 1) GO TO 2
C
C Forward elimination
C
      DO 1 I = 2,N
        IM1 = I-1
        IF (DIAGO(IM1) .EQ. 0.) RETURN
        DIAGO(IM1) = 1./DIAGO(IM1)
        SUBDO(I) = SUBDI(I)*DIAGO(IM1)
    1   DIAGO(I) = DIAGI(I)-SUBDO(I)*SUPD(IM1)
    2 IF (DIAGO(N) .EQ. 0.) RETURN
      DIAGO(N) = 1./DIAGO(N)
      IERR = 0
      RETURN
C
C N less than 1
C
    3 IERR = 1
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE TRISOL (N,SUBD,DIAG,SUPD,RHS,MRHS,NUMRHS,
     *                   INCRHS)
C
      INTEGER N,MRHS,NUMRHS,INCRHS
      REAL SUBD(N),DIAG(N),SUPD(N)
      REAL RHS(1+INCRHS*(N-1)+MRHS*(NUMRHS-1))
C
C                            From FITPACK -- August 31, 1981
C                                 Coded by Alan Kaylor Cline
C                            Department of Computer Sciences
C                              University of Texas at Austin
C                               Revised -- December 31, 1992
C                                            by Ludek Klimes
C                                   Institute of Geotechnics
C                               Czechosl. Acad. Sci., Prague
C
C This subroutine solves tridiagonal systems of linear
C equations with multiple right hand sides. The right hand
C sides may be stored row-wise or column-wise. The
C subroutine TRIDEC should be called earlier to determine a
C factorization of the tridiagonal matrix. The solution
C vectors over-write the right hand sides.
C
C On input--
C
C   N contains the order of the matrix (N .GE. 1).
C
C   SUBD, DIAG, and SUPD are arrays of length N containing
C   the subdiagonal, diagonal, and superdiagonal of the
C   factorization, respectively.
C
C   RHS is an array containing the right hand sides of the
C   tridiagonal system.
C
C   MRHS is the increment between the first components of
C   each of the right hand side vectors in storage (MRHS
C   .GE. 1).
C
C   NUMRHS is the number of right hand sides to be solved.
C
C And
C
C   INCRHS is the increment between components within each
C   of the right hand side vectors in storage (INCRHS .GE.
C   1).
C
C The parameters N, SUBD, DIAG, and SUPD may be input as the
C parameters N, SUBDO, DIAGO, and SUPD output by subroutine
C TRIDEC, respectively.
C
C On output--
C
C   RHS contains the solution vectors in the same storage
C   structure as for the right hand sides.
C
C And
C
C   N, SUBD, DIAG, SUPD, MRHS, NUMRHS, and INCRHS are
C   unaltered.
C
C-----------------------------------------------------------
C
      NP1 = N+1
C
C Loop on right hand sides
C
      DO 4 K = 1,NUMRHS
C
C Forward elimination
C
        IRHS = 1+MRHS*(K-1)
        IF (N .EQ. 1) GO TO 2
        DO 1 I = 2,N
          IM1RHS = IRHS
          IRHS = IRHS+INCRHS
    1     RHS(IRHS) = RHS(IRHS)-SUBD(I)*RHS(IM1RHS)
C
C Back substitution
C
    2   RHS(IRHS) = DIAG(N)*RHS(IRHS)
        IF (N .EQ. 1) GO TO 4
        DO 3 IBAK = 2,N
          I = NP1-IBAK
          RHS(IM1RHS) = DIAG(I)*(RHS(IM1RHS)-SUPD(I)
     *                                           *RHS(IRHS))
          IRHS = IM1RHS
    3     IM1RHS = IM1RHS-INCRHS
    4   CONTINUE
      RETURN
      END
C
C=======================================================================
C
C Part 2:
C
C=======================================================================
C
C     
C
      SUBROUTINE CURV2D (T,YY,YX,YXX,N,X,Y,YP,SIGMA)
C
      INTEGER N
      REAL T,YY,YX,YXX,X(N),Y(N),YP(N),SIGMA
C
C                            From FITPACK -- August 31, 1981
C                                 Coded by Alan Kaylor Cline
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This subroutine determines function value, first, and
C second derivatives of a curve at a given point using a
C spline under tension. The subroutine CURV1 should be
C called earlier to determine certain necessary parameters.
C
C On input--
C
C   T contains a real value at which the function and
C   derivatives are to be evaluated.
C
C   N contains the number of points which were specified to
C   determine the curve.
C
C   X and Y are arrays containing the abscissae and
C   ordinates, respectively, of the specified points.
C
C   YP is an array of second derivative values of the curve
C   at the nodes.
C
C And
C
C   SIGMA contains the tension factor (its sign is ignored).
C
C The parameters N, X, Y, YP, and SIGMA should be input
C unaltered from the output of CURV1.
C
C On output--
C
C   YY, YX, and YXX contain the function value, first and
C   second derivatives, respectively.
C
C None of the input parameters are altered.
C
C This subroutine references package modules INTRVL and
C SNHCSH.
C
C-----------------------------------------------------------
C
C Determine interval
C
      IM1 = INTRVL(T,X,N)
      I = IM1+1
C
C Denormalize tension factor
C
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
C
C Set up and perform interpolation
C
      DEL1 = T-X(IM1)
      DEL2 = X(I)-T
      DELS = X(I)-X(IM1)
      YY = (Y(I)*DEL1+Y(IM1)*DEL2)/DELS
      YX = (Y(I)-Y(IM1))/DELS
      IF (SIGMAP .NE. 0.) GO TO 1
      YY = YY-DEL1*DEL2*(YP(I)*(DEL1+DELS)+YP(IM1)*
     *        (DEL2+DELS))/(6.*DELS)
      YX = YX+(YP(I)*(2.*DEL1*DEL1-DEL2*(DEL1+DELS))-
     *         YP(IM1)*(2.*DEL2*DEL2-DEL1*(DEL2+DELS)))
     *         /(6.*DELS)
      YXX = (YP(I)*DEL1+YP(IM1)*DEL2)/DELS
      RETURN
    1 DELP1 = SIGMAP*(DEL1+DELS)/2.
      DELP2 = SIGMAP*(DEL2+DELS)/2.
      CALL SNHCSH (SINHM1,COSHM1,SIGMAP*DEL1,0)
      CALL SNHCSH (SINHM2,COSHM2,SIGMAP*DEL2,0)
      CALL SNHCSH (SINHMS,DUMMY,SIGMAP*DELS,-1)
      CALL SNHCSH (SINHP1,DUMMY,SIGMAP*DEL1/2.,-1)
      CALL SNHCSH (SINHP2,DUMMY,SIGMAP*DEL2/2.,-1)
      CALL SNHCSH (DUMMY,COSHP1,DELP1,1)
      CALL SNHCSH (DUMMY,COSHP2,DELP2,1)
      YY = YY+(YP(I)*(SINHM1*DEL2-DEL1*(2.*(COSHP1+1.)*
     *        SINHP2+SIGMAP*COSHP1*DEL2))+YP(IM1)*(SINHM2*
     *        DEL1-DEL2*(2.*(COSHP2+1.)*SINHP1+SIGMAP*
     *        COSHP2*DEL1)))/(SIGMAP*SIGMAP*DELS*(SINHMS+
     *        SIGMAP*DELS))
      YX = YX+(YP(I)*(DELS*SIGMAP*COSHM1-SINHMS)-
     *     YP(IM1)*(DELS*SIGMAP*COSHM2-SINHMS))/
     *     (SIGMAP*SIGMAP*DELS*(SINHMS+SIGMAP*DELS))
      YXX = (YP(I)*(SINHM1+SIGMAP*DEL1)+YP(IM1)*(SINHM2+
     *       SIGMAP*DEL2))/(SINHMS+SIGMAP*DELS)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE CURVBD (XX,W,WX,WXX,NX,X,C,VX,SIGMA)
C
      INTEGER NX
      REAL XX,W,WX,WXX,X(NX),VX(5,NX),C(NX),SIGMA
C
C                                      Complement to FITPACK
C                                       by Alan Kaylor Cline
C                                   Coded -- October 9, 1986
C                                            by Ludek Klimes
C                                      Inst. Geol. Geotechn.
C                               Czechosl. Acad. Sci., Prague
C
C This subroutine evaluates the function value, the
C first partial derivative, and the second partial
C derivative of a spline under tension in one variable.
C
C On input--
C
C   XX contains the X-coordinate of the point
C   at which the interpolation is to be performed
C
C   NX is the number of grid points
C
C   X is array containing the X-grid values.
C
C   C is an array of coefficients describing the function in
C   terms of a B-spline under tension basis. In the
C   expansion of the function, for I = 1,...,NX ,
C   the coefficient multiplying the basis
C   function I is stored in C(I).
C
C   VX is the array of length 5*NX
C   containing the B-spline basis data
C
C   SIGMA contains the tension factor (its sign is ignored).
C
C The parameters NX, X, C, VX, and SIGMA
C should be input unaltered from the output of CURVB1.
C
C On output--
C
C   W contains the interpolated function value.
C
C   WX contains the first derivative .
C
C   WXX contains the second derivative .
C
C And
C
C   None of the input parameters are altered.
C
C This subroutine references package modules DSPLNZ, INTRVL,
C and SNHCSH.
C
C--------------------------------------------------------------
C
      REAL BX(3,4)
C
C Evaluate basis functions at XX
C
      CALL DSPLNZ (XX,NX,X,VX,SIGMA,ISTART,BX)
C
C Accumulate basis functions
C
      SUM = 0.
      SUMX = 0.
      SUMXX = 0.
      DO 1 I = 1,4
        II = ISTART+I-1
        IF (II .EQ. 0 .OR. II .GT. NX) GO TO 1
        BX1I = BX(1,I)
        CI  = C(II)
        SUM = SUM+CI*BX1I
        SUMX = SUMX+CI*BX(2,I)
        SUMXX = SUMXX+CI*BX(3,I)
*V      CALL VAR2(II,BX1I,BX(2,I),0.,0.)
    1   CONTINUE
      W = SUM
      WX = SUMX
      WXX = SUMXX
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SURFBD (XX,YY,W,WX,WY,WXX,WXY,WYY,NX,NY,X,
     *                   Y,C,VX,VY,SIGMA)
C
      INTEGER NX,NY
      REAL XX,YY,W,WX,WY,WXX,WXY,WYY,X(NX),Y(NY),VX(5,NX),
     *     VY(5,NY),C(NX,NY),SIGMA
C
C                            From FITPACK -- August 31, 1981
C                                 Coded by Alan Kaylor Cline
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This subroutine evaluates the function value, the two
C first partial derivatives, and the six second partial
C derivatives of a tensor product spline under tension in
C two variables.
C
C On input--
C
C   XX and YY contain the X- and Y-coordinates of the point
C   at which the interpolation is to be performed.
C
C   NX and NY are the number of grid lines in the X- and Y-
C   directions, respectively, of the rectangular grid which
C   specified the function.
C
C   X and Y are arrays containing the X- and Y-grid values,
C   respectively.
C
C   C is an array of coefficients describing the function in
C   terms of a B-spline under tension basis. In the
C   expansion of the function, for I = 1,...,NX and J = 1,
C   ...,NY, the coefficient multiplying the product of basis
C   function I in X and basis function J in Y is stored in
C   C(I,J).
C
C   VX and VY VZ are arrays of length 5*NX and 5*NY,
C   respectively, containing the B-spline basis data for the
C   X- and Y-grids.
C
C And
C
C   SIGMA contains the tension factor (its sign is ignored).
C
C The parameters NX, NY, X, Y, Z, C, VX, VY, and SIGMA
C should be input unaltered from the output of SURFB1.
C
C On output--
C
C   W contains the interpolated function value.
C
C   WX and WY contain the X- and Y-partial derivatives,
C   respectively.
C
C   WXX, WXY, and WYY contain the XX-, XY-, and YY-partial
C   derivatives, respectively.
C
C And
C
C   None of the input parameters are altered.
C
C This subroutine references package modules DSPLNZ, INTRVL,
C and SNHCSH.
C
C--------------------------------------------------------- ----
C
      REAL BX(3,4),BY(3,4)
C
C Evaluate basis functions at XX and YY
C
      CALL DSPLNZ (XX,NX,X,VX,SIGMA,ISTART,BX)
      CALL DSPLNZ (YY,NY,Y,VY,SIGMA,JSTART,BY)
C
C Accumulate tensor products
C
      SUM = 0.
      SUMX = 0.
      SUMY = 0.
      SUMXX = 0.
      SUMXY = 0.
      SUMYY = 0.
      DO 2 J = 1,4
        JJ = JSTART+J-1
        IF (JJ .EQ. 0 .OR. JJ .GT. NY) GO TO 2
        BY1J = BY(1,J)
        BY2J = BY(2,J)
        BY3J = BY(3,J)
        DO 1 I = 1,4
          II = ISTART+I-1
          IF (II .EQ. 0 .OR. II .GT. NX) GO TO 1
          BX1I = BX(1,I)
          BX2I =  BX(2,I)
          CIJ = C(II,JJ)
          SUM = SUM+CIJ*BX1I*BY1J
          SUMX = SUMX+CIJ*BX2I*BY1J
          SUMY = SUMY+CIJ*BX1I*BY2J
          SUMXX = SUMXX+CIJ*BX(3,I)*BY1J
          SUMXY = SUMXY+CIJ*BX2I*BY2J
          SUMYY = SUMYY+CIJ*BX1I*BY3J
*V        CALL VAR2(II+NX*(JJ-1),BX1I*BY1J,BX2I*BY1J,BX1I*BY2J,0.)
    1     CONTINUE
    2   CONTINUE
      W = SUM
      WX = SUMX
      WY = SUMY
      WXX = SUMXX
      WXY = SUMXY
      WYY = SUMYY
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE VAL3BD (XX,YY,ZZ,W,WX,WY,WZ,WXX,WXY,WYY,
     *                   WYZ,WZZ,WXZ,NX,NY,NZ,X,Y,Z,C,VX,VY,
     *                   VZ,SIGMA)
C
      INTEGER NX,NY,NZ
      REAL XX,YY,ZZ,W,WX,WY,WZ,WXX,WXY,WYY,WYZ,WZZ,WXZ,
     *     X(NX),Y(NY),Z(NZ),VX(5,NX),VY(5,NY),VZ(5,NZ),
     *     C(NX,NY,NZ),SIGMA
C
C                            From FITPACK -- August 31, 1981
C                                 Coded by Alan Kaylor Cline
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This subroutine evaluates the function value, the three
C first partial derivatives, and the six second partial
C derivatives of a tensor product spline under tension in
C three variables.
C
C On input--
C
C   XX, YY, and ZZ contain the X-, Y-, and Z-coordinates of
C   the point at which the interpolation is to be performed.
C
C   NX, NY, and NZ are the number of grid lines in the X-,
C   Y-, and Z-directions, respectively, of the rectangular
C   grid which specified the function.
C
C   X, Y, and Z are arrays containing the X-, Y-, and Z-grid
C   values, respectively.
C
C   C is an array of coefficients describing the function in
C   terms of a B-spline under tension basis. In the
C   expansion of the function, for I = 1,...,NX, J = 1,...,
C   NY, AND K = 1,...,NZ, the coefficient multiplying the
C   product of basis function I in X, basis function J in Y,
C   and basis function K in Z is stored in C(I,J,K).
C
C   VX, VY, and VZ are arrays of length 5*NX, 5*NY, and
C   5*NZ, respectively, containing the B-spline basis data
C   for the X-, Y-, and Z-grids.
C
C And
C
C   SIGMA contains the tension factor (its sign is ignored).
C
C The parameters NX, NY, NZ, X, Y, Z, C, VX, VY, VZ, and
C SIGMA should be input unaltered from the output of
C VAL3B1.
C
C On output--
C
C   W contains the interpolated function value.
C
C   WX, WY, and WZ contain the X-, Y-, and Z-partial
C   derivatives, respectively.
C
C   WXX, WXY, WYY, WYZ, WZZ, and WXZ contain the XX-, XY-
C   YY-, YZ-, ZZ-, and XZ-partial derivatives, respectively.
C
C And
C
C   None of the input parameters are altered.
C
C This subroutine references package modules DSPLNZ, INTRVL,
C and SNHCSH.
C
C--------------------------------------------------------------
C
      REAL BX(3,4),BY(3,4),BZ(3,4)
C
C Evaluate basis functions at XX, YY, and ZZ
C
      CALL DSPLNZ (XX,NX,X,VX,SIGMA,ISTART,BX)
      CALL DSPLNZ (YY,NY,Y,VY,SIGMA,JSTART,BY)
      CALL DSPLNZ (ZZ,NZ,Z,VZ,SIGMA,KSTART,BZ)
C
C Accumulate tensor products
C
      SUM = 0.
      SUMX = 0.
      SUMY = 0.
      SUMZ = 0.
      SUMXX = 0.
      SUMXY = 0.
      SUMYY = 0.
      SUMYZ = 0.
      SUMZZ = 0.
      SUMXZ = 0.
      DO 3 K = 1,4
        KK = KSTART+K-1
        IF (KK .EQ. 0 .OR. KK .GT. NZ) GO TO 3
        BZ1K = BZ(1,K)
        BZ2K = BZ(2,K)
        BZ3K = BZ(3,K)
        DO 2 J = 1,4
          JJ = JSTART+J-1
          IF (JJ .EQ. 0 .OR. JJ .GT. NY) GO TO 2
          BY1J = BY(1,J)
          BY2J = BY(2,J)
          BY3J = BY(3,J)
          DO 1 I = 1,4
            II = ISTART+I-1
            IF (II .EQ. 0 .OR. II .GT. NX) GO TO 1
            BX1I = BX(1,I)
            BX2I =  BX(2,I)
            CIJK = C(II,JJ,KK)
            SUM = SUM+CIJK*BX1I*BY1J*BZ1K
            SUMX = SUMX+CIJK*BX2I*BY1J*BZ1K
            SUMY = SUMY+CIJK*BX1I*BY2J*BZ1K
            SUMZ = SUMZ+CIJK*BX1I*BY1J*BZ2K
            SUMXX = SUMXX+CIJK*BX(3,I)*BY1J*BZ1K
            SUMXY = SUMXY+CIJK*BX2I*BY2J*BZ1K
            SUMYY = SUMYY+CIJK*BX1I*BY3J*BZ1K
            SUMYZ = SUMYZ+CIJK*BX1I*BY2J*BZ2K
            SUMZZ = SUMZZ+CIJK*BX1I*BY1J*BZ3K
            SUMXZ = SUMXZ+CIJK*BX2I*BY1J*BZ2K
*V          CALL VAR2(II+NX*(JJ-1+NY*(KK-1)),BX1I*BY1J*BZ1K,
*V   *                     BX2I*BY1J*BZ1K,BX1I*BY2J*BZ1K,BX1I*BY1J*BZ2K)
    1       CONTINUE
    2     CONTINUE
    3   CONTINUE
      W = SUM
      WX = SUMX
      WY = SUMY
      WZ = SUMZ
      WXX = SUMXX
      WXY = SUMXY
      WYY = SUMYY
      WYZ = SUMYZ
      WZZ = SUMZZ
      WXZ = SUMXZ
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE DSPLNZ (T,N,X,V,SIGMA,ISTART,B)
C
      INTEGER N,ISTART
      REAL T,X(N),V(5,N),SIGMA,B(3,4)
C
C                            From FITPACK -- August 31, 1981
C                                 Coded by Alan Kaylor Cline
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This subroutine evaluates at a given point the four non-
C zero basis functions of a B-spline under tension basis and
C their first and second derivatives. The index of the first
C non-zero basis function is also determined. (the sense of
C the word non-zero is extended to include the special case
C where the given point coincides with a knot in which case
C the last of the four returned function values may be zero.
C ) the subroutine VGEN should be called earlier to
C determine certain necessary coefficients.
C
C On input--
C
C   T contains a real value at which the basis functions are
C   to be evaluated.
C
C   N contains the number of knots defining the basis.
C
C   X contains the array of knots.
C
C   V contains the array of coefficients determined by VGEN
C   for calculation of basis functions.
C
C   SIGMA contains the tension factor (its sign is ignored).
C
C   ISTART is an integer variable.
C
C And
C
C   B is a real array with 3 rows and 4 columns.
C
C The parameters N, X, V, and SIGMA should be input
C unaltered from the output of VGEN.
C
C On output--
C
C   ISTART contains the index of the first non-zero basis
C   function. Thus 0 .LE. ISTART .LE. N-2 and the non-zero
C   basis functions have indices ISTART, ... , ISTART+3.
C
C   B contains the values at T of basis functions ISTART,
C   ... , ISTART+3 in B(1,1), ... , B(1,4), respectively.
C   First and second derivatives of the corresponding
C   functions are contained in B(2,1), ... , B(2,4), and
C   B(3,1), ... , B(3,4), respectively.
C
C   T, N, X, V, and SIGMA are unaltered.
C
C This subroutine references package modules INTRVL and
C SNHCSH.
C
C-----------------------------------------------------------
C
C Denormalize tension factor
C
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
C
C Determine index of first non-zero basis function
C
      I = INTRVL (T,X,N)-1
C
C Compute distances to adjacent knots and lagrangian
C weights
C
      DEL1 = T-X(I+1)
      DEL2 = X(I+2)-T
      DELS = X(I+2)-X(I+1)
      C10 = DEL2/DELS
      C20 = DEL1/DELS
      C11 = -1./DELS
      C21 = 1./DELS
      IF (SIGMAP .NE. 0.) GO TO 1
      FAC = -DEL1*DEL2/(6.*DELS)
      CP10 = FAC*(DEL2+DELS)
      CP20 = FAC*(DEL1+DELS)
      CP11 = -(2.*DEL2*DEL2-DEL1*(DEL2+DELS))/(6.*DELS)
      CP21 = (2.*DEL1*DEL1-DEL2*(DEL1+DELS))/(6.*DELS)
      CP12 = C10
      CP22 = C20
      GO TO 2
    1 DELP1 = SIGMAP*(DEL1+DELS)/2.
      DELP2 = SIGMAP*(DEL2+DELS)/2.
      CALL SNHCSH (SINHM1,COSHM1,SIGMAP*DEL1,0)
      CALL SNHCSH (SINHM2,COSHM2,SIGMAP*DEL2,0)
      CALL SNHCSH (SINHMS,DUMMY,SIGMAP*DELS,-1)
      CALL SNHCSH (SINHP1,DUMMY,SIGMAP*DEL1/2.,-1)
      CALL SNHCSH (SINHP2,DUMMY,SIGMAP*DEL2/2.,-1)
      CALL SNHCSH (DUMMY,COSHP1,DELP1,1)
      CALL SNHCSH (DUMMY,COSHP2,DELP2,1)
      SINHS = SINHMS+SIGMAP*DELS
      DENOM = SIGMAP*SIGMAP*DELS*SINHS
      CP10 = (SINHM2*DEL1-DEL2*(2.*(COSHP2+1.)*SINHP1
     *       +SIGMAP*COSHP2*DEL1))/DENOM
      CP20 = (SINHM1*DEL2-DEL1*(2.*(COSHP1+1.)*SINHP2
     *       +SIGMAP*COSHP1*DEL2))/DENOM
      CP11 = -(DELS*SIGMAP*COSHM2-SINHMS)/DENOM
      CP21 = (DELS*SIGMAP*COSHM1-SINHMS)/DENOM
      CP12 = (SINHM2+SIGMAP*DEL2)/SINHS
      CP22 = (SINHM1+SIGMAP*DEL1)/SINHS
C
C Compute basis function values
C
    2 II = I
      IF (II .EQ. 0) II = N
      IIP1 = I+1
      IIP2 = I+2
      IIP3 = I+3
      IF (IIP2 .EQ. N) IIP3 = 1
      B(1,1) = C10*V(5,II)+CP10*V(3,II)
      B(1,2) = C10+C20*V(5,IIP1)+CP10*V(2,IIP1)+
     *         CP20*V(3,IIP1)
      B(1,3) = C10*V(4,IIP2)+C20+CP10*V(1,IIP2)+
     *         CP20*V(2,IIP2)
      B(1,4) = C20*V(4,IIP3)+CP20*V(1,IIP3)
      B(2,1) = C11*V(5,II)+CP11*V(3,II)
      B(2,2) = C11+C21*V(5,IIP1)+CP11*V(2,IIP1)+
     *         CP21*V(3,IIP1)
      B(2,3) = C11*V(4,IIP2)+C21+CP11*V(1,IIP2)+
     *         CP21*V(2,IIP2)
      B(2,4) = C21*V(4,IIP3)+CP21*V(1,IIP3)
      B(3,1) = CP12*V(3,II)
      B(3,2) = CP12*V(2,IIP1)+CP22*V(3,IIP1)
      B(3,3) = CP12*V(1,IIP2)+CP22*V(2,IIP2)
      B(3,4) = CP22*V(1,IIP3)
      ISTART = I
      RETURN
      END
C
C=======================================================================
C
C     
C
      FUNCTION INTRVL (T,X,N)
C
      INTEGER N
      REAL T,X(N)
C
C                            From FITPACK -- August 31, 1981
C                       Coded by A. K. Cline and R. J. Renka
C                            Department of Computer Sciences
C                              University of Texas at Austin
C
C This function determines the index of the interval
C (determined by a given increasing sequence) in which
C a given value lies.
C
C On input--
C
C   T is the given value.
C
C   X is a vector of strictly increasing values.
C
C And
C
C   N is the length of X (N .GE. 2).
C
C On output--
C
C   INTRVL returns an integer I such that
C
C          I = 1       if             T .LT. X(2)  ,
C          I = N-1     if X(N-1) .LE. T            ,
C          otherwise       X(I)  .LE. T .LT. X(I+1),
C
C None of the input parameters are altered.
C
C-----------------------------------------------------------
C
      TT = T
      IF (TT .LT. X(2)) GO TO 4
      IF (TT .GE. X(N-1)) GO TO 5
      IL = 2
      IH = N-1
C
C Linear interpolation
C
    1 I = MIN0(IL+IFIX(FLOAT(IH-IL)*(TT-X(IL))/(X(IH)-X(IL))),
     *         IH-1)
      IF (TT .LT. X(I)) GO TO 2
      IF (TT .LT. X(I+1)) GO TO 3
C
C Too high
C
      IL = I+1
      GO TO 1
C
C Too low
C
    2 IH = I
      GO TO 1
    3 INTRVL = I
      RETURN
C
C Left end
C
    4 INTRVL = 1
      RETURN
C
C Right end
C
    5 INTRVL = N-1
      RETURN
      END
C
C=======================================================================
C