      SUBROUTINE SPSP(NAX,NAY,NAZ,NBX,NBY,NBZ,NX,NY,NZ,
     *                X,Y,Z,VX,VY,VZ,SIGMA,WEIGHT,B,BX,BY,BZ)
C
      INTEGER NAX,NAY,NAZ,NBX,NBY,NBZ,NX,NY,NZ
      REAL    X(NX),Y(NY),Z(NZ),VX(5,NX),VY(5,NY),VZ(5,NZ)
      REAL    SIGMA,WEIGHT,B(*),BX(*),BY(*),BZ(*)
C
C                                      COMPLEMENT TO FITPACK
C                                       BY ALAN KAYLOR CLINE
C                                  CODED -- JANUARY 23, 1994
C                                            BY LUDEK KLIMES
C                                   DEPARTMENT OF GEOPHYSICS
C                                 CHARLES UNIVERSITY, PRAGUE
C
C THIS SUBROUTINE EVALUATES THE SOBOLEV SCALAR PRODUCTS
C OF SPLINE UNDER TENSION BASIS FUNCTIONS IN THREE VARIABLES
C (THE SOBOLEV SCALAR PRODUCT CONSISTS OF INTEGRALS OF THE
C PRODUCTS OF PARTIAL DERIVATIVES OF THE TWO ARGUMENT FUNCTIONS)
C
C ON INPUT--
C
C   NXA, NYA, NZA ARE THE ORDERS OF PARTIAL DERIVATIVES OF
C   THE FIRST ARGUMENT FUNCTION IN THE SCALAR PRODUCT
C
C   NXB, NYB, NZB ARE THE ORDERS OF PARTIAL DERIVATIVES OF
C   THE SECOND ARGUMENT FUNCTION IN THE SCALAR PRODUCT
C
C   NX, NY, NZ ARE THE NUMBERS OF GRID POINTS IN THE
C   X-, Y-, Z-DIRECTIONS, RESPECTIVELY. (NX, NY, NZ
C   SHOULD BE AT LEAST 1)
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   VX, VY,VZ ARE ARRAYS OF LENGTHS 5*NX, 5*NY, 5*NZ,
C   RESPECTIVELY, CONTAINING THE B-SPLINE BASIS DATA FOR THE
C   X-, Y- AND Z-GRIDS.  THEY CONTAIN CERTAIN COEFFICIENTS
C   TO BE USED 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. FUCTION 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   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   WEIGHT IS THE WEIGHT OF THE PRODUCT OF NXA,NYA,NZA-PARTIAL
C   DERIVATIVE OF THE FIRST ARGUMENT AND NXB,NYB,NZB-PARTIAL
C   DERIVATIVE OF THE SECOND ARGUMENT, IN THE SOBOLEV SCALAR
C   PRODUCT. THE INTEGRAL OF THE PRODUCT OF THE PARTIAL
C   DERIVATIVES MULTIPLIED BY WEIGHT IS ADDED TO MATRIX B.
C
C   B IS THE ARRAY CONTAINING NN*NN MATRIX B (NN=NX*NY*NZ),
C   STORED AS A SYMMETRIC MATRIX ( NN*(NN+1)/2 STORAGE
C   LOCATIONS ) IF NAX.EQ.NBX AND NAY.EQ.NBY AND NAZ.EQ.NBZ,
C   ELSE STORED AS A GENERAL MATRIX ( NN*NN STORAGE
C   LOCATIONS ).  THE II,JJ-ELEMENT OF THE MATRIX B
C   WILL BE INCREASED BY THE INTEGRAL OF THE PRODUCT 0F
C   NXA-,NYA-,NZA-PARTIAL DERIVATIVE OF THE II-TH BASIS
C   FUNCTION AND NXB-,NYB-,NZB-PARTIAL DERIVATIVE OF THE
C   JJ-TH BASIS FUNCTION, MULTIPLIED BY WEIGHT.
C   HERE THE BASIS FUNCTION IX,IY,IZ (1.LE.IX.LE.NX,
C   1.LE.IY.LE.NY, 1.LE.IZ.LE.NZ) IS INDEXED BY
C   II=IX+NX*(IY+NY*IZ).
C
C   BX IS AN AUXILIARY ARRAY OF AT LEAST NX*(NX+1)/2
C   LOCATIONS FOR NXA.EQ.NXB, OR OF AT LEAST NX*NX LOCATIONS
C   FOR NXA.NE.NXB.  IT IS USED FOR SCRATCH STORAGE.
C
C   BY IS AN AUXILIARY ARRAY OF AT LEAST NY*(NY+1)/2
C   LOCATIONS FOR NYA.EQ.NYB, OR OF AT LEAST NY*NY LOCATIONS
C   FOR NYA.NE.NYB.  IT IS USED FOR SCRATCH STORAGE.
C
C   BZ IS AN AUXILIARY ARRAY OF AT LEAST NZ*(NZ+1)/2
C   LOCATIONS FOR NZA.EQ.NZB, OR OF AT LEAST NZ*NZ LOCATIONS
C   FOR NZA.NE.NZB.  IT IS USED FOR SCRATCH STORAGE.
C
C AND
C
C   NONE OF THE INPUT PARAMETERS, EXCEPT B, BX, BY, BZ, ARE
C   ALTERED
C
C THE PARAMETERS NX, NY, NZ, X, Y, Z, VX, VY, VZ AND SIGMA
C SHOULD BE INPUT UNALTERED FROM THE OUTPUT OF VAL3B1
C (SURFB1, CURVB1).
C
C ON OUTPUT--
C
C   B IS THE INPUT ARRAY INCREASED BY THE INTEGRALS OF THE
C   PRODUCTS OF NXA-,NYA-,NZA-PARTIAL DERIVATIVES AND
C   NXB-,NYB-,NZB-PARTIAL DERIVATIVES OF THE SPLINE UNDER
C   TENSION BASIS FUNCTIONS, MULTIPLIED BY WEIGHT.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES QSPL, QINT,
C AND SNHCSH.
C
      EXTERNAL QSPL
C
C--------------------------------------------------------------
C
C OTHER VARIABLES USED INSIDE THE SUBROUTINE QSPL
C
      INTEGER IX,JX,KX,MX,IY,JY,KY,MY,IZ,JZ,KZ,MZ,II,JJ,KK,MM
C
C   THE MATRIX ELEMENT B(II,JJ) IS LOCATED IN THE ARRAY ELEMENT
C   B(KK), WHERE
C     FOR SYMMETRIC MATRIX B, II.LE.JJ :
C       KK= (JJ-1)*JJ/2+II
C     FOR SYMMETRIC MATRIX B, II.GT.JJ :
C       KK= (II-1)*II/2+JJ
C     FOR NONSYMMETRIC MATRIX B :
C       KK= (JJ-1)*NN+II
C     WITH NN=NX*NY*NZ BEING THE DIMENSION OF THE MATRIX B.
C
C   THE MATRIX ELEMENT BX(IX,JX) IS LOCATED IN THE ARRAY ELEMENT
C   BX(KX). THE MEANING OF IX,JX,KX IS SIMILAR AS THE MEANING
C   OF II,JJ,KK IN THE CASE OF MATRIX B.
C
C   THE MATRIX ELEMENT BY(IY,JY) IS LOCATED IN THE ARRAY ELEMENT
C   BZ(KY). THE MEANING OF IY,JY,KY IS SIMILAR AS THE MEANING
C   OF II,JJ,KK IN THE CASE OF MATRIX B.
C
C   THE MATRIX ELEMENT BZ(IZ,JZ) IS LOCATED IN THE ARRAY ELEMENT
C   BZ(KZ). THE MEANING OF IZ,JZ,KZ IS SIMILAR AS THE MEANING
C   OF II,JJ,KK IN THE CASE OF MATRIX B.
C
C   MM, MX, MY, MZ ARE AUXILIARY VARIABLES CONSIDERING THE
C   SYMMETRY OF THE MATRICES B, BX, BY, BZ.
C
C---------------------------------------------------------------
C
C     SCALAR PRODUCTS OF B-SPLINES IN X-DIRECTION
      KX= 0
      MX= NX
      DO 12 JX=1,NX
C       IS BX SYMMETRIC MATRIX ?
        IF(NAX.EQ.NBX) MX=JX
        DO 11 IX=1,MX
          KX= KX+1
          CALL QSPL(NAX,NBX,IX,JX,NX,X,VX,SIGMA,BX(KX))
   11   CONTINUE
   12 CONTINUE
C
C     SCALAR PRODUCTS OF B-SPLINES IN Y-DIRECTION
      KY= 0
      MY= NY
      DO 14 JY=1,NY
C       IS BY SYMMETRIC MATRIX ?
        IF(NAY.EQ.NBY) MY=JY
        DO 13 IY=1,MY
          KY= KY+1
          CALL QSPL(NAY,NBY,IY,JY,NY,Y,VY,SIGMA,BY(KY))
   13   CONTINUE
   14 CONTINUE
C
C     SCALAR PRODUCTS OF B-SPLINES IN Z-DIRECTION
      KZ= 0
      MZ= NZ
      DO 16 JZ=1,NZ
C       IS BZ SYMMETRIC MATRIX ?
        IF(NAZ.EQ.NBZ) MZ=JZ
        DO 15 IZ=1,MZ
          KZ= KZ+1
          CALL QSPL(NAZ,NBZ,IZ,JZ,NZ,Z,VZ,SIGMA,BZ(KZ))
   15   CONTINUE
   16 CONTINUE
C
C     SCALAR PRODUCTS OF 3-D B-SPLINES
C     IS B SYMMETRIC MATRIX ?
        IF(NAX.EQ.NBX.AND.NAY.EQ.NBY.AND.NAZ.EQ.NBZ) THEN
          MM= 1
        ELSE
          MM= 0
        END IF
      KK= 0
      JJ= 0
      DO 27 JZ=1,NZ
        DO 26 JY=1,NY
          DO 25 JX=1,NX
            JJ= JJ+1
            II= 0
C           IS BZ SYMMETRIC MATRIX ?
            IF(NAZ.EQ.NBZ) THEN
              KZ= (JZ-1)*JZ/2
            ELSE
              KZ= (JZ-1)*NZ
            END IF
            DO 23 IZ=1,NZ
              KZ= KZ+1
C             SUBDIAGONAL ELEMENT OF MATRIX BZ
              IF(NAZ.EQ.NBZ.AND.IZ.GT.JZ) KZ=KZ+IZ-2
C             IS BY SYMMETRIC MATRIX ?
              IF(NAY.EQ.NBY) THEN
                KY= (JY-1)*JY/2
              ELSE
                KY= (JY-1)*NY
              END IF
              DO 22 IY=1,NY
                KY= KY+1
C               SUBDIAGONAL ELEMENT OF MATRIX BY
                IF(NAY.EQ.NBY.AND.IY.GT.JY) KY=KY+IY-2
C               IS BX SYMMETRIC MATRIX ?
                IF(NAX.EQ.NBX) THEN
                  KX= (JX-1)*JX/2
                ELSE
                  KX= (JX-1)*NX
                END IF
                DO 21 IX=1,NX
                  KX= KX+1
C                 SUBDIAGONAL ELEMENT OF MATRIX BX
                  IF(NAX.EQ.NBX.AND.IX.GT.JX) KX=KX+IX-2
                  KK= KK+1
                  B(KK)= B(KK)+WEIGHT*BX(KX)*BY(KY)*BZ(KZ)
                  II= II+1
                  IF(MM*II.GE.JJ) GO TO 24
   21           CONTINUE
   22         CONTINUE
   23       CONTINUE
   24       CONTINUE
   25     CONTINUE
   26   CONTINUE
   27 CONTINUE
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE QSPL(NA,NB,IA,IB,N,X,V,SIGMA,Q)
C
      INTEGER NA,NB,IA,IB,N
      REAL    X(N),V(5,N),SIGMA,Q
C
C                                      COMPLEMENT TO FITPACK
C                                       BY ALAN KAYLOR CLINE
C                                  CODED -- JANUARY 23, 1994
C                                            BY LUDEK KLIMES
C                                   DEPARTMENT OF GEOPHYSICS
C                                 CHARLES UNIVERSITY, PRAGUE
C
C THIS SUBROUTINE EVALUATES THE SOBOLEV SCALAR PRODUCT
C OF SPLINE UNDER TENSION BASIS FUNCTIONS IN ONE VARIABLE
C (THE SOBOLEV SCALAR PRODUCT CONSISTS OF INTEGRALS OF THE
C PRODUCTS OF PARTIAL DERIVATIVES OF THE TWO ARGUMENT FUNCTIONS)
C
C ON INPUT--
C
C   NA IS THE ORDER OF THE PARTIAL DERIVATIVE OF
C   THE FIRST ARGUMENT FUNCTION IN THE SCALAR PRODUCT.
C
C   NB IS THE ORDER OF THE PARTIAL DERIVATIVE OF
C   THE SECOND ARGUMENT FUNCTION IN THE SCALAR PRODUCT.
C
C   IA IS THE INDEX OF THE FIRST ARGUMENT FUNCTION
C   (1.LE.IA.LE.N).
C
C   IB IS THE INDEX OF THE SECOND ARGUMENT FUNCTION
C   (1.LE.IB.LE.N).
C
C   N IS THE NUMBER OF GRID POINTS.
C   (N SHOULD BE AT LEAST 1)
C
C   X IS THE ARRAY OF THE N COORDINATES OF GRID POINTS.
C   THESE SHOULD BE STRICTLY INCREASING.
C
C   V IS THE ARRAY OF LENGTHS 5*N,
C   CONTAINING CERTAIN COEFFICIENTS TO BE USED
C   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. FUCTION 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   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   NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C THE PARAMETERS N, X, V, AND SIGMA
C SHOULD BE INPUT UNALTERED FROM THE OUTPUT OF VAL3B1
C (SURFB1, CURVB1).
C
C ON OUTPUT--
C
C   Q IS THE INTEGRAL OF THE PRODUCT OF NA-TH PARTIAL
C   DERIVATIVE OF THE IA-TH BASIS FUNCTION  AND
C   NB-TH PARTIAL DERIVATIVE OF THE IB-TH SPLINE UNDER
C   TENSION BASIS FUNCTION.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES QINT, SNHCSH.
C
      EXTERNAL QINT
C
C---------------------------------------------------------------
C
C OTHER VARIABLES USED INSIDE THE SUBROUTINE QSPL:
C
      INTEGER I,J
      REAL    SIGMAP,V1A,V2A,V3A,V4A,V5A,V1B,V2B,V3B,V4B,V5B
C
C   I...INDEX OF THE INTERVAL.
C   J...POSITION OF THE SECOND B-SPLINE WITH RESPECT TO THE
C             INTERVAL I.
C   SIGMAP...DENORMALIZED TENSION FACTOR.
C   V1A,V2A,V3A,V4A,V5A,V1B,V2B,V3B,V4B,V5B...AUXILIARY
C             STORAGE LOCATIONS FOR V(1,IA),...,V(5,IB).
C
C---------------------------------------------------------------
C
      IF(N.GT.1) GO TO 10
        Q = 1.
        IF(NA.NE.0.OR.NB.NE.0) Q=0.
        GO TO 90
C
   10 SIGMAP= ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
      V1A= V(1,IA)
      V2A= V(2,IA)
      V3A= V(3,IA)
      V4A= V(4,IA)
      V5A= V(5,IA)
      V1B= V(1,IB)
      V2B= V(2,IB)
      V3B= V(3,IB)
      V4B= V(4,IB)
      V5B= V(5,IB)
      Q = 0.
C
      I = IA-2
      IF(I.LT.1) GO TO 20
      J = I-IB+3
      IF(J.LT.1) GO TO 20
      IF(J.GT.4) GO TO 90
      GO TO (11,12,13,14),J
   11 CALL QINT(X(I),X(I+1),0. ,0. ,V4A,V1A,NA,
     *                      0. ,0. ,V4B,V1B,NB,SIGMAP,Q)
        GO TO 20
   12 CALL QINT(X(I),X(I+1),0. ,0. ,V4A,V1A,NA,
     *                      V4B,V1B,1. ,V2B,NB,SIGMAP,Q)
        GO TO 20
   13 CALL QINT(X(I),X(I+1),0. ,0. ,V4A,V1A,NA,
     *                      1. ,V2B,V5B,V3B,NB,SIGMAP,Q)
        GO TO 20
   14 CALL QINT(X(I),X(I+1),0. ,0. ,V4A,V1A,NA,
     *                      V5B,V3B,0. ,0. ,NB,SIGMAP,Q)
C
   20 I = IA-1
      IF(I.LT.1) GO TO 30
      J = I-IB+3
      IF(J.LT.1) GO TO 30
      IF(J.GT.4) GO TO 90
      GO TO (21,22,23,24),J
   21 CALL QINT(X(I),X(I+1),V4A,V1A,1. ,V2A,NA,
     *                      0. ,0. ,V4B,V1B,NB,SIGMAP,Q)
        GO TO 30
   22 CALL QINT(X(I),X(I+1),V4A,V1A,1. ,V2A,NA,
     *                      V4B,V1B,1. ,V2B,NB,SIGMAP,Q)
        GO TO 30
   23 CALL QINT(X(I),X(I+1),V4A,V1A,1. ,V2A,NA,
     *                      1. ,V2B,V5B,V3B,NB,SIGMAP,Q)
        GO TO 30
   24 CALL QINT(X(I),X(I+1),V4A,V1A,1. ,V2A,NA,
     *                      V5B,V3B,0. ,0. ,NB,SIGMAP,Q)
C
   30 I = IA
      IF(I.GE.N) GO TO 90
      J = I-IB+3
      IF(J.LT.1) GO TO 40
      IF(J.GT.4) GO TO 90
      GO TO (31,32,33,34),J
   31 CALL QINT(X(I),X(I+1),1. ,V2A,V5A,V3A,NA,
     *                      0. ,0. ,V4B,V1B,NB,SIGMAP,Q)
        GO TO 40
   32 CALL QINT(X(I),X(I+1),1. ,V2A,V5A,V3A,NA,
     *                      V4B,V1B,1. ,V2B,NB,SIGMAP,Q)
        GO TO 40
   33 CALL QINT(X(I),X(I+1),1. ,V2A,V5A,V3A,NA,
     *                      1. ,V2B,V5B,V3B,NB,SIGMAP,Q)
        GO TO 40
   34 CALL QINT(X(I),X(I+1),1. ,V2A,V5A,V3A,NA,
     *                      V5B,V3B,0. ,0. ,NB,SIGMAP,Q)
C
   40 I = IA+1
      IF(I.GE.N) GO TO 90
      J = I-IB+3
      IF(J.LT.1) GO TO 90
      IF(J.GT.4) GO TO 90
      GO TO (41,42,43,44),J
   41 CALL QINT(X(I),X(I+1),V5A,V3A,0. ,0. ,NA,
     *                      0. ,0. ,V4B,V1B,NB,SIGMAP,Q)
        GO TO 90
   42 CALL QINT(X(I),X(I+1),V5A,V3A,0. ,0. ,NA,
     *                      V4B,V1B,1. ,V2B,NB,SIGMAP,Q)
        GO TO 90
   43 CALL QINT(X(I),X(I+1),V5A,V3A,0. ,0. ,NA,
     *                      1. ,V2B,V5B,V3B,NB,SIGMAP,Q)
        GO TO 90
   44 CALL QINT(X(I),X(I+1),V5A,V3A,0. ,0. ,NA,
     *                      V5B,V3B,0. ,0. ,NB,SIGMAP,Q)
C
   90 CONTINUE
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE QINT(X1,X2,FA1,DA1,FA2,DA2,NA,
     *                      FB1,DB1,FB2,DB2,NB,SIGMAP,Q)
C
      INTEGER NA,NB
      REAL    X1,X2,FA1,DA1,FA2,DA2,FB1,DB1,FB2,DB2,SIGMAP,Q
C
C                                      COMPLEMENT TO FITPACK
C                                       BY ALAN KAYLOR CLINE
C                                  CODED -- JANUARY 23, 1994
C                                            BY LUDEK KLIMES
C                                   DEPARTMENT OF GEOPHYSICS
C                                 CHARLES UNIVERSITY, PRAGUE
C
C THIS SUBROUTINE EVALUATES THE INTEGRAL OF THE PRODUCT
C OF THE GIVEN DERIVATIVES OF THE TWO GIVEN CUBIC FUNCTIONS
C OR SPLINE UNDER TENSION BASIS FUNCTIONS IN ONE VARIABLE,
C OVER A SINGLE SPECIFIED INTERVAL.
C
C ON INPUT--
C
C   X1, X2 ENDPOINTS OF THE GIVEN INTERVAL.
C
C   FA1, DA1  FUNCTION VALUE AND SECOND DERIVATIVE OF THE
C   FIRST GIVEN FUNCTION AT X1.
C
C   FA2, DA2  FUNCTION VALUE AND SECOND DERIVATIVE OF THE
C   FIRST GIVEN FUNCTION AT X2.
C
C   NA IS THE ORDER OF THE PARTIAL DERIVATIVE OF
C   THE FIRST ARGUMENT FUNCTION IN THE SCALAR PRODUCT.
C
C   FB1, DB1, FB2, DB2  THE SAME AS FA1, DA1, FA2, DA2, BUT
C   FOR THE SECOND GIVEN FUNCTION.
C
C   NB IS THE ORDER OF THE PARTIAL DERIVATIVE OF
C   THE SECOND ARGUMENT FUNCTION IN THE SCALAR PRODUCT.
C
C   SIGMAP  IS THE DENORMALIZED TENSION FACTOR.
C
C AND
C
C   NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C ON OUTPUT--
C
C   Q IS THE INTEGRAL OF THE PRODUCT OF NA-TH PARTIAL
C   DERIVATIVE OF THE FIRST FUNCTION  AND
C   NB-TH PARTIAL DERIVATIVE OF THE SECOND FUNCTION,
C   OVER THE INTERVAL X1,X2.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULE SNHCSH.
C
      EXTERNAL SNHCSH
C
C---------------------------------------------------------------
C
C OTHER VARIABLES USED INSIDE THE SUBROUTINE QINT:
C
      INTEGER MA,MB,M
      REAL    QQ,H,SH,CH,SH1,CH1,SIGMA2
      REAL    A1,A2,A3,A4,B1,B2,B3,B4,AB11,AB21,AB12,AB22
C
C---------------------------------------------------------------
C
      MA= MOD(NA,2)
      MB= MOD(NB,2)
      M = MA+MA+MB+1
      QQ= 0.
C
      IF(SIGMAP.NE.0.) GO TO 40
C
C     NO TENSION:
      H = X2-X1
      IF(NA.LE.3.AND.NB.LE.3) GO TO 1
        GO TO 91
    1 IF(NA.LE.1) GO TO 3
C       COEFFICIENTS OF LINEAR FUNCTION
        A3= DA2/H
        A4=-DA1/H
        IF(NB.LE.1) GO TO 2
C         COEFFICIENTS OF LINEAR FUNCTION
          B3= DB2/H
          B4=-DB1/H
          GO TO 80
    2   CONTINUE
C         COEFFICIENTS OF CUBIC AND LINEAR FUNCTIONS
          B1= DB2/H
          B2=-DB1/H
          B3= FB2/H-DB2*H/6.
          B4=-FB1/H+DB1*H/6.
          GO TO 30
    3 CONTINUE
C       COEFFICIENTS OF CUBIC AND LINEAR FUNCTIONS
        A1= DA2/H
        A2=-DA1/H
        A3= FA2/H-DA2*H/6.
        A4=-FA1/H+DA1*H/6.
        IF(NB.LE.1) GO TO 4
C         COEFFICIENTS OF LINEAR FUNCTION
          B3= DB2/H
          B4=-DB1/H
          GO TO 20
    4   CONTINUE
C         COEFFICIENTS OF CUBIC AND LINEAR FUNCTIONS
          B1= DB2/H
          B2=-DB1/H
          B3= FB2/H-DB2*H/6.
          B4=-FB1/H+DB1*H/6.
C
C     INTEGRALS OF (CUBIC FUNCTION)*(CUBIC FUNCTION):
      GO TO (11,12,13,14),M
C     (EVEN DERIVATIVE)*(EVEN DERIVATIVE)
   11   AB11= (H**7)/252.
        AB21=-(H**7)/5040.
        AB12= AB21
        AB22= AB11
        GO TO 15
C     (EVEN DERIVATIVE)*(ODD DERIVATIVE)
   12   AB11= (H**6)/72.
        AB21=-(H**6)/720.
        AB12=-AB21
        AB22=-AB11
        GO TO 15
C     (ODD DERIVATIVE)*(EVEN DERIVATIVE)
   13   AB11= (H**6)/72.
        AB21= (H**6)/720.
        AB12=-AB21
        AB22=-AB11
        GO TO 15
C     (ODD DERIVATIVE)*(ODD DERIVATIVE)
   14   AB11= (H**5)/20.
        AB21= (H**5)/120.
        AB12= AB21
        AB22= AB11
C     ACCUMULATION OF THE COMPUTED INTEGRAL:
   15 QQ=QQ+A1*(AB11*B1+AB12*B2)+A2*(AB21*B1+AB22*B2)
C
C     INTEGRALS OF (CUBIC FUNCTION)*(LINEAR FUNCTION):
   20 GO TO (21,22,23,24),M
C     (EVEN DERIVATIVE)*(EVEN DERIVATIVE)
   21   AB11= (H**5)/30.
        AB21=-(H**5)/120.
        AB12= AB21
        AB22= AB11
        GO TO 25
C     (EVEN DERIVATIVE)*(ODD DERIVATIVE)
   22   AB11= (H**4)/24.
        AB21=-(H**4)/24.
        AB12=-AB21
        AB22=-AB11
        GO TO 25
C     (ODD DERIVATIVE)*(EVEN DERIVATIVE)
   23   AB11= (H**4)/8.
        AB21= (H**4)/24.
        AB12=-AB21
        AB22=-AB11
        GO TO 25
C     (ODD DERIVATIVE)*(ODD DERIVATIVE)
   24   AB11= (H**3)/6.
        AB21= (H**3)/6.
        AB12= AB21
        AB22= AB11
C     ACCUMULATION OF THE COMPUTED INTEGRAL:
   25 QQ=QQ+A1*(AB11*B3+AB12*B4)+A2*(AB21*B3+AB22*B4)
C
      IF(NB.GT.1) GO TO 80
C
C     INTEGRALS OF (LINEAR FUNCTION)*(CUBIC FUNCTION):
   30 GO TO (31,32,33,34),M
C     (EVEN DERIVATIVE)*(EVEN DERIVATIVE)
   31   AB11= (H**5)/30.
        AB21=-(H**5)/120.
        AB12= AB21
        AB22= AB11
        GO TO 35
C     (EVEN DERIVATIVE)*(ODD DERIVATIVE)
   32   AB11= (H**4)/8.
        AB21=-(H**4)/24.
        AB12=-AB21
        AB22=-AB11
        GO TO 35
C     (ODD DERIVATIVE)*(EVEN DERIVATIVE)
   33   AB11= (H**4)/24.
        AB21= (H**4)/24.
        AB12=-AB21
        AB22=-AB11
        GO TO 35
C     (ODD DERIVATIVE)*(ODD DERIVATIVE)
   34   AB11= (H**3)/6.
        AB21= (H**3)/6.
        AB12= AB21
        AB22= AB11
C     ACCUMULATION OF THE COMPUTED INTEGRAL:
   35 QQ=QQ+A3*(AB11*B1+AB12*B2)+A4*(AB21*B1+AB22*B2)
      GO TO 80
C
C
C     NONZERO TENSION:
   40 H = SIGMAP*(X2-X1)
      CALL SNHCSH(SH1,CH1,H,0)
      SH= SH1+H
      CH= CH1+1.
      SIGMA2= SIGMAP*SIGMAP
C
C     COEFFICIENTS OF HYPERBOLIC FUNCTIONS (MULTIPLIED BY SH)
      A1= DA2/SIGMA2
      A2=-DA1/SIGMA2
      B1= DB2/SIGMA2
      B2=-DB1/SIGMA2
C
C     DOUBLED
C     INTEGRALS OF (HYPERBOLIC FUNCTION)*(HYPERBOLIC FUNCTION):
      GO TO (51,52,53,54),M
C     (EVEN DERIVATIVE)*(EVEN DERIVATIVE)
   51   AB11= CH*SH1+H*CH1
        AB21=    SH1-H*CH1
        AB12= AB21
        AB22= AB11
        GO TO 55
C     (EVEN DERIVATIVE)*(ODD DERIVATIVE)
   52   AB11= SH*SH
        AB21= -H*SH
        AB12=-AB21
        AB22=-AB11
        GO TO 55
C     (ODD DERIVATIVE)*(EVEN DERIVATIVE)
   53   AB11= SH*SH
        AB21=  H*SH
        AB12=-AB21
        AB22=-AB11
        GO TO 55
C     (ODD DERIVATIVE)*(ODD DERIVATIVE)
   54   AB11= SH*CH+H
        AB21= SH+H*CH
        AB12= AB21
        AB22= AB11
C     ACCUMULATION OF THE COMPUTED INTEGRAL:
   55 QQ=QQ+(A1*(AB11*B1+AB12*B2)+A2*(AB21*B1+AB22*B2))/(2.*SH*SH)
C
      IF(NB.GT.1) GO TO 70
C
C     COEFFICIENTS OF LINEAR FUNCTION
      B3= ( FB2-B1)/H
      B4= (-FB1-B2)/H
C
C     INTEGRALS OF (HYPERBOLIC FUNCTION)*(LINEAR FUNCTION):
      GO TO (61,62,63,64),M
C     (EVEN DERIVATIVE)*(EVEN DERIVATIVE)
   61   AB11= H*CH1-SH1
        AB21=      -SH1
        AB12= AB21
        AB22= AB11
        GO TO 65
C     (EVEN DERIVATIVE)*(ODD DERIVATIVE)
   62   AB11= CH1
        AB21=-CH1
        AB12=-AB21
        AB22=-AB11
        GO TO 65
C     (ODD DERIVATIVE)*(EVEN DERIVATIVE)
   63   AB11= H*SH-CH1
        AB21=      CH1
        AB12=-AB21
        AB22=-AB11
        GO TO 65
C     (ODD DERIVATIVE)*(ODD DERIVATIVE)
   64   AB11= SH
        AB21= SH
        AB12= AB21
        AB22= AB11
C     ACCUMULATION OF THE COMPUTED INTEGRAL:
   65 QQ=QQ+(A1*(AB11*B3+AB12*B4)+A2*(AB21*B3+AB22*B4))/SH
C
   70 IF(NA.GT.1) GO TO 90
C
C     COEFFICIENTS OF LINEAR FUNCTION
      A3= ( FA2-A1)/H
      A4= (-FA1-A2)/H
C
C     INTEGRALS OF (LINEAR FUNCTION)*(HYPERBOLIC FUNCTION):
      GO TO (71,72,73,74),M
C     (EVEN DERIVATIVE)*(EVEN DERIVATIVE)
   71   AB11= H*CH1-SH1
        AB21=      -SH1
        AB12= AB21
        AB22= AB11
        GO TO 75
C     (EVEN DERIVATIVE)*(ODD DERIVATIVE)
   72   AB11= H*SH-CH1
        AB21=     -CH1
        AB12=-AB21
        AB22=-AB11
        GO TO 75
C     (ODD DERIVATIVE)*(EVEN DERIVATIVE)
   73   AB11= CH1
        AB21= CH1
        AB12=-AB21
        AB22=-AB11
        GO TO 75
C     (ODD DERIVATIVE)*(ODD DERIVATIVE)
   74   AB11= SH
        AB21= SH
        AB12= AB21
        AB22= AB11
C     ACCUMULATION OF THE COMPUTED INTEGRAL:
   75 QQ=QQ+(A3*(AB11*B1+AB12*B2)+A4*(AB21*B1+AB22*B2))/SH
C
      IF(NB.GT.1) GO TO 90
C
C     INTEGRALS OF (LINEAR FUNCTION)*(LINEAR FUNCTION):
   80 GO TO (81,82,83,84),M
C     (EVEN DERIVATIVE)*(EVEN DERIVATIVE)
   81   AB11= (H**3)/3.
        AB21=-(H**3)/6.
        AB12= AB21
        AB22= AB11
        GO TO 85
C     (EVEN DERIVATIVE)*(ODD DERIVATIVE)
   82   AB11= (H**2)/2.
        AB21=-AB11
        AB12=-AB21
        AB22=-AB11
        GO TO 85
C     (ODD DERIVATIVE)*(EVEN DERIVATIVE)
   83   AB11= (H**2)/2.
        AB21= AB11
        AB12=-AB21
        AB22=-AB11
        GO TO 85
C     (ODD DERIVATIVE)*(ODD DERIVATIVE)
   84   AB11= H
        AB21= H
        AB12= AB21
        AB22= AB11
C     ACCUMULATION OF THE COMPUTED INTEGRAL:
   85 QQ=QQ+A3*(AB11*B3+AB12*B4)+A4*(AB21*B3+AB22*B4)
C
C     TRANSFORMATION FROM INDEPENDENT VARIABLE SIGMAP*X TO X
   90 IF(SIGMAP.NE.0.) QQ=QQ*SIGMAP**(NA+NB-1)
   91 Q= Q+QQ
      RETURN
      END
C
C=======================================================================
C
