C
      SUBROUTINE MODEL (MTEXT,LU)
C
C     APPROXIMATION OF INTERFACES AND VELOCITY DISTRIBUTION
C     IN INDIVIDUAL LAYERS (B-SPLINE APPROXIMATION)
C
      CHARACTER*80 MTEXT
      INTEGER NX,NY,NZ,IERR
      REAL X1(10),Y1(10),Z1(10),WG1(10),WG2(10,10),W(10,10,10),C(1000),
     * VX(5,10),VY(5,10),VZ(5,10),TEMP(100),SIGMA
      DIMENSION AUX(12),DEP(6),Y(2),IPR(101)
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      COMMON /AUXX/  MMX(20),MMY(20),MMXY(20)
      COMMON /EPAR/ aij(21,100),aijx(5,21,100),aijy(5,21,100),
     * aijz(5,21,100),XX(200),YY(200),ZZ(200),SSIGMA(20),
     * NXX(20),NYY(20),NZZ(20),NXXYZ(20)
      COMMON /DENS/ RHO(20)
      COMMON /INTRF/ Z(1000),SX(350),SY(350),NNX(20),NNY(20),BRD(4),
     1 NINT,XINTA
C
      NW1=10
      NW2=10
      READ(LU,101) MPRINT,NINT
      IF(MPRINT.GE.0)WRITE(LOU,101) MPRINT,NINT
      NLAY=NINT-1
C
C  INPUT FOR 3D INTERFACES
C
      MX2=0
      MY2=0
      MXY2=0
      DO 2 I=1,NINT
      MX1=MX2+1
      MY1=MY2+1
      MXY1=MXY2+1
      READ(LU,101) MX,MY
      IF(MPRINT.GE.0) WRITE(LOU,101) MX,MY
      NNX(I)=MX
      NNY(I)=MY
      MX2=MX1+MX-1
      MY2=MY1+MY-1
      MXY2=MXY1+MX*MY-1
      READ(LU,102)(SX(J),J=MX1,MX2)
      READ(LU,102)(SY(J),J=MY1,MY2)
      IF(MPRINT.GE.0) WRITE(LOU,102)(SX(J),J=MX1,MX2)
      IF(MPRINT.GE.0) WRITE(LOU,102)(SY(J),J=MY1,MY2)
      M1=MXY1
      DO 1 L=1,MX
      M2=M1+MY-1
      READ(LU,102)(Z(J),J=M1,M2)
      IF(MPRINT.GE.0) WRITE(LOU,102)(Z(J),J=M1,M2)
    1 M1=M2+1
C
C  DETERMINATION OF COEFFICIENTS OF BICUBIC SPLINES
C  APPROXIMATING INTERFACES
C
      CALL BIAP(MX1,MX,MY1,MY,MXY1)
      MMX(I)=MX1
      MMY(I)=MY1
      MMXY(I)=MXY1
    2 CONTINUE
      NB1=MMX(1)
      NB2=MMX(2)-1
      BRD(1)=SX(NB1)
      BRD(2)=SX(NB2)
      NB1=MMY(1)
      NB2=MMY(2)-1
      BRD(3)=SY(NB1)
      BRD(4)=SY(NB2)
C
C  INPUT DATA FOR PRINTER PLOT OF 3-D INTERFACES
C
      DO 6 INTR=1,NINT
      IF(MPRINT.GE.1) WRITE(LOU,109) INTR
      READ(LU,102) ZMIN,ZMAX
      IF(MPRINT.GE.0) WRITE(LOU,102) ZMIN,ZMAX
      IF(MPRINT.EQ.0)GO TO 6
C
C  NUMERICAL FORM OF 3-D INTERFACES
C
      ZMM=ZMAX-ZMIN
      ZMMM=ZMM/10.
      IF(MPRINT.GE.0) WRITE(LOU,102) ZMIN,ZMAX,ZMMM
      DY=(BRD(4)-BRD(3))/50.
      DX=(BRD(2)-BRD(1))/5.
      Y(2)=BRD(3)-DY
      K1=1
      AUX(1)=BRD(1)
      DO 3 I=2,6
      AUX(I)=AUX(I-1)+DX
    3 CONTINUE
      IF(MPRINT.GE.1) WRITE(LOU,103)(AUX(I),I=1,6)
      DX=(BRD(2)-BRD(1))/100.
      MAUX=0
      NDER=0
      DO 5 K=1,51
      Y(2)=Y(2)+DY
      Y(1)=BRD(1)-DX
      DO 4 L=1,101
      Y(1)=Y(1)+DX
      CALL DISC(Y,DEP)
      AUX1=10.*(DEP(1)-ZMIN)/ZMM
      IPR(L)=IFIX(AUX1)
      IF(AUX1.LT.0.0.OR.AUX1.GT.10) IPR(L)=11
    4 CONTINUE
C
C  PRINTER PLOT OF 3-D INTERFACES
C
      IF(K1.EQ.K.AND.MPRINT.GE.1) WRITE(LOU,104)Y(2),IPR
      IF(K1.NE.K.AND.MPRINT.GE.1) WRITE(LOU,105)IPR
      IF(K1.EQ.K)K1=K1+10
    5 CONTINUE
C
C  END OF LOOP OVER ALL INTERFACES
C
    6 CONTINUE
C
C  INPUT OF ELASTIC PARAMETERS
C
      READ(LU,101)ISQRT,IRHO
      IF(MPRINT.EQ.0)WRITE(LOU,101)ISQRT,IRHO
      IF(ISQRT.EQ.0.AND.MPRINT.GT.0) WRITE(LOU,114)
      IF(ISQRT.NE.0.AND.MPRINT.GT.0) WRITE(LOU,113)
      if(irho.ne.0)then
        READ(LU,102)(RHO(L),L=1,NLAY)
        IF(MPRINT.GE.0)WRITE(LOU,102)(RHO(L),L=1,NLAY)
        DO 7 L=1,NLAY
        IF(ABS(RHO(L)).LT..00001)RHO(L)=1.
    7   CONTINUE
      end if
C
C  INPUT OF MATRIX OF ELASTIC PARAMETERS IN COMPRESSED FORM.
C  ELEMENTS OF THE MATRIX HAVE TO BE GIVEN IN (KM**2/SEC**2)
C  WHICH CORRESPONDS TO ELASTIC PARAMETERS DIVIDED BY DENSITY
C
      nxxu=0
      nyyu=0
      nzzu=0
      nxyzu=0
      DO 19 L=1,NLAY
      READ(LU,'(I10,f10.5)')IANI(L),sigma
      IF(MPRINT.EQ.0)WRITE(LOU,'(I10,f10.5)')IANI(L),sigma
      ssigma(l)=sigma
      READ(LU,100)nx,ny,nz
      WRITE(LOU,100)nx,ny,nz
      READ(LU,102)(x1(i),i=1,nx)
      WRITE(LOU,102)(x1(i),i=1,nx)
      READ(LU,102)(y1(i),i=1,ny)
      WRITE(LOU,102)(y1(i),i=1,ny)
      READ(LU,102)(z1(i),i=1,nz)
      WRITE(LOU,102)(z1(i),i=1,nz)
      nxyz=nx*ny*nz
      do 8 i=1,nx
        ii=nxxu+i
        xx(ii)=x1(i)
    8 continue
      do 9 i=1,ny
        ii=nyyu+i
        yy(ii)=y1(i)
    9 continue
      do 10 i=1,nz
        ii=nzzu+i
        zz(ii)=z1(i)
   10 continue
      IF(IANI(L).NE.0) THEN
        IF(MPRINT.GE.1)WRITE(LOU,115) L
        iii=1
      else
        IF(MPRINT.GE.1)WRITE(LOU,118) L
        iii=15
      end if
      do 18 ii=1,21,iii
        do 12 i=1,nx
          do 11 j=1,ny
            READ(LU,102)(w(i,j,k),k=1,nz)
            WRITE(LOU,102)(w(i,j,k),k=1,nz)
            if(isqrt.ne.0.and.iani(l).eq.0)then
              do 50 k=1,nz
                w(i,j,k)=sqrt(w(i,j,k))
   50         continue
            end if
            if(nx.eq.1.and.ny.eq.1.and.nz.eq.1)go to 11
            if(nx.eq.1)then
              if(ny.eq.1)then
                do 51 k=1,nz
                  wg1(k)=w(1,1,k)
   51           continue
              end if
              if(nz.eq.1)then
                wg1(j)=w(1,j,1)
              end if
              if(ny.ne.1.and.nz.ne.1)then
                do 52 k=1,nz
                  wg2(j,k)=w(1,j,k)
   52           continue
              end if
            end if
            if(ny.eq.1)then
              if(nz.eq.1)then
                wg1(i)=w(i,1,1)
              end if
              if(nx.ne.1.and.nz.ne.1)then
                do 53 k=1,nz
                  wg2(i,k)=w(i,1,k)
   53           continue
              end if
            end if
            if(nz.eq.1)then
              if(nx.ne.1.and.ny.ne.1)then
                wg2(i,j)=w(i,j,1)
              end if
            end if
   11     continue
   12   continue
        if(nx.eq.1.and.ny.eq.1.and.nz.eq.1)then
          aij(ii,nxyzu+1)=w(1,1,1)
          do 54 m=1,5
            aijx(m,ii,nxxu+1)=0.
            aijy(m,ii,nyyu+1)=0.
            aijz(m,ii,nzzu+1)=0.
   54     continue
          go to 18
        end if
        if(nx.eq.1)then
          if(ny.eq.1)then
            call CURVB1 (NZ,Z1,WG1,C,VZ,TEMP,SIGMA,IERR)
          end if
          if(nz.eq.1)then
            call CURVB1 (NY,Y1,WG1,C,VY,TEMP,SIGMA,IERR)
          end if
          if(ny.ne.1.and.nz.ne.1)then
            call SURFB1 (NY,NZ,Y1,Z1,WG2,NW1,C,VY,VZ,TEMP,SIGMA,IERR)
          end if
        end if
        if(ny.eq.1)then
          if(nz.eq.1)then
            call CURVB1 (NX,X1,WG1,C,VX,TEMP,SIGMA,IERR)
          end if
          if(nx.ne.1.and.nz.ne.1)then
            call SURFB1 (NX,NZ,X1,Z1,WG2,NW1,C,VX,VZ,TEMP,SIGMA,IERR)
          end if
        end if
        if(nz.eq.1)then
          if(nx.ne.1.and.ny.ne.1)then
            call SURFB1 (NX,NY,X1,Y1,WG2,NW1,C,VX,VY,TEMP,SIGMA,IERR)
          end if
        end if
        if(nx.ne.1.and.ny.ne.1.and.nz.ne.1)then
          call VAL3B1 (NX,NY,NZ,X1,Y1,Z1,W,NW1,NW2,C,VX,VY,
     *                 VZ,TEMP,SIGMA,IERR)
        end if
        do 13 i=1,nxyz
          kk=nxyzu+i
          aij(ii,kk)=c(i)
   13   continue
        do 17 m=1,5
          do 14 i=1,nx
            kk=nxxu+i
            aijx(m,ii,kk)=vx(m,i)
            if(nx.eq.1)aijx(m,ii,kk)=0.
   14     continue
          do 15 i=1,ny
            kk=nyyu+i
            aijy(m,ii,kk)=vy(m,i)
            if(ny.eq.1)aijy(m,ii,kk)=0.
   15     continue
          do 16 i=1,nz
            kk=nzzu+i
            aijz(m,ii,kk)=vz(m,i)
            if(nz.eq.1)aijz(m,ii,kk)=0.
   16     continue
   17   continue
   18 continue
      nxxu=nxxu+nx
      nyyu=nyyu+ny
      nzzu=nzzu+nz
      nxyzu=nxyzu+nxyz
      nxx(l)=nxxu
      nyy(l)=nyyu
      nzz(l)=nzzu
      nxxyz(l)=nxyzu
   19 CONTINUE
C
      WRITE(LOU,107)MTEXT
C
C  FORMATS
C
  100 FORMAT(26I3)
  101 FORMAT(14I5)
  102 FORMAT(8F10.5)
  103 FORMAT(5X,5(F7.3,13X),F7.3)
  104 FORMAT(F7.3,1X,101I1)
  105 FORMAT(8X,101I1)
  107 FORMAT(////,3X,A///)
  109 FORMAT(///,' INTERFACE NUMBER ',I5)
  111 FORMAT(6F10.5)
  113 FORMAT(//' INTERPOLATION IN SQUARE ROOTS OF ELASTIC PARAMETERES'/)
  114 FORMAT(//' INTERPOLATION IN VALUES OF ELASTIC PARAMETERES'/)
  115 FORMAT(/' LAYER',I4,' IS ANISOTROPIC ',/,' MATRIX OF
     1 ELASTIC PARAMETERS GIVEN IN (KM/SEC)**2:
     2  IT CONTAINS ELASTIC PARAMETERS/DENSITY'/)
  118 FORMAT(/' LAYER',I4,' IS ISOTROPIC'/)
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE PARDIS(Y,IAY)
C
      save X1,Y1,Z1,SIGMA,NX,NY,NZ,nxyz,nxxl,nyyl,nzzl,nxyzl,llay
      DIMENSION X1(10),Y1(10),Z1(10),C(1000),
     * VX(5,10),VY(5,10),VZ(5,10)
      DIMENSION Y(18),E(21),EX(21),EY(21),EZ(21)
      COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66,
     1               DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23,
     1               DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44,
     1               DXA45,DXA46,DXA55,DXA56,DXA66,
     1               DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23,
     1               DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44,
     1               DYA45,DYA46,DYA55,DYA56,DYA66,
     1               DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23,
     1               DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44,
     1               DZA45,DZA46,DZA55,DZA56,DZA66,
     1               A2546,A1266,A1355,A1456,A3645,A2344
      COMMON /APROX1/ D(21),DX(21),DY(21),DZ(21),DXX(21),
     1          DXY(21),DXZ(21),DYY(21),DYZ(21),DZZ(21)
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMMON /EPAR/ aij(21,100),aijx(5,21,100),aijy(5,21,100),
     * aijz(5,21,100),XX(200),YY(200),ZZ(200),SSIGMA(20),
     * NXX(20),NYY(20),NZZ(20),NXXYZ(20)
      COMPLEX PS
      COMMON /RAY/   AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),N,IREF,IND,IND1
      COMMON /RAY2/  DRY(3,400)
C
      EQUIVALENCE(E(1),A11),(E(2),A12),(E(3),A13),(E(4),A14),(E(5),A15)
     1           ,(E(6),A16),(E(7),A22),(E(8),A23),(E(9),A24)
     2     ,(E(10),A25),(E(11),A26),(E(12),A33),(E(13),A34),(E(14),A35)
     2     ,(E(15),A36),(E(16),A44),(E(17),A45),(E(18),A46),(E(19),A55)
     2     ,(E(20),A56),(E(21),A66)
     1          ,(EX(1),DXA11),(EX(2),DXA12),(EX(3),DXA13),(EX(4),DXA14)
     1          ,(EX(5),DXA15),(EX(6),DXA16),(EX(7),DXA22),(EX(8),DXA23)
     1       ,(EX(9),DXA24),(EX(10),DXA25),(EX(11),DXA26),(EX(12),DXA33)
     1      ,(EX(13),DXA34),(EX(14),DXA35),(EX(15),DXA36),(EX(16),DXA44)
     1      ,(EX(17),DXA45),(EX(18),DXA46),(EX(19),DXA55),(EX(20),DXA56)
     1      ,(EX(21),DXA66)
      EQUIVALENCE
     1           (EY(1),DYA11),(EY(2),DYA12),(EY(3),DYA13),(EY(4),DYA14)
     1          ,(EY(5),DYA15),(EY(6),DYA16),(EY(7),DYA22),(EY(8),DYA23)
     1       ,(EY(9),DYA24),(EY(10),DYA25),(EY(11),DYA26),(EY(12),DYA33)
     1      ,(EY(13),DYA34),(EY(14),DYA35),(EY(15),DYA36),(EY(16),DYA44)
     1      ,(EY(17),DYA45),(EY(18),DYA46),(EY(19),DYA55),(EY(20),DYA56)
     1      ,(EY(21),DYA66)
     1          ,(EZ(1),DZA11),(EZ(2),DZA12),(EZ(3),DZA13),(EZ(4),DZA14)
     1          ,(EZ(5),DZA15),(EZ(6),DZA16),(EZ(7),DZA22),(EZ(8),DZA23)
     1       ,(EZ(9),DZA24),(EZ(10),DZA25),(EZ(11),DZA26),(EZ(12),DZA33)
     1      ,(EZ(13),DZA34),(EZ(14),DZA35),(EZ(15),DZA36),(EZ(16),DZA44)
     1      ,(EZ(17),DZA45),(EZ(18),DZA46),(EZ(19),DZA55),(EZ(20),DZA56)
     1      ,(EZ(21),DZA66)
C
      if(lay.eq.llay)go to 4
      sigma=ssigma(lay)
      if(lay.gt.1)then
        nxxl=nxx(lay-1)
        nyyl=nyy(lay-1)
        nzzl=nzz(lay-1)
        nxyzl=nxxyz(lay-1)
      else
        nxxl=0
        nyyl=0
        nzzl=0
        nxyzl=0
      end if
      nxxu=nxx(lay)
      nyyu=nyy(lay)
      nzzu=nzz(lay)
      nxyzu=nxxyz(lay)
      NX=nxxu-nxxl
      do 1 i=1,nx
        ii=nxxl+i
        x1(i)=xx(ii)
    1 continue
      NY=nyyu-nyyl
      do 2 i=1,ny
        ii=nyyl+i
        y1(i)=yy(ii)
    2 continue
      NZ=nzzu-nzzl
      do 3 i=1,nz
        ii=nzzl+i
        z1(i)=zz(ii)
    3 continue
      nxyz=nxyzu-nxyzl
    4 continue
      IF(IANI(Lay).NE.0) THEN
        iii=1
      else
        iii=15
      end if
      do 10 ii=1,21,iii
        WX=0.
        WY=0.
        WZ=0.
        WXX=0.
        WXY=0.
        WXZ=0.
        WYY=0.
        WYZ=0.
        WZZ=0.
        do 5 i=1,nxyz
          kk=nxyzl+i
          c(i)=aij(ii,kk)
    5   continue
        do 9 m=1,5
          do 6 i=1,nx
            kk=nxxl+i
            vx(m,i)=aijx(m,ii,kk)
    6     continue
          do 7 i=1,ny
            kk=nyyl+i
            vy(m,i)=aijy(m,ii,kk)
    7     continue
          do 8 i=1,nz
            kk=nzzl+i
            vz(m,i)=aijz(m,ii,kk)
    8     continue
    9   continue
        if(nx.eq.1.and.ny.eq.1.and.nz.eq.1)then
          w=c(1)
          go to 20
        end if
        if(nx.eq.1)then
          if(ny.eq.1)then
            call CURVBD (Y(3),W,WZ,WZZ,NZ,Z1,C,VZ,SIGMA)
          end if
          if(nz.eq.1)then
            call CURVBD (Y(2),W,WY,WYY,NY,Y1,C,VY,SIGMA)
          end if
          if(ny.ne.1.and.nz.ne.1)then
            call SURFBD (Y(2),Y(3),W,WY,WZ,WYY,WYZ,WZZ,NY,NZ,Y1,
     *                   Z1,C,VY,VZ,SIGMA)
          end if
        end if
        if(ny.eq.1)then
          if(nz.eq.1)then
            call CURVBD (Y(1),W,WX,WXX,NX,X1,C,VX,SIGMA)
          end if
          if(nx.ne.1.and.nz.ne.1)then
            call SURFBD (Y(1),Y(3),W,WX,WZ,WXX,WXZ,WZZ,NX,NZ,X1,
     *                   Z1,C,VX,VZ,SIGMA)
          end if
        end if
        if(nz.eq.1)then
          if(nx.ne.1.and.ny.ne.1)then
            call SURFBD (Y(1),Y(2),W,WX,WY,WXX,WXY,WYY,NX,NY,X1,
     *                   Y1,C,VX,VY,SIGMA)
          end if
        end if
        if(nx.ne.1.and.ny.ne.1.and.nz.ne.1)then
          call VAL3BD(Y(1),Y(2),Y(3),W,WX,WY,WZ,WXX,WXY,WYY,WYZ,WZZ,WXZ,
     *                NX,NY,NZ,X1,Y1,Z1,C,VX,VY,VZ,SIGMA)
        end if
   20   E(ii)=W
        EX(ii)=WX
        EY(ii)=WY
        EZ(ii)=WZ
        D(ii)=W
        DX(ii)=WX
        DY(ii)=WY
        DZ(ii)=WZ
        DXX(ii)=WXX
        DXY(ii)=WXY
        DXZ(ii)=WXZ
        DYY(ii)=WYY
        DYZ(ii)=WYZ
        DZZ(ii)=WZZ
        if(isqrt.ne.0.and.iani(lay).eq.0)then
          EE=2.*W
          E(ii)=W*W
          EX(ii)=WX*EE
          EY(ii)=WY*EE
          EZ(ii)=WZ*EE
          D(ii)=W*W
          DX(ii)=WX*EE
          DY(ii)=WY*EE
          DZ(ii)=WZ*EE
          DXX(ii)=WXX*EE+2.*WX*WX
          DXY(ii)=WXY*EE+2.*WX*WY
          DXZ(ii)=WXZ*EE+2.*WX*WZ
          DYY(ii)=WYY*EE+2.*WY*WY
          DYZ(ii)=WYZ*EE+2.*WY*WZ
          DZZ(ii)=WZZ*EE+2.*WZ*WZ
        end if
   10 continue
      llay=lay
      IF(IANI(LAY).EQ.0) GOTO 12
      A2546=A25+A46
      A1266=A12+A66
      A1355=A13+A55
      A1456=A14+A56
      A3645=A36+A45
      A2344=A23+A44
      IF(IAY.EQ.0)RETURN
      DO 11 I=1,21
   11 AY(I+7,N)=E(I)
      return
   12 if(iay.eq.0)return
      ay(8,n)=a11
      AY(9,N)=DXA11
      AY(10,N)=DYA11
      AY(11,N)=DZA11
      ay(12,n)=a44
      AY(13,N)=DXA44
      AY(14,N)=DYA44
      AY(15,N)=DZA44
      RETURN
      END
C
C     *********************************************************
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         COMMON TO ALL THE FOLLOWING PARTS.
C     (1) THE SUBROUTINES PREPARING THE PARAMETERS NECESSARY TO COMPUTE
C         AN INTERPOLATORY 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     (2) THE SUBROUTINES EVALUATING THE VALUE, FIRST AND SECOND PARTIAL
C         DERIVATIVES OF THE INTERPOLATORY FUNCTION AT A GIVEN POINT:
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
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     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 APPEARENCES OF THE STATEMENTS 'CALL VAR2' ARE DENOTED BY
C     'V' IN THE FIRST COLUMN.  THE THREE LINES SHOULD BE REMOVED OR
C     MODIFIED BEFORE COMPILATION.
C
C=======================================================================
C
C PART 0:
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
      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) 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
      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
      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
      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. 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   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
      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
      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
      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
      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)
*       CALL VAR2(II,BX1I,BX(2,I),0.,0.)
    1   CONTINUE
      W = SUM
      WX = SUMX
      WXX = SUMXX
      RETURN
      END
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-COORDINES OF THE POINT
C   AT WHICH THE INTERPOLATION IS TO BE PERFORME .D.
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
*         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
      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
*           CALL VAR2(II+NX*(JJ-1+NY*(KK-1)),BX1I*BY1J*BZ1K,
*    *                     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
      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 N0N-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
      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