      SUBROUTINE NSURF2 (DXMIN,DXMAX,MD,DYMIN,DYMAX,ND,DZ,
     *                   IDZ,M,N,X,Y,Z,IZ,ZP,WORK,SIGMA)
C
      INTEGER MD,ND,IDZ,M,N,IZ
      REAL DXMIN,DXMAX,DYMIN,DYMAX,DZ(IDZ,ND),X(M),Y(N),
     *     Z(IZ,N),ZP(M,N,3),WORK(4,MD),SIGMA
C
C                      FROM THE SPLINE UNDER TENSION PACKAGE
C                                 CODED BY ALAN KAYLOR CLINE
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE MAPS VALUES ONTO A SURFACE AT EVERY POINT
C OF A GRID EQUALLY SPACED IN BOTH X AND Y COORDINATES. THE
C SURFACE INTERPOLATION IS PERFORMED USING A BI-SPLINE
C UNDER TENSION. THE SUBROUTINE SURF1 OR NSURF1 SHOULD BE
C CALLED EARLIER TO DETERMINE CERTAIN NECESSARY PARAMETERS.
C
C ON INPUT--
C
C   DXMIN AND DXMAX CONTAIN THE LOWER AND UPPER LIMITS,
C   RESPECTIVELY, OF THE X-COORDINATES OF THE SECOND GRID.
C
C   MD CONTAINS THE NUMBER OF GRID LINES IN THE X DIRECTION
C   OF THE SECOND GRID (MD .GE. 1).
C
C   DYMIN AND DYMAX CONTAIN THE LOWER AND UPPER LIMITS,
C   RESPECTIVELY, OF THE Y-COORDINATES OF THE SECOND GRID.
C
C   ND CONTAINS THE NUMBER OF GRID LINES IN THE Y DIRECTION
C   OF THE SECOND GRID (ND .GE. 1).
C
C   IDZ CONTAINS THE ROW DIMENSION OF THE ARRAY DZ AS
C   DECLARED IN THE CALLING PROGRAM.
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   WORK IS AN ARRAY OF 4*MD LOCATIONS TO BE USED INTERNALLY
C   FOR WORKSPACE.
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 OR NSURF1.
C
C ON OUTPUT--
C
C   DZ CONTAINS THE MD BY ND ARRAY OF SURFACE VALUES
C   INTERPOLATED AT THE POINTS OF THE SECOND GRID.
C
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C THIS FUNCTION REFERENCES PACKAGE MODULE SNHCSH.
C
C-----------------------------------------------------------
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 FIND INTERVALS OF SECOND X GRID WITH RESPECT TO ORIGINAL X
C GRID
C
      DELTDX = 0.
      IF (MD .GE. 2) DELTDX = (DXMAX-DXMIN)/FLOAT(MD-1)
      LASTI = 1
      DO 3 II = 1,MD
        XII = DXMIN+FLOAT(II-1)*DELTDX
        I = LASTI
    1   I = I+1
        IF (XII .GT. X(I) .AND. I .LT. M) GO TO 1
        IM1 = I-1
        LASTI = IM1
        DEL1 = XII-X(IM1)
        DEL2 = X(I)-XII
        DELS = X(I)-X(IM1)
        WORK(1,II) = DEL2/DELS
        WORK(2,II) = DEL1/DELS
        IF (SIGMAX .NE. 0.) GO TO 2
        TEMP = -DEL1*DEL2/(6.*DELS)
        WORK(3,II) = TEMP*(DEL2+DELS)
        WORK(4,II) = TEMP*(DEL1+DELS)
        GO TO 3
    2   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)
        TEMP = SIGMAX*SIGMAX*DELS*(SINHMS+SIGMAX*DELS)
        WORK(3,II) = (SINHM2*DEL1-DEL2*(2.*(COSHP2+1.)*
     *                SINHP1+SIGMAX*COSHP2*DEL1))/TEMP
        WORK(4,II) = (SINHM1*DEL2-DEL1*(2.*(COSHP1+1.)*
     *                SINHP2+SIGMAX*COSHP1*DEL2))/TEMP
    3   CONTINUE
C
C FIND INTERVALS OF SECOND Y GRID WITH RESPECT TO ORIGINAL Y
C GRID AND PERFORM INTRPOLATION
C
      DELTDY = 0.
      IF (ND .GE. 2) DELTDY = (DYMAX-DYMIN)/FLOAT(ND-1)
      LASTJ = 1
      DO 8 JJ=1,ND
        YJJ = DYMIN+FLOAT(JJ-1)*DELTDY
        J = LASTJ
    4   J = J+1
        IF (YJJ .GT. Y(J) .AND. J .LT. N) GO TO 4
        JM1 = J-1
        LASTJ = JM1
        DEL1 = YJJ-Y(JM1)
        DEL2 = Y(J)-YJJ
        DELS = Y(J)-Y(JM1)
        C1 = DEL2/DELS
        C2 = DEL1/DELS
        IF (SIGMAY .NE. 0.) GO TO 5
        TEMP = -DEL1*DEL2/(6.*DELS)
        C3 = TEMP*(DEL2+DELS)
        C4 = TEMP*(DEL1+DELS)
        GO TO 6
    5   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)
        TEMP = SIGMAY*SIGMAY*DELS*(SINHMS+SIGMAY*DELS)
        C3 = (SINHM2*DEL1-DEL2*(2.*(COSHP2+1.)*
     *        SINHP1+SIGMAY*COSHP2*DEL1))/TEMP
        C4 = (SINHM1*DEL2-DEL1*(2.*(COSHP1+1.)*
     *        SINHP2+SIGMAY*COSHP1*DEL2))/TEMP
    6   LASTI = 0
        DO 8 II=1,MD
          XII = DXMIN+FLOAT(II-1)*DELTDX
          I = MAX0 (1, LASTI)
    7     I = I+1
          IF (XII .GT. X(I) .AND. I .LT. M) GO TO 7
          IM1 = I-1
          IF (IM1 .EQ. LASTI) GO TO 8
          LASTI = IM1
          ZIM1 = C1*Z(IM1,JM1)+C2*Z(IM1,J)
     *           +C3*ZP(IM1,JM1,1)+C4*ZP(IM1,J,1)
          ZI = C1*Z(I,JM1)+C2*Z(I,J)
     *         +C3*ZP(I,JM1,1)+C4*ZP(I,J,1)
          ZXXIM1 = C1*ZP(IM1,JM1,2)+C2*ZP(IM1,J,2)
     *             +C3*ZP(IM1,JM1,3)+C4*ZP(IM1,J,3)
          ZXXI = C1*ZP(I,JM1,2)+C2*ZP(I,J,2)
     *           +C3*ZP(I,JM1,3)+C4*ZP(I,J,3)
    8     DZ(II,JJ) = WORK(1,II)*ZIM1+WORK(2,II)*ZI
     *                +WORK(3,II)*ZXXIM1+WORK(4,II)*ZXXI
      RETURN
      END
