C     ALGORITHM 624 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 4,
C     DEC., 1984, P. 453.
      SUBROUTINE ADNODE (KK,X,Y, IADJ,IEND, IER)
      INTEGER KK, IADJ(1), IEND(KK), IER
      REAL    X(KK), Y(KK)
      LOGICAL SWPTST
      EXTERNAL INDEX
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE ADDS NODE KK TO A TRIANGULATION OF A SET
C OF POINTS IN THE PLANE PRODUCING A NEW TRIANGULATION.  A
C SEQUENCE OF EDGE SWAPS IS THEN APPLIED TO THE MESH,
C RESULTING IN AN OPTIMAL TRIANGULATION.  ADNODE IS PART
C OF AN INTERPOLATION PACKAGE WHICH ALSO PROVIDES ROUTINES
C TO INITIALIZE THE DATA STRUCTURE, PLOT THE MESH, AND
C DELETE ARCS.
C
C INPUT PARAMETERS -   KK - INDEX OF THE NODE TO BE ADDED
C                           TO THE MESH.  KK .GE. 4.
C
C                     X,Y - VECTORS OF COORDINATES OF THE
C                           NODES IN THE MESH.  (X(I),Y(I))
C                           DEFINES NODE I FOR I = 1,..,KK.
C
C                    IADJ - SET OF ADJACENCY LISTS OF NODES
C                           1,..,KK-1.
C
C                    IEND - POINTERS TO THE ENDS OF
C                           ADJACENCY LISTS IN IADJ FOR
C                           EACH NODE IN THE MESH.
C
C IADJ AND IEND MAY BE CREATED BY TRMESH.
C
C KK, X, AND Y ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ADDITION
C                                 OF NODE KK AS THE LAST
C                                 ENTRY.
C
C                           IER - ERROR INDICATOR
C                                 IER = 0 IF NO ERRORS
C                                         WERE ENCOUNTERED.
C                                 IER = 1 IF ALL NODES
C                                         (INCLUDING KK) ARE
C                                         COLLINEAR.
C
C MODULES REFERENCED BY ADNODE - TRFIND, INTADD, BDYADD,
C                                SHIFTD, INDEX, SWPTST,
C                                SWAP
C
C***********************************************************
C
      INTEGER K, KM1, I1, I2, I3, INDKF, INDKL, NABOR1,
     .        IO1, IO2, IN1, INDK1, IND2F, IND21
      REAL    XK, YK
C
C LOCAL PARAMETERS -
C
C K =        LOCAL COPY OF KK
C KM1 =      K - 1
C I1,I2,I3 = VERTICES OF A TRIANGLE CONTAINING K
C INDKF =    IADJ INDEX OF THE FIRST NEIGHBOR OF K
C INDKL =    IADJ INDEX OF THE LAST NEIGHBOR OF K
C NABOR1 =   FIRST NEIGHBOR OF K BEFORE ANY SWAPS OCCUR
C IO1,IO2 =  ADJACENT NEIGHBORS OF K DEFINING AN ARC TO
C              BE TESTED FOR A SWAP
C IN1 =      VERTEX OPPOSITE K -- FIRST NEIGHBOR OF IO2
C              WHICH PRECEDES IO1.  IN1,IO1,IO2 ARE IN
C              COUNTERCLOCKWISE ORDER.
C INDK1 =    INDEX OF IO1 IN THE ADJACENCY LIST FOR K
C IND2F =    INDEX OF THE FIRST NEIGHBOR OF IO2
C IND21 =    INDEX OF IO1 IN THE ADJACENCY LIST FOR IO2
C XK,YK =    X(K), Y(K)
C
      IER = 0
      K = KK
C
C INITIALIZATION
C
      KM1 = K - 1
      XK = X(K)
      YK = Y(K)
C
C ADD NODE K TO THE MESH
C
      CALL TRFIND(KM1,XK,YK,X,Y,IADJ,IEND, I1,I2,I3)
      IF (I1 .EQ. 0) GO TO 5
      IF (I3 .EQ. 0) CALL BDYADD(K,I1,I2, IADJ,IEND )
      IF (I3 .NE. 0) CALL INTADD(K,I1,I2,I3, IADJ,IEND )
C
C INITIALIZE VARIABLES FOR OPTIMIZATION OF THE MESH
C
      INDKF = IEND(KM1) + 1
      INDKL = IEND(K)
      NABOR1 = IADJ(INDKF)
      IO2 = NABOR1
      INDK1 = INDKF + 1
      IO1 = IADJ(INDK1)
C
C BEGIN LOOP -- FIND THE VERTEX OPPOSITE K
C
    1 IND2F = 1
      IF (IO2 .NE. 1) IND2F = IEND(IO2-1) + 1
      IND21 = INDEX(IO2,IO1,IADJ,IEND)
      IF (IND2F .EQ. IND21) GO TO 2
      IN1 = IADJ(IND21-1)
      GO TO 3
C
C IN1 IS THE LAST NEIGHBOR OF IO2
C
    2 IND21 = IEND(IO2)
      IN1 = IADJ(IND21)
      IF (IN1 .EQ. 0) GO TO 4
C
C SWAP TEST -- IF A SWAP OCCURS, TWO NEW ARCS ARE OPPOSITE K
C              AND MUST BE TESTED.  INDK1 AND INDKF MUST BE
C              DECREMENTED.
C
    3 IF ( .NOT. SWPTST(IN1,K,IO1,IO2,X,Y) ) GO TO 4
      CALL SWAP(IN1,K,IO1,IO2, IADJ,IEND )
      IO1 = IN1
      INDK1 = INDK1 - 1
      INDKF = INDKF - 1
      GO TO 1
C
C NO SWAP OCCURRED.  RESET IO2 AND IO1, AND TEST FOR
C   TERMINATION.
C
    4 IF (IO1 .EQ. NABOR1) RETURN
      IO2 = IO1
      INDK1 = INDK1 + 1
      IF (INDK1 .GT. INDKL) INDK1 = INDKF
      IO1 = IADJ(INDK1)
      IF (IO1 .NE. 0) GO TO 1
      RETURN
C
C ALL NODES ARE COLLINEAR
C
    5 IER = 1
      RETURN
      END
      FUNCTION AREA (X,Y,NB,NODES)
      INTEGER NB, NODES(NB)
      REAL    X(1), Y(1)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A SEQUENCE OF NB POINTS IN THE PLANE, THIS
C FUNCTION COMPUTES THE AREA BOUNDED BY THE CLOSED POLY-
C GONAL CURVE WHICH PASSES THROUGH THE POINTS IN THE
C SPECIFIED ORDER.  EACH SIMPLE CLOSED CURVE IS POSITIVELY
C ORIENTED (BOUNDS POSITIVE AREA) IF AND ONLY IF THE POINTS
C ARE SPECIFIED IN COUNTERCLOCKWISE ORDER.  THE LAST POINT
C OF THE CURVE IS TAKEN TO BE THE FIRST POINT SPECIFIED, AND
C THUS THIS POINT NEED NOT BE SPECIFIED TWICE.  HOWEVER, ANY
C POINT MAY BE SPECIFIED MORE THAN ONCE IN ORDER TO DEFINE A
C MULTIPLY CONNECTED DOMAIN.
C   THE AREA OF A TRIANGULATION MAY BE COMPUTED BY CALLING
C AREA WITH VALUES OF NB AND NODES DETERMINED BY SUBROUTINE
C BNODES.
C
C INPUT PARAMETERS -   X,Y - N-VECTORS OF COORDINATES OF
C                            POINTS IN THE PLANE FOR N .GE.
C                            NB.  NODE I HAS COORDINATES
C                            (X(I),Y(I)) FOR I = 1, 2, ...,
C                            N.
C
C                       NB - LENGTH OF NODES.
C
C                    NODES - VECTOR OF NODE INDICES IN THE
C                            RANGE 1 TO N DEFINING THE
C                            POLYGONAL CURVE.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION.
C
C OUTPUT PARAMETER -  AREA - SIGNED AREA BOUNDED BY THE
C                            POLYGONAL CURVE DEFINED
C                            ABOVE.
C
C MODULES REFERENCED BY AREA - NONE
C
C***********************************************************
C
      INTEGER NNB, ND, I
      REAL    A, X0, Y0, DX1, DY1, DX2, DY2
C
C LOCAL PARAMETERS -
C
C NNB =     LOCAL COPY OF NB
C ND =      ELEMENT OF NODES
C I =       DO-LOOP AND NODES INDEX
C A =       PARTIAL SUM OF SIGNED (AND DOUBLED) TRIANGLE
C             AREAS
C X0,Y0 =   X(NODES(1)), Y(NODES(1))
C DX1,DY1 = COMPONENTS OF THE VECTOR NODES(1)-NODES(I) FOR
C             I = 2, 3, ..., NB-1
C DX2,DY2 = COMPONENTS OF THE VECTOR NODES(1)-NODES(I) FOR
C             I = 3, 4, ..., NB
C
      NNB = NB
      A = 0.
      IF (NNB .LT. 3) GO TO 2
C
C INITIALIZATION
C
      ND = NODES(1)
      X0 = X(ND)
      Y0 = Y(ND)
      ND = NODES(2)
      DX1 = X(ND) - X0
      DY1 = Y(ND) - Y0
C
C LOOP ON TRIANGLES (NODES(1),NODES(I),NODES(I+1)),
C   I = 2, 3, ..., NB-1, ADDING TWICE THEIR SIGNED
C   AREAS TO A
C
      DO 1 I = 3,NNB
        ND = NODES(I)
        DX2 = X(ND) - X0
        DY2 = Y(ND) - Y0
        A = A + DX1*DY2 - DX2*DY1
        DX1 = DX2
        DY1 = DY2
    1   CONTINUE
C
C A CONTAINS TWICE THE SIGNED AREA OF THE REGION
C
    2 AREA = A/2.
      RETURN
      END
      SUBROUTINE BDYADD (KK,I1,I2, IADJ,IEND )
      INTEGER KK, I1, I2, IADJ(1), IEND(KK)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE ADDS A BOUNDARY NODE TO A TRIANGULATION
C OF A SET OF KK-1 POINTS IN THE PLANE.  IADJ AND IEND ARE
C UPDATED WITH THE INSERTION OF NODE KK.
C
C INPUT PARAMETERS -   KK - INDEX OF AN EXTERIOR NODE TO BE
C                           ADDED.  KK .GE. 4.
C
C                      I1 - FIRST (RIGHTMOST AS VIEWED FROM
C                           KK) BOUNDARY NODE IN THE MESH
C                           WHICH IS VISIBLE FROM KK - THE
C                           LINE SEGMENT KK-I1 INTERSECTS
C                           NO ARCS.
C
C                      I2 - LAST (LEFTMOST) BOUNDARY NODE
C                           WHICH IS VISIBLE FROM KK.
C
C                    IADJ - SET OF ADJACENCY LISTS OF NODES
C                           IN THE MESH.
C
C                    IEND - POINTERS TO THE ENDS OF
C                           ADJACENCY LISTS IN IADJ FOR
C                           EACH NODE IN THE MESH.
C
C   IADJ AND IEND MAY BE CREATED BY TRMESH AND MUST CONTAIN
C THE VERTICES I1 AND I2.  I1 AND I2 MAY BE DETERMINED BY
C TRFIND.
C
C KK, I1, AND I2 ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ADDITION
C                                 OF NODE KK AS THE LAST
C                                 ENTRY.  NODE KK WILL BE
C                                 CONNECTED TO I1, I2, AND
C                                 ALL BOUNDARY NODES BETWEEN
C                                 THEM.  NO OPTIMIZATION OF
C                                 THE MESH IS PERFORMED.
C
C MODULE REFERENCED BY BDYADD - SHIFTD
C
C INTRINSIC FUNCTIONS CALLED BY BDYADD - MIN0, MAX0
C
C***********************************************************
C
      INTEGER K, KM1, NRIGHT, NLEFT, NF, NL, N1, N2, I,
     .        IMIN, IMAX, KEND, NEXT, INDX
C
C LOCAL PARAMETERS -
C
C K =            LOCAL COPY OF KK
C KM1 =          K - 1
C NRIGHT,NLEFT = LOCAL COPIES OF I1, I2
C NF,NL =        INDICES OF IADJ BOUNDING THE PORTION OF THE
C                  ARRAY TO BE SHIFTED
C N1 =           IADJ INDEX OF THE FIRST NEIGHBOR OF NLEFT
C N2 =           IADJ INDEX OF THE LAST NEIGHBOR OF NRIGHT
C I =            DO-LOOP INDEX
C IMIN,IMAX =    BOUNDS ON DO-LOOP INDEX -- FIRST AND LAST
C                  ELEMENTS OF IEND TO BE INCREMENTED
C KEND =         POINTER TO THE LAST NEIGHBOR OF K IN IADJ
C NEXT =         NEXT BOUNDARY NODE TO BE CONNECTED TO KK
C INDX =         INDEX FOR IADJ
C
      K = KK
      KM1 = K - 1
      NRIGHT = I1
      NLEFT = I2
C
C INITIALIZE VARIABLES
C
      NL = IEND(KM1)
      N1 = 1
      IF (NLEFT .NE. 1) N1 = IEND(NLEFT-1) + 1
      N2 = IEND(NRIGHT)
      NF = MAX0(N1,N2)
C
C INSERT K AS A NEIGHBOR OF MAX(NRIGHT,NLEFT)
C
      CALL SHIFTD(NF,NL,2, IADJ )
      IADJ(NF+1) = K
      IMIN = MAX0(NRIGHT,NLEFT)
      DO 1 I = IMIN,KM1
        IEND(I) = IEND(I) + 2
    1   CONTINUE
C
C INITIALIZE KEND AND INSERT K AS A NEIGHBOR OF
C   MIN(NRIGHT,NLEFT)
C
      KEND = NL + 3
      NL = NF - 1
      NF = MIN0(N1,N2)
      CALL SHIFTD(NF,NL,1, IADJ )
      IADJ(NF) = K
      IMAX = IMIN - 1
      IMIN = MIN0(NRIGHT,NLEFT)
      DO 2 I = IMIN,IMAX
        IEND(I) = IEND(I) + 1
    2   CONTINUE
C
C INSERT NRIGHT AS THE FIRST NEIGHBOR OF K
C
      IADJ(KEND) = NRIGHT
C
C INITIALIZE INDX FOR LOOP ON BOUNDARY NODES BETWEEN NRIGHT
C   AND NLEFT
C
      INDX = IEND(NRIGHT) - 2
    3 NEXT = IADJ(INDX)
      IF (NEXT .EQ. NLEFT) GO TO 4
C
C CONNECT NEXT AND K
C
      KEND = KEND + 1
      IADJ(KEND) = NEXT
      INDX = IEND(NEXT)
      IADJ(INDX) = K
      INDX = INDX - 1
      GO TO 3
C
C INSERT NLEFT AND 0 AS THE LAST NEIGHBORS OF K
C
    4 IADJ(KEND+1) = NLEFT
      KEND = KEND + 2
      IADJ(KEND) = 0
      IEND(K) = KEND
      RETURN
      END
      SUBROUTINE BNODES (N,IADJ,IEND, NB,NA,NT,NODES)
      INTEGER N, IADJ(1), IEND(N), NB, NA, NT, NODES(1)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A TRIANGULATION OF N POINTS IN THE PLANE, THIS
C ROUTINE RETURNS A VECTOR CONTAINING THE INDICES, IN
C COUNTERCLOCKWISE ORDER, OF THE NODES ON THE BOUNDARY OF
C THE CONVEX HULL OF THE SET OF POINTS.
C
C INPUT PARAMETERS -     N - NUMBER OF NODES IN THE MESH.
C
C                     IADJ - SET OF ADJACENCY LISTS OF
C                            NODES IN THE MESH.
C
C                     IEND - POINTERS TO THE ENDS OF
C                            ADJACENCY LISTS IN IADJ FOR
C                            EACH NODE IN THE MESH.
C
C                    NODES - VECTOR OF LENGTH .GE. NB.
C                            (NB .LE. N).
C
C   IADJ AND IEND MAY BE CREATED BY TRMESH AND ARE NOT
C ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETERS -   NB - NUMBER OF BOUNDARY NODES.
C
C                    NA,NT - NUMBER OF ARCS AND TRIANGLES,
C                            RESPECTIVELY, IN THE MESH.
C
C                    NODES - VECTOR OF NB BOUNDARY NODE
C                            INDICES RANGING FROM 1 TO N.
C
C MODULES REFERENCED BY BNODES - NONE
C
C***********************************************************
C
      INTEGER NST, INDL, K, N0, INDF
C
C LOCAL PARAMETERS -
C
C NST =  FIRST ELEMENT OF NODES -- ARBITRARILY CHOSEN
C INDL = IADJ INDEX OF THE LAST NEIGHBOR OF NST
C K =    NODES INDEX
C N0 =   BOUNDARY NODE TO BE ADDED TO NODES
C INDF = IADJ INDEX OF THE FIRST NEIGHBOR OF N0
C
C SET NST TO THE FIRST BOUNDARY NODE ENCOUNTERED
C
      NST = 1
    1 INDL = IEND(NST)
      IF (IADJ(INDL) .EQ. 0) GO TO 2
      NST = NST + 1
      GO TO 1
C
C INITIALIZATION
C
    2 NODES(1) = NST
      K = 1
      N0 = NST
C
C TRAVERSE THE BOUNDARY IN COUNTERCLOCKWISE ORDER
C
    3 INDF = 1
      IF (N0 .GT. 1) INDF = IEND(N0-1) + 1
      N0 = IADJ(INDF)
      IF (N0 .EQ. NST) GO TO 4
      K = K + 1
      NODES(K) = N0
      GO TO 3
C
C TERMINATION
C
    4 NB = K
      NT = 2*N - NB - 2
      NA = NT + N - 1
      RETURN
      END
      SUBROUTINE DELETE (NN,NOUT1,NOUT2, IADJ,IEND, IER)
      INTEGER NN, NOUT1, NOUT2, IADJ(1), IEND(NN), IER
      EXTERNAL INDEX
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE DELETES A BOUNDARY EDGE FROM A TRIANGU-
C LATION OF A SET OF POINTS IN THE PLANE.  IT MAY BE NEC-
C ESSARY TO FORCE CERTAIN EDGES TO BE PRESENT BEFORE CALL-
C ING DELETE (SEE SUBROUTINE EDGE).  NOTE THAT SUBROUTINES
C EDGE, TRFIND, AND THE ROUTINES WHICH CALL TRFIND (ADNODE,
C UNIF, INTRC1, AND INTRC0) SHOULD NOT BE CALLED FOLLOWING
C A DELETION.
C
C INPUT PARAMETERS -    NN - NUMBER OF NODES IN THE TRIAN-
C                            GULATION.
C
C              NOUT1,NOUT2 - PAIR OF ADJACENT NODES ON THE
C                            BOUNDARY DEFINING THE ARC TO
C                            BE REMOVED.  NOUT2 MUST BE THE
C                            LAST NONZERO NEIGHBOR OF NOUT1.
C
C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C                IADJ,IEND - DATA STRUCTURE DEFINING THE
C                            TRIANGULATION (SEE SUBROUTINE
C                            TRMESH).
C
C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE REMOVAL
C                                 OF THE ARC NOUT1-NOUT2
C                                 IF IER .EQ. 0.
C
C                           IER - ERROR INDICATOR
C                                 IER = 0 IF NO ERRORS WERE
C                                         ENCOUNTERED.
C                                 IER = 1 IF NOUT1 OR NOUT2
C                                         IS NOT ON THE
C                                         BOUNDARY.
C                                 IER = 2 IF NOUT1 OR NOUT2
C                                         HAS ONLY 2 NONZERO
C                                         NEIGHBORS.
C                                 IER = 3 IF NOUT2 IS NOT
C                                         THE LAST NEIGHBOR
C                                         OF NOUT1.
C                                 IER = 4 IF A DELETION
C                                         WOULD DIVIDE THE
C                                         MESH INTO TWO
C                                         REGIONS.
C
C MODULES REFERENCED BY DELETE - SHIFTD, INDEX
C
C***********************************************************
C
      INTEGER N, IOUT1, IOUT2, IO1, IO2, IND12, IND21,
     .        ITEMP, IND1F, IND1L, IND2F, IND2L, NEWBD,
     .        INDNF, INDNL, INDN0, INDFP2, INDLM3, NF, NL,
     .        I, IMAX
C
C LOCAL PARAMETERS -
C
C N =           LOCAL COPY OF NN
C IOUT1,IOUT2 = LOCAL COPIES OF NOUT1 AND NOUT2
C IO1,IO2 =     NOUT1,NOUT2 IN ORDER OF INCREASING MAGNITUDE
C IND12 =       INDEX OF IO2 IN THE ADJACENCY LIST FOR IO1
C IND21 =       INDEX OF IO1 IN THE ADJACENCY LIST FOR IO2
C ITEMP =       TEMPORARY STORAGE LOCATION FOR PERMUTATIONS
C IND1F =       IADJ INDEX OF THE FIRST NEIGHBOR OF IO1
C IND1L =       IADJ INDEX OF THE LAST NEIGHBOR OF IO1
C IND2F =       IADJ INDEX OF THE FIRST NEIGHBOR OF IO2
C IND2L =       IADJ INDEX OF THE LAST NEIGHBOR OF IO2
C NEWBD =       THE NEIGHBOR COMMON TO NOUT1 AND NOUT2
C INDNF =       IADJ INDEX OF THE FIRST NEIGHBOR OF NEWBD
C INDNL =       IADJ INDEX OF THE LAST NEIGHBOR OF NEWBD
C INDN0 =       INDEX OF 0 IN THE ADJACENCY LIST FOR NEWBD
C                 BEFORE PERMUTING THE NEIGHBORS
C INDFP2 =      INDNF + 2
C INDLM3 =      INDNL - 3
C NF,NL =       BOUNDS ON THE PORTION OF IADJ TO BE SHIFTED
C I =           DO-LOOP INDEX
C IMAX =        UPPER BOUND ON DO-LOOP FOR SHIFTING IEND
C
      N = NN
      IOUT1 = NOUT1
      IOUT2 = NOUT2
C
C INITIALIZE INDICES
C
      IND1F = 1
      IF (IOUT1 .GT. 1) IND1F = IEND(IOUT1-1) + 1
      IND1L = IEND(IOUT1)
      IND2F = 1
      IF (IOUT2 .GT. 1) IND2F = IEND(IOUT2-1) + 1
      IND2L = IEND(IOUT2)
      NEWBD = IADJ(IND1L-2)
      INDN0 = INDEX(NEWBD,IOUT2,IADJ,IEND)
      INDNL = IEND(NEWBD)
C
C ORDER VERTICES SUCH THAT THE NEIGHBORS OF IO1 PRECEDE
C   THOSE OF IO2
C
      IF (IOUT1 .GT. IOUT2) GO TO 1
      IO1 = IOUT1
      IO2 = IOUT2
      IND12 = IND1L - 1
      IND21 = IND2F
      GO TO 2
    1 IO1 = IOUT2
      IO2 = IOUT1
      IND12 = IND2F
      IND21 = IND1L - 1
C
C CHECK FOR ERRORS
C
    2 IF ( (IADJ(IND1L) .NE. 0) .OR. (IADJ(IND2L) .NE. 0) )
     .   GO TO 21
      IF ( (IND1L-IND1F .LE. 2) .OR. (IND2L-IND2F .LE. 2) )
     .   GO TO 22
      IF (IADJ(IND1L-1) .NE. IOUT2) GO TO 23
      IF (IADJ(INDNL) .EQ. 0) GO TO 24
C
C DELETE THE EDGE IO1-IO2 AND MAKE NEWBD A BOUNDARY NODE
C
      IF (NEWBD .LT. IO1) GO TO 8
      IF (NEWBD .LT. IO2) GO TO 6
C
C THE VERTICES ARE ORDERED IO1, IO2, NEWBD.
C DELETE IO2 AS A NEIGHBOR OF IO1.
C
      NF = IND12 + 1
      NL = IND21 - 1
      CALL SHIFTD(NF,NL,-1, IADJ )
      IMAX = IO2 - 1
      DO 3 I = IO1,IMAX
        IEND(I) = IEND(I) - 1
    3   CONTINUE
C
C DELETE IO1 AS A NEIGHBOR OF IO2
C
      NF = NL + 2
      NL = INDN0
      CALL SHIFTD(NF,NL,-2, IADJ )
      IMAX = NEWBD - 1
      DO 4 I = IO2,IMAX
        IEND(I) = IEND(I) - 2
    4   CONTINUE
C
C SHIFT THE BOTTOM OF IADJ UP 1 LEAVING ROOM FOR 0 AS A
C   NEIGHBOR OF NEWBD
C
      INDN0 = INDN0 - 1
      NF = NL + 1
      NL = IEND(N)
      IF (NF .LE. NL) CALL SHIFTD(NF,NL,-1, IADJ )
      DO 5 I = NEWBD,N
        IEND(I) = IEND(I) - 1
    5   CONTINUE
      GO TO 12
C
C THE VERTICES ARE ORDERED IO1, NEWBD, IO2.
C DELETE IO2 AS A NEIGHBOR OF IO1 LEAVING ROOM FOR 0 AS A
C   NEIGHBOR OF NEWBD.
C
    6 NF = IND12 + 1
      NL = INDN0
      CALL SHIFTD(NF,NL,-1, IADJ )
      IMAX = NEWBD - 1
      DO 7 I = IO1,IMAX
        IEND(I) = IEND(I) - 1
    7   CONTINUE
      GO TO 10
C
C THE VERTICES ARE ORDERED NEWBD, IO1, IO2.
C DELETE IO2 AS A NEIGHBOR OF IO1 LEAVING ROOM FOR 0 AS A
C   NEIGHBOR OF NEWBD.
C
    8 INDN0 = INDN0 + 1
      NF = INDN0
      NL = IND12 - 1
      IF (NF .LE. NL) CALL SHIFTD(NF,NL,1, IADJ )
      IMAX = IO1 - 1
      DO 9 I = NEWBD,IMAX
        IEND(I) = IEND(I) + 1
    9   CONTINUE
C
C DELETE IO1 AS A NEIGHBOR OF IO2
C
   10 NF = IND21 + 1
      NL = IEND(N)
      CALL SHIFTD(NF,NL,-1, IADJ )
      DO 11 I = IO2,N
        IEND(I) = IEND(I) - 1
   11   CONTINUE
C
C PERMUTE THE NEIGHBORS OF NEWBD WITH END-AROUND SHIFTS SO
C   THAT 0 IS THE LAST NEIGHBOR
C
   12 INDNF = 1
      IF (NEWBD .GT. 1) INDNF = IEND(NEWBD-1) + 1
      INDNL = IEND(NEWBD)
      IF (INDN0-INDNF .GE. INDNL-INDN0) GO TO 16
C
C SHIFT UPWARD
C
      IF (INDN0 .GT. INDNF) GO TO 13
      CALL SHIFTD(INDNF+1,INDNL,-1, IADJ )
      GO TO 20
   13 INDFP2 = INDNF + 2
      IF (INDN0 .LT. INDFP2) GO TO 15
      DO 14 I = INDFP2,INDN0
        ITEMP = IADJ(INDNF)
        CALL SHIFTD(INDNF+1,INDNL,-1, IADJ )
        IADJ(INDNL) = ITEMP
   14   CONTINUE
C
C THE LAST SHIFT IS BY 2
C
   15 ITEMP = IADJ(INDNF)
      CALL SHIFTD(INDFP2,INDNL,-2, IADJ )
      IADJ(INDNL-1) = ITEMP
      GO TO 20
C
C SHIFT DOWNWARD
C
   16 IF (INDN0 .EQ. INDNL) GO TO 20
      IF (INDN0 .LT. INDNL-1) GO TO 17
      CALL SHIFTD(INDNF,INDNL-2,1, IADJ )
      IADJ(INDNF) = IADJ(INDNL)
      GO TO 20
   17 INDLM3 = INDNL - 3
      IF (INDN0 .GT. INDLM3) GO TO 19
      DO 18 I = INDN0,INDLM3
        ITEMP = IADJ(INDNL)
        CALL SHIFTD(INDNF,INDNL-1,1, IADJ )
        IADJ(INDNF) = ITEMP
   18   CONTINUE
C
C THE LAST SHIFT IS BY 2
C
   19 ITEMP = IADJ(INDNL-1)
      CALL SHIFTD(INDNF,INDLM3,2, IADJ )
      IADJ(INDNF+1) = IADJ(INDNL)
      IADJ(INDNF) = ITEMP
C
C INSERT 0 AS THE LAST NEIGHBOR OF NEWBD
C
   20 IADJ(INDNL) = 0
      IER = 0
      RETURN
C
C ONE OF THE VERTICES IS NOT ON THE BOUNDARY
C
   21 IER = 1
      RETURN
C
C ONE OF THE VERTICES HAS ONLY TWO NONZERO NEIGHBORS.  THE
C   TRIANGULATION WOULD BE DESTROYED BY A DELETION
C
   22 IER = 2
      RETURN
C
C NOUT2 IS NOT THE LAST NONZERO NEIGHBOR OF NOUT1
C
   23 IER = 3
      RETURN
C
C A DELETION WOULD DIVIDE THE MESH INTO TWO REGIONS
C   CONNECTED AT A SINGLE NODE
C
   24 IER = 4
      RETURN
      END
      SUBROUTINE EDGE (IN1,IN2,X,Y, LWK,IWK,IADJ,IEND, IER)
      LOGICAL SWPTST
      INTEGER IN1, IN2, LWK, IWK(2,1), IADJ(1), IEND(1), IER
      REAL    X(1), Y(1)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A TRIANGULATION OF N NODES AND A PAIR OF NODAL
C INDICES IN1 AND IN2, THIS ROUTINE SWAPS ARCS AS NECESSARY
C TO FORCE IN1 AND IN2 TO BE ADJACENT.  ONLY ARCS WHICH
C INTERSECT IN1-IN2 ARE SWAPPED OUT.  IF A THIESSEN TRIANGU-
C LATION IS INPUT, THE RESULTING TRIANGULATION IS AS CLOSE
C AS POSSIBLE TO A THIESSEN TRIANGULATION IN THE SENSE THAT
C ALL ARCS OTHER THAN IN1-IN2 ARE LOCALLY OPTIMAL.
C   A SEQUENCE OF CALLS TO EDGE MAY BE USED TO FORCE THE
C PRESENCE OF A SET OF EDGES DEFINING THE BOUNDARY OF A NON-
C CONVEX REGION.  SUBSEQUENT DELETION OF EDGES OUTSIDE THIS
C REGION (BY SUBROUTINE DELETE) RESULTS IN A NONCONVEX TRI-
C ANGULATION WHICH MAY SERVE AS A FINITE ELEMENT GRID.
C (EDGE SHOULD NOT BE CALLED AFTER A CALL TO DELETE.)  IF,
C ON THE OTHER HAND, INTERPOLATION IS TO BE PERFORMED IN THE
C NONCONVEX REGION, EDGES MUST NOT BE DELETED, BUT IT IS
C STILL ADVANTAGEOUS TO HAVE THE NONCONVEX BOUNDARY PRESENT
C IF IT IS DESIRABLE THAT INTERPOLATED VALUES BE INFLUENCED
C BY THE GEOMETRY.  NOTE THAT SUBROUTINE GETNP WHICH IS USED
C TO SELECT THE NODES ENTERING INTO LOCAL DERIVATIVE ESTI-
C MATES WILL NOT NECESSARILY RETURN CLOSEST NODES IF THE
C TRIANGULATION HAS BEEN RENDERED NONOPTIMAL BY A CALL TO
C EDGE.  HOWEVER, THE EFFECT WILL BE MERELY TO FURTHER EN-
C HANCE THE INFLUENCE OF THE NONCONVEX GEOMETRY ON INTERPO-
C LATED VALUES.
C
C INPUT PARAMETERS - IN1,IN2 - INDICES (OF X AND Y) IN THE
C                              RANGE 1,...,N DEFINING A PAIR
C                              OF NODES TO BE CONNECTED BY
C                              AN ARC.
C
C                        X,Y - N-VECTORS CONTAINING CARTE-
C                              SIAN COORDINATES OF THE
C                              NODES.
C
C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C                        LWK - NUMBER OF COLUMNS RESERVED
C                              FOR IWK.  THIS MUST BE AT
C                              LEAST NI -- THE NUMBER OF
C                              ARCS WHICH INTERSECT IN1-IN2.
C                              (NI IS BOUNDED BY N-3).
C
C                        IWK - INTEGER WORK ARRAY DIMENSION-
C                              ED 2 BY LWK (OR VECTOR OF
C                              LENGTH .GE. 2*LWK).
C
C                  IADJ,IEND - DATA STRUCTURE DEFINING THE
C                              TRIANGULATION.  SEE SUBROU-
C                              TINE TRMESH.
C
C OUTPUT PARAMETERS - LWK - NUMBER OF IWK COLUMNS REQUIRED
C                           IF IER = 0 OR IER = 2.  LWK = 0
C                           IFF IN1 AND IN2 WERE ADJACENT
C                           ON INPUT.
C
C                     IWK - CONTAINS THE INDICES OF THE END-
C                           POINTS OF THE NEW ARCS OTHER
C                           THAN IN1-IN2 UNLESS IER .GT. 0
C                           OR LWK = 0.  NEW ARCS TO THE
C                           LEFT OF IN1->IN2 ARE STORED IN
C                           THE FIRST K-1 COLUMNS (LEFT POR-
C                           TION OF IWK), COLUMN K CONTAINS
C                           ZEROS, AND NEW ARCS TO THE RIGHT
C                           OF IN1->IN2 OCCUPY COLUMNS K+1,
C                           ...,LWK.  (K CAN BE DETERMINED
C                           BY SEARCHING IWK FOR THE ZEROS.)
C
C               IADJ,IEND - UPDATED IF NECESSARY TO REFLECT
C                           THE PRESENCE OF AN ARC CONNECT-
C                           ING IN1 AND IN2, UNALTERED IF
C                           IER .NE. 0.
C
C                     IER - ERROR INDICATOR
C                           IER = 0 IF NO ERRORS WERE EN-
C                                   COUNTERED.
C                           IER = 1 IF IN1 .LT. 1, IN2 .LT.
C                                   1, IN1 = IN2, OR LWK
C                                   .LT. 0 ON INPUT.
C                           IER = 2 IF MORE SPACE IS REQUIR-
C                                   ED IN IWK.  SEE LWK.
C                           IER = 3 IF IN1 AND IN2 COULD NOT
C                                   BE CONNECTED DUE TO AN
C                                   INVALID DATA STRUCTURE.
C
C MODULES REFERENCED BY EDGE - SWAP, INDEX, SHIFTD, SWPTST
C
C***********************************************************
C
      INTEGER N1, N2, IWEND, IWL, INDF, INDX, N1LST, NL, NR,
     .        NEXT, IWF, LFT, N0, IWC, IWCP1, IWCM1, I, IO1,
     .        IO2, INDL
      REAL    X1, Y1, X2, Y2, X0, Y0
      LOGICAL SWP, LEFT
C
C LOCAL PARAMETERS -
C
C N1,N2 =   LOCAL COPIES OF IN1 AND IN2 OR NODES OPPOSITE AN
C             ARC IO1-IO2 TO BE TESTED FOR A SWAP IN THE
C             OPTIMIZATION LOOPS
C IWEND =   INPUT OR OUTPUT VALUE OF LWK
C IWL =     IWK (COLUMN) INDEX OF THE LAST (RIGHTMOST) ARC
C             WHICH INTERSECTS IN1->IN2
C INDF =    IADJ INDEX OF THE FIRST NEIGHBOR OF IN1 OR IO1
C INDX =    IADJ INDEX OF A NEIGHBOR OF IN1, NL, OR IO1
C N1LST =   LAST NEIGHBOR OF IN1
C NL,NR =   ENDPOINTS OF AN ARC WHICH INTERSECTS IN1-IN2
C             WITH NL LEFT IN1->IN2
C NEXT =    NODE OPPOSITE NL->NR
C IWF =     IWK (COLUMN) INDEX OF THE FIRST (LEFTMOST) ARC
C             WHICH INTERSECTS IN1->IN2
C LFT =     FLAG USED TO DETERMINE IF A SWAP RESULTS IN THE
C             NEW ARC INTERSECTING IN1-IN2 -- LFT = 0 IFF
C             N0 = IN1, LFT = -1 IMPLIES N0 LEFT IN1->IN2,
C             AND LFT = 1 IMPLIES N0 LEFT IN2->IN1
C N0 =      NODE OPPOSITE NR->NL
C IWC =     IWK INDEX BETWEEN IWF AND IWL -- NL->NR IS
C             STORED IN IWK(1,IWC)->IWK(2,IWC)
C IWCP1 =   IWC + 1
C IWCM1 =   IWC - 1
C I =       DO-LOOP INDEX AND COLUMN INDEX FOR IWK
C IO1,IO2 = ENDPOINTS OF AN ARC TO BE TESTED FOR A SWAP IN
C             THE OPTIMIZATION LOOPS
C INDL =    IADJ INDEX OF THE LAST NEIGHBOR OF IO1
C X1,Y1 =   COORDINATES OF IN1
C X2,Y2 =   COORDINATES OF IN2
C X0,Y0 =   COORDINATES OF N0
C SWP =     FLAG SET TO .TRUE. IFF A SWAP OCCURS IN AN OPTI-
C             MIZATION LOOP
C LEFT =    STATEMENT FUNCTION WHICH RETURNS THE VALUE
C             .TRUE. IFF (XP,YP) IS ON OR TO THE LEFT OF THE
C             VECTOR (XA,YA)->(XB,YB)
C
      LEFT(XA,YA,XB,YB,XP,YP) = (XB-XA)*(YP-YA) .GE.
     .                          (XP-XA)*(YB-YA)
C
C STORE IN1, IN2, AND LWK IN LOCAL VARIABLES AND CHECK FOR
C   ERRORS.
C
      N1 = IN1
      N2 = IN2
      IWEND = LWK
      IF (N1 .LT. 1  .OR.  N2 .LT. 1  .OR.  N1 .EQ. N2  .OR.
     .    IWEND .LT. 0) GO TO 35
C
C STORE THE COORDINATES OF N1 AND N2 AND INITIALIZE IWL.
C
      X1 = X(N1)
      Y1 = Y(N1)
      X2 = X(N2)
      Y2 = Y(N2)
      IWL = 0
C
C SET NR AND NL TO ADJACENT NEIGHBORS OF N1 SUCH THAT
C   NR LEFT N2->N1 AND NL LEFT N1->N2.
C
C   SET INDF AND INDX TO THE INDICES OF THE FIRST AND LAST
C     NEIGHBORS OF N1 AND SET N1LST TO THE LAST NEIGHBOR.
C
      INDF = 1
      IF (N1 .GT. 1) INDF = IEND(N1-1) + 1
      INDX = IEND(N1)
      N1LST = IADJ(INDX)
      IF (N1LST .EQ. 0) INDX = INDX - 1
      IF (N1LST .EQ. 0) GO TO 2
C
C   N1 IS AN INTERIOR NODE.  LOOP THROUGH THE NEIGHBORS NL
C     IN REVERSE ORDER UNTIL NL LEFT N1->N2.
C
      NL = N1LST
    1 IF ( LEFT(X1,Y1,X2,Y2,X(NL),Y(NL)) ) GO TO 2
      INDX = INDX - 1
      NL = IADJ(INDX)
      IF (INDX .GT. INDF) GO TO 1
C
C   NL IS THE FIRST NEIGHBOR OF N1.  SET NR TO THE LAST
C     NEIGHBOR AND TEST FOR AN ARC N1-N2.
C
      NR = N1LST
      IF (NL .EQ. N2) GO TO 34
      GO TO 4
C
C   NL = IADJ(INDX) LEFT N1->N2 AND INDX .GT. INDF.  SET
C     NR TO THE PRECEDING NEIGHBOR OF N1.
C
    2 INDX = INDX - 1
      NR = IADJ(INDX)
      IF ( LEFT(X2,Y2,X1,Y1,X(NR),Y(NR)) ) GO TO 3
      IF (INDX .GT. INDF) GO TO 2
C
C   SET NL AND NR TO THE FIRST AND LAST NEIGHBORS OF N1 AND
C     TEST FOR AN INVALID DATA STRUCTURE (N1 CANNOT BE A
C     BOUNDARY NODE AND CANNOT BE ADJACENT TO N2).
C
      NL = NR
      NR = N1LST
      IF (NR .EQ. 0  .OR.  NR .EQ. N2) GO TO 37
      GO TO 4
C
C   SET NL TO THE NEIGHBOR FOLLOWING NR AND TEST FOR AN ARC
C     N1-N2.
C
    3 NL = IADJ(INDX+1)
      IF (NL .EQ. N2  .OR.  NR .EQ. N2) GO TO 34
C
C STORE THE ORDERED SEQUENCE OF INTERSECTING EDGES NL->NR IN
C   IWK(1,IWL)->IWK(2,IWL).
C
    4 IWL = IWL + 1
      IF (IWL .LE. IWEND) IWK(1,IWL) = NL
      IF (IWL .LE. IWEND) IWK(2,IWL) = NR
C
C   SET NEXT TO THE NEIGHBOR OF NL WHICH FOLLOWS NR.
C
      INDX = IEND(NL)
      IF (IADJ(INDX) .NE. NR) GO TO 5
C
C   NR IS THE LAST NEIGHBOR OF NL.  SET NEXT TO THE FIRST
C     NEIGHBOR.
C
      INDX = 0
      IF (NL .NE. 1) INDX = IEND(NL-1)
      GO TO 6
C
C   NR IS NOT THE LAST NEIGHBOR OF NL.  LOOP THROUGH THE
C     NEIGHBORS IN REVERSE ORDER.
C
    5 INDX = INDX - 1
      IF (IADJ(INDX) .NE. NR) GO TO 5
C
C   STORE NEXT, TEST FOR AN INVALID TRIANGULATION (NL->NR
C     CANNOT BE A BOUNDARY EDGE), AND TEST FOR TERMINATION
C     OF THE LOOP.
C
    6 NEXT = IADJ(INDX+1)
      IF (NEXT .EQ. 0) GO TO 37
      IF (NEXT .EQ. N2) GO TO 8
C
C   SET NL OR NR TO NEXT.
C
      IF ( LEFT(X1,Y1,X2,Y2,X(NEXT),Y(NEXT)) ) GO TO 7
      NR = NEXT
      GO TO 4
    7 NL = NEXT
      GO TO 4
C
C IWL IS THE NUMBER OF ARCS WHICH INTERSECT N1-N2.  STORE
C   LWK AND TEST FOR SUFFICIENT SPACE.
C
    8 LWK = IWL
      IF (IWL .GT. IWEND) GO TO 36
      IWEND = IWL
C
C INITIALIZE FOR EDGE SWAPPING LOOP -- ALL POSSIBLE SWAPS
C   ARE APPLIED (EVEN IF THE NEW ARC AGAIN INTERSECTS
C   N1-N2), ARCS TO THE LEFT OF N1->N2 ARE STORED IN THE
C   LEFT PORTION OF IWK, AND ARCS TO THE RIGHT ARE STORED IN
C   THE RIGHT PORTION.  IWF AND IWL INDEX THE FIRST AND LAST
C   INTERSECTING ARCS.
C
      IER = 0
      IWF = 1
C
C TOP OF LOOP -- SET N0 TO N1 AND NL->NR TO THE FIRST EDGE.
C   IWC POINTS TO THE ARC CURRENTLY BEING PROCESSED.  LFT
C   .LE. 0 IFF N0 LEFT N1->N2.
C
    9 LFT = 0
      N0 = N1
      X0 = X1
      Y0 = Y1
      NL = IWK(1,IWF)
      NR = IWK(2,IWF)
      IWC = IWF
C
C   SET NEXT TO THE NODE OPPOSITE NL->NR UNLESS IWC IS THE
C     LAST ARC.
C
   10 IF (IWC .EQ. IWL) GO TO 21
      IWCP1 = IWC + 1
      NEXT = IWK(1,IWCP1)
      IF (NEXT .NE. NL) GO TO 15
      NEXT = IWK(2,IWCP1)
C
C   NEXT RIGHT N1->N2 AND IWC .LT. IWL.  TEST FOR A POSSIBLE
C     SWAP.
C
      IF ( .NOT. LEFT(X0,Y0,X(NR),Y(NR),X(NEXT),Y(NEXT)) )
     .   GO TO 13
      IF (LFT .GE. 0) GO TO 11
      IF ( .NOT. LEFT(X(NL),Y(NL),X0,Y0,X(NEXT),Y(NEXT)) )
     .   GO TO 13
C
C   REPLACE NL->NR WITH N0->NEXT.
C
      CALL SWAP(NEXT,N0,NL,NR, IADJ,IEND )
      IWK(1,IWC) = N0
      IWK(2,IWC) = NEXT
      GO TO 14
C
C   SWAP NL-NR FOR N0-NEXT, SHIFT COLUMNS IWC+1,...,IWL TO
C     THE LEFT, AND STORE N0-NEXT IN THE RIGHT PORTION OF
C     IWK.
C
   11 CALL SWAP(NEXT,N0,NL,NR, IADJ,IEND )
      DO 12 I = IWCP1,IWL
        IWK(1,I-1) = IWK(1,I)
   12   IWK(2,I-1) = IWK(2,I)
      IWK(1,IWL) = N0
      IWK(2,IWL) = NEXT
      IWL = IWL - 1
      NR = NEXT
      GO TO 10
C
C   A SWAP IS NOT POSSIBLE.  SET N0 TO NR.
C
   13 N0 = NR
      X0 = X(N0)
      Y0 = Y(N0)
      LFT = 1
C
C   ADVANCE TO THE NEXT ARC.
C
   14 NR = NEXT
      IWC = IWC + 1
      GO TO 10
C
C   NEXT LEFT N1->N2, NEXT .NE. N2, AND IWC .LT. IWL.
C     TEST FOR A POSSIBLE SWAP.
C
   15 IF ( .NOT. LEFT(X(NL),Y(NL),X0,Y0,X(NEXT),Y(NEXT)) )
     .   GO TO 19
      IF (LFT .LE. 0) GO TO 16
      IF ( .NOT. LEFT(X0,Y0,X(NR),Y(NR),X(NEXT),Y(NEXT)) )
     .   GO TO 19
C
C   REPLACE NL->NR WITH NEXT->N0.
C
      CALL SWAP(NEXT,N0,NL,NR, IADJ,IEND )
      IWK(1,IWC) = NEXT
      IWK(2,IWC) = N0
      GO TO 20
C
C   SWAP NL-NR FOR N0-NEXT, SHIFT COLUMNS IWF,...,IWC-1 TO
C     THE RIGHT, AND STORE N0-NEXT IN THE LEFT PORTION OF
C     IWK.
C
   16 CALL SWAP(NEXT,N0,NL,NR, IADJ,IEND )
      I = IWC
   17 IF (I .EQ. IWF) GO TO 18
      IWK(1,I) = IWK(1,I-1)
      IWK(2,I) = IWK(2,I-1)
      I = I - 1
      GO TO 17
   18 IWK(1,IWF) = N0
      IWK(2,IWF) = NEXT
      IWF = IWF + 1
      GO TO 20
C
C   A SWAP IS NOT POSSIBLE.  SET N0 TO NL.
C
   19 N0 = NL
      X0 = X(N0)
      Y0 = Y(N0)
      LFT = -1
C
C   ADVANCE TO THE NEXT ARC.
C
   20 NL = NEXT
      IWC = IWC + 1
      GO TO 10
C
C   N2 IS OPPOSITE NL->NR (IWC = IWL).
C
   21 IF (N0 .EQ. N1) GO TO 24
      IF (LFT .LT. 0) GO TO 22
C
C   N0 RIGHT N1->N2.  TEST FOR A POSSIBLE SWAP.
C
      IF ( .NOT. LEFT(X0,Y0,X(NR),Y(NR),X2,Y2) ) GO TO 9
C
C   SWAP NL-NR FOR N0-N2 AND STORE N0-N2 IN THE RIGHT
C     PORTION OF IWK.
C
      CALL SWAP(N2,N0,NL,NR, IADJ,IEND )
      IWK(1,IWL) = N0
      IWK(2,IWL) = N2
      IWL = IWL - 1
      GO TO 9
C
C   N0 LEFT N1->N2.  TEST FOR A POSSIBLE SWAP.
C
   22 IF ( .NOT. LEFT(X(NL),Y(NL),X0,Y0,X2,Y2) ) GO TO 9
C
C   SWAP NL-NR FOR N0-N2, SHIFT COLUMNS IWF,...,IWL-1 TO THE
C     RIGHT, AND STORE N0-N2 IN THE LEFT PORTION OF IWK.
C
      CALL SWAP(N2,N0,NL,NR, IADJ,IEND )
      I = IWL
   23 IWK(1,I) = IWK(1,I-1)
      IWK(2,I) = IWK(2,I-1)
      I = I - 1
      IF (I .GT. IWF) GO TO 23
      IWK(1,IWF) = N0
      IWK(2,IWF) = N2
      IWF = IWF + 1
      GO TO 9
C
C IWF = IWC = IWL.  SWAP OUT THE LAST ARC FOR N1-N2 AND
C   STORE ZEROS IN IWK.
C
   24 CALL SWAP(N2,N1,NL,NR, IADJ,IEND )
      IWK(1,IWC) = 0
      IWK(2,IWC) = 0
      IF (IWC .EQ. 1) GO TO 29
C
C OPTIMIZATION LOOPS -- OPTIMIZE THE SET OF NEW ARCS TO THE
C   LEFT OF IN1->IN2.  THE LOOP IS REPEATED UNTIL NO SWAPS
C   ARE PERFORMED.
C
      IWCM1 = IWC - 1
   25 SWP = .FALSE.
      DO 28 I = 1,IWCM1
        IO1 = IWK(1,I)
        IO2 = IWK(2,I)
C
C   SET N1 TO THE NEIGHBOR OF IO1 WHICH FOLLOWS IO2 AND SET
C     N2 TO THE NEIGHBOR OF IO1 WHICH PRECEDES IO2.
C
        INDF = 1
        IF (IO1 .GT. 1) INDF = IEND(IO1-1) + 1
        INDL = IEND(IO1)
        INDX = INDL
        IF (IADJ(INDX) .NE. IO2) GO TO 26
C
C   IO2 IS THE LAST NEIGHBOR OF IO1.
C
        N1 = IADJ(INDF)
        N2 = IADJ(INDX-1)
        GO TO 27
C
C   IO2 IS NOT THE LAST NEIGHBOR OF IO1.  LOOP THROUGH THE
C     NEIGHBORS IN REVERSE ORDER.
C
   26   INDX = INDX - 1
        IF (IADJ(INDX) .NE. IO2) GO TO 26
        N1 = IADJ(INDX+1)
        IF (INDX .NE. INDF) N2 = IADJ(INDX-1)
        IF (INDX .EQ. INDF) N2 = IADJ(INDL)
C
C   TEST IO1-IO2 FOR A SWAP.
C
   27   IF ( .NOT. SWPTST(N1,N2,IO1,IO2,X,Y) ) GO TO 28
        SWP = .TRUE.
        CALL SWAP(N1,N2,IO1,IO2, IADJ,IEND )
        IWK(1,I) = N1
        IWK(2,I) = N2
   28   CONTINUE
      IF (SWP) GO TO 25
C
C TEST FOR TERMINATION.
C
   29 IF (IWC .EQ. IWEND) RETURN
      IWCP1 = IWC + 1
C
C OPTIMIZE THE SET OF NEW ARCS TO THE RIGHT OF IN1->IN2.
C
   30 SWP = .FALSE.
      DO 33 I = IWCP1,IWEND
        IO1 = IWK(1,I)
        IO2 = IWK(2,I)
C
C   SET N1 AND N2 TO THE NODES OPPOSITE IO1->IO2 AND
C     IO2->IO1, RESPECTIVELY.
C
        INDF = 1
        IF (IO1 .GT. 1) INDF = IEND(IO1-1) + 1
        INDL = IEND(IO1)
        INDX = INDL
        IF (IADJ(INDX) .NE. IO2) GO TO 31
C
        N1 = IADJ(INDF)
        N2 = IADJ(INDX-1)
        GO TO 32
C
   31   INDX = INDX - 1
        IF (IADJ(INDX) .NE. IO2) GO TO 31
        N1 = IADJ(INDX+1)
        IF (INDX .NE. INDF) N2 = IADJ(INDX-1)
        IF (INDX .EQ. INDF) N2 = IADJ(INDL)
C
   32   IF ( .NOT. SWPTST(N1,N2,IO1,IO2,X,Y) ) GO TO 33
        SWP = .TRUE.
        CALL SWAP(N1,N2,IO1,IO2, IADJ,IEND )
        IWK(1,I) = N1
        IWK(2,I) = N2
   33   CONTINUE
      IF (SWP) GO TO 30
      RETURN
C
C IN1 AND IN2 WERE ADJACENT ON INPUT.
C
   34 IER = 0
      LWK = 0
      RETURN
C
C PARAMETER OUT OF RANGE
C
   35 IER = 1
      RETURN
C
C INSUFFICIENT SPACE IN IWK
C
   36 IER = 2
      RETURN
C
C INVALID TRIANGULATION DATA STRUCTURE
C
   37 IER = 3
      RETURN
      END
      SUBROUTINE GETNP (X,Y,IADJ,IEND,L, NPTS, DS,IER)
      INTEGER IADJ(1), IEND(1), L, NPTS(L), IER
      REAL    X(1), Y(1), DS
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A THIESSEN TRIANGULATION OF N NODES AND AN ARRAY
C NPTS CONTAINING THE INDICES OF L-1 NODES ORDERED BY
C EUCLIDEAN DISTANCE FROM NPTS(1), THIS SUBROUTINE SETS
C NPTS(L) TO THE INDEX OF THE NEXT NODE IN THE SEQUENCE --
C THE NODE, OTHER THAN NPTS(1),...,NPTS(L-1), WHICH IS
C CLOSEST TO NPTS(1).  THUS, THE ORDERED SEQUENCE OF K
C CLOSEST NODES TO N1 (INCLUDING N1) MAY BE DETERMINED BY
C K-1 CALLS TO GETNP WITH NPTS(1) = N1 AND L = 2,3,...,K
C FOR K .GE. 2.
C   THE ALGORITHM USES THE FACT THAT, IN A THIESSEN TRIAN-
C GULATION, THE K-TH CLOSEST NODE TO A GIVEN NODE N1 IS A
C NEIGHBOR OF ONE OF THE K-1 CLOSEST NODES TO N1.
C
C INPUT PARAMETERS - X,Y - VECTORS OF LENGTH N CONTAINING
C                          THE CARTESIAN COORDINATES OF THE
C                          NODES.
C
C                   IADJ - SET OF ADJACENCY LISTS OF NODES
C                          IN THE TRIANGULATION.
C
C                   IEND - POINTERS TO THE ENDS OF ADJACENCY
C                          LISTS FOR EACH NODE IN THE TRI-
C                          ANGULATION.
C
C                      L - NUMBER OF NODES IN THE SEQUENCE
C                          ON OUTPUT.  2 .LE. L .LE. N.
C
C                   NPTS - ARRAY OF LENGTH .GE. L CONTAIN-
C                          ING THE INDICES OF THE L-1 CLOS-
C                          EST NODES TO NPTS(1) IN THE FIRST
C                          L-1 LOCATIONS.
C
C IADJ AND IEND MAY BE CREATED BY SUBROUTINE TRMESH.
C
C INPUT PARAMETERS OTHER THAN NPTS ARE NOT ALTERED BY THIS
C   ROUTINE.
C
C OUTPUT PARAMETERS - NPTS - UPDATED WITH THE INDEX OF THE
C                            L-TH CLOSEST NODE TO NPTS(1) IN
C                            POSITION L UNLESS IER = 1.
C
C                       DS - SQUARED EUCLIDEAN DISTANCE BE-
C                            TWEEN NPTS(1) AND NPTS(L)
C                            UNLESS IER = 1.
C
C                      IER - ERROR INDICATOR
C                            IER = 0 IF NO ERRORS WERE EN-
C                                    COUNTERED.
C                            IER = 1 IF L IS OUT OF RANGE.
C
C MODULES REFERENCED BY GETNP - NONE
C
C INTRINSIC FUNCTION CALLED BY GETNP - IABS
C
C***********************************************************
C
      INTEGER LM1, N1, I, NI, NP, INDF, INDL, INDX, NB
      REAL    X1, Y1, DNP, DNB
C
C LOCAL PARAMETERS -
C
C LM1 =     L - 1
C N1 =      NPTS(1)
C I =       NPTS INDEX AND DO-LOOP INDEX
C NI =      NPTS(I)
C NP =      CANDIDATE FOR NPTS(L)
C INDF =    IADJ INDEX OF THE FIRST NEIGHBOR OF NI
C INDL =    IADJ INDEX OF THE LAST NEIGHBOR OF NI
C INDX =    IADJ INDEX IN THE RANGE INDF,...,INDL
C NB =      NEIGHBOR OF NI AND CANDIDATE FOR NP
C X1,Y1 =   COORDINATES OF N1
C DNP,DNB = SQUARED DISTANCES FROM N1 TO NP AND NB,
C             RESPECTIVELY
C
      LM1 = L - 1
      IF (LM1 .LT. 1) GO TO 4
      IER = 0
      N1 = NPTS(1)
      X1 = X(N1)
      Y1 = Y(N1)
C
C MARK THE ELEMENTS OF NPTS
C
      DO 1 I = 1,LM1
        NI = NPTS(I)
        IEND(NI) = -IEND(NI)
    1   CONTINUE
C
C CANDIDATES FOR NP = NPTS(L) ARE THE UNMARKED NEIGHBORS
C   OF NODES IN NPTS.  NP=0 IS A FLAG TO SET NP TO THE
C   FIRST CANDIDATE ENCOUNTERED.
C
      NP = 0
      DNP = 0.
C
C LOOP ON NODES NI IN NPTS
C
      DO 2 I = 1,LM1
        NI = NPTS(I)
        INDF = 1
        IF (NI .GT. 1) INDF = IABS(IEND(NI-1)) + 1
        INDL = -IEND(NI)
C
C LOOP ON NEIGHBORS NB OF NI
C
        DO 2 INDX = INDF,INDL
          NB = IADJ(INDX)
          IF (NB .EQ. 0  .OR.  IEND(NB) .LT. 0) GO TO 2
C
C NB IS AN UNMARKED NEIGHBOR OF NI.  REPLACE NP IF NB IS
C   CLOSER TO N1 OR IS THE FIRST CANDIDATE ENCOUNTERED.
C
          DNB = (X(NB)-X1)**2 + (Y(NB)-Y1)**2
          IF (NP .NE. 0  .AND.  DNB .GE. DNP) GO TO 2
          NP = NB
          DNP = DNB
    2     CONTINUE
      NPTS(L) = NP
      DS = DNP
C
C UNMARK THE ELEMENTS OF NPTS
C
      DO 3 I = 1,LM1
        NI = NPTS(I)
        IEND(NI) = -IEND(NI)
    3   CONTINUE
      RETURN
C
C L IS OUT OF RANGE
C
    4 IER = 1
      RETURN
      END
      INTEGER FUNCTION INDEX (NVERTX,NABOR,IADJ,IEND)
      INTEGER NVERTX, NABOR, IADJ(1), IEND(1)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS FUNCTION RETURNS THE INDEX OF NABOR IN THE
C ADJACENCY LIST FOR NVERTX.
C
C INPUT PARAMETERS - NVERTX - NODE WHOSE ADJACENCY LIST IS
C                             TO BE SEARCHED.
C
C                     NABOR - NODE WHOSE INDEX IS TO BE
C                             RETURNED.  NABOR MUST BE
C                             CONNECTED TO NVERTX.
C
C                      IADJ - SET OF ADJACENCY LISTS.
C
C                      IEND - POINTERS TO THE ENDS OF
C                             ADJACENCY LISTS IN IADJ.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION.
C
C OUTPUT PARAMETER -  INDEX - IADJ(INDEX) = NABOR.
C
C MODULES REFERENCED BY INDEX - NONE
C
C***********************************************************
C
      INTEGER NB, INDX
C
C LOCAL PARAMETERS -
C
C NB =   LOCAL COPY OF NABOR
C INDX = INDEX FOR IADJ
C
      NB = NABOR
C
C INITIALIZATION
C
      INDX = IEND(NVERTX) + 1
C
C SEARCH THE LIST OF NVERTX NEIGHBORS FOR NB
C
    1 INDX = INDX - 1
      IF (IADJ(INDX) .NE. NB) GO TO 1
C
      INDEX = INDX
      RETURN
      END
      SUBROUTINE INTADD (KK,I1,I2,I3, IADJ,IEND )
      INTEGER KK, I1, I2, I3, IADJ(1), IEND(KK)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE ADDS AN INTERIOR NODE TO A TRIANGULATION
C OF A SET OF KK-1 POINTS IN THE PLANE.  IADJ AND IEND ARE
C UPDATED WITH THE INSERTION OF NODE KK IN THE TRIANGLE
C WHOSE VERTICES ARE I1, I2, AND I3.
C
C INPUT PARAMETERS -        KK - INDEX OF NODE TO BE
C                                INSERTED.  KK .GE. 4.
C
C                     I1,I2,I3 - INDICES OF THE VERTICES OF
C                                A TRIANGLE CONTAINING NODE
C                                KK -- IN COUNTERCLOCKWISE
C                                ORDER.
C
C                         IADJ - SET OF ADJACENCY LISTS
C                                OF NODES IN THE MESH.
C
C                         IEND - POINTERS TO THE ENDS OF
C                                ADJACENCY LISTS IN IADJ FOR
C                                EACH NODE IN THE MESH.
C
C   IADJ AND IEND MAY BE CREATED BY TRMESH AND MUST CONTAIN
C THE VERTICES I1, I2, AND I3.  I1,I2,I3 MAY BE DETERMINED
C BY TRFIND.
C
C KK, I1, I2, AND I3 ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ADDITION
C                                 OF NODE KK AS THE LAST
C                                 ENTRY.  NODE KK WILL BE
C                                 CONNECTED TO NODES I1, I2,
C                                 AND I3.  NO OPTIMIZATION
C                                 OF THE MESH IS PERFORMED.
C
C MODULE REFERENCED BY INTADD - SHIFTD
C
C INTRINSIC FUNCTION CALLED BY INTADD - MOD
C
C***********************************************************
C
      INTEGER K, KM1, N(3), NFT(3), IP1, IP2, IP3, INDX, NF,
     .        NL, N1, N2, IMIN, IMAX, I, ITEMP
C
C LOCAL PARAMETERS -
C
C K =           LOCAL COPY OF KK
C KM1 =         K - 1
C N =           VECTOR CONTAINING I1, I2, I3
C NFT =         POINTERS TO THE TOPS OF THE 3 SETS OF IADJ
C                 ELEMENTS TO BE SHIFTED DOWNWARD
C IP1,IP2,IP3 = PERMUTATION INDICES FOR N AND NFT
C INDX =        INDEX FOR IADJ AND N
C NF,NL =       INDICES OF FIRST AND LAST ENTRIES IN IADJ
C                 TO BE SHIFTED DOWN
C N1,N2 =       FIRST 2 VERTICES OF A NEW TRIANGLE --
C                 (N1,N2,KK)
C IMIN,IMAX =   BOUNDS ON DO-LOOP INDEX -- FIRST AND LAST
C                 ELEMENTS OF IEND TO BE INCREMENTED
C I =           DO-LOOP INDEX
C ITEMP =       TEMPORARY STORAGE LOCATION
C
      K = KK
C
C INITIALIZATION
C
      N(1) = I1
      N(2) = I2
      N(3) = I3
C
C SET UP NFT
C
      DO 2 I = 1,3
        N1 = N(I)
        INDX = MOD(I,3) + 1
        N2 = N(INDX)
        INDX = IEND(N1) + 1
C
C FIND THE INDEX OF N2 AS A NEIGHBOR OF N1
C
    1   INDX = INDX - 1
        IF (IADJ(INDX) .NE. N2) GO TO 1
        NFT(I) = INDX + 1
    2   CONTINUE
C
C ORDER THE VERTICES BY DECREASING MAGNITUDE.
C   N(IP(I+1)) PRECEDES N(IP(I)) IN IEND FOR
C   I = 1,2.
C
      IP1 = 1
      IP2 = 2
      IP3 = 3
      IF ( N(2) .LE. N(1) ) GO TO 3
      IP1 = 2
      IP2 = 1
    3 IF ( N(3) .LE. N(IP1) ) GO TO 4
      IP3 = IP1
      IP1 = 3
    4 IF ( N(IP3) .LE. N(IP2) )  GO TO 5
      ITEMP = IP2
      IP2 = IP3
      IP3 = ITEMP
C
C ADD NODE K TO THE ADJACENCY LISTS OF EACH VERTEX AND
C   UPDATE IEND.  FOR EACH VERTEX, A SET OF IADJ ELEMENTS
C   IS SHIFTED DOWNWARD AND K IS INSERTED.  SHIFTING STARTS
C   AT THE END OF THE ARRAY.
C
    5 KM1 = K - 1
      NL = IEND(KM1)
      NF = NFT(IP1)
      IF (NF .LE. NL) CALL SHIFTD(NF,NL,3, IADJ )
      IADJ(NF+2) = K
      IMIN = N(IP1)
      IMAX = KM1
      DO 6 I = IMIN,IMAX
        IEND(I) = IEND(I) + 3
    6   CONTINUE
C
      NL = NF - 1
      NF = NFT(IP2)
      CALL SHIFTD(NF,NL,2, IADJ )
      IADJ(NF+1) = K
      IMAX = IMIN - 1
      IMIN = N(IP2)
      DO 7 I = IMIN,IMAX
        IEND(I) = IEND(I) + 2
    7   CONTINUE
C
      NL = NF - 1
      NF = NFT(IP3)
      CALL SHIFTD(NF,NL,1, IADJ )
      IADJ(NF) = K
      IMAX = IMIN - 1
      IMIN = N(IP3)
      DO 8 I = IMIN,IMAX
        IEND(I) = IEND(I) + 1
    8   CONTINUE
C
C ADD NODE K TO IEND AND ITS NEIGHBORS TO IADJ
C
      INDX = IEND(KM1)
      IEND(K) = INDX + 3
      DO 9 I = 1,3
        INDX = INDX + 1
        IADJ(INDX) = N(I)
    9   CONTINUE
      RETURN
      END
      SUBROUTINE PERMUT (NN,IP, A )
      INTEGER NN, IP(NN)
      REAL    A(NN)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE APPLIES A SET OF PERMUTATIONS TO A VECTOR.
C
C INPUT PARAMETERS - NN - LENGTH OF A AND IP.
C
C                    IP - VECTOR CONTAINING THE SEQUENCE OF
C                         INTEGERS 1,...,NN PERMUTED IN THE
C                         SAME FASHION THAT A IS TO BE PER-
C                         MUTED.
C
C                     A - VECTOR TO BE PERMUTED.
C
C NN AND IP ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETERS - A - REORDERED VECTOR REFLECTING THE
C                         PERMUTATIONS DEFINED BY IP.
C
C MODULES REFERENCED BY PERMUT - NONE
C
C***********************************************************
C
      INTEGER N, K, J, IPJ
      REAL    TEMP
C
C LOCAL PARAMETERS -
C
C N =    LOCAL COPY OF NN
C K =    INDEX FOR IP AND FOR THE FIRST ELEMENT OF A IN A
C          PERMUTATION
C J =    INDEX FOR IP AND A, J .GE. K
C IPJ =  IP(J)
C TEMP = TEMPORARY STORAGE FOR A(K)
C
      N = NN
      IF (N .LT. 2) RETURN
      K = 1
C
C LOOP ON PERMUTATIONS
C
    1 J = K
      TEMP = A(K)
C
C APPLY PERMUTATION TO A.  IP(J) IS MARKED (MADE NEGATIVE)
C   AS BEING INCLUDED IN THE PERMUTATION.
C
    2 IPJ = IP(J)
      IP(J) = -IPJ
      IF (IPJ .EQ. K) GO TO 3
      A(J) = A(IPJ)
      J = IPJ
      GO TO 2
    3 A(J) = TEMP
C
C SEARCH FOR AN UNMARKED ELEMENT OF IP
C
    4 K = K + 1
      IF (K .GT. N) GO TO 5
      IF (IP(K) .GT. 0) GO TO 1
      GO TO 4
C
C ALL PERMUTATIONS HAVE BEEN APPLIED.  UNMARK IP.
C
    5 DO 6 K = 1,N
        IP(K) = -IP(K)
    6   CONTINUE
      RETURN
      END
      SUBROUTINE QSORT (N,X, IND)
      INTEGER N, IND(N)
      REAL    X(N)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS SUBROUTINE USES AN ORDER N*LOG(N) QUICK SORT TO
C SORT THE REAL ARRAY X INTO INCREASING ORDER.  THE ALGOR-
C ITHM IS AS FOLLOWS.  IND IS INITIALIZED TO THE ORDERED
C SEQUENCE OF INDICES 1,...,N, AND ALL INTERCHANGES ARE
C APPLIED TO IND.  X IS DIVIDED INTO TWO PORTIONS BY PICKING
C A CENTRAL ELEMENT T.  THE FIRST AND LAST ELEMENTS ARE COM-
C PARED WITH T, AND INTERCHANGES ARE APPLIED AS NECESSARY SO
C THAT THE THREE VALUES ARE IN ASCENDING ORDER.  INTER-
C CHANGES ARE THEN APPLIED SO THAT ALL ELEMENTS GREATER THAN
C T ARE IN THE UPPER PORTION OF THE ARRAY AND ALL ELEMENTS
C LESS THAN T ARE IN THE LOWER PORTION.  THE UPPER AND LOWER
C INDICES OF ONE OF THE PORTIONS ARE SAVED IN LOCAL ARRAYS,
C AND THE PROCESS IS REPEATED ITERATIVELY ON THE OTHER
C PORTION.  WHEN A PORTION IS COMPLETELY SORTED, THE PROCESS
C BEGINS AGAIN BY RETRIEVING THE INDICES BOUNDING ANOTHER
C UNSORTED PORTION.
C
C INPUT PARAMETERS -   N - LENGTH OF THE ARRAY X.
C
C                      X - VECTOR OF LENGTH N TO BE SORTED.
C
C                    IND - VECTOR OF LENGTH .GE. N.
C
C N AND X ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETER - IND - SEQUENCE OF INDICES 1,...,N
C                          PERMUTED IN THE SAME FASHION AS X
C                          WOULD BE.  THUS, THE ORDERING ON
C                          X IS DEFINED BY Y(I) = X(IND(I)).
C
C MODULES REFERENCED BY QSORT - NONE
C
C INTRINSIC FUNCTIONS CALLED BY QSORT - IFIX, FLOAT
C
C***********************************************************
C
C NOTE -- IU AND IL MUST BE DIMENSIONED .GE. LOG(N) WHERE
C         LOG HAS BASE 2.
C
C***********************************************************
C
      INTEGER IU(21), IL(21)
      INTEGER M, I, J, K, L, IJ, IT, ITT, INDX
      REAL    R, T
C
C LOCAL PARAMETERS -
C
C IU,IL =  TEMPORARY STORAGE FOR THE UPPER AND LOWER
C            INDICES OF PORTIONS OF THE ARRAY X
C M =      INDEX FOR IU AND IL
C I,J =    LOWER AND UPPER INDICES OF A PORTION OF X
C K,L =    INDICES IN THE RANGE I,...,J
C IJ =     RANDOMLY CHOSEN INDEX BETWEEN I AND J
C IT,ITT = TEMPORARY STORAGE FOR INTERCHANGES IN IND
C INDX =   TEMPORARY INDEX FOR X
C R =      PSEUDO RANDOM NUMBER FOR GENERATING IJ
C T =      CENTRAL ELEMENT OF X
C
      IF (N .LE. 0) RETURN
C
C INITIALIZE IND, M, I, J, AND R
C
      DO 1 I = 1,N
        IND(I) = I
    1   CONTINUE
      M = 1
      I = 1
      J = N
      R = .375
C
C TOP OF LOOP
C
    2 IF (I .GE. J) GO TO 10
      IF (R .GT. .5898437) GO TO 3
      R = R + .0390625
      GO TO 4
    3 R = R - .21875
C
C INITIALIZE K
C
    4 K = I
C
C SELECT A CENTRAL ELEMENT OF X AND SAVE IT IN T
C
      IJ = I + IFIX(R*FLOAT(J-I))
      IT = IND(IJ)
      T = X(IT)
C
C IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T,
C   INTERCHANGE IT WITH T
C
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 5
      IND(IJ) = INDX
      IND(I) = IT
      IT = INDX
      T = X(IT)
C
C INITIALIZE L
C
    5 L = J
C
C IF THE LAST ELEMENT OF THE ARRAY IS LESS THAN T,
C   INTERCHANGE IT WITH T
C
      INDX = IND(J)
      IF (X(INDX) .GE. T) GO TO 7
      IND(IJ) = INDX
      IND(J) = IT
      IT = INDX
      T = X(IT)
C
C IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T,
C   INTERCHANGE IT WITH T
C
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 7
      IND(IJ) = INDX
      IND(I) = IT
      IT = INDX
      T = X(IT)
      GO TO 7
C
C INTERCHANGE ELEMENTS K AND L
C
    6 ITT = IND(L)
      IND(L) = IND(K)
      IND(K) = ITT
C
C FIND AN ELEMENT IN THE UPPER PART OF THE ARRAY WHICH IS
C   NOT LARGER THAN T
C
    7 L = L - 1
      INDX = IND(L)
      IF (X(INDX) .GT. T) GO TO 7
C
C FIND AN ELEMENT IN THE LOWER PART OF THE ARRAY WHCIH IS
C   NOT SMALLER THAN T
C
    8 K = K + 1
      INDX = IND(K)
      IF (X(INDX) .LT. T) GO TO 8
C
C IF K .LE. L, INTERCHANGE ELEMENTS K AND L
C
      IF (K .LE. L) GO TO 6
C
C SAVE THE UPPER AND LOWER SUBSCRIPTS OF THE PORTION OF THE
C   ARRAY YET TO BE SORTED
C
      IF (L-I .LE. J-K) GO TO 9
      IL(M) = I
      IU(M) = L
      I = K
      M = M + 1
      GO TO 11
C
    9 IL(M) = K
      IU(M) = J
      J = L
      M = M + 1
      GO TO 11
C
C BEGIN AGAIN ON ANOTHER UNSORTED PORTION OF THE ARRAY
C
   10 M = M - 1
      IF (M .EQ. 0) RETURN
      I = IL(M)
      J = IU(M)
C
   11 IF (J-I .GE. 11) GO TO 4
      IF (I .EQ. 1) GO TO 2
      I = I - 1
C
C SORT ELEMENTS I+1,...,J.  NOTE THAT 1 .LE. I .LT. J AND
C   J-I .LT. 11.
C
   12 I = I + 1
      IF (I .EQ. J) GO TO 10
      INDX = IND(I+1)
      T = X(INDX)
      IT = INDX
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 12
      K = I
C
   13 IND(K+1) = IND(K)
      K = K - 1
      INDX = IND(K)
      IF (T .LT. X(INDX)) GO TO 13
      IND(K+1) = IT
      GO TO 12
      END
      SUBROUTINE REORDR (N,IFLAG, A,B,C, IND)
      INTEGER N, IFLAG, IND(N)
      REAL    A(N), B(N), C(N)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS SUBROUTINE USES AN ORDER N*LOG(N) QUICK SORT TO
C REORDER THE REAL ARRAY A INTO INCREASING ORDER.  A RECORD
C OF THE PERMUTATIONS APPLIED TO A IS STORED IN IND, AND
C THESE PERMUTATIONS MAY BE APPLIED TO ONE OR TWO ADDITIONAL
C VECTORS BY THIS ROUTINE.  ANY OTHER VECTOR V MAY BE PER-
C MUTED IN THE SAME FASHION BY CALLING SUBROUTINE PERMUT
C WITH N, IND, AND V AS PARAMETERS.
C   A SET OF NODES (X(I),Y(I)) AND DATA VALUES Z(I) MAY BE
C PREPROCESSED BY REORDR FOR INCREASED EFFICIENCY IN THE
C TRIANGULATION ROUTINE TRMESH.  EFFICIENCY IS INCREASED BY
C A FACTOR OF APPROXIMATELY SQRT(N)/6 FOR RANDOMLY DISTRIB-
C UTED NODES, AND THE PREPROCESSING IS ALSO USEFUL FOR
C DETECTING DUPLICATE NODES.  EITHER X OR Y MAY BE USED AS
C THE SORT KEY (ASSOCIATED WITH A).
C
C INPUT PARAMETERS - N - NUMBER OF NODES.
C
C                IFLAG - NUMBER OF VECTORS TO BE PERMUTED.
C                        IFLAG .LE. 0 IF A, B, AND C ARE TO
C                                     REMAIN UNALTERED.
C                        IFLAG .EQ. 1 IF ONLY A IS TO BE
C                                     PERMUTED.
C                        IFLAG .EQ. 2 IF A AND B ARE TO BE
C                                     PERMUTED.
C                        IFLAG .GE. 3 IF A, B, AND C ARE TO
C                                     BE PERMUTED.
C
C                A,B,C - VECTORS OF LENGTH N TO BE SORTED
C                        (ON THE COMPONENTS OF A), OR DUMMY
C                        PARAMETERS, DEPENDING ON IFLAG.
C
C                  IND - VECTOR OF LENGTH .GE. N.
C
C N, IFLAG, AND ANY DUMMY PARAMETERS ARE NOT ALTERED BY THIS
C   ROUTINE.
C
C OUTPUT PARAMETERS - A,B,C - SORTED OR UNALTERED VECTORS.
C
C                       IND - SEQUENCE OF INDICES 1,...,N
C                             PERMUTED IN THE SAME FASHION
C                             AS THE REAL VECTORS.  THUS,
C                             THE ORDERING MAY BE APPLIED TO
C                             A REAL VECTOR V AND STORED IN
C                             W BY SETTING W(I) = V(IND(I)),
C                             OR V MAY BE OVERWRITTEN WITH
C                             THE ORDERING BY A CALL TO PER-
C                             MUT.
C
C MODULES REFERENCED BY REORDR - QSORT, PERMUT
C
C***********************************************************
C
      INTEGER NN, NV
C
C LOCAL PARAMETERS -
C
C NN = LOCAL COPY OF N
C NV = LOCAL COPY OF IFLAG
C
      NN = N
      NV = IFLAG
      CALL QSORT(NN,A, IND)
      IF (NV .LE. 0) RETURN
      CALL PERMUT(NN,IND, A )
      IF (NV .EQ. 1) RETURN
      CALL PERMUT(NN,IND, B )
      IF (NV .EQ. 2) RETURN
      CALL PERMUT(NN,IND, C )
      RETURN
      END
      SUBROUTINE SHIFTD (NFRST,NLAST,KK, IARR )
      INTEGER NFRST, NLAST, KK, IARR(1)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE SHIFTS A SET OF CONTIGUOUS ELEMENTS OF AN
C INTEGER ARRAY KK POSITIONS DOWNWARD (UPWARD IF KK .LT. 0).
C THE LOOPS ARE UNROLLED IN ORDER TO INCREASE EFFICIENCY.
C
C INPUT PARAMETERS - NFRST,NLAST - BOUNDS ON THE PORTION OF
C                                  IARR TO BE SHIFTED.  ALL
C                                  ELEMENTS BETWEEN AND
C                                  INCLUDING THE BOUNDS ARE
C                                  SHIFTED UNLESS NFRST .GT.
C                                  NLAST, IN WHICH CASE NO
C                                  SHIFT OCCURS.
C
C                             KK - NUMBER OF POSITIONS EACH
C                                  ELEMENT IS TO BE SHIFTED.
C                                  IF KK .LT. 0 SHIFT UP.
C                                  IF KK .GT. 0 SHIFT DOWN.
C
C                           IARR - INTEGER ARRAY OF LENGTH
C                                  .GE. NLAST + MAX(KK,0).
C
C NFRST, NLAST, AND KK ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETER -        IARR - SHIFTED ARRAY.
C
C MODULES REFERENCED BY SHIFTD - NONE
C
C***********************************************************
C
      INTEGER INC, K, NF, NL, NLP1, NS, NSL, I, IBAK, INDX,
     .        IMAX
      DATA    INC/5/
C
C LOCAL PARAMETERS -
C
C INC =  DO-LOOP INCREMENT (UNROLLING FACTOR) -- IF INC IS
C          CHANGED, STATEMENTS MUST BE ADDED TO OR DELETED
C          FROM THE DO-LOOPS
C K =    LOCAL COPY OF KK
C NF =   LOCAL COPY OF NFRST
C NL =   LOCAL COPY OF NLAST
C NLP1 = NL + 1
C NS =   NUMBER OF SHIFTS
C NSL =  NUMBER OF SHIFTS DONE IN UNROLLED DO-LOOP (MULTIPLE
C          OF INC)
C I =    DO-LOOP INDEX AND INDEX FOR IARR
C IBAK = INDEX FOR DOWNWARD SHIFT OF IARR
C INDX = INDEX FOR IARR
C IMAX = BOUND ON DO-LOOP INDEX
C
      K = KK
      NF = NFRST
      NL = NLAST
      IF (NF .GT. NL  .OR.  K .EQ. 0) RETURN
      NLP1 = NL + 1
      NS = NLP1 - NF
      NSL = INC*(NS/INC)
      IF ( K .LT. 0) GO TO 4
C
C SHIFT DOWNWARD STARTING FROM THE BOTTOM
C
      IF (NSL .LE. 0) GO TO 2
      DO 1 I = 1,NSL,INC
        IBAK = NLP1 - I
        INDX = IBAK + K
        IARR(INDX) = IARR(IBAK)
        IARR(INDX-1) = IARR(IBAK-1)
        IARR(INDX-2) = IARR(IBAK-2)
        IARR(INDX-3) = IARR(IBAK-3)
        IARR(INDX-4) = IARR(IBAK-4)
    1   CONTINUE
C
C PERFORM THE REMAINING NS-NSL SHIFTS ONE AT A TIME
C
    2 IBAK = NLP1 - NSL
    3 IF (IBAK .LE. NF) RETURN
      IBAK = IBAK - 1
      INDX = IBAK + K
      IARR(INDX) = IARR(IBAK)
      GO TO 3
C
C SHIFT UPWARD STARTING FROM THE TOP
C
    4 IF (NSL .LE. 0) GO TO 6
      IMAX = NLP1 - INC
      DO 5 I = NF,IMAX,INC
        INDX = I + K
        IARR(INDX) = IARR(I)
        IARR(INDX+1) = IARR(I+1)
        IARR(INDX+2) = IARR(I+2)
        IARR(INDX+3) = IARR(I+3)
        IARR(INDX+4) = IARR(I+4)
    5   CONTINUE
C
C PERFORM THE REMAINING NS-NSL SHIFTS ONE AT A TIME
C
    6 I = NSL + NF
    7 IF (I .GT. NL) RETURN
      INDX = I + K
      IARR(INDX) = IARR(I)
      I = I + 1
      GO TO 7
      END
      SUBROUTINE SWAP (NIN1,NIN2,NOUT1,NOUT2, IADJ,IEND )
      INTEGER NIN1, NIN2, NOUT1, NOUT2, IADJ(1), IEND(1)
      EXTERNAL INDEX
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS SUBROUTINE SWAPS THE DIAGONALS IN A CONVEX QUADRI-
C LATERAL.
C
C INPUT PARAMETERS -  NIN1,NIN2,NOUT1,NOUT2 - NODAL INDICES
C                            OF A PAIR OF ADJACENT TRIANGLES
C                            WHICH FORM A CONVEX QUADRILAT-
C                            ERAL.  NOUT1 AND NOUT2 ARE CON-
C                            NECTED BY AN ARC WHICH IS TO BE
C                            REPLACED BY THE ARC NIN1-NIN2.
C                            (NIN1,NOUT1,NOUT2) MUST BE TRI-
C                            ANGLE VERTICES IN COUNTERCLOCK-
C                            WISE ORDER.
C
C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C                IADJ,IEND - TRIANGULATION DATA STRUCTURE
C                            (SEE SUBROUTINE TRMESH).
C
C OUTPUT PARAMETERS - IADJ,IEND - UPDATED WITH THE ARC
C                                 REPLACEMENT.
C
C MODULES REFERENCED BY SWAP - INDEX, SHIFTD
C
C***********************************************************
C
      INTEGER IN(2), IO(2), IP1, IP2, J, K, NF, NL, I,
     .        IMIN, IMAX
C
C LOCAL PARAMETERS -
C
C IN =        NIN1 AND NIN2 ORDERED BY INCREASING MAGNITUDE
C               (THE NEIGHBORS OF IN(1) PRECEDE THOSE OF
C               IN(2) IN IADJ)
C IO =        NOUT1 AND NOUT2 IN INCREASING ORDER
C IP1,IP2 =   PERMUTATION OF (1,2) SUCH THAT IO(IP1)
C               PRECEDES IO(IP2) AS A NEIGHBOR OF IN(1)
C J,K =       PERMUTATION OF (1,2) USED AS INDICES OF IN
C               AND IO
C NF,NL =     IADJ INDICES BOUNDARY A PORTION OF THE ARRAY
C               TO BE SHIFTED
C I =         IEND INDEX
C IMIN,IMAX = BOUNDS ON THE PORTION OF IEND TO BE INCRE-
C               MENTED OR DECREMENTED
C
      IN(1) = NIN1
      IN(2) = NIN2
      IO(1) = NOUT1
      IO(2) = NOUT2
      IP1 = 1
C
C ORDER THE INDICES SO THAT IN(1) .LT. IN(2) AND IO(1) .LT.
C   IO(2), AND CHOOSE IP1 AND IP2 SUCH THAT (IN(1),IO(IP1),
C   IO(IP2)) FORMS A TRIANGLE.
C
      IF (IN(1) .LT. IN(2)) GO TO 1
      IN(1) = IN(2)
      IN(2) = NIN1
      IP1 = 2
    1 IF (IO(1) .LT. IO(2)) GO TO 2
      IO(1) = IO(2)
      IO(2) = NOUT1
      IP1 = 3 - IP1
    2 IP2 = 3 - IP1
      IF (IO(2) .LT. IN(1)) GO TO 8
      IF (IN(2) .LT. IO(1)) GO TO 12
C
C IN(1) AND IO(1) PRECEDE IN(2) AND IO(2).  FOR (J,K) =
C   (1,2) AND (2,1), DELETE IO(K) AS A NEIGHBOR OF IO(J)
C   BY SHIFTING A PORTION OF IADJ EITHER UP OR DOWN AND
C   AND INSERT IN(K) AS A NEIGHBOR OF IN(J).
C
      DO 7 J = 1,2
        K = 3 - J
        IF (IN(J) .GT. IO(J)) GO TO 4
C
C   THE NEIGHBORS OF IN(J) PRECEDE THOSE OF IO(J) -- SHIFT
C     DOWN BY 1
C
        NF = 1 + INDEX(IN(J),IO(IP1),IADJ,IEND)
        NL = -1 + INDEX(IO(J),IO(K),IADJ,IEND)
        IF (NF .LE. NL) CALL SHIFTD(NF,NL,1, IADJ )
        IADJ(NF) = IN(K)
        IMIN = IN(J)
        IMAX = IO(J)-1
        DO 3 I = IMIN,IMAX
    3     IEND(I) = IEND(I) + 1
        GO TO 6
C
C   THE NEIGHBORS OF IO(J) PRECEDE THOSE OF IN(J) -- SHIFT
C     UP BY 1
C
    4   NF = 1 + INDEX(IO(J),IO(K),IADJ,IEND)
        NL = -1 + INDEX(IN(J),IO(IP2),IADJ,IEND)
        IF (NF .LE. NL) CALL SHIFTD(NF,NL,-1, IADJ )
        IADJ(NL) = IN(K)
        IMIN = IO(J)
        IMAX = IN(J) - 1
        DO 5 I = IMIN,IMAX
    5     IEND(I) = IEND(I) - 1
C
C   REVERSE (IP1,IP2) FOR (J,K) = (2,1)
C
    6   IP1 = IP2
        IP2 = 3 - IP1
    7   CONTINUE
      RETURN
C
C THE VERTICES ARE ORDERED (IO(1),IO(2),IN(1),IN(2)).
C   DELETE IO(2) BY SHIFTING UP BY 1
C
    8 NF = 1 + INDEX(IO(1),IO(2),IADJ,IEND)
      NL = -1 + INDEX(IO(2),IO(1),IADJ,IEND)
      IF (NF .LE. NL) CALL SHIFTD(NF,NL,-1, IADJ )
      IMIN = IO(1)
      IMAX = IO(2)-1
      DO 9 I = IMIN,IMAX
    9   IEND(I) = IEND(I) - 1
C
C   DELETE IO(1) BY SHIFTING UP BY 2 AND INSERT IN(2)
C
      NF = NL + 2
      NL = -1 + INDEX(IN(1),IO(IP2),IADJ,IEND)
      IF (NF .LE. NL) CALL SHIFTD(NF,NL,-2, IADJ )
      IADJ(NL-1) = IN(2)
      IMIN = IO(2)
      IMAX = IN(1)-1
      DO 10 I = IMIN,IMAX
   10   IEND(I) = IEND(I) - 2
C
C   SHIFT UP BY 1 AND INSERT IN(1)
C
      NF = NL + 1
      NL = -1 + INDEX(IN(2),IO(IP1),IADJ,IEND)
      CALL SHIFTD(NF,NL,-1, IADJ )
      IADJ(NL) = IN(1)
      IMIN = IN(1)
      IMAX = IN(2)-1
      DO 11 I = IMIN,IMAX
   11   IEND(I) = IEND(I) - 1
      RETURN
C
C THE VERTICES ARE ORDERED (IN(1),IN(2),IO(1),IO(2)).
C   DELETE IO(1) BY SHIFTING DOWN BY 1
C
   12 NF = 1 + INDEX(IO(1),IO(2),IADJ,IEND)
      NL = -1 + INDEX(IO(2),IO(1),IADJ,IEND)
      IF (NF .LE. NL) CALL SHIFTD(NF,NL,1, IADJ )
      IMIN = IO(1)
      IMAX = IO(2) - 1
      DO 13 I = IMIN,IMAX
   13   IEND(I) = IEND(I) + 1
C
C   DELETE IO(2) BY SHIFTING DOWN BY 2 AND INSERT IN(1)
C
      NL = NF - 2
      NF = 1 + INDEX(IN(2),IO(IP2),IADJ,IEND)
      IF (NF .LE. NL) CALL SHIFTD(NF,NL,2, IADJ )
      IADJ(NF+1) = IN(1)
      IMIN = IN(2)
      IMAX = IO(1) - 1
      DO 14 I = IMIN,IMAX
   14   IEND(I) = IEND(I) + 2
C
C   SHIFT DOWN BY 1 AND INSERT IN(2)
C
      NL = NF - 1
      NF = 1 + INDEX(IN(1),IO(IP1),IADJ,IEND)
      CALL SHIFTD(NF,NL,1, IADJ )
      IADJ(NF) = IN(2)
      IMIN = IN(1)
      IMAX = IN(2) - 1
      DO 15 I = IMIN,IMAX
   15   IEND(I) = IEND(I) + 1
      RETURN
      END
      LOGICAL FUNCTION SWPTST (IN1,IN2,IO1,IO2,X,Y)
      INTEGER IN1, IN2, IO1, IO2
      REAL    X(1), Y(1)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS FUNCTION DECIDES WHETHER OR NOT TO REPLACE A
C DIAGONAL ARC IN A QUADRILATERAL WITH THE OTHER DIAGONAL.
C THE DETERMINATION IS BASED ON THE SIZES OF THE ANGLES
C CONTAINED IN THE 2 TRIANGLES DEFINED BY THE DIAGONAL.
C THE DIAGONAL IS CHOSEN TO MAXIMIZE THE SMALLEST OF THE
C SIX ANGLES OVER THE TWO PAIRS OF TRIANGLES.
C
C INPUT PARAMETERS -  IN1,IN2,IO1,IO2 - NODE INDICES OF THE
C                              FOUR POINTS DEFINING THE
C                              QUADRILATERAL.  IO1 AND IO2
C                              ARE CURRENTLY CONNECTED BY A
C                              DIAGONAL ARC.  THIS ARC
C                              SHOULD BE REPLACED BY AN ARC
C                              CONNECTING IN1, IN2 IF THE
C                              DECISION IS MADE TO SWAP.
C                              IN1,IO1,IO2 MUST BE IN
C                              COUNTERCLOCKWISE ORDER.
C
C                        X,Y - VECTORS OF NODAL COORDINATES.
C                              (X(I),Y(I)) ARE THE COORD-
C                              INATES OF NODE I FOR I = IN1,
C                              IN2, IO1, OR IO2.
C
C NONE OF THE INPUT PARAMETERS ARE ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETER -  SWPTST - .TRUE. IFF THE ARC CONNECTING
C                              IO1 AND IO2 IS TO BE REPLACED
C
C MODULES REFERENCED BY SWPTST - NONE
C
C***********************************************************
C
      REAL DX11, DX12, DX22, DX21, DY11, DY12, DY22, DY21,
     .     SIN1, SIN2, COS1, COS2, SIN12
C
C LOCAL PARAMETERS -
C
C DX11,DY11 = X,Y COORDINATES OF THE VECTOR IN1-IO1
C DX12,DY12 = X,Y COORDINATES OF THE VECTOR IN1-IO2
C DX22,DY22 = X,Y COORDINATES OF THE VECTOR IN2-IO2
C DX21,DY21 = X,Y COORDINATES OF THE VECTOR IN2-IO1
C SIN1 =      CROSS PRODUCT OF THE VECTORS IN1-IO1 AND
C               IN1-IO2 -- PROPORTIONAL TO SIN(T1) WHERE T1
C               IS THE ANGLE AT IN1 FORMED BY THE VECTORS
C COS1 =      INNER PRODUCT OF THE VECTORS IN1-IO1 AND
C               IN1-IO2 -- PROPORTIONAL TO COS(T1)
C SIN2 =      CROSS PRODUCT OF THE VECTORS IN2-IO2 AND
C               IN2-IO1 -- PROPORTIONAL TO SIN(T2) WHERE T2
C               IS THE ANGLE AT IN2 FORMED BY THE VECTORS
C COS2 =      INNER PRODUCT OF THE VECTORS IN2-IO2 AND
C               IN2-IO1 -- PROPORTIONAL TO COS(T2)
C SIN12 =     SIN1*COS2 + COS1*SIN2 -- PROPORTIONAL TO
C               SIN(T1+T2)
C
      SWPTST = .FALSE.
C
C COMPUTE THE VECTORS CONTAINING THE ANGLES T1, T2
C
      DX11 = X(IO1) - X(IN1)
      DX12 = X(IO2) - X(IN1)
      DX22 = X(IO2) - X(IN2)
      DX21 = X(IO1) - X(IN2)
C
      DY11 = Y(IO1) - Y(IN1)
      DY12 = Y(IO2) - Y(IN1)
      DY22 = Y(IO2) - Y(IN2)
      DY21 = Y(IO1) - Y(IN2)
C
C COMPUTE INNER PRODUCTS
C
      COS1 = DX11*DX12 + DY11*DY12
      COS2 = DX22*DX21 + DY22*DY21
C
C THE DIAGONALS SHOULD BE SWAPPED IFF (T1+T2) .GT. 180
C   DEGREES.  THE FOLLOWING TWO TESTS INSURE NUMERICAL
C   STABILITY.
C
      IF (COS1 .GE. 0.  .AND.  COS2 .GE. 0.) RETURN
      IF (COS1 .LT. 0.  .AND.  COS2 .LT. 0.) GO TO 1
C
C COMPUTE VECTOR CROSS PRODUCTS
C
      SIN1 = DX11*DY12 - DX12*DY11
      SIN2 = DX22*DY21 - DX21*DY22
      SIN12 = SIN1*COS2 + COS1*SIN2
      IF (SIN12 .GE. 0.) RETURN
    1 SWPTST = .TRUE.
      RETURN
      END
      SUBROUTINE TRFIND (NST,PX,PY,X,Y,IADJ,IEND, I1,I2,I3)
      INTEGER NST, IADJ(1), IEND(1), I1, I2, I3
      REAL    PX, PY, X(1), Y(1)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE LOCATES A POINT P IN A THIESSEN TRIANGU-
C LATION, RETURNING THE VERTEX INDICES OF A TRIANGLE WHICH
C CONTAINS P.  TRFIND IS PART OF AN INTERPOLATION PACKAGE
C WHICH PROVIDES SUBROUTINES FOR CREATING THE MESH.
C
C INPUT PARAMETERS -    NST - INDEX OF NODE AT WHICH TRFIND
C                             BEGINS SEARCH.  SEARCH TIME
C                             DEPENDS ON THE PROXIMITY OF
C                             NST TO P.
C
C                     PX,PY - X AND Y-COORDINATES OF THE
C                             POINT TO BE LOCATED.
C
C                       X,Y - VECTORS OF COORDINATES OF
C                             NODES IN THE MESH.  (X(I),Y(I))
C                             DEFINES NODE I FOR I = 1,...,N
C                             WHERE N .GE. 3.
C
C                      IADJ - SET OF ADJACENCY LISTS OF
C                             NODES IN THE MESH.
C
C                      IEND - POINTERS TO THE ENDS OF
C                             ADJACENCY LISTS IN IADJ FOR
C                             EACH NODE IN THE MESH.
C
C IADJ AND IEND MAY BE CREATED BY TRMESH.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETERS - I1,I2,I3 - VERTEX INDICES IN COUNTER-
C                                CLOCKWISE ORDER - VERTICES
C                                OF A TRIANGLE CONTAINING P
C                                IF P IS AN INTERIOR NODE.
C                                IF P IS OUTSIDE OF THE
C                                BOUNDARY OF THE MESH, I1
C                                AND I2 ARE THE FIRST (RIGHT
C                                -MOST) AND LAST (LEFTMOST)
C                                NODES WHICH ARE VISIBLE
C                                FROM P, AND I3 = 0.  IF P
C                                AND ALL OF THE NODES LIE ON
C                                A SINGLE LINE THEN I1 = I2
C                                = I3 = 0.
C
C MODULES REFERENCED BY TRFIND - NONE
C
C INTRINSIC FUNCTION CALLED BY TRFIND - MAX0
C
C***********************************************************
C
      INTEGER N0, N1, N2, N3, N4, INDX, IND, NF,
     .        NL, NEXT
      REAL    XP, YP
      LOGICAL LEFT
C
C LOCAL PARAMETERS -
C
C XP,YP =     LOCAL VARIABLES CONTAINING PX AND PY
C N0,N1,N2 =  NODES IN COUNTERCLOCKWISE ORDER DEFINING A
C               CONE (WITH VERTEX N0) CONTAINING P
C N3,N4 =     NODES OPPOSITE N1-N2 AND N2-N1, RESPECTIVELY
C INDX,IND =  INDICES FOR IADJ
C NF,NL =     FIRST AND LAST NEIGHBORS OF N0 IN IADJ, OR
C               FIRST (RIGHTMOST) AND LAST (LEFTMOST) NODES
C               VISIBLE FROM P WHEN P IS OUTSIDE THE
C               BOUNDARY
C NEXT =      CANDIDATE FOR I1 OR I2 WHEN P IS OUTSIDE OF
C               THE BOUNDARY
C LEFT =      STATEMENT FUNCTION WHICH COMPUTES THE SIGN OF
C               A CROSS PRODUCT (Z-COMPONENT).  LEFT(X1,...,
C               Y0) = .TRUE. IFF (X0,Y0) IS ON OR TO THE
C               LEFT OF THE VECTOR FROM (X1,Y1) TO (X2,Y2).
C
      LEFT(X1,Y1,X2,Y2,X0,Y0) = (X2-X1)*(Y0-Y1) .GE.
     .                          (X0-X1)*(Y2-Y1)
      XP = PX
      YP = PY
C
C INITIALIZE VARIABLES AND FIND A CONE CONTAINING P
C
      N0 = MAX0(NST,1)
    1 INDX = IEND(N0)
      NL = IADJ(INDX)
      INDX = 1
      IF (N0 .NE. 1) INDX = IEND(N0-1) + 1
      NF = IADJ(INDX)
      N1 = NF
      IF (NL .NE. 0) GO TO 3
C
C N0 IS A BOUNDARY NODE.  SET NL TO THE LAST NONZERO
C   NEIGHBOR OF N0.
C
      IND = IEND(N0) - 1
      NL = IADJ(IND)
      IF ( LEFT(X(N0),Y(N0),X(NF),Y(NF),XP,YP) ) GO TO 2
C
C P IS OUTSIDE THE BOUNDARY
C
      NL = N0
      GO TO 16
    2 IF ( LEFT(X(NL),Y(NL),X(N0),Y(N0),XP,YP) ) GO TO 4
C
C P IS OUTSIDE THE BOUNDARY AND N0 IS THE RIGHTMOST
C   VISIBLE BOUNDARY NODE
C
      I1 = N0
      GO TO 18
C
C N0 IS AN INTERIOR NODE.  FIND N1.
C
    3 IF ( LEFT(X(N0),Y(N0),X(N1),Y(N1),XP,YP) ) GO TO 4
      INDX = INDX + 1
      N1 = IADJ(INDX)
      IF (N1 .EQ. NL) GO TO 7
      GO TO 3
C
C P IS TO THE LEFT OF ARC N0-N1.  INITIALIZE N2 TO THE NEXT
C   NEIGHBOR OF N0.
C
    4 INDX = INDX + 1
      N2 = IADJ(INDX)
      IF ( .NOT. LEFT(X(N0),Y(N0),X(N2),Y(N2),XP,YP) )
     .   GO TO 8
      N1 = N2
      IF (N1 .NE. NL) GO TO 4
      IF ( .NOT. LEFT(X(N0),Y(N0),X(NF),Y(NF),XP,YP) )
     .   GO TO 7
      IF (XP .EQ. X(N0) .AND. YP .EQ. Y(N0)) GO TO 6
C
C P IS LEFT OF OR ON ARCS N0-NB FOR ALL NEIGHBORS NB
C   OF N0.
C ALL POINTS ARE COLLINEAR IFF P IS LEFT OF NB-N0 FOR
C   ALL NEIGHBORS NB OF N0.  SEARCH THE NEIGHBORS OF N0
C   IN REVERSE ORDER.  NOTE -- N1 = NL AND INDX POINTS TO
C   NL.
C
    5 IF ( .NOT. LEFT(X(N1),Y(N1),X(N0),Y(N0),XP,YP) )
     .   GO TO 6
      IF (N1 .EQ. NF) GO TO 20
      INDX = INDX - 1
      N1 = IADJ(INDX)
      GO TO 5
C
C P IS TO THE RIGHT OF N1-N0, OR P=N0.  SET N0 TO N1 AND
C   START OVER.
C
    6 N0 = N1
      GO TO 1
C
C P IS BETWEEN ARCS N0-N1 AND N0-NF
C
    7 N2 = NF
C
C P IS CONTAINED IN A CONE DEFINED BY LINE SEGMENTS N0-N1
C   AND N0-N2 WHERE N1 IS ADJACENT TO N2
C
    8 N3 = N0
    9 IF ( LEFT(X(N1),Y(N1),X(N2),Y(N2),XP,YP) ) GO TO 13
C
C SET N4 TO THE FIRST NEIGHBOR OF N2 FOLLOWING N1
C
      INDX = IEND(N2)
      IF (IADJ(INDX) .NE. N1) GO TO 10
C
C N1 IS THE LAST NEIGHBOR OF N2.
C SET N4 TO THE FIRST NEIGHBOR.
C
      INDX = 1
      IF (N2 .NE. 1) INDX = IEND(N2-1) + 1
      N4 = IADJ(INDX)
      GO TO 11
C
C N1 IS NOT THE LAST NEIGHBOR OF N2
C
   10 INDX = INDX-1
      IF (IADJ(INDX) .NE. N1) GO TO 10
      N4 = IADJ(INDX+1)
      IF (N4 .NE. 0) GO TO 11
C
C P IS OUTSIDE THE BOUNDARY
C
      NF = N2
      NL = N1
      GO TO 16
C
C DEFINE A NEW ARC N1-N2 WHICH INTERSECTS THE LINE
C   SEGMENT N0-P
C
   11 IF ( LEFT(X(N0),Y(N0),X(N4),Y(N4),XP,YP) ) GO TO 12
      N3 = N2
      N2 = N4
      GO TO 9
   12 N3 = N1
      N1 = N4
      GO TO 9
C
C P IS IN THE TRIANGLE (N1,N2,N3) AND NOT ON N2-N3.  IF
C   N3-N1 OR N1-N2 IS A BOUNDARY ARC CONTAINING P, TREAT P
C   AS EXTERIOR.
C
   13 INDX = IEND(N1)
      IF (IADJ(INDX) .NE. 0) GO TO 15
C
C N1 IS A BOUNDARY NODE.  N3-N1 IS A BOUNDARY ARC IFF N3
C   IS THE LAST NONZERO NEIGHBOR OF N1.
C
      IF (N3 .NE. IADJ(INDX-1)) GO TO 14
C
C N3-N1 IS A BOUNDARY ARC
C
      IF ( .NOT. LEFT(X(N1),Y(N1),X(N3),Y(N3),XP,YP) )
     .   GO TO 14
C
C P LIES ON N1-N3
C
      I1 = N1
      I2 = N3
      I3 = 0
      RETURN
C
C N3-N1 IS NOT A BOUNDARY ARC CONTAINING P.  N1-N2 IS A
C   BOUNDARY ARC IFF N2 IS THE FIRST NEIGHBOR OF N1.
C
   14 INDX = 1
      IF (N1 .NE. 1) INDX = IEND(N1-1) + 1
      IF (N2 .NE. IADJ(INDX)) GO TO 15
C
C N1-N2 IS A BOUNDARY ARC
C
      IF ( .NOT. LEFT(X(N2),Y(N2),X(N1),Y(N1),XP,YP) )
     .   GO TO 15
C
C P LIES ON N1-N2
C
      I1 = N2
      I2 = N1
      I3 = 0
      RETURN
C
C P DOES NOT LIE ON A BOUNDARY ARC.
C
   15 I1 = N1
      I2 = N2
      I3 = N3
      RETURN
C
C NF AND NL ARE ADJACENT BOUNDARY NODES WHICH ARE VISIBLE
C   FROM P.  FIND THE FIRST VISIBLE BOUNDARY NODE.
C SET NEXT TO THE FIRST NEIGHBOR OF NF.
C
   16 INDX = 1
      IF (NF .NE. 1) INDX = IEND(NF-1) + 1
      NEXT = IADJ(INDX)
      IF ( LEFT(X(NF),Y(NF),X(NEXT),Y(NEXT),XP,YP) )
     .   GO TO 17
      NF = NEXT
      GO TO 16
C
C NF IS THE FIRST (RIGHTMOST) VISIBLE BOUNDARY NODE
C
   17 I1 = NF
C
C FIND THE LAST VISIBLE BOUNDARY NODE.  NL IS THE FIRST
C   CANDIDATE FOR I2.
C SET NEXT TO THE LAST NEIGHBOR OF NL.
C
   18 INDX = IEND(NL) - 1
      NEXT = IADJ(INDX)
      IF ( LEFT(X(NEXT),Y(NEXT),X(NL),Y(NL),XP,YP) )
     .   GO TO 19
      NL = NEXT
      GO TO 18
C
C NL IS THE LAST (LEFTMOST) VISIBLE BOUNDARY NODE
C
   19 I2 = NL
      I3 = 0
      RETURN
C
C ALL POINTS ARE COLLINEAR
C
   20 I1 = 0
      I2 = 0
      I3 = 0
      RETURN
      END
      SUBROUTINE TRMESH (N,X,Y, IADJ,IEND,IER)
      INTEGER N, IADJ(1), IEND(N), IER
      REAL    X(N), Y(N)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE CREATES A THIESSEN TRIANGULATION OF N
C ARBITRARILY SPACED POINTS IN THE PLANE REFERRED TO AS
C NODES.  THE TRIANGULATION IS OPTIMAL IN THE SENSE THAT IT
C IS AS NEARLY EQUIANGULAR AS POSSIBLE.  TRMESH IS PART OF
C AN INTERPOLATION PACKAGE WHICH ALSO PROVIDES SUBROUTINES
C TO REORDER THE NODES, ADD A NEW NODE, DELETE AN ARC, PLOT
C THE MESH, AND PRINT THE DATA STRUCTURE.
C   UNLESS THE NODES ARE ALREADY ORDERED IN SOME REASONABLE
C FASHION, THEY SHOULD BE REORDERED BY SUBROUTINE REORDR FOR
C INCREASED EFFICIENCY BEFORE CALLING TRMESH.
C
C INPUT PARAMETERS -     N - NUMBER OF NODES IN THE MESH.
C                            N .GE. 3.
C
C                      X,Y - N-VECTORS OF COORDINATES.
C                            (X(I),Y(I)) DEFINES NODE I.
C
C                     IADJ - VECTOR OF LENGTH .GE. 6*N-9.
C
C                     IEND - VECTOR OF LENGTH .GE. N.
C
C N, X, AND Y ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETERS - IADJ - ADJACENCY LISTS OF NEIGHBORS IN
C                            COUNTERCLOCKWISE ORDER.  THE
C                            LIST FOR NODE I+1 FOLLOWS THAT
C                            FOR NODE I WHERE X AND Y DEFINE
C                            THE ORDER.  THE VALUE 0 DENOTES
C                            THE BOUNDARY (OR A PSEUDO-NODE
C                            AT INFINITY) AND IS ALWAYS THE
C                            LAST NEIGHBOR OF A BOUNDARY
C                            NODE.  IADJ IS UNCHANGED IF IER
C                            .NE. 0.
C
C                     IEND - POINTERS TO THE ENDS OF
C                            ADJACENCY LISTS (SETS OF
C                            NEIGHBORS) IN IADJ.  THE
C                            NEIGHBORS OF NODE 1 BEGIN IN
C                            IADJ(1).  FOR K .GT. 1, THE
C                            NEIGHBORS OF NODE K BEGIN IN
C                            IADJ(IEND(K-1)+1) AND K HAS
C                            IEND(K) - IEND(K-1) NEIGHBORS
C                            INCLUDING (POSSIBLY) THE
C                            BOUNDARY.  IADJ(IEND(K)) .EQ. 0
C                            IFF NODE K IS ON THE BOUNDARY.
C                            IEND IS UNCHANGED IF IER = 1.
C                            IF IER = 2 IEND CONTAINS THE
C                            INDICES OF A SEQUENCE OF N
C                            NODES ORDERED FROM LEFT TO
C                            RIGHT WHERE LEFT AND RIGHT ARE
C                            DEFINED BY ASSUMING NODE 1 IS
C                            TO THE LEFT OF NODE 2.
C
C                      IER - ERROR INDICATOR
C                            IER = 0 IF NO ERRORS WERE
C                                    ENCOUNTERED.
C                            IER = 1 IF N .LT. 3.
C                            IER = 2 IF N .GE. 3 AND ALL
C                                    NODES ARE COLLINEAR.
C
C MODULES REFERENCED BY TRMESH - SHIFTD, ADNODE, TRFIND,
C                                INTADD, BDYADD, SWPTST,
C                                SWAP, INDEX
C
C***********************************************************
C
      INTEGER NN, K, KM1, NL, NR, IND, INDX, N0, ITEMP,
     .        IERR, KM1D2, KMI, I, KMIN
      REAL    XL, YL, XR, YR, DXR, DYR, XK, YK, DXK, DYK,
     .        CPROD, SPROD
C
C LOCAL PARAMETERS -
C
C NN =          LOCAL COPY OF N
C K =           NODE (INDEX) TO BE INSERTED INTO IEND
C KM1 =         K-1 - (VARIABLE) LENGTH OF IEND
C NL,NR =       IEND(1), IEND(KM1) -- LEFTMOST AND RIGHTMOST
C                 NODES IN IEND AS VIEWED FROM THE RIGHT OF
C                 1-2 WHEN IEND CONTAINS THE INITIAL ORDERED
C                 SET OF NODAL INDICES
C XL,YL,XR,YR = X AND Y COORDINATES OF NL AND NR
C DXR,DYR =     XR-XL, YR-YL
C XK,YK =       X AND Y COORDINATES OF NODE K
C DXK,DYK =     XK-XL, YK-YL
C CPROD =       VECTOR CROSS PRODUCT OF NL-NR AND NL-K --
C                 USED TO DETERMINE THE POSITION OF NODE K
C                 WITH RESPECT TO THE LINE DEFINED BY THE
C                 NODES IN IEND
C SPROD =       SCALAR PRODUCT USED TO DETERMINE THE
C                 INTERVAL CONTAINING NODE K WHEN K IS ON
C                 THE LINE DEFINED BY THE NODES IN IEND
C IND,INDX =    INDICES FOR IEND AND IADJ, RESPECTIVELY
C N0,ITEMP =    TEMPORARY NODES (INDICES)
C IERR =        DUMMY PARAMETER FOR CALL TO ADNODE
C KM1D2,KMI,I = KM1/2, K-I, DO-LOOP INDEX -- USED IN IEND
C                 REORDERING LOOP
C KMIN =        FIRST NODE INDEX SENT TO ADNODE
C
      NN = N
      IER = 1
      IF (NN .LT. 3) RETURN
      IER = 0
C
C INITIALIZE IEND, NL, NR, AND K
C
      IEND(1) = 1
      IEND(2) = 2
      XL = X(1)
      YL = Y(1)
      XR = X(2)
      YR = Y(2)
      K = 2
C
C BEGIN LOOP ON NODES 3,4,...
C
    1 DXR = XR-XL
      DYR = YR-YL
C
C NEXT LOOP BEGINS HERE IF NL AND NR ARE UNCHANGED
C
    2 IF (K .EQ. NN) GO TO 13
      KM1 = K
      K = KM1 + 1
      XK = X(K)
      YK = Y(K)
      DXK = XK-XL
      DYK = YK-YL
      CPROD = DXR*DYK - DXK*DYR
      IF (CPROD .GT. 0.) GO TO 6
      IF (CPROD .LT. 0.) GO TO 8
C
C NODE K LIES ON THE LINE CONTAINING NODES 1,2,...,K-1.
C   SET SPROD TO (NL-NR,NL-K).
C
      SPROD = DXR*DXK + DYR*DYK
      IF (SPROD .GT. 0.) GO TO 3
C
C NODE K IS TO THE LEFT OF NL.  INSERT K AS THE FIRST
C   (LEFTMOST) NODE IN IEND AND SET NL TO K.
C
      CALL SHIFTD(1,KM1,1, IEND )
      IEND(1) = K
      XL = XK
      YL = YK
      GO TO 1
C
C NODE K IS TO THE RIGHT OF NL.  FIND THE LEFTMOST NODE
C   N0 WHICH LIES TO THE RIGHT OF K.
C   SET SPROD TO (N0-NL,N0-K).
C
    3 DO 4 IND = 2,KM1
        N0 = IEND(IND)
        SPROD = (XL-X(N0))*(XK-X(N0)) +
     .          (YL-Y(N0))*(YK-Y(N0))
        IF (SPROD .GE. 0.) GO TO 5
    4   CONTINUE
C
C NODE K IS TO THE RIGHT OF NR.  INSERT K AS THE LAST
C   (RIGHTMOST) NODE IN IEND AND SET NR TO K.
C
      IEND(K) = K
      XR = XK
      YR = YK
      GO TO 1
C
C NODE K LIES BETWEEN IEND(IND-1) AND IEND(IND).  INSERT K
C   IN IEND.
C
    5 CALL SHIFTD(IND,KM1,1, IEND )
      IEND(IND) = K
      GO TO 2
C
C NODE K IS TO THE LEFT OF NL-NR.  REORDER IEND SO THAT NL
C   IS THE LEFTMOST NODE AS VIEWED FROM K.
C
    6 KM1D2 = KM1/2
      DO 7 I = 1,KM1D2
        KMI = K-I
        ITEMP = IEND(I)
        IEND(I) = IEND(KMI)
        IEND(KMI) = ITEMP
    7   CONTINUE
C
C NODE K IS TO THE RIGHT OF NL-NR.  CREATE A TRIANGULATION
C   CONSISTING OF NODES 1,2,...,K.
C
    8 NL = IEND(1)
      NR = IEND(KM1)
C
C CREATE THE ADJACENCY LISTS FOR THE FIRST K-1 NODES.
C   INSERT NEIGHBORS IN REVERSE ORDER.  EACH NODE HAS FOUR
C   NEIGHBORS EXCEPT NL AND NR WHICH HAVE THREE.
C
      DO 9 IND = 1,KM1
        N0 = IEND(IND)
        INDX = 4*N0
        IF (N0 .GE. NL) INDX = INDX-1
        IF (N0 .GE. NR) INDX = INDX-1
        IADJ(INDX) = 0
        INDX = INDX-1
        IF (IND .LT. KM1) IADJ(INDX) = IEND(IND+1)
        IF (IND .LT. KM1) INDX = INDX-1
        IADJ(INDX) = K
        IF (IND .EQ. 1) GO TO 9
        IADJ(INDX-1) = IEND(IND-1)
    9   CONTINUE
C
C CREATE THE ADJACENCY LIST FOR NODE K
C
      INDX = 5*KM1 - 1
      IADJ(INDX) = 0
      DO 10 IND = 1,KM1
        INDX = INDX-1
        IADJ(INDX) = IEND(IND)
   10   CONTINUE
C
C REPLACE IEND ELEMENTS WITH POINTERS TO IADJ
C
      INDX = 0
      DO 11 IND = 1,KM1
        INDX = INDX + 4
        IF (IND .EQ. NL  .OR.  IND .EQ. NR) INDX = INDX-1
        IEND(IND) = INDX
   11   CONTINUE
      INDX = INDX + K
      IEND(K) = INDX
C
C ADD THE REMAINING NODES TO THE TRIANGULATION
C
      IF (K .EQ. NN) RETURN
      KMIN = K+1
      DO 12 K = KMIN,NN
        CALL ADNODE(K,X,Y, IADJ,IEND, IERR)
   12   CONTINUE
      RETURN
C
C ALL NODES ARE COLLINEAR
C
   13 IER = 2
      RETURN
      END
      SUBROUTINE TRPLOT (N,X,Y,IADJ,IEND,ITITLE,NC,
     .                   NUMBR, IER)
      INTEGER N, IADJ(1), IEND(N), ITITLE(1), NUMBR, IER
      REAL    X(N), Y(N)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS SUBROUTINE PLOTS THE ARCS OF A TRIANGULATION OF N
C NODES IN THE PLANE.  CARDS WITH C* IN THE FIRST TWO COL-
C UMNS MUST BE REPLACED WITH CALLS TO USER-SUPPLIED GRAPHICS
C SUBROUTINES IN ORDER TO GET ANY USE OUT OF THIS ROUTINE.
C
C INPUT PARAMETERS - N - NUMBER OF NODES IN THE TRIANGULA-
C                        TION.  N .GE. 3.
C
C                  X,Y - CARTESIAN COORDINATES OF THE NODES.
C
C            IADJ,IEND - TRIANGULATION DATA STRUCTURE (SEE
C                        SUBROUTINE TRMESH).
C
C               ITITLE - INTEGER ARRAY CONTAINING A LINE OF
C                        TEXT TO BE CENTERED ABOVE THE PLOT
C                        IF NC .GT. 0.  ITITLE MUST BE INIT-
C                        IALIZED WITH HOLLERITH CONSTANTS OR
C                        READ WITH AN A-FORMAT.  ITS DIMEN-
C                        SION DEPENDS ON NC AND THE NUMBER
C                        OF CHARACTERS STORED IN A COMPUTER
C                        WORD.
C
C                   NC - NUMBER OF CHARACTERS IN ITITLE.
C                        0 .LE. NC .LE. 40.  NO TITLE IS
C                        PLOTTED IF NC = 0.
C
C                NUMBR - OPTION INDICATOR.  IF NUMBR .NE. 0,
C                        THE NODAL INDICES ARE PLOTTTED NEXT
C                        TO THE NODES.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETER - IER - ERROR INDICATOR
C                          IER = 0 IF NO ERRORS WERE ENCOUN-
C                                  TERED.
C                          IER = 1 IF N OR NC IS OUT OF
C                                  RANGE.
C                          IER = 2 IF THE NODES LIE ON A
C                                  HORIZONTAL OR VERTICAL
C                                  LINE.
C
C MODULES REFERENCED BY TRPLOT - NONE
C
C INTRINSIC FUNCTIONS CALLED BY TRPLOT - AMIN1, AMAX1, IABS
C
C***********************************************************
C
      INTEGER NN, I, N0, INDF, INDL, INDX, N1, NEW
      REAL    XMIN, XMAX, YMIN, YMAX, DX, DY, X0, Y0
      LOGICAL ST0
C
C LOCAL PARAMETERS -
C
C NN =        LOCAL COPY OF N
C I =         DO-LOOP INDEX
C N0 =        NODE WHICH IS TO BE CONNECTED TO ITS NEIGHBORS
C               WITH LINE SEGMENTS
C INDF =      IADJ INDEX OF THE FIRST NEIGHBOR OF N0
C INDL =      IADJ INDEX OF THE LAST NEIGHBOR OF N0
C INDX =      IADJ INDEX IN THE RANGE INDF,...,INDL
C N1 =        NEIGHBOR OF N0
C NEW =       NEIGHBOR OF N0 AND CANDIDATE FOR NEXT VALUE
C               OF N0
C XMIN,XMAX = MINIMUM AND MAXIMUM X COORDINATES
C YMIN,YMAX = MINIMUM AND MAXIMUM Y COORDINATES
C DX,DY =     XMAX-XMIN AND YMAX-YMIN (DATA SPACE DIMEN-
C               SIONS)
C X0,Y0 =     COORDINATES OF N0
C ST0 =       SWITCH USED TO ALTERNATE DIRECTION OF PEN
C               MOVEMENT -- TRUE IFF PEN STARTS AT N0
C
      NN = N
C
C CHECK FOR INVALID PARAMETERS
C
      IF (NN .LT. 3  .OR.  NC .LT. 0  .OR.  NC .GT. 40)
     .   GO TO 13
      IER = 0
C
C COMPUTE DATA SPACE DIMENSIONS
C
      XMIN = X(1)
      YMIN = Y(1)
      XMAX = XMIN
      YMAX = YMIN
      DO 1 I = 1,NN
        XMIN = AMIN1(X(I),XMIN)
        YMIN = AMIN1(Y(I),YMIN)
        XMAX = AMAX1(X(I),XMAX)
    1   YMAX = AMAX1(Y(I),YMAX)
      DX = XMAX - XMIN
      DY = YMAX - YMIN
      IF (DX .EQ. 0.  .OR.  DY .EQ. 0.) GO TO 14
C*  COMMANDS WHICH PERFORM THE FOLLOWING FUNCTIONS SHOULD
C*  BE INSERTED HERE --
C*    INITIALIZE THE PLOTTING ENVIRONMENT IF NECESSARY,
C*    COMPUTE PLOTTER SPACE DIMENSIONS,
C*    ESTABLISH A LINEAR MAPPING FROM THE DATA SPACE TO
C*      THE PLOTTER SPACE,
C*    OPTIONALLY, DRAW AND LABEL AXES, AND
C*    PLOT THE TITLE IF NC .NE. 0.
C
C INITIALIZE FOR LOOP ON NODES.  EACH NODE N0 IS TO BE CON-
C   NECTED TO ITS NEIGHBORS WHICH HAVE LARGER INDICES.  N0
C   IS THEN MARKED BY MAKING THE CORRESPONDING IEND ENTRY
C   NEGATIVE, AND THE SEARCH FOR THE NEXT UNMARKED NODE BE-
C   GINS WITH THE NEIGHBORS OF N0.
C
      N0 = 1
      INDF = 1
C
C TOP OF LOOP -- SET INDL, X0, AND Y0, AND INITIALIZE ST0
C   AND INDX
C
    2 INDL = IEND(N0)
      X0 = X(N0)
      Y0 = Y(N0)
      ST0 = .TRUE.
      INDX = INDL
      IF (IADJ(INDX) .EQ. 0) INDX = INDX - 1
C
C   LOOP ON NEIGHBORS OF N0 IN REVERSE ORDER
C
    3 N1 = IADJ(INDX)
      IF (N1 .LT. N0) GO TO 4
C
C   CONNECT N0 AND N1 -- THE DIRECTION OF PEN MOVEMENT
C     ALTERNATES BETWEEN AWAY FROM N0 AND TOWARD N0 FOR
C     REDUCED PEN-UP TIME.
C
      IF (ST0)       print 1111,X0/1e1,X(N1)/1e1,Y0,Y(N1)
      IF (.NOT. ST0) print 1111,X(N1)/1e1,X0/1e1,Y(N1),Y0
      ST0 = .NOT. ST0
C
C   TEST FOR TERMINATION OF LOOP ON NEIGHBORS
C
    4 IF (INDX .EQ. INDF) GO TO 5
      INDX = INDX - 1
      GO TO 3
C
C   MARK N0 AS HAVING BEEN PROCESSED
C
    5 IEND(N0) = -INDL
C
C   SEARCH THE NEIGHBORS OF N0 FOR AN UNMARKED NODE
C     STARTING WITH IADJ(INDX) = IADJ(INDF)
C
    6 NEW = IADJ(INDX)
      IF (NEW .EQ. 0) GO TO 8
      IF (IEND(NEW) .LT. 0) GO TO 7
      N0 = NEW
      INDF = IABS(IEND(N0-1)) + 1
      GO TO 2
C
C   TEST FOR TERMINATION
C
    7 IF (INDX .EQ. INDL) GO TO 8
      INDX = INDX + 1
      GO TO 6
C
C   ALL NEIGHBORS OF N0 ARE MARKED.  SEARCH IEND FOR AN
C     UNMARKED NODE.
C
    8 DO 9 N0 = 2,NN
        IF (IEND(N0) .LT. 0) GO TO 9
        INDF = -IEND(N0-1) + 1
        GO TO 2
    9   CONTINUE
C
C ALL NODES HAVE BEEN MARKED.  RESTORE IEND.
C
      DO 10 N0 = 1,NN
   10   IEND(N0) = -IEND(N0)
      IF (NUMBR .EQ. 0) GO TO 12
C
C PLOT THE NODAL INDICES
C
C*    THE NUMBR OPTION MAY BE IMPLEMENTED BY INSERTING A
C*    LOOP ON THE NODAL COORDINATES WHICH WRITES INDICES
C*    NEXT TO THE NODES.
C
C TERMINATE PLOTTING -- MOVE TO A NEW FRAME.
C
   12 CONTINUE
C*    CALL FRAME
      RETURN
C
C N OR NC IS OUT OF RANGE
C
   13 IER = 1
      RETURN
C
C NODES ARE COLLINEAR
C
   14 IER = 2
      RETURN
 1111 format(' oplot,[',f6.4,',',f6.4,'],[',f6.4,',',f6.4,']')
      END
      SUBROUTINE TRPRNT (N,LUNIT,X,Y,IADJ,IEND,IFLAG)
      INTEGER N, LUNIT, IADJ(1), IEND(N), IFLAG
      REAL    X(N), Y(N)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A TRIANGULATION OF A SET OF POINTS IN THE PLANE,
C THIS ROUTINE PRINTS THE ADJACENCY LISTS AND, OPTIONALLY,
C THE NODAL COORDINATES.  THE NUMBERS OF BOUNDARY NODES,
C TRIANGLES, AND ARCS ARE ALSO PRINTED.
C
C INPUT PARAMETERS -     N - NUMBER OF NODES IN THE MESH.
C                            3 .LE. N .LE. 9999.
C
C                    LUNIT - LOGICAL UNIT FOR OUTPUT.  1
C                            .LE. LUNIT .LE. 99.  OUTPUT IS
C                            PRINTED ON UNIT 6 IF LUNIT IS
C                            OUT OF RANGE.
C
C                      X,Y - VECTORS OF COORDINATES OF THE
C                            NODES IN THE MESH.  NOT USED
C                            UNLESS IFLAG = 0.
C
C                     IADJ - SET OF ADJACENCY LISTS OF NODES
C                            IN THE MESH.
C
C                     IEND - POINTERS TO THE ENDS OF
C                            ADJACENCY LISTS IN IADJ FOR
C                            EACH NODE IN THE MESH.
C
C                    IFLAG - OPTION INDICATOR
C                            IFLAG = 0 IF X AND Y ARE TO BE
C                                      PRINTED (TO 6 DECIMAL
C                                      PLACES).
C                            IFLAG = 1 IF X AND Y ARE NOT
C                                      TO BE PRINTED.
C
C IADJ AND IEND MAY BE CREATED BY TRMESH.
C
C NONE OF THE PARAMETERS ARE ALTERED BY THIS ROUTINE.
C
C MODULES REFERENCED BY TRPRNT - NONE
C
C***********************************************************
C
      INTEGER NN, NMAX, LUN, NODE, INDF, INDL, NL, NLMAX,
     .        INC, I, NB, NT, NA
      DATA    NMAX/9999/, NLMAX/60/
C
C LOCAL PARAMETERS -
C
C NN =        LOCAL COPY OF N
C NMAX =      UPPER BOUND ON N
C LUN =       LOCAL COPY OF LUNIT
C NODE =      INDEX OF A NODE
C INDF,INDL = IADJ INDICES OF THE FIRST AND LAST NEIGHBORS
C               OF NODE
C NL =        NUMBER OF LINES PRINTED ON A PAGE
C NLMAX =     MAXIMUM NUMBER OF PRINT LINES PER PAGE EXCEPT
C               FOR THE LAST PAGE WHICH HAS 3 ADDITIONAL
C               LINES
C INC =       INCREMENT FOR NL
C I =         IADJ INDEX FOR IMPLIED DO-LOOP
C NB =        NUMBER OF BOUNDARY NODES
C NT =        NUMBER OF TRIANGLES
C NA =        NUMBER OF ARCS (UNDIRECTED EDGES)
C
      NN = N
      LUN = LUNIT
      IF (LUN .LT. 1  .OR.  LUN .GT. 99) LUN = 6
C
C PRINT HEADING AND INITIALIZE NL
C
      WRITE (LUN,100) NN
      IF (NN .LT. 3  .OR.  NN .GT. NMAX) GO TO 5
      NL = 6
      IF (IFLAG .EQ. 0) GO TO 2
C
C PRINT IADJ ONLY
C
      WRITE (LUN,101)
      NB = 0
      INDF = 1
      DO 1 NODE = 1,NN
        INDL = IEND(NODE)
        IF (IADJ(INDL) .EQ. 0) NB = NB + 1
        INC = (INDL - INDF)/14 + 2
        NL = NL + INC
        IF (NL .GT. NLMAX) WRITE (LUN,106)
        IF (NL .GT. NLMAX) NL = INC
        WRITE (LUN,103) NODE, (IADJ(I), I = INDF,INDL)
        IF (INDL-INDF .NE. 13) WRITE (LUN,105)
        INDF = INDL + 1
    1   CONTINUE
      GO TO 4
C
C PRINT X, Y, AND IADJ
C
    2 WRITE (LUN,102)
      NB = 0
      INDF = 1
      DO 3 NODE = 1,NN
        INDL = IEND(NODE)
        IF (IADJ(INDL) .EQ. 0) NB = NB + 1
        INC = (INDL - INDF)/8 + 2
        NL = NL + INC
        IF (NL .GT. NLMAX) WRITE (LUN,106)
        IF (NL .GT. NLMAX) NL = INC
        WRITE (LUN,104) NODE, X(NODE), Y(NODE),
     .                  (IADJ(I), I = INDF,INDL)
        IF (INDL-INDF .NE. 7) WRITE (LUN,105)
        INDF = INDL + 1
    3   CONTINUE
C
C PRINT NB, NA, AND NT
C
    4 NT = 2*NN - NB - 2
      NA = NT + NN - 1
      WRITE (LUN,107) NB, NA, NT
      RETURN
C
C N IS OUT OF RANGE
C
    5 WRITE (LUN,108)
      RETURN
C
C PRINT FORMATS
C
  100 FORMAT (1H1,26X,23HADJACENCY SETS,    N = ,I5//)
  101 FORMAT (1H ,4HNODE,32X,17HNEIGHBORS OF NODE//)
  102 FORMAT (1H ,4HNODE,5X,7HX(NODE),8X,7HY(NODE),
     .        20X,17HNEIGHBORS OF NODE//)
  103 FORMAT (1H ,I4,5X,14I5/(1H ,9X,14I5))
  104 FORMAT (1H ,I4,2E15.6,5X,8I5/(1H ,39X,8I5))
  105 FORMAT (1H )
  106 FORMAT (1H1)
  107 FORMAT (/1H ,5HNB = ,I4,15H BOUNDARY NODES,10X,
     .        5HNA = ,I5,5H ARCS,10X,5HNT = ,I5,
     .        10H TRIANGLES)
  108 FORMAT (1H ,10X,25H*** N IS OUT OF RANGE ***)
      END
      SUBROUTINE COORDS (X,Y,X1,X2,X3,Y1,Y2,Y3, R,IER)
      INTEGER IER
      REAL    X, Y, X1, X2, X3, Y1, Y2, Y3, R(3)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE COMPUTES THE THREE BARYCENTRIC COORDINATES
C OF A POINT IN THE PLANE FOR A GIVEN TRIANGLE.
C
C INPUT PARAMETERS - X,Y - X AND Y COORDINATES OF THE POINT
C                          WHOSE BARYCENTRIC COORDINATES ARE
C                          DESIRED.
C
C      X1,X2,X3,Y1,Y2,Y3 - COORDINATES OF THE VERTICES OF
C                          THE TRIANGLE.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETERS -  R - 3-VECTOR OF BARYCENTRIC COORDI-
C                          NATES UNLESS IER = 1.  NOTE THAT
C                          R(I) .LT. 0. IFF (X,Y) IS TO THE
C                          RIGHT OF THE VECTOR FROM VERTEX
C                          I+1 TO VERTEX I+2 (CYCLICAL
C                          ARITHMETIC).
C
C                    IER - ERROR INDICATOR
C                          IER = 0 IF NO ERRORS WERE
C                                  ENCOUNTERED.
C                          IER = 1 IF THE VERTICES OF THE
C                                  TRIANGLE ARE COLLINEAR.
C
C MODULES REFERENCED BY COORDS - NONE
C
C***********************************************************
C
      REAL    U(3), V(3), AREA, XP, YP
C
C LOCAL PARAMETERS -
C
C U(K),V(K) = X AND Y COMPONENTS OF THE VECTOR REPRESENTING
C               THE SIDE OPPOSITE VERTEX K FOR K = 1,2,3.
C AREA =      TWICE THE AREA OF THE TRIANGLE.
C XP,YP =     X-X1, Y-Y1
C
      U(1) = X3-X2
      U(2) = X1-X3
      U(3) = X2-X1
C
      V(1) = Y3-Y2
      V(2) = Y1-Y3
      V(3) = Y2-Y1
C
C AREA = 3-1 X 3-2
C
      AREA = U(1)*V(2) - U(2)*V(1)
      IF (AREA .EQ. 0.) GO TO 1
C
C R(1) = (2-3 X 2-(X,Y))/AREA, R(2) = (1-(X,Y) X 1-3)/AREA,
C   R(3) = (1-2 X 1-(X,Y))/AREA
C
      R(1) = (U(1)*(Y-Y2) - V(1)*(X-X2))/AREA
      XP = X - X1
      YP = Y - Y1
      R(2) = (U(2)*YP - V(2)*XP)/AREA
      R(3) = (U(3)*YP - V(3)*XP)/AREA
      IER = 0
      RETURN
C
C VERTICES ARE COLLINEAR
C
    1 IER = 1
      RETURN
      END
      SUBROUTINE GIVENS ( A,B, C,S)
      REAL A, B, C, S
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE CONSTRUCTS THE GIVENS PLANE ROTATION --
C     ( C  S)
C G = (     ) WHERE C*C + S*S = 1 -- WHICH ZEROS THE SECOND
C     (-S  C)
C ENTRY OF THE 2-VECTOR (A B)-TRANSPOSE.  A CALL TO GIVENS
C IS NORMALLY FOLLOWED BY A CALL TO ROTATE WHICH APPLIES
C THE TRANSFORMATION TO A 2 BY N MATRIX.  THIS ROUTINE WAS
C TAKEN FROM LINPACK.
C
C INPUT PARAMETERS - A,B - COMPONENTS OF THE 2-VECTOR TO BE
C                          ROTATED.
C
C OUTPUT PARAMETERS -  A - OVERWRITTEN BY R = +/-SQRT(A*A
C                          + B*B)
C
C                      B - OVERWRITTEN BY A VALUE Z WHICH
C                          ALLOWS C AND S TO BE RECOVERED
C                          AS FOLLOWS -
C                          C = SQRT(1-Z*Z), S=Z IF ABS(Z)
C                              .LE. 1.
C                          C = 1/Z, S = SQRT(1-C*C) IF
C                              ABS(Z) .GT. 1.
C
C                      C - +/-(A/R)
C
C                      S - +/-(B/R)
C
C MODULES REFERENCED BY GIVENS - NONE
C
C INTRINSIC FUNCTIONS CALLED BY GIVENS - ABS, SQRT
C
C***********************************************************
C
      REAL AA, BB, R, U, V
C
C LOCAL PARAMETERS -
C
C AA,BB = LOCAL COPIES OF A AND B
C R =     C*A + S*B = +/-SQRT(A*A+B*B)
C U,V =   VARIABLES USED TO SCALE A AND B FOR COMPUTING R
C
      AA = A
      BB = B
      IF (ABS(AA) .LE. ABS(BB)) GO TO 1
C
C ABS(A) .GT. ABS(B)
C
      U = AA + AA
      V = BB/U
      R = SQRT(.25 + V*V) * U
      C = AA/R
      S = V * (C + C)
C
C NOTE THAT R HAS THE SIGN OF A, C .GT. 0, AND S HAS
C   SIGN(A)*SIGN(B)
C
      B = S
      A = R
      RETURN
C
C ABS(A) .LE. ABS(B)
C
    1 IF (BB .EQ. 0.) GO TO 2
      U = BB + BB
      V = AA/U
C
C STORE R IN A
C
      A = SQRT(.25 + V*V) * U
      S = BB/A
      C = V * (S + S)
C
C NOTE THAT R HAS THE SIGN OF B, S .GT. 0, AND C HAS
C   SIGN(A)*SIGN(B)
C
      B = 1.
      IF (C .NE. 0.) B = 1./C
      RETURN
C
C A = B = 0.
C
    2 C = 1.
      S = 0.
      RETURN
      END
      SUBROUTINE GRADG (N,X,Y,Z,IADJ,IEND,EPS, NIT,
     .                  ZXZY, IER)
      INTEGER N, IADJ(1), IEND(N), NIT, IER
      REAL    X(N), Y(N), Z(N), EPS, ZXZY(2,N)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A TRIANGULATION OF N NODES IN THE PLANE WITH
C ASSOCIATED DATA VALUES, THIS ROUTINE USES A GLOBAL METHOD
C TO COMPUTE ESTIMATED GRADIENTS AT THE NODES.  THE METHOD
C CONSISTS OF MINIMIZING A QUADRATIC FUNCTIONAL Q(G) OVER
C THE N-VECTOR G OF GRADIENTS WHERE Q APPROXIMATES THE LIN-
C EARIZED CURVATURE OF AN INTERPOLANT F OVER THE TRIANGULA-
C TION.  THE RESTRICTION OF F TO AN ARC OF THE TRIANGULATION
C IS TAKEN TO BE THE HERMITE CUBIC INTERPOLANT OF THE DATA
C VALUES AND TANGENTIAL GRADIENT COMPONENTS AT THE END-
C POINTS OF THE ARC, AND Q IS THE SUM OF THE LINEARIZED
C CURVATURES OF F ALONG THE ARCS -- THE INTEGRALS OVER THE
C ARCS OF D2F(T)**2 WHERE D2F(T) IS THE SECOND DERIVATIVE
C OF F WITH RESPECT TO DISTANCE T ALONG THE ARC.  THIS MIN-
C IMIZATION PROBLEM CORRESPONDS TO AN ORDER 2N SYMMETRIC
C POSITIVE-DEFINITE SPARSE LINEAR SYSTEM WHICH IS SOLVED FOR
C THE X AND Y PARTIAL DERIVATIVES BY THE BLOCK GAUSS-SEIDEL
C METHOD WITH 2 BY 2 BLOCKS.
C   AN ALTERNATIVE METHOD, SUBROUTINE GRADL, COMPUTES A
C LOCAL APPROXIMATION TO THE PARTIALS AT A SINGLE NODE AND
C MAY BE MORE ACCURATE, DEPENDING ON THE DATA VALUES AND
C DISTRIBUTION OF NODES (NEITHER METHOD EMERGED AS SUPERIOR
C IN TESTS FOR ACCURACY).  HOWEVER, IN TESTS RUN ON AN IBM
C 370, GRADG WAS FOUND TO BE ABOUT 3.6 TIMES AS FAST FOR
C NIT = 4.
C
C INPUT PARAMETERS - N - NUMBER OF NODES.  N .GE. 3.
C
C                  X,Y - CARTESIAN COORDINATES OF THE NODES.
C
C                    Z - DATA VALUES AT THE NODES.  Z(I) IS
C                        ASSOCIATED WITH (X(I),Y(I)).
C
C            IADJ,IEND - DATA STRUCTURE DEFINING THE TRIAN-
C                        GULATION.  SEE SUBROUTINE TRMESH.
C
C                  EPS - NONNEGATIVE CONVERGENCE CRITERION.
C                        THE METHOD IS TERMINATED WHEN THE
C                        MAXIMUM CHANGE IN A GRADIENT COMPO-
C                        NENT BETWEEN ITERATIONS IS AT MOST
C                        EPS.  EPS = 1.E-2 IS SUFFICIENT FOR
C                        EFFECTIVE CONVERGENCE.
C
C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C                  NIT - MAXIMUM NUMBER OF GAUSS-SEIDEL
C                        ITERATIONS TO BE APPLIED.  THIS
C                        MAXIMUM WILL LIKELY BE ACHIEVED IF
C                        EPS IS SMALLER THAN THE MACHINE
C                        PRECISION.  OPTIMAL EFFICIENCY WAS
C                        ACHIEVED IN TESTING WITH EPS = 0
C                        AND NIT = 3 OR 4.
C
C                 ZXZY - 2 BY N ARRAY WHOSE COLUMNS CONTAIN
C                        INITIAL ESTIMATES OF THE PARTIAL
C                        DERIVATIVES (ZERO VECTORS ARE
C                        SUFFICIENT).
C
C OUTPUT PARAMETERS - NIT - NUMBER OF GAUSS-SEIDEL ITERA-
C                           TIONS EMPLOYED.
C
C                    ZXZY - ESTIMATED X AND Y PARTIAL DERIV-
C                           ATIVES AT THE NODES WITH X PAR-
C                           TIALS IN THE FIRST ROW.  ZXZY IS
C                           NOT CHANGED IF IER = 2.
C
C                     IER - ERROR INDICATOR
C                           IER = 0 IF THE CONVERGENCE CRI-
C                                   TERION WAS ACHIEVED.
C                           IER = 1 IF CONVERGENCE WAS NOT
C                                   ACHIEVED WITHIN NIT
C                                   ITERATIONS.
C                           IER = 2 IF N OR EPS IS OUT OF
C                                   RANGE OR NIT .LT. 0 ON
C                                   INPUT.
C
C MODULES REFERENCED BY GRADG - NONE
C
C INTRINSIC FUNCTIONS CALLED BY GRADG - SQRT, AMAX1, ABS
C
C***********************************************************
C
      INTEGER NN, MAXIT, ITER, K, INDF, INDL, INDX, NB
      REAL    TOL, DGMAX, XK, YK, ZK, ZXK, ZYK, A11, A12,
     .        A22, R1, R2, DELX, DELY, DELXS, DELYS, DSQ,
     .        DCUB, T, DZX, DZY
C
C LOCAL PARAMETERS -
C
C NN =          LOCAL COPY OF N
C MAXIT =       INPUT VALUE OF NIT
C ITER =        NUMBER OF ITERATIONS USED
C K =           DO-LOOP AND NODE INDEX
C INDF,INDL =   IADJ INDICES OF THE FIRST AND LAST NEIGHBORS
C                 OF K
C INDX =        IADJ INDEX IN THE RANGE INDF,...,INDL
C NB =          NEIGHBOR OF K
C TOL =         LOCAL COPY OF EPS
C DGMAX =       MAXIMUM CHANGE IN A GRADIENT COMPONENT BE-
C                 TWEEN ITERATIONS
C XK,YK,ZK =    X(K), Y(K), Z(K)
C ZXK,ZYK =     INITIAL VALUES OF ZXZY(1,K) AND ZXZY(2,K)
C A11,A12,A22 = MATRIX COMPONENTS OF THE 2 BY 2 BLOCK A*DG
C                 = R WHERE A IS SYMMETRIC, DG = (DZX,DZY)
C                 IS THE CHANGE IN THE GRADIENT AT K, AND R
C                 IS THE RESIDUAL
C R1,R2 =       COMPONENTS OF THE RESIDUAL -- DERIVATIVES OF
C                 Q WITH RESPECT TO THE COMPONENTS OF THE
C                 GRADIENT AT NODE K
C DELX,DELY =   COMPONENTS OF THE ARC NB-K
C DELXS,DELYS = DELX**2, DELY**2
C DSQ =         SQUARE OF THE DISTANCE D BETWEEN K AND NB
C DCUB =        D**3
C T =           FACTOR OF R1 AND R2
C DZX,DZY =     SOLUTION OF THE 2 BY 2 SYSTEM -- CHANGE IN
C                 DERIVATIVES AT K FROM THE PREVIOUS ITERATE
C
      NN = N
      TOL = EPS
      MAXIT = NIT
C
C ERROR CHECKS AND INITIALIZATION
C
      IF (NN .LT. 3  .OR.  TOL .LT. 0.  .OR.  MAXIT .LT. 0)
     .   GO TO 5
      ITER = 0
C
C TOP OF ITERATION LOOP
C
    1 IF (ITER .EQ. MAXIT) GO TO 4
      DGMAX = 0.
      INDL = 0
      DO 3 K = 1,NN
        XK = X(K)
        YK = Y(K)
        ZK = Z(K)
        ZXK = ZXZY(1,K)
        ZYK = ZXZY(2,K)
C
C   INITIALIZE COMPONENTS OF THE 2 BY 2 SYSTEM
C
        A11 = 0.
        A12 = 0.
        A22 = 0.
        R1 = 0.
        R2 = 0.
C
C   LOOP ON NEIGHBORS NB OF K
C
        INDF = INDL + 1
        INDL = IEND(K)
        DO 2 INDX = INDF,INDL
          NB = IADJ(INDX)
          IF (NB .EQ. 0) GO TO 2
C
C   COMPUTE THE COMPONENTS OF ARC NB-K
C
          DELX = X(NB) - XK
          DELY = Y(NB) - YK
          DELXS = DELX*DELX
          DELYS = DELY*DELY
          DSQ = DELXS + DELYS
          DCUB = DSQ*SQRT(DSQ)
C
C   UPDATE THE SYSTEM COMPONENTS FOR NODE NB
C
          A11 = A11 + DELXS/DCUB
          A12 = A12 + DELX*DELY/DCUB
          A22 = A22 + DELYS/DCUB
          T = ( 1.5*(Z(NB)-ZK) - ((ZXZY(1,NB)/2.+ZXK)*DELX +
     .          (ZXZY(2,NB)/2.+ZYK)*DELY) )/DCUB
          R1 = R1 + T*DELX
          R2 = R2 + T*DELY
    2     CONTINUE
C
C   SOLVE THE 2 BY 2 SYSTEM AND UPDATE DGMAX
C
        DZY = (A11*R2 - A12*R1)/(A11*A22 - A12*A12)
        DZX = (R1 - A12*DZY)/A11
        DGMAX = AMAX1(DGMAX,ABS(DZX),ABS(DZY))
C
C   UPDATE THE PARTIALS AT NODE K
C
        ZXZY(1,K) = ZXK + DZX
    3   ZXZY(2,K) = ZYK + DZY
C
C   INCREMENT ITER AND TEST FOR CONVERGENCE
C
      ITER = ITER + 1
      IF (DGMAX .GT. TOL) GO TO 1
C
C METHOD CONVERGED
C
      NIT = ITER
      IER = 0
      RETURN
C
C METHOD FAILED TO CONVERGE WITHIN NIT ITERATIONS
C
    4 IER = 1
      RETURN
C
C PARAMETER OUT OF RANGE
C
    5 NIT = 0
      IER = 2
      RETURN
      END
      SUBROUTINE GRADL (N,K,X,Y,Z,IADJ,IEND, DX,DY,IER)
      INTEGER N, K, IADJ(1), IEND(N), IER
      REAL    X(N), Y(N), Z(N), DX, DY
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A THIESSEN TRIANGULATION OF N POINTS IN THE PLANE
C WITH ASSOCIATED DATA VALUES Z, THIS SUBROUTINE ESTIMATES
C X AND Y PARTIAL DERIVATIVES AT NODE K.  THE DERIVATIVES
C ARE TAKEN TO BE THE PARTIALS AT K OF A QUADRATIC FUNCTION
C WHICH INTERPOLATES Z(K) AND FITS THE DATA VALUES AT A SET
C OF NEARBY NODES IN A WEIGHTED LEAST SQUARES SENSE. A MAR-
C QUARDT STABILIZATION FACTOR IS USED IF NECESSARY TO ENSURE
C A WELL-CONDITIONED SYSTEM AND A LINEAR FITTING FUNCTION IS
C USED IF N .LT. 6.  THUS, A UNIQUE SOLUTION EXISTS UNLESS
C THE NODES ARE COLLINEAR.
C   AN ALTERNATIVE ROUTINE, GRADG, EMPLOYS A GLOBAL METHOD
C TO COMPUTE THE PARTIAL DERIVATIVES AT ALL OF THE NODES AT
C ONCE.  THAT METHOD IS MORE EFFICIENT (WHEN ALL PARTIALS
C ARE NEEDED) AND MAY BE MORE ACCURATE, DEPENDING ON THE
C DATA.
C
C INPUT PARAMETERS - N - NUMBER OF NODES IN THE TRIANGULA-
C                        TION.  N .GE. 3.
C
C                    K - NODE AT WHICH DERIVATIVES ARE
C                        SOUGHT.  1 .LE. K .LE. N.
C
C                  X,Y - N-VECTORS CONTAINING THE CARTESIAN
C                        COORDINATES OF THE NODES.
C
C                    Z - N-VECTOR CONTAINING THE DATA VALUES
C                        ASSOCIATED WITH THE NODES.
C
C                 IADJ - SET OF ADJACENCY LISTS.
C
C                 IEND - POINTERS TO THE ENDS OF ADJACENCY
C                        LISTS FOR EACH NODE.
C
C IADJ AND IEND MAY BE CREATED BY SUBROUTINE TRMESH.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETERS - DX,DY - ESTIMATED PARTIAL DERIVATIVES
C                             AT NODE K UNLESS IER .LT. 0.
C
C                       IER - ERROR INDICATOR
C                             IER .GT. 0 IF NO ERRORS WERE
C                                      ENCOUNTERED.  IER
C                                      CONTAINS THE NUMBER
C                                      OF NODES (INCLUDING
C                                      K) USED IN THE FIT.
C                                      IER = 3, 4, OR 5 IM-
C                                      PLIES A LINEAR FIT.
C                             IER = -1 IF N OR K IS OUT OF
C                                      RANGE.
C                             IER = -2 IF ALL NODES ARE
C                                      COLLINEAR.
C
C MODULES REFERENCED BY GRADL - GETNP, SETUP, GIVENS,
C                               ROTATE
C
C INTRINSIC FUNCTIONS CALLED BY GRADL - MIN0, FLOAT, SQRT,
C                                       AMIN1, ABS
C
C***********************************************************
C
      INTEGER NN, KK, LMN, LMX, LMIN, LMAX, LM1, LNP,
     .        NPTS(30), IERR, NP, I, J, IM1, JP1, IP1, L
      REAL    SUM, DS, R, RS, RTOL, AVSQ, AV, XK, YK, ZK,
     .        A(6,6), C, S, DMIN, DTOL, SF
      DATA    LMN/10/
      DATA    LMX/30/, RTOL/1.E-5/, DTOL/.01/, SF/1./
C
C LOCAL PARAMETERS -
C
C NN,KK =     LOCAL COPIES OF N AND K
C LMN,LMX =   MINIMUM AND MAXIMUM VALUES OF LNP FOR N
C               SUFFICIENTLY LARGE.  IN MOST CASES LMN-1
C               NODES ARE USED IN THE FIT.  4 .LE. LMN .LE.
C               LMX.
C LMIN,LMAX = MIN(LMN,N), MIN(LMX,N)
C LM1 =       LMIN-1 OR LNP-1
C LNP =       LENGTH OF NPTS
C NPTS =      ARRAY CONTAINING THE INDICES OF A SEQUENCE OF
C               NODES ORDERED BY DISTANCE FROM K.  NPTS(1)=K
C               AND THE FIRST LNP-1 ELEMENTS OF NPTS ARE
C               USED IN THE LEAST SQUARES FIT.  UNLESS LNP
C               EXCEEDS LMAX, NPTS(LNP) DETERMINES R.
C IERR =      ERROR FLAG FOR CALLS TO GETNP (NOT CHECKED)
C NP =        ELEMENT OF NPTS TO BE ADDED TO THE SYSTEM
C I,J =       DO-LOOP INDICES
C IM1,JP1 =   I-1, J+1
C IP1 =       I+1
C L =         NUMBER OF COLUMNS OF A**T TO WHICH A ROTATION
C               IS APPLIED
C SUM =       SUM OF SQUARED EUCLIDEAN DISTANCES BETWEEN
C               NODE K AND THE NODES USED IN THE LEAST
C               SQUARES FIT
C DS =        SQUARED DISTANCE BETWEEN NODE K AND AN ELE-
C               MENT OF NPTS
C R =         DISTANCE BETWEEN NODE K AND NPTS(LNP) OR SOME
C               POINT FURTHER FROM K THAN NPTS(LMAX) IF
C               NPTS(LMAX) IS USED IN THE FIT.  R IS A
C               RADIUS OF INFLUENCE WHICH ENTERS INTO THE
C               WEIGHTS (SEE SUBROUTINE SETUP).
C RS =        R*R
C RTOL =      TOLERANCE FOR DETERMINING R.  IF THE RELATIVE
C               CHANGE IN DS BETWEEN TWO ELEMENTS OF NPTS IS
C               NOT GREATER THAN RTOL THEY ARE TREATED AS
C               BEING THE SAME DISTANCE FROM NODE K
C AVSQ =      AV*AV
C AV =        ROOT-MEAN-SQUARE DISTANCE BETWEEN K AND THE
C               NODES (OTHER THAN K) IN THE LEAST SQUARES
C               FIT.  THE FIRST 3 COLUMNS OF THE SYSTEM ARE
C               SCALED BY 1/AVSQ, THE NEXT 2 BY 1/AV.
C XK,YK,ZK =  COORDINATES AND DATA VALUE ASSOCIATED WITH K
C A =         TRANSPOSE OF THE AUGMENTED REGRESSION MATRIX
C C,S =       COMPONENTS OF THE PLANE ROTATION DETERMINED
C               BY SUBROUTINE GIVENS
C DMIN =      MINIMUM OF THE MAGNITUDES OF THE DIAGONAL
C               ELEMENTS OF THE REGRESSION MATRIX AFTER
C               ZEROS ARE INTRODUCED BELOW THE DIAGONAL
C DTOL =      TOLERANCE FOR DETECTING AN ILL-CONDITIONED
C               SYSTEM.  THE SYSTEM IS ACCEPTED WHEN DMIN
C               .GE. DTOL
C SF =        MARQUARDT STABILIZATION FACTOR USED TO DAMP
C               OUT THE FIRST 3 SOLUTION COMPONENTS (SECOND
C               PARTIALS OF THE QUADRATIC) WHEN THE SYSTEM
C               IS ILL-CONDITIONED.  AS SF INCREASES, THE
C               FITTING FUNCTION APPROACHES A LINEAR
C
      NN = N
      KK = K
C
C CHECK FOR ERRORS AND INITIALIZE LMIN, LMAX
C
      IF (NN .LT. 3  .OR.  KK .LT. 1  .OR.  KK .GT. NN)
     .   GO TO 16
      LMIN = MIN0(LMN,NN)
      LMAX = MIN0(LMX,NN)
C
C COMPUTE NPTS, LNP, AVSQ, AV, AND R.
C   SET NPTS TO THE CLOSEST LMIN-1 NODES TO K.
C
      SUM = 0.
      NPTS(1) = KK
      LM1 = LMIN - 1
      DO 1 LNP = 2,LM1
        CALL GETNP (X,Y,IADJ,IEND,LNP, NPTS, DS,IERR)
    1   SUM = SUM + DS
C
C ADD ADDITIONAL NODES TO NPTS UNTIL THE RELATIVE INCREASE
C   IN DS IS AT LEAST RTOL.
C
      DO 2 LNP = LMIN,LMAX
        CALL GETNP (X,Y,IADJ,IEND,LNP, NPTS, RS,IERR)
        IF ((RS-DS)/DS .LE. RTOL) GO TO 2
        IF (LNP .GT. 6) GO TO 3
    2   SUM = SUM + RS
C
C USE ALL LMAX NODES IN THE LEAST SQUARES FIT.  RS IS
C   ARBITRARILY INCREASED BY 10 PER CENT.
C
      RS = 1.1*RS
      LNP = LMAX + 1
C
C THERE ARE LNP-2 EQUATIONS CORRESPONDING TO NODES NPTS(2),
C   ...,NPTS(LNP-1).
C
    3 AVSQ = SUM/FLOAT(LNP-2)
      AV = SQRT(AVSQ)
      R = SQRT(RS)
      XK = X(KK)
      YK = Y(KK)
      ZK = Z(KK)
      IF (LNP .LT. 7) GO TO 12
C
C SET UP THE FIRST 5 EQUATIONS OF THE AUGMENTED REGRESSION
C   MATRIX (TRANSPOSED) AS THE COLUMNS OF A, AND ZERO OUT
C   THE LOWER TRIANGLE (UPPER TRIANGLE OF A) WITH GIVENS
C   ROTATIONS
C
      DO 5 I = 1,5
        NP = NPTS(I+1)
        CALL SETUP (XK,YK,ZK,X(NP),Y(NP),Z(NP),AV,AVSQ,
     .              R, A(1,I))
        IF (I .EQ. 1) GO TO 5
        IM1 = I - 1
        DO 4 J = 1,IM1
          JP1 = J + 1
          L = 6 - J
          CALL GIVENS (A(J,J),A(J,I),C,S)
    4     CALL ROTATE (L,C,S,A(JP1,J),A(JP1,I))
    5   CONTINUE
C
C ADD THE ADDITIONAL EQUATIONS TO THE SYSTEM USING
C   THE LAST COLUMN OF A -- I .LE. LNP.
C
      I = 7
    6   IF (I .EQ. LNP) GO TO 8
        NP = NPTS(I)
        CALL SETUP (XK,YK,ZK,X(NP),Y(NP),Z(NP),AV,AVSQ,
     .              R, A(1,6))
        DO 7 J = 1,5
          JP1 = J + 1
          L = 6 - J
          CALL GIVENS (A(J,J),A(J,6),C,S)
    7     CALL ROTATE (L,C,S,A(JP1,J),A(JP1,6))
        I = I + 1
        GO TO 6
C
C TEST THE SYSTEM FOR ILL-CONDITIONING
C
    8 DMIN = AMIN1( ABS(A(1,1)),ABS(A(2,2)),ABS(A(3,3)),
     .              ABS(A(4,4)),ABS(A(5,5)) )
      IF (DMIN .GE. DTOL) GO TO 15
      IF (LNP .GT. LMAX) GO TO 9
C
C ADD ANOTHER NODE TO THE SYSTEM AND INCREASE R --
C   I .EQ. LNP
C
      LNP = LNP + 1
      IF (LNP .LE. LMAX) CALL GETNP (X,Y,IADJ,IEND,LNP,
     .                               NPTS, RS,IERR)
      R = SQRT(1.1*RS)
      GO TO 6
C
C STABILIZE THE SYSTEM BY DAMPING SECOND PARTIALS --ADD
C   MULTIPLES OF THE FIRST THREE UNIT VECTORS TO THE FIRST
C   THREE EQUATIONS.
C
    9 DO 11 I = 1,3
        A(I,6) = SF
        IP1 = I + 1
        DO 10 J = IP1,6
   10     A(J,6) = 0.
        DO 11 J = I,5
          JP1 = J + 1
          L = 6 - J
          CALL GIVENS (A(J,J),A(J,6),C,S)
   11     CALL ROTATE (L,C,S,A(JP1,J),A(JP1,6))
      GO TO 14
C
C 4 .LE. LNP .LE. 6 (2, 3, OR 4 EQUATIONS) -- FIT A PLANE TO
C   THE DATA USING THE LAST 3 COLUMNS OF A.
C
   12 NP = NPTS(2)
      CALL SETUP (XK,YK,ZK,X(NP),Y(NP),Z(NP),AV,AVSQ,
     .            R, A(1,4))
      NP = NPTS(3)
      CALL SETUP (XK,YK,ZK,X(NP),Y(NP),Z(NP),AV,AVSQ,
     .            R, A(1,5))
      CALL GIVENS (A(4,4),A(4,5),C,S)
      CALL ROTATE (2,C,S,A(5,4),A(5,5))
      IF (LNP .EQ. 4) GO TO 14
C
      LM1 = LNP - 1
      DO 13 I = 4,LM1
        NP = NPTS(I)
        CALL SETUP (XK,YK,ZK,X(NP),Y(NP),Z(NP),AV,AVSQ,
     .              R, A(1,6))
        CALL GIVENS (A(4,4),A(4,6),C,S)
        CALL ROTATE (2,C,S,A(5,4),A(5,6))
        CALL GIVENS (A(5,5),A(5,6),C,S)
   13   CALL ROTATE (1,C,S,A(6,5),A(6,6))
C
C TEST THE LINEAR FIT FOR ILL-CONDITIONING
C
   14 DMIN = AMIN1( ABS(A(4,4)),ABS(A(5,5)) )
      IF (DMIN .LT. DTOL) GO TO 17
C
C SOLVE THE 2 BY 2 TRIANGULAR SYSTEM FOR THE DERIVATIVES
C
   15 DY = A(6,5)/A(5,5)
      DX = (A(6,4) - A(5,4)*DY)/A(4,4)/AV
      DY = DY/AV
      IER = LNP - 1
      RETURN
C
C N OR K IS OUT OF RANGE
C
   16 IER = -1
      RETURN
C
C NO UNIQUE SOLUTION DUE TO COLLINEAR NODES
C
   17 IER = -2
      RETURN
      END
      SUBROUTINE INTRC0 (N,PX,PY,X,Y,Z,IADJ,IEND, IST, PZ,
     .                   IER)
      INTEGER N, IADJ(1), IEND(N), IST, IER
      REAL    PX, PY, X(N), Y(N), Z(N), PZ
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A TRIANGULATION OF A SET OF POINTS IN THE PLANE,
C THIS ROUTINE COMPUTES THE VALUE AT (PX,PY) OF A PIECEWISE
C LINEAR SURFACE WHICH INTERPOLATES DATA VALUES AT THE
C VERTICES OF THE TRIANGLES.  THE SURFACE IS EXTENDED IN A
C CONTINUOUS FASHION BEYOND THE BOUNDARY OF THE TRIANGULAR
C MESH, ALLOWING EXTRAPOLATION.  INTRC0 IS PART OF AN
C INTERPOLATION PACKAGE WHICH PROVIDES ROUTINES TO GENERATE,
C UPDATE, AND PLOT THE MESH.
C
C INPUT PARAMETERS -     N - NUMBER OF NODES IN THE MESH.
C                            N .GE. 3.
C
C                    PX,PY - POINT AT WHICH THE INTERPOLATED
C                            VALUE IS DESIRED.
C
C                      X,Y - VECTORS OF COORDINATES OF THE
C                            NODES IN THE MESH.
C
C                        Z - VECTOR OF DATA VALUES AT THE
C                            NODES.
C
C                     IADJ - SET OF ADJACENCY LISTS OF NODES
C                            IN THE MESH.
C
C                     IEND - POINTERS TO THE ENDS OF
C                            ADJACENCY LISTS IN IADJ FOR
C                            EACH NODE IN THE MESH.
C
C                      IST - INDEX OF THE STARTING NODE IN
C                            THE SEARCH FOR A TRIANGLE CON-
C                            TAINING (PX,PY).  1 .LE. IST
C                            .LE. N.  THE OUTPUT VALUE OF
C                            IST FROM A PREVIOUS CALL MAY
C                            BE A GOOD CHOICE.
C
C IADJ AND IEND MAY BE CREATED BY TRMESH.
C
C INPUT PARAMETERS OTHER THAN IST ARE NOT ALTERED BY THIS
C   ROUTINE.
C
C OUTPUT PARAMETERS -  IST - INDEX OF ONE OF THE VERTICES OF
C                            THE TRIANGLE CONTAINING (PX,PY)
C                            UNLESS IER .LT. 0.
C
C                       PZ - VALUE OF THE INTERPOLATORY
C                            SURFACE AT (PX,PY) OR ZERO
C                            IF IER .LT. 0.
C
C                      IER - ERROR INDICATOR
C                            IER = 0 IF NO ERRORS WERE
C                                    ENCOUNTERED.
C                            IER = 1 IF NO ERRORS WERE EN-
C                                    COUNTERED AND EXTRAPO-
C                                    LATION WAS PERFORMED.
C                            IER = -1 IF N OR IST IS OUT OF
C                                     RANGE.
C                            IER = -2 IF THE NODES ARE COL-
C                                     LINEAR.
C
C MODULES REFERENCED BY INTRC0 - TRFIND, COORDS
C
C***********************************************************
C
      INTEGER I1, I2, I3, N1, N2, INDX
      REAL    XP, YP, R(3), X1, Y1, X2, Y2, DP
C
C LOCAL PARAMETERS -
C
C I1,I2,I3 = VERTEX INDICES RETURNED BY TRFIND
C N1,N2 =    ENDPOINTS OF THE CLOSEST BOUNDARY EDGE TO P
C              WHEN P IS OUTSIDE OF THE MESH BOUNDARY
C INDX =     IADJ INDEX OF N1 AS A NEIGHBOR OF N2
C XP,YP =    LOCAL COPIES OF THE COORDINATES OF P=(PX,PY)
C R =        BARYCENTRIC COORDINATES
C X1,Y1 =    X,Y COORDINATES OF N1
C X2,Y2 =    X,Y COORDINATES OF N2
C DP =       INNER PRODUCT OF N1-N2 AND P-N2
C
      IF (N .LT. 3  .OR.  IST .LT. 1  .OR.  IST .GT. N)
     .   GO TO 5
      XP = PX
      YP = PY
C
C FIND A TRIANGLE CONTAINING P IF P IS WITHIN THE MESH
C   BOUNDARY
C
      CALL TRFIND(IST,XP,YP,X,Y,IADJ,IEND, I1,I2,I3)
      IF (I1 .EQ. 0) GO TO 6
      IST = I1
      IF (I3 .EQ. 0) GO TO 1
C
C COMPUTE BARYCENTRIC COORDINATES
C
      CALL COORDS(XP,YP,X(I1),X(I2),X(I3),Y(I1),Y(I2),
     .            Y(I3), R,IER)
      IF (IER .NE. 0) GO TO 6
      PZ = R(1)*Z(I1) + R(2)*Z(I2) + R(3)*Z(I3)
      RETURN
C
C P IS OUTSIDE OF THE MESH BOUNDARY.  EXTRAPOLATE TO P BY
C   EXTENDING THE INTERPOLATORY SURFACE AS A CONSTANT
C   BEYOND THE BOUNDARY.  THUS PZ IS THE SURFACE FUNCTION
C   VALUE AT Q WHERE Q IS THE CLOSEST BOUNDARY POINT TO P.
C
C DETERMINE Q BY TRAVERSING THE BOUNDARY STARTING FROM THE
C   RIGHTMOST VISIBLE NODE I1.
C
    1 N2 = I1
C
C SET N1 TO THE LAST NONZERO NEIGHBOR OF N2 AND COMPUTE DP
C
    2 INDX = IEND(N2) - 1
      N1 = IADJ(INDX)
      X1 = X(N1)
      Y1 = Y(N1)
      X2 = X(N2)
      Y2 = Y(N2)
      DP = (X1-X2)*(XP-X2) + (Y1-Y2)*(YP-Y2)
      IF (DP .LE. 0.) GO TO 3
      IF ((XP-X1)*(X2-X1) + (YP-Y1)*(Y2-Y1) .GT. 0.) GO TO 4
      N2 = N1
      GO TO 2
C
C N2 IS THE CLOSEST BOUNDARY POINT TO P
C
    3 PZ = Z(N2)
      IER = 1
      RETURN
C
C THE CLOSEST BOUNDARY POINT TO P LIES ON N2-N1.  COMPUTE
C   ITS COORDINATES WITH RESPECT TO N2-N1.
C
    4 R(1) = DP/( (X2-X1)**2 + (Y2-Y1)**2 )
      R(2) = 1. - R(1)
      PZ = R(1)*Z(N1) + R(2)*Z(N2)
      IER = 1
      RETURN
C
C N .LT. 3 OR IST IS OUT OF RANGE
C
    5 PZ = 0.
      IER = -1
      RETURN
C
C NODES ARE COLLINEAR
C
    6 PZ = 0.
      IER = -2
      RETURN
      END
      SUBROUTINE INTRC1 (N,PX,PY,X,Y,Z,IADJ,IEND,IFLAG,
     .                   ZXZY, IST, PZ,IER)
      INTEGER N, IADJ(1), IEND(N), IFLAG, IST, IER
      REAL    PX, PY, X(N), Y(N), Z(N), ZXZY(2,N), PZ
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A TRIANGULATION OF A SET OF POINTS IN THE PLANE,
C THIS ROUTINE DETERMINES A PIECEWISE CUBIC FUNCTION F(X,Y)
C WHICH INTERPOLATES A SET OF DATA VALUES AND PARTIAL
C DERIVATIVES AT THE VERTICES.  F HAS CONTINUOUS FIRST
C DERIVATIVES OVER THE MESH AND EXTENDS BEYOND THE MESH
C BOUNDARY ALLOWING EXTRAPOLATION.  INTERPOLATION IS EXACT
C FOR QUADRATIC DATA.  THE VALUE OF F AT (PX,PY) IS
C RETURNED.  INTRC1 IS PART OF AN INTERPOLATION PACKAGE
C WHICH PROVIDES ROUTINES TO GENERATE, UPDATE AND PLOT THE
C MESH.
C
C INPUT PARAMETERS -     N - NUMBER OF NODES IN THE MESH.
C                            N .GE. 3.
C
C                    PX,PY - COORDINATES OF A POINT AT WHICH
C                            F IS TO BE EVALUATED.
C
C                      X,Y - VECTORS OF COORDINATES OF THE
C                            NODES IN THE MESH.
C
C                        Z - VECTOR OF DATA VALUES AT THE
C                            NODES.
C
C                     IADJ - SET OF ADJACENCY LISTS OF NODES
C                            IN THE MESH.
C
C                     IEND - POINTERS TO THE ENDS OF
C                            ADJACENCY LISTS IN IADJ FOR
C                            EACH NODE IN THE MESH.
C
C                    IFLAG - OPTION INDICATOR
C                            IFLAG = 0 IF INTRC1 IS TO
C                                      PROVIDE DERIVATIVE
C                                      ESTIMATES (FROM
C                                      GRADL).
C                            IFLAG = 1 IF DERIVATIVES ARE
C                                      USER PROVIDED.
C
C                     ZXZY - 2 BY N ARRAY WHOSE COLUMNS
C                            CONTAIN ESTIMATED PARTIAL DER-
C                            IVATIVES AT THE NODES (X PAR-
C                            TIALS IN THE FIRST ROW) IF
C                            IFLAG = 1, NOT USED IF IFLAG
C                            = 0.
C
C                      IST - INDEX OF THE STARTING NODE IN
C                            THE SEARCH FOR A TRIANGLE CON-
C                            TAINING (PX,PY).  1 .LE. IST
C                            .LE. N.  THE OUTPUT VALUE OF
C                            IST FROM A PREVIOUS CALL MAY
C                            BE A GOOD CHOICE.
C
C IADJ AND IEND MAY BE CREATED BY TRMESH AND DERIVATIVE
C   ESTIMATES MAY BE COMPUTED BY GRADL OR GRADG.
C
C INPUT PARAMETERS OTHER THAN IST ARE NOT ALTERED BY THIS
C   ROUTINE.
C
C OUTPUT PARAMETERS - IST - INDEX OF ONE OF THE VERTICES OF
C                           THE TRIANGLE CONTAINING (PX,PY)
C                           UNLESS IER .LT. 0.
C
C                      PZ - VALUE OF F AT (PX,PY), OR 0 IF
C                           IER .LT. 0.
C
C                     IER - ERROR INDICATOR
C                           IER = 0 IF NO ERRORS WERE
C                                   ENCOUNTERED.
C                           IER = 1 IF NO ERRORS WERE EN-
C                                   COUNTERED AND EXTRAPOLA-
C                                   TION WAS PERFORMED.
C                           IER = -1 IF N, IFLAG, OR IST IS
C                                    OUT OF RANGE.
C                           IER = -2 IF THE NODES ARE COL-
C                                    LINEAR.
C
C MODULES REFERENCED BY INTRC1 - TRFIND, TVAL,
C             (AND OPTIONALLY)   GRADL, GETNP, SETUP,
C                                GIVENS, ROTATE
C
C***********************************************************
C
      INTEGER NN, I1, I2, I3, IERR, N1, N2, INDX
      REAL    XP, YP, ZX1, ZY1, ZX2, ZY2, ZX3, ZY3, X1, Y1,
     .        X2, Y2, X3, Y3, Z1, Z2, Z3, DUM, DP, U, V, XQ,
     .        YQ, R1, R2, A1, A2, B1, B2, C1, C2, F1, F2
C
C LOCAL PARAMETERS -
C
C NN =                      LOCAL COPY OF N
C I1,I2,I3 =                VERTICES DETERMINED BY TRFIND
C IERR =                    ERROR FLAG FOR CALLS TO GRADL
C                             AND TVAL
C N1,N2 =                   ENDPOINTS OF THE CLOSEST BOUND-
C                             ARY EDGE TO P WHEN P IS OUT-
C                             SIDE OF THE MESH BOUNDARY
C INDX =                    IADJ INDEX OF N1 AS A NEIGHBOR
C                             OF N2
C XP,YP =                   LOCAL COPIES OF THE COORDINATES
C                             OF P=(PX,PY)
C ZX1,ZY1,ZX2,ZY2,ZX3,ZY3 = X AND Y DERIVATIVES AT THE
C                             VERTICES OF A TRIANGLE T WHICH
C                             CONTAINS P OR AT N1 AND N2
C X1,Y1,X2,Y2,X3,Y3 =       X,Y COORDINATES OF THE VERTICES
C                             OF T OR OF N1 AND N2
C Z1,Z2,Z3 =                DATA VALUES AT THE VERTICES OF T
C DUM =                     DUMMY VARIABLE FOR CALL TO TVAL
C DP =                      INNER PRODUCT OF N1-N2 AND P-N2
C U,V =                     X,Y COORDINATES OF THE VECTOR
C                             N2-N1
C XQ,YQ =                   X,Y COORDINATES OF THE CLOSEST
C                             BOUNDARY POINT TO P WHEN P IS
C                             OUTSIDE OF THE MESH BOUNDARY
C R1,R2 =                   BARYCENTRIC COORDINATES OF Q
C                             WITH RESPECT TO THE LINE SEG-
C                             MENT N2-N1 CONTAINING Q
C A1,A2,B1,B2,C1,C2 =       CARDINAL FUNCTIONS FOR EVALUAT-
C                             ING THE INTERPOLATORY SURFACE
C                             AT Q
C F1,F2 =                   CUBIC FACTORS USED TO COMPUTE
C                             THE CARDINAL FUNCTIONS
C
      NN = N
      PZ = 0.
      IF (NN .LT. 3  .OR.  IFLAG .LT. 0  .OR.  IFLAG .GT. 1
     .   .OR.  IST .LT. 1  .OR.  IST .GT. NN) GO TO 11
      XP = PX
      YP = PY
C
C FIND A TRIANGLE CONTAINING P IF P IS WITHIN THE MESH
C   BOUNDARY
C
      CALL TRFIND(IST,XP,YP,X,Y,IADJ,IEND, I1,I2,I3)
      IF (I1 .EQ. 0) GO TO 12
      IST = I1
      IF (I3 .EQ. 0) GO TO 3
      IF (IFLAG .NE. 1) GO TO 1
C
C DERIVATIVES ARE USER PROVIDED
C
      ZX1 = ZXZY(1,I1)
      ZX2 = ZXZY(1,I2)
      ZX3 = ZXZY(1,I3)
      ZY1 = ZXZY(2,I1)
      ZY2 = ZXZY(2,I2)
      ZY3 = ZXZY(2,I3)
      GO TO 2
C
C COMPUTE DERIVATIVE ESTIMATES AT THE VERTICES
C
    1 CALL GRADL(NN,I1,X,Y,Z,IADJ,IEND, ZX1,ZY1,IERR)
      CALL GRADL(NN,I2,X,Y,Z,IADJ,IEND, ZX2,ZY2,IERR)
      CALL GRADL(NN,I3,X,Y,Z,IADJ,IEND, ZX3,ZY3,IERR)
C
C SET LOCAL PARAMETERS FOR CALL TO TVAL
C
    2 X1 = X(I1)
      Y1 = Y(I1)
      X2 = X(I2)
      Y2 = Y(I2)
      X3 = X(I3)
      Y3 = Y(I3)
      Z1 = Z(I1)
      Z2 = Z(I2)
      Z3 = Z(I3)
      CALL TVAL(XP,YP,X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,ZX1,ZX2,
     .          ZX3,ZY1,ZY2,ZY3,0, PZ,DUM,DUM,IERR)
      IF (IERR .NE. 0) GO TO 12
      IER = 0
      RETURN
C
C P IS OUTSIDE OF THE MESH BOUNDARY.  EXTRAPOLATE TO P BY
C   PASSING A LINEAR FUNCTION OF ONE VARIABLE THROUGH THE
C   VALUE AND DIRECTIONAL DERIVATIVE (IN THE DIRECTION
C   P-Q) OF THE INTERPOLATORY SURFACE (TVAL) AT Q WHERE
C   Q IS THE CLOSEST BOUNDARY POINT TO P.
C
C DETERMINE Q BY TRAVERSING THE BOUNDARY STARTING FROM
C   THE RIGHTMOST VISIBLE NODE I1.
C
    3 N2 = I1
C
C SET N1 TO THE LAST NONZERO NEIGHBOR OF N2 AND COMPUTE DP
C
    4 INDX = IEND(N2) - 1
      N1 = IADJ(INDX)
      X1 = X(N1)
      Y1 = Y(N1)
      X2 = X(N2)
      Y2 = Y(N2)
      DP = (X1-X2)*(XP-X2) + (Y1-Y2)*(YP-Y2)
      IF (DP .LE. 0.) GO TO 5
      IF ((XP-X1)*(X2-X1) + (YP-Y1)*(Y2-Y1) .GT. 0.) GO TO 8
      N2 = N1
      GO TO 4
C
C N2 IS THE CLOSEST BOUNDARY POINT TO P.  COMPUTE PARTIAL
C   DERIVATIVES AT N2.
C
    5 IF (IFLAG .NE. 1) GO TO 6
      ZX2 = ZXZY(1,N2)
      ZY2 = ZXZY(2,N2)
      GO TO 7
    6 CALL GRADL(NN,N2,X,Y,Z,IADJ,IEND, ZX2,ZY2,IERR)
C
C COMPUTE EXTRAPOLATED VALUE AT P
C
    7 PZ = Z(N2) + ZX2*(XP-X2) + ZY2*(YP-Y2)
      IER = 1
      RETURN
C
C THE CLOSEST BOUNDARY POINT Q LIES ON N2-N1.  COMPUTE
C   PARTIALS AT N1 AND N2.
C
    8 IF (IFLAG .NE. 1) GO TO 9
      ZX1 = ZXZY(1,N1)
      ZY1 = ZXZY(2,N1)
      ZX2 = ZXZY(1,N2)
      ZY2 = ZXZY(2,N2)
      GO TO 10
    9 CALL GRADL(NN,N1,X,Y,Z,IADJ,IEND, ZX1,ZY1,IERR)
      CALL GRADL(NN,N2,X,Y,Z,IADJ,IEND, ZX2,ZY2,IERR)
C
C COMPUTE Q, ITS BARYCENTRIC COORDINATES, AND THE CARDINAL
C   FUNCTIONS FOR EXTRAPOLATION
C
   10 U = X2-X1
      V = Y2-Y1
      R1 = DP/(U**2 + V**2)
      R2 = 1. - R1
      XQ = R1*X1 + R2*X2
      YQ = R1*Y1 + R2*Y2
      F1 = R1*R1*R2
      F2 = R1*R2*R2
      A1 = R1 + (F1-F2)
      A2 = R2 - (F1-F2)
      B1 = U*F1
      B2 = -U*F2
      C1 = V*F1
      C2 = -V*F2
C
C COMPUTE THE VALUE OF THE INTERPOLATORY SURFACE (TVAL)
C   AT Q
C
      PZ = A1*Z(N1) + A2*Z(N2) + B1*ZX1 + B2*ZX2 +
     .     C1*ZY1 + C2*ZY2
C
C COMPUTE THE EXTRAPOLATED VALUE AT P
C
      PZ = PZ + (R1*ZX1 + R2*ZX2)*(XP-XQ) +
     .          (R1*ZY1 + R2*ZY2)*(YP-YQ)
      IER = 1
      RETURN
C
C N, IFLAG, OR IST OUT OF RANGE
C
   11 IER = -1
      RETURN
C
C NODES ARE COLLINEAR
C
   12 IER = -2
      RETURN
      END
      SUBROUTINE ROTATE (N,C,S, X,Y )
      INTEGER N
      REAL    C, S, X(N), Y(N)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C                                            ( C  S)
C   THIS ROUTINE APPLIES THE GIVENS ROTATION (     ) TO THE
C                                            (-S  C)
C               (X(1) ... X(N))
C 2 BY N MATRIX (             ).  THIS ROUTINE WAS TAKEN
C               (Y(1) ... Y(N))
C LINPACK.
C
C INPUT PARAMETERS -   N - NUMBER OF COLUMNS TO BE ROTATED.
C
C                    C,S - ELEMENTS OF THE GIVENS ROTATION.
C                          THESE MAY BE DETERMINED BY
C                          SUBROUTINE GIVENS.
C
C                    X,Y - VECTORS OF LENGTH .GE. N
C                          CONTAINING THE 2-VECTORS TO BE
C                          ROTATED.
C
C   THE PARAMETERS N, C, AND S ARE NOT ALTERED BY THIS
C ROUTINE.
C
C OUTPUT PARAMETERS - X,Y - ROTATED VECTORS
C
C MODULES REFERENCED BY ROTATE - NONE
C
C***********************************************************
C
      INTEGER I
      REAL    XI, YI
C
C LOCAL PARAMETERS -
C
C I =     DO-LOOP INDEX
C XI,YI = X(I), Y(I)
C
      IF (N .LE. 0 .OR. (C .EQ. 1. .AND. S .EQ. 0.)) RETURN
      DO 1 I = 1,N
        XI = X(I)
        YI = Y(I)
        X(I) = C*XI + S*YI
        Y(I) = -S*XI + C*YI
    1   CONTINUE
      RETURN
      END
      SUBROUTINE SETUP (XK,YK,ZK,XI,YI,ZI,S1,S2,R, ROW)
      REAL XK, YK, ZK, XI, YI, ZI, S1, S2, R, ROW(6)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE SETS UP THE I-TH ROW OF AN AUGMENTED RE-
C GRESSION MATRIX FOR A WEIGHTED LEAST-SQUARES FIT OF A
C QUADRATIC FUNCTION Q(X,Y) TO A SET OF DATA VALUES Z WHERE
C Q(XK,YK) = ZK.  THE FIRST 3 COLUMNS (QUADRATIC TERMS) ARE
C SCALED BY 1/S2 AND THE FOURTH AND FIFTH COLUMNS (LINEAR
C TERMS) ARE SCALED BY 1/S1.  THE WEIGHT IS (R-D)/(R*D) IF
C R .GT. D AND 0 IF R .LE. D, WHERE D IS THE DISTANCE
C BETWEEN NODES I AND K.
C
C INPUT PARAMETERS - XK,YK,ZK - COORDINATES AND DATA VALUE
C                               AT NODE K -- INTERPOLATED
C                               BY Q.
C
C                    XI,YI,ZI - COORDINATES AND DATA VALUE
C                               AT NODE I.
C
C                       S1,S2 - INVERSE SCALE FACTORS.
C
C                           R - RADIUS OF INFLUENCE ABOUT
C                               NODE K DEFINING THE WEIGHT.
C
C                         ROW - VECTOR OF LENGTH 6.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETER - ROW - VECTOR CONTAINING A ROW OF THE
C                          AUGMENTED REGRESSION MATRIX.
C
C MODULES REFERENCED BY SETUP - NONE
C
C INTRINSIC FUNCTION CALLED BY SETUP - SQRT
C
C***********************************************************
C
      INTEGER I
      REAL    DX, DY, DXSQ, DYSQ, D, W, W1, W2
C
C LOCAL PARAMETERS -
C
C I =    DO-LOOP INDEX
C DX =   XI - XK
C DY =   YI - YK
C DXSQ = DX*DX
C DYSQ = DY*DY
C D =    DISTANCE BETWEEN NODES K AND I
C W =    WEIGHT ASSOCIATED WITH THE ROW
C W1 =   W/S1
C W2 =   W/S2
C
      DX = XI - XK
      DY = YI - YK
      DXSQ = DX*DX
      DYSQ = DY*DY
      D = SQRT(DXSQ + DYSQ)
      IF (D .LE. 0.  .OR.  D .GE. R) GO TO 1
      W = (R-D)/R/D
      W1 = W/S1
      W2 = W/S2
      ROW(1) = DXSQ*W2
      ROW(2) = DX*DY*W2
      ROW(3) = DYSQ*W2
      ROW(4) = DX*W1
      ROW(5) = DY*W1
      ROW(6) = (ZI - ZK)*W
      RETURN
C
C NODES K AND I COINCIDE OR NODE I IS OUTSIDE OF THE RADIUS
C   OF INFLUENCE.  SET ROW TO THE ZERO VECTOR.
C
    1 DO 2 I = 1,6
    2   ROW(I) = 0.
      RETURN
      END
      FUNCTION TRVOL (X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3)
      REAL X1, X2, X3, Y1, Y2, Y3, Z1, Z2, Z3
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS FUNCTION COMPUTES THE INTEGRAL OVER A TRIANGLE OF
C THE PLANAR SURFACE WHICH INTERPOLATES DATA VALUES AT THE
C VERTICES.
C
C INPUT PARAMETERS - X1,X2,X3 - X COORDINATES OF THE
C                               VERTICES OF THE TRIANGLE IN
C                               COUNTERCLOCKWISE ORDER.
C
C                    Y1,Y2,Y3 - Y COORDINATES OF THE
C                               VERTICES OF THE TRIANGLE IN
C                               THE SAME ORDER AS THE X
C                               COORDINATES.
C
C                    Z1,Z2,Z3 - DATA VALUES AT THE VERTICES
C                               IN THE SAME ORDER AS THE
C                               COORDINATES.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION.
C
C OUTPUT PARAMETER -    TRVOL - VOLUME UNDER THE INTERPOLA-
C                               TORY SURFACE ABOVE THE
C                               TRIANGLE OR ZERO IF THE
C                               COORDINATES ARE INCORRECTLY
C                               ORDERED OR COLLINEAR.
C
C MODULES REFERENCED BY TRVOL - NONE
C
C***********************************************************
C
      T1 = X2*Y3 - X3*Y2
      T2 = X3*Y1 - X1*Y3
      T3 = X1*Y2 - X2*Y1
      AREA = T1 + T2 + T3
      IF (AREA .LT. 0.) AREA = 0.
C
C AREA IS TWICE THE AREA OF THE TRIANGLE.  TRVOL IS THE MEAN
C   OF THE DATA VALUES TIMES THE AREA OF THE TRIANGLE.
C
      TRVOL = (Z1 + Z2 + Z3)*AREA/6.
      RETURN
      END
      SUBROUTINE TVAL (X,Y,X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,ZX1,
     .                 ZX2,ZX3,ZY1,ZY2,ZY3,IFLAG, W,WX,WY,
     .                 IER)
      INTEGER IFLAG, IER
      REAL    X, Y, X1, X2, X3, Y1, Y2, Y3, Z1, Z2, Z3,
     .        ZX1, ZX2, ZX3, ZY1, ZY2, ZY3, W, WX, WY
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN FUNCTION VALUES AND FIRST PARTIAL DERIVATIVES AT
C THE THREE VERTICES OF A TRIANGLE, THIS ROUTINE DETERMINES
C A FUNCTION W WHICH AGREES WITH THE GIVEN DATA, RETURNING
C THE VALUE AND (OPTIONALLY) FIRST PARTIAL DERIVATIVES OF W
C AT A POINT (X,Y) IN THE TRIANGLE.  THE INTERPOLATION
C METHOD IS EXACT FOR QUADRATIC POLYNOMIAL DATA.  THE
C TRIANGLE IS PARTITIONED INTO THREE SUBTRIANGLES WITH
C EQUAL AREAS.  W IS CUBIC IN EACH SUBTRIANGLE AND ALONG
C THE EDGES, BUT HAS ONLY ONE CONTINUOUS DERIVATIVE ACROSS
C EDGES.  THE NORMAL DERIVATIVE OF W VARIES LINEARLY ALONG
C EACH OUTER EDGE.  THE VALUES AND PARTIAL DERIVATIVES OF W
C ALONG A TRIANGLE EDGE DEPEND ONLY ON THE DATA VALUES AT
C THE ENDPOINTS OF THE EDGE.  THUS THE METHOD YIELDS C-1
C CONTINUITY WHEN USED TO INTERPOLATE OVER A TRIANGULAR
C GRID.  THIS ALGORITHM IS DUE TO C. L. LAWSON.
C
C INPUT PARAMETERS -   X,Y - COORDINATES OF A POINT AT WHICH
C                            W IS TO BE EVALUATED.
C
C        X1,X2,X3,Y1,Y2,Y3 - COORDINATES OF THE VERTICES OF
C                            A TRIANGLE CONTAINING (X,Y).
C
C                 Z1,Z2,Z3 - FUNCTION VALUES AT THE VERTICES
C                            TO BE INTERPOLATED.
C
C              ZX1,ZX2,ZX3 - X-DERIVATIVE VALUES AT THE
C                            VERTICES.
C
C              ZY1,ZY2,ZY3 - Y-DERIVATIVE VALUES AT THE
C                            VERTICES.
C
C                    IFLAG - OPTION INDICATOR
C                            IFLAG = 0 IF ONLY W IS TO BE
C                                      COMPUTED.
C                            IFLAG = 1 IF W, WX, AND WY ARE
C                                      TO BE RETURNED.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETERS -   W - ESTIMATED VALUE OF THE INTERP-
C                           OLATORY FUNCTION AT (X,Y) IF
C                           IER = 0.  OTHERWISE W = 0.
C
C                   WX,WY - PARTIAL DERIVATIVES OF W AT
C                           (X,Y) IF IER = 0 AND IFLAG = 1,
C                           UNCHANGED IF IFLAG .NE. 1, ZERO
C                           IF IER .NE. 0 AND IFLAG = 1.
C
C                     IER - ERROR INDICATOR
C                           IER = 0 IF NO ERRORS WERE
C                                   ENCOUNTERED.
C                           IER = 1 IF THE VERTICES OF THE
C                                   TRIANGLE ARE COLLINEAR.
C
C MODULES REFERENCED BY TVAL - NONE
C
C INTRINSIC FUNCTION CALLED BY TVAL - AMIN1
C
C***********************************************************
C
      INTEGER I, IP1, IP2, IP3
      REAL    U(3), V(3), SL(3), AREA, XP, YP, R(3), RX(3),
     .        RY(3), PHI(3), PHIX(3), PHIY(3), RMIN, C1, C2,
     .        RO(3), ROX(3), ROY(3), F(3), G(3), GX(3),
     .        GY(3), P(3), PX(3), PY(3), Q(3), QX(3), QY(3),
     .        A(3), AX(3), AY(3), B(3), BX(3), BY(3), C(3),
     .        CX(3), CY(3)
C
C LOCAL PARAMETERS -
C
C I =               DO-LOOP INDEX
C IP1,IP2,IP3 =     PERMUTED INDICES FOR COMPUTING RO, ROX,
C                     AND ROY
C U(K) =            X-COMPONENT OF THE VECTOR REPRESENTING
C                     THE SIDE OPPOSITE VERTEX K
C V(K) =            Y-COMPONENT OF THE VECTOR REPRESENTING
C                     THE SIDE OPPOSITE VERTEX K
C SL(K) =           SQUARE OF THE LENGTH OF THE SIDE
C                     OPPOSITE VERTEX K
C AREA =            TWICE THE AREA OF THE TRIANGLE
C XP,YP =           X-X1, Y-Y1
C R(K) =            K-TH BARYCENTRIC COORDINATE
C RX(K),RY(K) =     X,Y PARTIAL DERIVATIVES OF R(K)
C PHI(K)            R(K-1)*R(K+1) -- QUADRATIC
C PHIX(K),PHIY(K) = X,Y PARTIALS OF PHI(K)
C RMIN =            MIN(R1,R2,R3)
C C1,C2 =           FACTORS FOR COMPUTING RO
C RO(K) =           FACTORS FOR COMPUTING G -- CUBIC
C                     CORRECTION TERMS
C ROX(K),ROY(K) =   X,Y PARTIALS OF RO(K)
C F(K) =            FACTORS FOR COMPUTING G, GX, AND GY --
C                     CONSTANT
C G(K) =            FACTORS FOR COMPUTING THE CARDINAL
C                     FUNCTIONS -- CUBIC
C GX(K),GY(K) =     X,Y PARTIALS OF G(K)
C P(K) =            G(K) + PHI(K)
C PX(K),PY(K) =     X,Y PARTIALS OF P(K)
C Q(K) =            G(K) - PHI(K)
C QX(K),QY(K) =     X,Y PARTIALS OF Q(K)
C A(K) =            CARDINAL FUNCTION WHOSE COEFFICIENT IS
C                     Z(K)
C AX(K),AY(K) =     X,Y PARTIALS OF A(K) -- CARDINAL
C                     FUNCTIONS FOR WX AND WY
C B(K) =            TWICE THE CARDINAL FUNCTION WHOSE
C                     COEFFICIENT IS ZX(K)
C BX(K),BY(K) =     X,Y PARTIALS OF B(K)
C C(K) =            TWICE THE CARDINAL FUNCTION WHOSE
C                     COEFFICIENT IS ZY(K)
C CX(K),CY(K) =     X,Y PARTIALS OF C(K)
C
      U(1) = X3 - X2
      U(2) = X1 - X3
      U(3) = X2 - X1
C
      V(1) = Y3 - Y2
      V(2) = Y1 - Y3
      V(3) = Y2 - Y1
C
      DO 1 I = 1,3
        SL(I) = U(I)*U(I) + V(I)*V(I)
    1   CONTINUE
C
C AREA = 3-1 X 3-2
C
      AREA = U(1)*V(2) - U(2)*V(1)
      IF (AREA .EQ. 0.) GO TO 9
C
C R(1) = (2-3 X 2-(X,Y))/AREA, R(2) = (1-(X,Y) X 1-3)/AREA,
C   R(3) = (1-2 X 1-(X,Y))/AREA
C
      R(1) = (U(1)*(Y-Y2) - V(1)*(X-X2))/AREA
      XP = X - X1
      YP = Y - Y1
      R(2) = (U(2)*YP - V(2)*XP)/AREA
      R(3) = (U(3)*YP - V(3)*XP)/AREA
      IER = 0
C
      PHI(1) = R(2)*R(3)
      PHI(2) = R(3)*R(1)
      PHI(3) = R(1)*R(2)
C
      RMIN = AMIN1(R(1),R(2),R(3))
      IF (RMIN .NE. R(1)) GO TO 3
      IP1 = 1
      IP2 = 2
      IP3 = 3
      GO TO 5
    3 IF (RMIN .NE. R(2)) GO TO 4
      IP1 = 2
      IP2 = 3
      IP3 = 1
      GO TO 5
    4 IP1 = 3
      IP2 = 1
      IP3 = 2
C
    5 C1 = RMIN*RMIN/2.
      C2 = RMIN/3.
      RO(IP1) = (PHI(IP1) + 5.*C1/3.)*R(IP1) - C1
      RO(IP2) = C1*(R(IP3) - C2)
      RO(IP3) = C1*(R(IP2) - C2)
C
      F(1) = 3.*(SL(2)-SL(3))/SL(1)
      F(2) = 3.*(SL(3)-SL(1))/SL(2)
      F(3) = 3.*(SL(1)-SL(2))/SL(3)
C
      G(1) = (R(2)-R(3))*PHI(1) + F(1)*RO(1) - RO(2) + RO(3)
      G(2) = (R(3)-R(1))*PHI(2) + F(2)*RO(2) - RO(3) + RO(1)
      G(3) = (R(1)-R(2))*PHI(3) + F(3)*RO(3) - RO(1) + RO(2)
C
      DO 6 I = 1,3
        P(I) = G(I) + PHI(I)
        Q(I) = G(I) - PHI(I)
    6   CONTINUE
C
      A(1) = R(1) + G(3) - G(2)
      A(2) = R(2) + G(1) - G(3)
      A(3) = R(3) + G(2) - G(1)
C
      B(1) = U(3)*P(3) + U(2)*Q(2)
      B(2) = U(1)*P(1) + U(3)*Q(3)
      B(3) = U(2)*P(2) + U(1)*Q(1)
C
      C(1) = V(3)*P(3) + V(2)*Q(2)
      C(2) = V(1)*P(1) + V(3)*Q(3)
      C(3) = V(2)*P(2) + V(1)*Q(1)
C
C W IS A LINEAR COMBINATION OF THE CARDINAL FUNCTIONS
C
      W = A(1)*Z1 + A(2)*Z2 + A(3)*Z3 + (B(1)*ZX1 + B(2)*ZX2
     .    + B(3)*ZX3 + C(1)*ZY1 + C(2)*ZY2 + C(3)*ZY3)/2.
      IF (IFLAG .NE. 1) RETURN
C
C COMPUTE WX AND WY
C
      DO 7 I = 1,3
        RX(I) = -V(I)/AREA
        RY(I) = U(I)/AREA
    7   CONTINUE
      PHIX(1) = R(2)*RX(3) + RX(2)*R(3)
      PHIY(1) = R(2)*RY(3) + RY(2)*R(3)
      PHIX(2) = R(3)*RX(1) + RX(3)*R(1)
      PHIY(2) = R(3)*RY(1) + RY(3)*R(1)
      PHIX(3) = R(1)*RX(2) + RX(1)*R(2)
      PHIY(3) = R(1)*RY(2) + RY(1)*R(2)
C
      ROX(IP1) = RX(IP1)*(PHI(IP1) + 5.*C1) +
     .           R(IP1)*(PHIX(IP1) - RX(IP1))
      ROY(IP1) = RY(IP1)*(PHI(IP1) + 5.*C1) +
     .           R(IP1)*(PHIY(IP1) - RY(IP1))
      ROX(IP2) = RX(IP1)*(PHI(IP2) - C1) + C1*RX(IP3)
      ROY(IP2) = RY(IP1)*(PHI(IP2) - C1) + C1*RY(IP3)
      ROX(IP3) = RX(IP1)*(PHI(IP3) - C1) + C1*RX(IP2)
      ROY(IP3) = RY(IP1)*(PHI(IP3) - C1) + C1*RY(IP2)
C
      GX(1) = (RX(2) - RX(3))*PHI(1) + (R(2) - R(3))*PHIX(1)
     .        + F(1)*ROX(1) - ROX(2) + ROX(3)
      GY(1) = (RY(2) - RY(3))*PHI(1) + (R(2) - R(3))*PHIY(1)
     .        + F(1)*ROY(1) - ROY(2) + ROY(3)
      GX(2) = (RX(3) - RX(1))*PHI(2) + (R(3) - R(1))*PHIX(2)
     .        + F(2)*ROX(2) - ROX(3) + ROX(1)
      GY(2) = (RY(3) - RY(1))*PHI(2) + (R(3) - R(1))*PHIY(2)
     .        + F(2)*ROY(2) - ROY(3) + ROY(1)
      GX(3) = (RX(1) - RX(2))*PHI(3) + (R(1) - R(2))*PHIX(3)
     .        + F(3)*ROX(3) - ROX(1) + ROX(2)
      GY(3) = (RY(1) - RY(2))*PHI(3) + (R(1) - R(2))*PHIY(3)
     .        + F(3)*ROY(3) - ROY(1) + ROY(2)
C
      DO 8 I = 1,3
        PX(I) = GX(I) + PHIX(I)
        PY(I) = GY(I) + PHIY(I)
        QX(I) = GX(I) - PHIX(I)
        QY(I) = GY(I) - PHIY(I)
    8   CONTINUE
C
      AX(1) = RX(1) + GX(3) - GX(2)
      AY(1) = RY(1) + GY(3) - GY(2)
      AX(2) = RX(2) + GX(1) - GX(3)
      AY(2) = RY(2) + GY(1) - GY(3)
      AX(3) = RX(3) + GX(2) - GX(1)
      AY(3) = RY(3) + GY(2) - GY(1)
C
      BX(1) = U(3)*PX(3) + U(2)*QX(2)
      BY(1) = U(3)*PY(3) + U(2)*QY(2)
      BX(2) = U(1)*PX(1) + U(3)*QX(3)
      BY(2) = U(1)*PY(1) + U(3)*QY(3)
      BX(3) = U(2)*PX(2) + U(1)*QX(1)
      BY(3) = U(2)*PY(2) + U(1)*QY(1)
C
      CX(1) = V(3)*PX(3) + V(2)*QX(2)
      CY(1) = V(3)*PY(3) + V(2)*QY(2)
      CX(2) = V(1)*PX(1) + V(3)*QX(3)
      CY(2) = V(1)*PY(1) + V(3)*QY(3)
      CX(3) = V(2)*PX(2) + V(1)*QX(1)
      CY(3) = V(2)*PY(2) + V(1)*QY(1)
C
C WX AND WY ARE LINEAR COMBINATIONS OF THE CARDINAL
C   FUNCTIONS
C
      WX = AX(1)*Z1 + AX(2)*Z2 + AX(3)*Z3 + (BX(1)*ZX1 +
     .     BX(2)*ZX2 + BX(3)*ZX3 + CX(1)*ZY1 + CX(2)*ZY2 +
     .     CX(3)*ZY3)/2.
      WY = AY(1)*Z1 + AY(2)*Z2 + AY(3)*Z3 + (BY(1)*ZX1 +
     .     BY(2)*ZX2 + BY(3)*ZX3 + CY(1)*ZY1 + CY(2)*ZY2 +
     .     CY(3)*ZY3)/2.
      RETURN
C
C VERTICES ARE COLLINEAR
C
    9 IER = 1
      W = 0.
      IF (IFLAG .NE. 1) RETURN
      WX = 0.
      WY = 0.
      RETURN
      END
      SUBROUTINE UNIF (N,X,Y,Z,IADJ,IEND,NROW,NX,NY,PX,PY,
     .                 IFLAG, ZXZY, ZZ,IER)
      INTEGER N, IADJ(1), IEND(N), NROW, NX, NY, IFLAG, IER
      REAL    X(N), Y(N), Z(N), PX(NX), PY(NY), ZXZY(2,N),
     .        ZZ(NROW,NY)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A THIESSEN TRIANGULATION OF A SET OF POINTS IN THE
C PLANE WITH CORRESPONDING DATA VALUES, THIS ROUTINE INTERP-
C OLATES THE DATA VALUES TO A SET OF RECTANGULAR GRID POINTS
C FOR SUCH APPLICATIONS AS CONTOURING.  THE INTERPOLANT IS
C ONCE CONTINUOUSLY DIFFERENTIABLE.  EXTRAPOLATION IS PER-
C FORMED AT GRID POINTS EXTERIOR TO THE TRIANGULATION.
C
C INPUT PARAMETERS - N - NUMBER OF NODES IN THE TRIANGULA-
C                        TION.  N .GE. 3
C
C                X,Y,Z - N-VECTORS OF NODAL COORDINATES AND
C                        DATA VALUES.
C
C            IADJ,IEND - TRIANGULATION DATA STRUCTURE -- MAY
C                        BE CREATED BY TRMESH.
C
C                 NROW - NUMBER OF ROWS IN THE DIMENSION
C                        STATEMENT OF ZZ.
C
C                NX,NY - NUMBER OF ROWS AND COLUMNS, RESPEC-
C                        TIVELY, IN THE RECTANGULAR GRID.
C                        1 .LE. NX .LE. NROW AND 1 .LE. NY.
C
C                PX,PY - VECTORS OF LENGTH NX AND NY, RE-
C                        SPECTIVELY, CONTAINING THE COORDI-
C                        NATES OF THE GRID LINES.
C
C                IFLAG - OPTION INDICATOR
C                        IFLAG = 0 IF DERIVATIVE ESTIMATES
C                                  AT THE VERTICES OF A
C                                  TRIANGLE ARE TO BE RECOM-
C                                  PUTED FOR EACH GRID POINT
C                                  IN THE TRIANGLE AND NOT
C                                  SAVED.
C                        IFLAG = 1 IF DERIVATIVE ESTIMATES
C                                  ARE INPUT IN ZXZY.
C                        IFLAG = 2 IF DERIVATIVE ESTIMATES
C                                  ARE TO BE COMPUTED ONCE
C                                  FOR EACH NODE (BY GRADL)
C                                  AND SAVED IN ZXZY.
C
C                 ZXZY - NOT USED IF IFLAG = 0, 2 BY N ARRAY
C                        WHOSE COLUMNS CONTAIN THE X AND Y
C                        PARTIAL DERIVATIVE ESTIMATES (X
C                        PARTIALS IN THE FIRST ROW) IF
C                        IFLAG = 1, OR 2 BY N ARRAY (OR WORK
C                        SPACE OF LENGTH .GE. 2*N) IF IFLAG
C                        = 2.
C
C DERIVATIVE ESTIMATES MAY BE COMPUTED BY GRADL OR GRADG.
C
C                   ZZ - NROW BY NCOL ARRAY WITH NROW .GE.
C                        NX AND NCOL .GE. NY.
C
C NONE OF THE INPUT PARAMETERS ARE ALTERED EXCEPT AS
C   NOTED BELOW.
C
C OUTPUT PARAMETERS - ZXZY - 2 BY N ARRAY WHOSE COLUMNS CON-
C                            TAIN X AND Y PARTIAL DERIVATIVE
C                            ESTIMATES AT THE NODES IF IFLAG
C                            = 2 AND IER .GE. 0, NOT ALTERED
C                            IF IFLAG .NE. 2.
C
C                       ZZ - INTERPOLATED VALUES AT THE GRID
C                            POINTS IF IER .GE. 0.
C                            ZZ(I,J) = F(PX(I),PY(J)) FOR
C                            I = 1,...,NX AND J = 1,...,NY.
C
C                      IER - ERROR INDICATOR
C                            IER .GE. 0 IF NO ERRORS WERE
C                                       ENCOUNTERED.  IER
C                                       CONTAINS THE NUMBER
C                                       OF GRID POINTS EXT-
C                                       ERIOR TO THE TRIAN-
C                                       GULATION BOUNDARY.
C                            IER  =  -1 IF N, NX, NY, OR
C                                       IFLAG IS OUT OF
C                                       RANGE.
C                            IER  =  -2 IF THE NODES ARE
C                                       COLLINEAR.
C
C MODULES REFERENCED BY UNIF - INTRC1, TRFIND, TVAL,
C           (AND OPTIONALLY)   GRADL, GETNP, SETUP, GIVENS,
C                                AND ROTATE
C
C***********************************************************
C
      INTEGER NST, IST, NN, NI, NJ, IFL, I, J, NIT, IERR,
     .        NEX
      DATA    NST/1/
C
C LOCAL PARAMETERS -
C
C IST =   PARAMETER FOR INTRC1
C NST =   INITIAL VALUE FOR IST
C NN =    LOCAL COPY OF N
C NI,NJ = LOCAL COPIES OF NX AND NY
C IFL =   LOCAL COPY OF IFLAG FOR INTRC1
C I,J =   DO-LOOP INDICES
C IERR =  ERROR FLAG FOR CALLS TO GRADL AND INTRC1
C NEX =   NUMBER OF GRID POINTS EXTERIOR TO THE TRIANGULA-
C           TION BOUNDARY (NUMBER OF EXTRAPOLATED VALUES)
C
      NN = N
      NI = NX
      NJ = NY
      IFL = IFLAG
      IF (NN .LT. 3  .OR.  NI .LT. 1  .OR.  NI .GT. NROW
     .   .OR.  NJ .LT. 1  .OR.  IFL .LT. 0  .OR.
     .   IFL .GT. 2) GO TO 4
      IST = NST
      IF (IFL .NE. 2) GO TO 2
C
C COMPUTE DERIVATIVE ESTIMATES AT THE NODES.
C
      IFL = 1
      DO 1 I = 1,NN
        CALL GRADL(NN,I,X,Y,Z,IADJ,IEND, ZXZY(1,I),
     .             ZXZY(2,I),IERR)
        IF (IERR .LT. 0) GO TO 5
    1   CONTINUE
C
C COMPUTE INTERPOLATED VALUES
C
    2 NEX = 0
      DO 3 J = 1,NJ
        DO 3 I = 1,NI
          CALL INTRC1(NN,PX(I),PY(J),X,Y,Z,IADJ,IEND,IFL,
     .                ZXZY, IST, ZZ(I,J),IERR)
          IF (IERR .LT. 0) GO TO 5
          NEX = NEX + IERR
    3     CONTINUE
      IER = NEX
      RETURN
C
C N, NX, NY, OR IFLAG IS OUT OF RANGE
C
    4 IER = -1
      RETURN
C
C TRIANGULATION NODES ARE COLLINEAR
C
    5 IER = -2
      RETURN
      END
      FUNCTION VOLUME (N,X,Y,Z,IADJ,IEND)
      INTEGER N, IADJ(1), IEND(N)
      REAL    X(N), Y(N), Z(N)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A SET OF N DATA POINTS (X(I),Y(I)) AND FUNCTION
C VALUES Z(I)=F(X(I),Y(I)) AND A TRIANGULATION COVERING THE
C CONVEX HULL H OF THE DATA POINTS, THIS FUNCTION APPROXI-
C MATES THE INTEGRAL OF F OVER H BY INTEGRATING THE PIECE-
C WISE LINEAR INTERPOLANT OF THE DATA VALUES.  VOLUME IS
C PART OF AN INTERPOLATION PACKAGE WHICH PROVIDES ROUTINES
C TO CREATE, UPDATE, AND PLOT THE TRIANGULAR MESH.
C
C INPUT PARAMETERS -      N - NUMBER OF NODES IN THE MESH.
C                             N .GE. 3.
C
C                       X,Y - VECTORS OF COORDINATES OF
C                             THE NODES IN THE MESH.
C
C                         Z - VECTOR OF DATA VALUES AT THE
C                             NODES.
C
C                      IADJ - SET OF ADJACENCY LISTS OF
C                             NODES IN THE MESH.
C
C                      IEND - POINTERS TO THE ENDS OF
C                             ADJACENCY LISTS IN IADJ FOR
C                             EACH NODE IN THE MESH.
C
C IADJ AND IEND MAY BE CREATED BY TRMESH.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION.
C
C OUTPUT PARAMETER - VOLUME - SUM OF THE VOLUMES OF THE
C                             LINEAR INTERPOLANTS ON THE
C                             TRIANGLES.
C
C MODULE REFERENCED BY VOLUME - TRVOL
C
C***********************************************************
C
      INTEGER NN, NM2, N1, N2, N3, INDF, INDL, INDX
      REAL    SUM, XN1, YN1, ZN1
C
C LOCAL PARAMETERS -
C
C NN =          LOCAL COPY OF N
C NM2 =         N-2
C N1,N2,N3 =    VERTICES OF A TRIANGLE IN COUNTERCLOCKWISE
C                 ORDER
C INDF =        IADJ INDEX OF THE FIRST NEIGHBOR OF N1
C INDL =        IADJ INDEX OF THE LAST NEIGHBOR OF N1
C INDX =        IADJ INDEX VARYING FROM INDF TO INDL
C SUM =         TEMPORARY STORAGE FOR ACCUMULATED VOLUME
C XN1,YN1,ZN1 = X(N1), Y(N1), Z(N1) -- STORED LOCALLY FOR
C                 EFFICIENCY
C
      NN = N
      IF (NN .LT. 3) GO TO 3
C
C INITIALIZATION
C
      NM2 = NN-2
      INDF = 1
      SUM = 0.
C
C LOOP ON TRIANGLES (N1,N2,N3) SUCH THAT N2 AND N3 ARE
C   ADJACENT NEIGHBORS OF N1 WHICH ARE BOTH LARGER THAN N1
C
      DO 2 N1 = 1,NM2
        XN1 = X(N1)
        YN1 = Y(N1)
        ZN1 = Z(N1)
        INDL = IEND(N1)
        DO 1 INDX = INDF,INDL
          N2 = IADJ(INDX)
          N3 = IADJ(INDX+1)
          IF (INDX .EQ. INDL) N3 = IADJ(INDF)
          IF (N2 .LT. N1  .OR.  N3 .LT. N1) GO TO 1
          SUM = SUM + TRVOL(XN1,X(N2),X(N3),YN1,Y(N2),Y(N3),
     .                      ZN1,Z(N2),Z(N3))
    1     CONTINUE
        INDF = INDL + 1
    2   CONTINUE
C
      VOLUME = SUM
      RETURN
C
C N IS OUT OF RANGE
C
    3 VOLUME = 0.
      RETURN
      END
C
      SUBROUTINE ARCTST (N,X,Y,IADJ,IEND, LL, LIST,IER)
      INTEGER N, IADJ(1), IEND(N), LL, LIST(2,1), IER
      REAL    X(N), Y(N)
      LOGICAL SWPTST
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE DETERMINES THE NUMBER OF ARCS IN A TRIANGU-
C LATION WHICH ARE NOT LOCALLY OPTIMAL AS DEFINED BY LOGICAL
C FUNCTION SWPTST.  A THIESSEN TRIANGULATION (SEE SUBROUTINE
C TRMESH) IS CHARACTERIZED BY THE PROPERTY THAT ALL ARCS ARE
C LOCALLY OPTIMAL.  THE ALGORITHM CONSISTS OF APPLYING THE
C SWAP TEST TO ALL INTERIOR ARCS.
C
C INPUT PARAMETERS -
C
C       N - NUMBER OF NODES IN THE TRIANGULATION.  N .GE. 3.
C
C       X,Y - COORDINATES OF THE NODES.
C
C       IADJ,IEND - DATA STRUCTURE CONTAINING THE TRIANGULA-
C                   TION.  SEE SUBROUTINE TRMESH.
C
C THE ABOVE PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C       LL - NUMBER OF COLUMNS RESERVED FOR LIST (SEE THE
C            OUTPUT VALUE OF LIST).  LL .GE. 0.
C
C       LIST - INTEGER ARRAY DIMENSIONED 2 BY LL (OR VECTOR
C              OF LENGTH .GE. 2*LL) IF LL .GT. 0.
C
C OUTPUT PARAMETERS -
C
C       LL - NUMBER OF ARCS WHICH ARE NOT LOCALLY OPTIMAL
C            UNLESS IER = 1.
C
C       LIST - COLUMNS CONTAIN THE ENDPOINT NODAL INDICES
C              (SMALLER INDEX IN THE FIRST ROW) OF THE FIRST
C              K NONOPTIMAL ARCS ENCOUNTERED, WHERE K IS THE
C              MINIMUM OF THE INPUT AND OUTPUT VALUES OF LL,
C              IF IER = 0 AND K .GT. 0.  THE NUMBER OF INTE-
C              RIOR ARCS IS 3N-2NB-3 .LE. 3(N-3) WHERE NB IS
C              THE NUMBER OF BOUNDARY NODES.  BOUNDARY ARCS
C              ARE OPTIMAL BY DEFINITION.
C
C       IER - ERROR INDICATOR
C             IER = 0 IF NO ERRORS WERE ENCOUNTERED.
C             IER = 1 IF N .LT. 3 OR LL .LT. 0.
C
C MODULES REQUIRED BY ARCTST - SWPTST
C
C***********************************************************
C
      NM1 = N - 1
      LMAX = LL
C
C TEST FOR ERRORS AND INITIALIZE FOR LOOP ON INTERIOR ARCS.
C
      IF (NM1 .LT. 2  .OR.  LMAX .LT. 0) GO TO 3
      IER = 0
      L = 0
      INDF = 1
C
C OUTER LOOP ON NODES IO1
C
      DO 2 IO1 = 1,NM1
        INDL = IEND(IO1)
        IN2 = IADJ(INDL)
C
C INNER LOOP ON NEIGHBORS IO2 OF IO1 SUCH THAT IO2 .GT. IO1
C   AND IO1-IO2 IS AN INTERIOR ARC -- (IO1,IO2,IN1) AND
C   (IO2,IO1,IN2) ARE TRIANGLES.
C
        DO 1 INDX = INDF,INDL
          IO2 = IADJ(INDX)
          IN1 = IADJ(INDX+1)
          IF (INDX .EQ. INDL) IN1 = IADJ(INDF)
          IF (IO2 .LT. IO1  .OR.  IN1 .EQ. 0  .OR.
     .        IN2 .EQ. 0) GO TO 1
C
C TEST FOR A SWAP.
C
          IF (.NOT. SWPTST(IN1,IN2,IO1,IO2,X,Y)) GO TO 1
          L = L + 1
          IF (L .GT. LMAX) GO TO 1
          LIST(1,L) = IO1
          LIST(2,L) = IO2
    1     IN2 = IO2
    2   INDF = INDL + 1
      LL = L
      RETURN
C
C N OR LL OUT OF RANGE
C
    3 IER = 1
      RETURN
      END
      SUBROUTINE CIRCUM (X1,X2,X3,Y1,Y2,Y3, CX,CY,IER)
      INTEGER IER
      REAL    X1, X2, X3, Y1, Y2, Y3, CX, CY
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS SUBROUTINE COMPUTES THE COORDINATES OF THE CENTER
C OF A CIRCLE DEFINED BY THREE POINTS IN THE PLANE.
C
C INPUT PARAMETERS -
C
C       X1,...,Y3 - X AND Y COORDINATES OF THREE POINTS IN
C                   THE PLANE.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETERS -
C
C       CX,CY - COORDINATES OF THE CENTER OF THE CIRCLE
C               UNLESS IER = 1.
C
C       IER - ERROR INDICATOR
C             IER = 0 IF NO ERRORS WERE ENCOUNTERED.
C             IER = 1 IF THE POINTS ARE COLLINEAR.
C
C MODULES REQUIRED BY CIRCUM - NONE
C
C***********************************************************
C
      REAL U(3), V(3), DS(3)
C
C SET U(K) AND V(K) TO THE X AND Y COMPONENTS OF THE EDGE
C   OPPOSITE VERTEX K, TREATING THE POINTS AS VERTICES OF
C   A TRIANGLE.
C
      U(1) = X3 - X2
      U(2) = X1 - X3
      U(3) = X2 - X1
      V(1) = Y3 - Y2
      V(2) = Y1 - Y3
      V(3) = Y2 - Y1
C
C SET A TO TWICE THE SIGNED AREA OF THE TRIANGLE.  A .GT. 0
C   IFF (X3,Y3) IS STRICTLY TO THE LEFT OF THE EDGE FROM
C   (X1,Y1) TO (X2,Y2).
C
      A = U(1)*V(2) - U(2)*V(1)
      IF (A .EQ. 0.) GO TO 2
C
C SET DS(K) TO THE SQUARED DISTANCE FROM THE ORIGIN TO
C   VERTEX K.
C
      DS(1) = X1**2 + Y1**2
      DS(2) = X2**2 + Y2**2
      DS(3) = X3**2 + Y3**2
C
C COMPUTE FACTORS OF CX AND CY.
C
      FX = 0.
      FY = 0.
      DO 1 I = 1,3
        FX = FX - DS(I)*V(I)
    1   FY = FY + DS(I)*U(I)
      CX = FX/2./A
      CY = FY/2./A
      IER = 0
      RETURN
C
C COLLINEAR POINTS
C
    2 IER = 1
      RETURN
      END
      INTEGER FUNCTION LOPTST (N1,N2,X,Y,IADJ,IEND)
      INTEGER N1, N2, IADJ(1), IEND(1)
      REAL    X(1), Y(1)
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   GIVEN A PAIR OF INDICES DEFINING A TRIANGULATION ARC,
C THIS FUNCTION DETERMINES WHETHER OR NOT THE ARC IS LOCALLY
C OPTIMAL AS DEFINED BY LOGICAL FUNCTION SWPTST.
C
C INPUT PARAMETERS -
C
C       N1,N2 - X,Y INDICES OF A PAIR OF ADJACENT NODES.
C
C       X,Y - NODAL COORDINATES.
C
C       IADJ,IEND - DATA STRUCTURE CONTAINING THE TRIANGULA-
C                   TION.  SEE SUBROUTINE TRMESH.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS FUNCTION.
C
C OUTPUT PARAMETER -
C
C       LOPTST = -2 IF N1 AND N2 ARE NOT ADJACENT,
C              = -1 IF N1-N2 IS AN INTERIOR ARC WHICH IS
C                   NOT LOCALLY OPTIMAL,
C              =  0 IF N1-N2 SATISFIES THE NEUTRAL CASE (THE
C                   VERTICES OF THE CORRESPONDING QUADRILAT-
C                   ERAL LIE ON A COMMON CIRCLE),
C              =  1 IF N1-N2 IS A LOCALLY OPTIMAL INTERIOR
C                   ARC,
C              =  2 IF N1-N2 IS A BOUNDARY ARC.
C      NOTE THAT N1-N2 IS LOCALLY OPTIMAL IFF LOPTST .GE. 0.
C
C MODULES REQUIRED BY LOPTST - NONE
C
C***********************************************************
C
      IO1 = N1
      IO2 = N2
C
C FIND THE INDEX OF IO2 AS A NEIGHBOR OF IO1
C
      INDF = 1
      IF (IO1 .GT. 1) INDF = IEND(IO1-1) + 1
      INDL = IEND(IO1)
    1 INDX = INDX - 1
      IF (IADJ(INDX) .EQ. IO2) GO TO 2
      IF (INDX .NE. INDF) GO TO 1
C
C N1 AND N2 ARE NOT ADJACENT
C
      LOPTST = -2
      RETURN
C
C DETERMINE IN1 AND IN2 SUCH THAT (IO1,IO2,IN1) AND
C   (IO2,IO1,IN2) ARE TRIANGLES.
C
    2 IF (INDX .NE. INDL) IN1 = IADJ(INDX+1)
      IF (INDX .EQ. INDL) IN1 = IADJ(INDF)
      IF (INDX .NE. INDF) IN2 = IADJ(INDX-1)
      IF (INDX .EQ. INDF) IN2 = IADJ(INDL)
      IF (IN1 .NE. 0  .AND.  IN2 .NE. 0) GO TO 3
C
C N1-N2 IS A BOUNDARY ARC
C
      LOPTST = 2
      RETURN
C
C COMPUTE COMPONENTS OF THE QUADRILATERAL SIDES.
C
    3 DX11 = X(IO1) - X(IN1)
      DX12 = X(IO2) - X(IN1)
      DX22 = X(IO2) - X(IN2)
      DX21 = X(IO1) - X(IN2)
C
      DY11 = Y(IO1) - Y(IN1)
      DY12 = Y(IO2) - Y(IN1)
      DY22 = Y(IO2) - Y(IN2)
      DY21 = Y(IO1) - Y(IN2)
C
C COMPUTE INNER PRODUCTS.
C
      COS1 = DX11*DX12 + DY11*DY12
      COS2 = DX22*DX21 + DY22*DY21
C
C IO1-IO2 IS LOCALLY OPTIMAL IFF A1+A2 .LE. 180 DEGREES
C   WHERE A1 AND A2 DENOTE THE ANGLES AT IN1 AND IN2.
C
      IF (COS1 .LT. 0.  .AND.  COS2 .LT. 0.) GO TO 4
      IF (COS1 .GT. 0.  .AND.  COS2 .GT. 0.) GO TO 5
C
C COMPUTE A QUANTITY WITH THE SIGN OF SIN(A1+A2).
C
      SIN12 = (DX11*DY12 - DX12*DY11)*COS2 +
     .        (DX22*DY21 - DX21*DY22)*COS1
      IF (SIN12 .LT. 0.) GO TO 4
      IF (SIN12 .GT. 0.) GO TO 5
C
C NEUTRAL CASE
C
      LOPTST = 0
      RETURN
C
C N1-N2 NOT LOCALLY OPTIMAL
C
    4 LOPTST = -1
      RETURN
C
C N1-N2 LOCALLY OPTIMAL
C
    5 LOPTST = 1
      RETURN
      END
*     FUNCTION STORE (X)
*     REAL X
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS FUNCTION FORCES ITS ARGMENT X TO BE STORED IN MAIN
C MEMORY.  THIS IS USEFUL FOR COMPUTING MACHINE DEPENDENT
C PARAMETERS (SUCH AS THE MACHINE PRECISION) WHERE IT IS
C NECESSARY TO AVOID COMPUTATION IN HIGH PRECISION REG-
C ISTERS.
C
C INPUT PARAMETER -
C
C       X - VALUE TO BE STORED.
C
C X IS NOT ALTERED BY THIS FUNCTION.
C
C OUTPUT PARAMETER -
C
C       STORE - VALUE OF X AFTER IT HAS BEEN STORED AND
C               (POSSIBLY) TRUNCATED OR ROUNDED TO THE
C               MACHINE WORD LENGTH.
C
C MODULES REQUIRED BY STORE - NONE
C
C***********************************************************
C
*     COMMON/STCOM/Y
*     Y = X
*     STORE = Y
*     RETURN
*     END
      SUBROUTINE TRMTST (N,X,Y,IADJ,IEND,TOL,LUN, IER)
      INTEGER N, IADJ(1), IEND(N), LUN, IER
      REAL    X(N), Y(N), TOL
C
C***********************************************************
C
C                                               ROBERT RENKA
C                                       OAK RIDGE NATL. LAB.
C                                             (615) 576-5139
C
C   THIS ROUTINE TESTS THE VALIDITY OF THE DATA STRUCTURE
C REPRESENTING A THIESSEN TRIANGULATION CREATED BY SUBROU-
C TINE TRMESH.  THE FOLLOWING PROPERTIES ARE TESTED --
C   1)  IEND(1) .GE. 3 AND IEND(K) .GE. IEND(K-1)+3 FOR K =
C       2,...,N (EACH NODE HAS AT LEAST THREE NEIGHBORS).
C   2)  0 .LE. IADJ(K) .LE. N FOR K = 1,...,IEND(N) (IADJ
C       ENTRIES ARE NODAL INDICES OR ZEROS REPRESENTING THE
C       BOUNDARY).
C   3)  NB .GE. 3, NT = 2N-NB-2, AND NA = 3N-NB-3 WHERE NB,
C       NT, AND NA ARE THE NUMBERS OF BOUNDARY NODES, TRI-
C       ANGLES, AND ARCS, RESPECTIVELY.
C   4)  EACH CIRCUMCIRCLE DEFINED BY THE VERTICES OF A TRI-
C       ANGLE CONTAINS NO NODES IN ITS INTERIOR.  THIS PROP-
C       ERTY DISTINGUISHES A THIESSEN TRIANGULATION FROM AN
C       ARBITRARY TRIANGULATION OF THE NODES.
C NOTE THAT NO TEST IS MADE FOR THE PROPERTY THAT A TRIANGU-
C LATION COVERS THE CONVEX HULL OF THE NODES, AND THUS A
C TEST ON A DATA STRUCTURE ALTERED BY SUBROUTINE DELETE
C SHOULD NOT RESULT IN AN ERROR.
C
C INPUT PARAMETERS -
C
C       N - NUMBER OF NODES.  N .GE. 3.
C
C       X,Y - NODAL COORDINATES.
C
C       IADJ,IEND - TRIANGULATION DATA STRUCTURE.  SEE SUB-
C                   ROUTINE TRMESH.
C
C       TOL - NONNEGATIVE TOLERANCE TO ALLOW FOR FLOATING-
C             POINT ERRORS IN THE CIRCUMCIRCLE TEST.  AN
C             ERROR SITUATION IS DEFINED AS R**2 - D**2 .GT.
C             TOL WHERE R IS THE RADIUS OF A CIRCUMCIRCLE
C             AND D IS THE DISTANCE FROM THE CIRCUMCENTER
C             TO THE NEAREST NODE.  A REASONABLE VALUE
C             FOR TOL IS 10*EPS WHERE EPS IS THE MACHINE
C             PRECISION.  THE TEST IS EFFECTIVELY BYPASSED
C             BY MAKING TOL LARGER THAN THE DIAMETER OF THE
C             CONVEX HULL OF THE NODES.
C
C       LUN - LOGICAL UNIT NUMBER FOR PRINTING ERROR MES-
C             SAGES.  IF LUN .LT. 1 OR LUN .GT. 99, NO MES-
C             SAGES ARE PRINTED.
C
C INPUT PARAMETERS ARE NOT ALTERED BY THIS ROUTINE.
C
C OUTPUT PARAMETER -
C
C       IER - ERROR INDICATOR
C             IER = -1 IF ONE OR MORE NULL TRIANGLES (AREA =
C                      0) ARE PRESENT BUT NO (OTHER) ERRORS
C                      WERE ENCOUNTERED.  A NULL TRIANGLE IS
C                      AN ERROR ONLY IF IT OCCURS IN THE
C                      THE INTERIOR.
C             IER = 0 IF NO ERRORS OR NULL TRIANGLES WERE
C                     ENCOUNTERED.
C             IER = 1 IF N .LT. 3 OR TOL .LT. 0.
C             IER = 2 IF AN IEND OR IADJ ENTRY IS OUT OF
C                     RANGE.
C             IER = 3 IF THE TRIANGULATION PARAMETERS (NB,
C                     NT, AND NA) ARE INCONSISTENT.
C             IER = 4 IF A TRIANGLE CONTAINS A NODE INTERIOR
C                     TO ITS CIRCUMCIRCLE.
C             THE ERROR SITUATIONS ARE TESTED IN THE ORDER
C               DEFINED BY THE (POSITIVE) IER VALUES.
C
C MODULE REQUIRED BY TRMTST - CIRCUM
C
C***********************************************************
C
      LOGICAL RITE
C
C SET LOCAL VARIABLES, TEST FOR ERRORS IN INPUT, AND
C   INITIALIZE COUNTS.
C
      NN = N
      RTOL = TOL
      RITE = LUN .GE. 1  .AND.  LUN .LE. 99
      IF (NN .LT. 3  .OR.  RTOL .LT. 0.) GO TO 5
      NB = 0
      NT = 0
      NULL = 0
      NFAIL = 0
C
C LOOP ON TRIANGLES (N1,N2,N3) SUCH THAT N2 AND N3 INDEX
C   ADJACENT NEIGHBORS OF N1 AND ARE BOTH LARGER THAN N1
C   (TRIANGLES ARE ASSOCIATED WITH THEIR SMALLEST INDEX).
C
      INDF = 1
      DO 4 N1 = 1,NN
        INDL = IEND(N1)
        IF (INDL .LT. INDF+2) GO TO 6
        IF (IADJ(INDL) .EQ. 0) NB = NB + 1
C
C   LOOP ON NEIGHBORS OF N1
C
        DO 3 INDX = INDF,INDL
          N2 = IADJ(INDX)
          IF (N2 .LT. 0  .OR.  N2 .GT. NN  .OR.
     .        (INDX .LT. INDL  .AND.  N2 .EQ. 0)) GO TO 7
          IF (INDX .LT. INDL) N3 = IADJ(INDX+1)
          IF (INDX .EQ. INDL) N3 = IADJ(INDF)
          IF (N2 .LT. N1  .OR.  N3 .LT. N1) GO TO 3
          NT = NT + 1
C
C   COMPUTE THE COORDINATES OF THE CIRCUMCENTER OF
C     (N1,N2,N3).
C
          CALL CIRCUM (X(N1),X(N2),X(N3),Y(N1),Y(N2),Y(N3),
     .                 CX,CY,IERR)
          IF (IERR .NE. 0) NULL = NULL + 1
          IF (IERR .NE. 0) GO TO 3
C
C   TEST FOR NODES WITHIN THE CIRCUMCIRCLE.
C
          R = (CX-X(N1))**2 + (CY-Y(N1))**2 - RTOL
          IF (R .LE. 0.) GO TO 3
          DO 1 I = 1,NN
            IF (I .EQ. N1  .OR.  I .EQ. N2  .OR.
     .          I .EQ. N3) GO TO 1
            IF ((CX-X(I))**2 + (CY-Y(I))**2 .LT. R) GO TO 2
    1       CONTINUE
          GO TO 3
C
C   NODE I IS INTERIOR TO THE CIRCUMCIRCLE OF (N1,N2,N3).
C
    2     NFAIL = NFAIL + 1
    3     CONTINUE
    4   INDF = INDL + 1
C
C CHECK PARAMETERS FOR CONSISTENCY AND TEST FOR NFAIL = 0.
C
      NA = (IEND(NN)-NB)/2
      IF (NB .LT. 3  .OR.  NT .NE. 2*NN-NB-2  .OR.
     .    NA .NE. 3*NN-NB-3) GO TO 8
      IF (NFAIL .NE. 0) GO TO 9
C
C NO ERRORS WERE ENCOUNTERED.
C
      IER = 0
      IF (NULL .EQ. 0) RETURN
      IER = -1
      IF (RITE) WRITE (LUN,100) NULL
  100 FORMAT (1H ,10HTRMTST -- ,I5,16H NULL TRIANGLES ,
     .        11HARE PRESENT/1H ,10X,16H(NULL TRIANGLES ,
     .        32HON THE BOUNDARY ARE UNAVOIDABLE)//)
      RETURN
C
C N OR TOL IS OUT OF RANGE.
C
    5 IER = 1
      IF (RITE) WRITE (LUN,110) N, RTOL
  110 FORMAT (1H ,33HTRMTST -- INVALID INPUT PARAMETER/
     .        1H ,10X,4HN = ,I5,8H, TOL = ,E8.1)
      RETURN
C
C IEND ENTRY OUT OF RANGE
C
    6 IER = 2
      IF (RITE) WRITE (LUN,120) N1
  120 FORMAT (1H ,15HTRMTST -- NODE ,I5,
     .        26H HAS LESS THAN 3 NEIGHBORS)
      RETURN
C
C IADJ ENTRY OUT OF RANGE
C
    7 IER = 2
      IF (RITE) WRITE (LUN,130) INDX
  130 FORMAT (1H ,33HTRMTST -- IADJ(K) IS NOT A VALID ,
     .        14HINDEX FOR K = ,I5)
      RETURN
C
C INCONSISTENT TRIANGULATION PARAMETERS
C
    8 IER = 3
      IF (RITE) WRITE (LUN,140) NB, NT, NA
  140 FORMAT (1H ,33HTRMTST -- INCONSISTENT PARAMETERS/
     .        1H ,10X,I5,15H BOUNDARY NODES,3X,I5,
     .        10H TRIANGLES,3X,I5,5H ARCS)
      RETURN
C
C CIRCUMCIRCLE TEST FAILURE
C
    9 IER = 4
      IF (RITE) WRITE (LUN,150) NFAIL
  150 FORMAT (1H ,10HTRMTST -- ,I5,15H CIRCUMCIRCLES ,
     .        32HCONTAIN NODES IN THEIR INTERIORS)
      RETURN
      END
