      FUNCTION SURF2 (XX,YY,M,N,X,Y,Z,IZ,ZP,SIGMA)
C
      INTEGER M,N,IZ
      REAL XX,YY,X(M),Y(N),Z(IZ,N),ZP(M,N,3),SIGMA
C
C                      FROM THE SPLINE UNDER TENSION PACKAGE
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 INTERPOLATES A SURFACE AT A GIVEN COORDINATE
C PAIR USING A BI-SPLINE UNDER TENSION. THE SUBROUTINE SURF1
C SHOULD BE CALLED EARLIER TO DETERMINE CERTAIN NECESSARY
C PARAMETERS.
C
C ON INPUT--
C
C   XX AND YY CONTAIN THE X- AND Y-COORDINATES OF THE POINT
C   TO BE MAPPED ONTO THE INTERPOLATING SURFACE.
C
C   M AND N CONTAIN THE NUMBER OF GRID LINES IN THE X- AND
C   Y-DIRECTIONS, RESPECTIVELY, OF THE RECTANGULAR GRID
C   WHICH SPECIFIED THE SURFACE.
C
C   X AND Y ARE ARRAYS CONTAINING THE X- AND Y-GRID VALUES,
C   RESPECTIVELY, EACH IN INCREASING ORDER.
C
C   Z IS A MATRIX CONTAINING THE M * N FUNCTIONAL VALUES
C   CORRESPONDING TO THE GRID VALUES (I. E. Z(I,J) IS THE
C   SURFACE VALUE AT THE POINT (X(I),Y(J)) FOR I = 1,...,M
C   AND J = 1,...,N).
C
C   IZ CONTAINS THE ROW DIMENSION OF THE ARRAY Z AS DECLARED
C   IN THE CALLING PROGRAM.
C
C   ZP IS AN ARRAY OF 3*M*N LOCATIONS STORED WITH THE
C   VARIOUS SURFACE DERIVATIVE INFORMATION DETERMINED BY
C   SURF1.
C
C AND
C
C   SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).
C
C THE PARAMETERS M, N, X, Y, Z, IZ, ZP, AND SIGMA SHOULD BE
C INPUT UNALTERED FROM THE OUTPUT OF SURF1.
C
C ON OUTPUT--
C
C   SURF2 CONTAINS THE INTERPOLATED SURFACE VALUE.
C
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C THIS FUNCTION REFERENCES PACKAGE MODULES INTRVL AND
C SNHCSH.
C
C-----------------------------------------------------------
C
C INLINE ONE DIMENSIONAL CUBIC SPLINE INTERPOLATION
C
      HERMZ (F1,F2,FP1,FP2) = (F2*DEL1+F1*DEL2)/DELS-DEL1*
     *                        DEL2*(FP2*(DEL1+DELS)+
     *                              FP1*(DEL2+DELS))/
     *                        (6.*DELS)
C
C INLINE ONE DIMENSIONAL SPLINE UNDER TENSION INTERPOLATION
C
      HERMNZ (F1,F2,FP1,FP2,SIGMAP) = (F2*DEL1+F1*DEL2)/DELS
     *          +(FP2*(SINHM1*DEL2-DEL1*(2.*(COSHP1+1.)*
     *                           SINHP2+SIGMAP*COSHP1*DEL2))
     *           +FP1*(SINHM2*DEL1-DEL2*(2.*(COSHP2+1.)*
     *                           SINHP1+SIGMAP*COSHP2*DEL1))
     *          )/(SIGMAP*SIGMAP*DELS*(SINHMS+SIGMAP*DELS))
C
C DENORMALIZE TENSION FACTOR IN X AND Y DIRECTION
C
      SIGMAX = ABS(SIGMA)*FLOAT(M-1)/(X(M)-X(1))
      SIGMAY = ABS(SIGMA)*FLOAT(N-1)/(Y(N)-Y(1))
C
C DETERMINE Y INTERVAL
C
      JM1 = INTRVL (YY,Y,N)
      J = JM1+1
C
C DETERMINE X INTERVAL
C
      IM1 = INTRVL (XX,X,M)
      I = IM1+1
      DEL1 = YY-Y(JM1)
      DEL2 = Y(J)-YY
      DELS = Y(J)-Y(JM1)
      IF (SIGMAY .NE. 0.) GO TO 1
C
C PERFORM FOUR INTERPOLATIONS IN Y-DIRECTION
C
      ZIM1 = HERMZ(Z(I-1,J-1),Z(I-1,J),ZP(I-1,J-1,1),
     *                                  ZP(I-1,J,1))
      ZI = HERMZ(Z(I,J-1),Z(I,J),ZP(I,J-1,1),ZP(I,J,1))
      ZXXIM1 = HERMZ(ZP(I-1,J-1,2),ZP(I-1,J,2),
     *                ZP(I-1,J-1,3),ZP(I-1,J,3))
      ZXXI = HERMZ(ZP(I,J-1,2),ZP(I,J,2),
     *              ZP(I,J-1,3),ZP(I,J,3))
      GO TO 2
    1 DELP1 = (DEL1+DELS)/2.
      DELP2 = (DEL2+DELS)/2.
      CALL SNHCSH (SINHM1,DUMMY,SIGMAY*DEL1,-1)
      CALL SNHCSH (SINHM2,DUMMY,SIGMAY*DEL2,-1)
      CALL SNHCSH (SINHMS,DUMMY,SIGMAY*DELS,-1)
      CALL SNHCSH (SINHP1,DUMMY,SIGMAY*DEL1/2.,-1)
      CALL SNHCSH (SINHP2,DUMMY,SIGMAY*DEL2/2.,-1)
      CALL SNHCSH (DUMMY,COSHP1,SIGMAY*DELP1,1)
      CALL SNHCSH (DUMMY,COSHP2,SIGMAY*DELP2,1)
      ZIM1 = HERMNZ(Z(I-1,J-1),Z(I-1,J),ZP(I-1,J-1,1),
     *               ZP(I-1,J,1),SIGMAY)
      ZI = HERMNZ(Z(I,J-1),Z(I,J),ZP(I,J-1,1),ZP(I,J,1),
     *             SIGMAY)
      ZXXIM1 = HERMNZ(ZP(I-1,J-1,2),ZP(I-1,J,2),
     *                 ZP(I-1,J-1,3),ZP(I-1,J,3),SIGMAY)
      ZXXI = HERMNZ(ZP(I,J-1,2),ZP(I,J,2),
     *               ZP(I,J-1,3),ZP(I,J,3),SIGMAY)
C
C PERFORM FINAL INTERPOLATION IN X-DIRECTION
C
    2 DEL1 = XX-X(IM1)
      DEL2 = X(I)-XX
      DELS = X(I)-X(IM1)
      IF (SIGMAX .NE. 0.) GO TO 3
      SURF2 = HERMZ(ZIM1,ZI,ZXXIM1,ZXXI)
      RETURN
    3 DELP1 = (DEL1+DELS)/2.
      DELP2 = (DEL2+DELS)/2.
      CALL SNHCSH (SINHM1,DUMMY,SIGMAX*DEL1,-1)
      CALL SNHCSH (SINHM2,DUMMY,SIGMAX*DEL2,-1)
      CALL SNHCSH (SINHMS,DUMMY,SIGMAX*DELS,-1)
      CALL SNHCSH (SINHP1,DUMMY,SIGMAX*DEL1/2.,-1)
      CALL SNHCSH (SINHP2,DUMMY,SIGMAX*DEL2/2.,-1)
      CALL SNHCSH (DUMMY,COSHP1,SIGMAX*DELP1,1)
      CALL SNHCSH (DUMMY,COSHP2,SIGMAX*DELP2,1)
      SURF2 = HERMNZ(ZIM1,ZI,ZXXIM1,ZXXI,SIGMAX)
      RETURN
      END
