      SUBROUTINE DTGGRD(X,Y,NP,IP,W,TRIANG,MT,B,MB,NT,INFO)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c  File: dtggrd.[f|for]
c  contains dtggrd, dtgang, dtgadj, dtgupd
c
c>> 1997-07-01 DTGGRD Krogh Reversed subscripts in B (CLL suggestion)
c>> 1997-06-14 DTGGRD Krogh Removed implicit none, added external stmt.
c>> 1997-06-14 DTGGRD CLL Argument INFO() replaces NBUSED, MODE, & OK1.
C>> 1997-06-09 DTGGRD  Krogh  Converted SFTRAN code to Fortran.
C>> 1995-09-26 DTGGRD CLL Editing for inclusion into MATH77.
c>> 1991-11-20 DTGGRD CLL Using Math77 subroutine _SORTP for sorting.
C
C         THIS SUBR CONSTRUCTS A LIST OF POINTERS IN THE ARRAY
C    TRIANG() DEFINING A TRIANGULAR GRID HAVING THE GIVEN (X,Y)
C    DATA AS VERTICES.
C         THE GRID WILL BE OPTIMAL IN THE SENSE OF THE MAX-MIN ANGLE
C    CRITERION.  IT HAS BEEN SHOWN, INDEPENDENTLY BY LAWSON AND BY
C    R. SIBSON, THAT THIS IS EQUIVALENT TO SOME OTHER OPTIMALLITY
C    CRITERIA.  THE GRID PRODUCED IS A DELAUNAY TRIANGULATION AND
C    IS THE GRAPH THEORETIC DUAL OF A TESSELLATION MADE UP OF
C    PROXIMITY REGIONS ASSOCIATED WITH THE NAMES DIRICHLET,
C    VORONOI, AND THIESSEN.
C
C    C. L. LAWSON, JPL, 1976 NOV 28
C     C.L.L., 1979 MAR 3. CONVERTED TO SFTRAN3.
C     C.L.L., 1979 JUL 22. ADDED CALL TO DTGSIZ.
C
C  X(1:NP),Y(1:NP) [in]    (X,Y) DATA POINTS.
C  NP [in]                    NO. OF DATA POINTS.  REQUIRE NP .GE. 3.
C
C  IP() [scratch]   WORK SPACE OF LENGTH at least NP.
C  W()  [scratch]   WORK SPACE OF LENGTH at least NP.
C  TRIANG() [out]  ARRAY IN WHICH THIS SUBR WILL BUILD A LIST DEFINING
C             A triangular grid.  access to this list is via the
C             four subrs DTGGET, DTGPUT, DTGSIZ, and DTGSET.
C  MT [in]   DIMENSION OF ARRAY TRIANG().  NO. OF TRIANGLES IT
C            CAN HOLD, NTLIM, IS DETERMINED BY USE OF SUBROUTINE
C            DTGSIZ.    NTLIM WILL BE MT/6 IF POINTERS ARE NOT
C            PACKED.  NTLIM = (3*MT)/6 IF POINTERS ARE PACKED
C            THREE PER WORD.
c  Bdry(1:4, MB) [integer, out]  Array containing pointers defining the
c      boundary of the (convex) triangular grid.
C         Bdry(1, K) = FWD POINTER. Points to next vertex in
c                     counterclockwise order.
C         Bdry(2, K) = BACKWARD POINTER. Points to next vertex in
c                     clockwise order.
C         Bdry(3, K) = A BOUNDARY POINT
C         Bdry(4, K) = A BOUNDARY TRIANGLE
c      The triangle Bdry(4, K) has a boundary edge that
c      connects the points Bdry(3, K) and Bdry(3, K+1).
c
c      On return from this subroutine, some elements of the array Bdry()
c      may not be members of the linked list defining the boundary. The
c      entry with K = 1 will always be a member of the boundary list.
c      Thus for another subroutine to scan the Bdry() list made by this
c      subroutine it should start with K = 1 and follow the forward or
c      backward pointers.
c
C  MB [in]   FIRST DIMENSION OF B(,).  SETTING MB = NP+1 WILL
C            ALWAYS BE ADEQUATE, HOWEVER A MUCH SMALLER VALUE
C            WILL OFTEN SUFFICE.  SUGGEST TRYING MB = 6 * (cube root
C            of NP) WHEN NP .GE. 64.
C  NT [out]  OUTPUT..  NO. OF TRIANGLES
c
c  INFO(1:3) [out, integer]  Termination status information.
c
c     INFO(1) INDICATES STATUS ON TERMINATION.
C              0  NORMAL TERMINATION.  The triangular grid is complete
c                 and successfully optimized.
C              1  ERROR. The triangular grid is complete but not
c                 successfully optimized.  Apparent looping in DTGADJ.
C                 GRID MAY NOT BE OPTIMAL.
C              2  ERROR.  ALL GIVEN POINTS ARE COLINEAR.  NO TRIANGLES
C                 CONSTRUCTED.  (OR ELSE THERE ARE SOME DUPLICATE PTS)
C              3  ERROR.  DUPLICATE POINTS.  TRIANGULAR GRID NOT
C                 COMPLETED.
C              4  ERROR.  NOT ENOUGH SPACE IN B(,).  USER MUST
C                 INCREASE MB.
C              5  NOT ENOUGH SPACE IN TRIANG().  USER MUST MAKE
C                 MT BIGGER.
C     INFO(2)  NO. OF BOUNDARY POINTS.
C     INFO(3)  SMALLEST SETTING OF MB THAT WOULD SUFFICE FOR THIS DATA.
c
c        [Relations to internal variables or previously used variable
c         names:  INFO(1) = MODE, INFO(2) = NB, INFO(3) = MBUSED,
c         The value of INFO(1) = MODE will be  set to 0 when
c         OK1 = .true. and 1 when OK1 = .false.]
C
C                  USE SUBR DTGSIZ TO COMPUTE NTLIM, THE MAX NO. OF
C                 TRIANGLES, GIVEN THE DIMENSION PARAMETER MT.
C
c     ------------------------------------------------------------------
c--D replaces "?": ?TGGRD, ?TGANG, ?TGSET, ?TGGET, ?TGPUT, ?TGSIZ,
c--&               ?TGADJ, ?SORTP, ?TGUPD
c     ------------------------------------------------------------------
      INTEGER NP, mb, mt
      INTEGER AVAIL, B(4, MB), ext, TRI(6)
      integer i, i1, i2, i3, i3prev, index, INFO(3), ip(np), isave
      integer jhi, jlo, jnew
      integer k, kb1, kbhi, kblo, kbnew, kbz
      integer mode
      integer nb, nt, ntlim, mbused,  TRIANG(MT)
      double precision aref, cross
      external DTGANG
      double precision DTGANG
      double precision dx12, dy12, eight, pangle
      DOUBLE PRECISION W(np), X(NP), xmean, Y(NP), ymean, zero
        parameter(zero = 0.0d0, eight = 8.0d0)
      LOGICAL  OK, OK1
c     ------------------------------------------------------------------
      call DTGSIZ( MT, NTLIM )
      MODE = 0
      OK1 = .TRUE.
C
C        FIND AN EXTREME POINT (X(EXT),Y(EXT))
C
      EXT=1
      do 10 I = 1, NP
         if ( X(I) .lt. X(EXT)
     *     .or. (X(I) .eq. X(EXT)  .and.  Y(I) .lt. Y(EXT) )) EXT=I
   10 continue
C
C             COMPUTE SQUARED DISTANCES FROM (X(EXT),Y(EXT)).
C             ALSO INITIALIZE IP().
      do 20 I = 1, NP
         W(I) = (X(I)-X(EXT))**2 + (Y(I)-Y(EXT))**2
   20 continue
c
c        Use _SORTP from the JPL MATH77 library to
c        determine a sorted order for the contents of W() and set
c        indices indicating the sorted order in IP().
c        IP() will be set so that W(IP(1)) is the smallest number in
c        W(), W(IP(2)) is the next larger number, etc.  The contents
c        of IP() will be used as indices into the X() and Y() arrays.
c        The contents of W() will not be reordered by _SORTP, however
c        this does not matter here since the next time we reference W()
c        will be to store other things into it.
c
      call DSORTP(W, 1, NP, IP)
C
C             BUILD THE FIRST TRIANGLE USING POINTS IP(1), IP(2), AND
C             THE NEXT POINT THAT IS NOT COLINEAR WITH THESE TWO.
C
      I1=IP(1)
      I2=IP(2)
      DX12= X(I2)-X(I1)
      DY12= Y(I2)-Y(I1)
      I3PREV = 0
      do 40 I = 3, NP
         I3 = IP(I)
         IP(I) = I3PREV
         CROSS =-(X(I3)-X(I2))*DY12 + (Y(I3)-Y(I2))*DX12
         if (CROSS .ne. 0.) then
            if (CROSS .lt. 0.) then
C                                 SWAP I2 AND I3
               ISAVE=I2
               I2   =I3
               I3   =ISAVE
            end if
            go to 50
         end if
         I3PREV=I3
   40 continue
C
C             ERROR IF DROP THRU TO HERE.  ALL NP POINTS ARE COLINEAR
C             OR ELSE POINTS I1 AND I2 ARE IDENTICAL.
C
      MODE = 2
      go to 999
C
C
   50 continue
C                  HERE THE VERTICES OF THE FIRST TRIANGLE ARE I1,I2,I3
C
C                  COMPUTE MEAN OF VERTICES OF FIRST TRIANGLE.
C
      XMEAN = (X(I1)+X(I2)+X(I3))/3.
      YMEAN = (Y(I1)+Y(I2)+Y(I3))/3.
C
C             COMPUTE THE PSEUDOANGLE OF VECTOR FROM (XMEAN,YMEAN)
C             TO I1.
C
      AREF = DTGANG(X(I1)-XMEAN, Y(I1)-YMEAN, zero)
C
C        INITIALIZE BOUNDARY POINTER STRUCTURE B(,).
C        B(1, K) = FWD POINTER
C        B(2, K) = BACKWARD POINTER
C        B(3, K) = A BOUNDARY POINT
C        B(4, K) = A BOUNDARY TRIANGLE
C        W(K)   = PSEUDOANGLE OF POINT B(3, K)
C
      do 70 K = 1, MB
         B(1, K) = K+1
   70 continue
      B(1, MB) = 0
      B(1, 4)  = 0
      AVAIL   = 5
      do 80 K = 1, 4
         B(2, K) = K-1
         B(4, K) = 1
   80 continue
      B(4, 4) = 0
C
      B(3, 1) = I1
      B(3, 2) = I2
      B(3, 3) = I3
      B(3, 4) = I1
C
      W(1) = zero
      W(2) = DTGANG(X(I2)-XMEAN,Y(I2)-YMEAN,AREF)
      W(3) = DTGANG(X(I3)-XMEAN,Y(I3)-YMEAN,AREF)
      W(4) = eight
      NB = 3
      MBUSED = 3
C
C                  BUILD FIRST TRIANGLE IN TRIANG().
C
      if (NTLIM .le. 0) then
         MODE = 5
         go to 999
      end if
      call DTGSET(1,0,0,0,I1,I2,I3,TRIANG,MT)
      NT=1
C    ******************************************************************
C
C                  BEGIN MAIN LOOP.  ADD ONE NEW POINT AT A TIME INTO
C                  THE TRIANGULAR GRID STRUCTURE.
C
C    ******************************************************************
      KBNEW=2
      do 200 INDEX = 4, NP
         JNEW = IP(INDEX)
         PANGLE = DTGANG(X(JNEW)-XMEAN, Y(JNEW)-YMEAN, AREF)
C
C                  USE  PANGLE  AS SEARCH KEY IN THE TABLE OF
C                  PSEUDO ANGLES OF BOUNDARY POINTS.
C                  START SEARCH AT PREVIOUS NEW BOUNDARY POINT.
C
         KB1= KBNEW
         if ( PANGLE .gt. W(KB1)) then
            KBHI = KB1
  100       KBLO = KBHI
            KBHI = B(1, KBLO)
            if ( PANGLE .gt. W(KBHI) ) go to 100
         else
C                            HERE WE HAVE  PANGLE .le. W(KB1)
            KBLO = KB1
  110       KBHI = KBLO
            KBLO = B(2, KBHI)
            if ( PANGLE .lt. W(KBLO) ) go to 110
         end if
C
         JLO = B(3, KBLO)
         JHI = B(3, KBHI)
C
C                            TEST FOR DUPLICATE POINTS
C
         if ( PANGLE .eq. W(KBLO) ) then
            if ( X(JNEW) .eq. X(JLO) .and. Y(JNEW).eq.Y(JLO)) then
               MODE = 3
               go to 999
            end if
         else if (PANGLE .eq. W(KBHI) ) then
            if (X(JNEW).eq.X(JHI) .and. Y(JNEW).eq.Y(JHI)) then
               MODE = 3
               go to 999
            end if
         end if
C
C                            Attach point JNEW to JHI and JLOW.
         if (AVAIL .eq. 0) then
C                                 ERROR..  INSUFFICIENT STORAGE IN B()
C                                 TO REPRESENT THE BOUNDARY.
            MODE = 4
              go to 999
           endif
         if (NT   .ge. NTLIM) then
C                                 ERROR..  INSUFFICIENT STORAGE IN
C                                 TRIANG()
            MODE=5
            go to 999
         endif
C                                 UPDATE TRIANG().
            NT=NT+1
            call DTGUPD(B(3, KBLO), B(4, KBLO), TRI, TRIANG, MT, NT)
            call DTGSET(NT,B(4, KBLO),0,0,JHI,JLO,JNEW,TRIANG,MT)
C
C                                 UPDATE B() AND W().
            NB = NB+1
            if (NB .gt. MBUSED) MBUSED = NB
            KBNEW = AVAIL
            AVAIL = B(1, AVAIL)
            B(1, KBNEW) = KBHI
            B(2, KBNEW) = KBLO
            B(3, KBNEW) = JNEW
            B(4, KBNEW) = NT
            W(KBNEW)   = PANGLE
            B(1, KBLO)  = KBNEW
            B(2, KBHI)  = KBNEW
            B(4, KBLO)  = NT
C
            call DTGADJ(NT,1,X,Y,NP,TRIANG,MT,B,MB, KBNEW, OK)
            OK1 = OK1 .and. OK

C  Done attaching point JNEW to JHI and JLOW.
C  We have just assigned JNEW the index KBNEW in the boundary list.
C
C             LOOP BACKWARD ALONG BOUNDARY FROM KBNEW TO FIND OTHER
C             POINTS TO WHICH POINT JNEW CAN BE CONNECTED.
C
         KBHI= KBLO
         KBLO = B(2, KBHI)
  130    if (KBLO .ne. 0) then
            JHI = B(3, KBHI)
            JLO = B(3, KBLO)
            CROSS = -(Y(JHI)-Y(JNEW)) * (X(JLO)-X(JHI))
     *              +(X(JHI)-X(JNEW)) * (Y(JLO)-Y(JHI))
            if (CROSS .gt. 0.) then
C                  CONNECT POINT JNEW TO JLO
               if (NT .ge. NTLIM) then
C                    ERROR..  INSUFFICIENT SPACE IN TRIANG().
C
                  MODE=5
                    go to 999
               endif
C                              UPDATE TRIANG()
                  NT=NT+1
                  call DTGUPD(B(3,KBLO), B(4,KBLO), TRI, TRIANG, MT, NT)
                  call DTGUPD(B(3,KBHI), B(4,KBHI), TRI, TRIANG, MT, NT)
                  call DTGSET(NT, B(4, KBLO), 0, B(4, KBHI), JHI, JLO,
     1               JNEW, TRIANG, MT)
C
C                              UPDATE B()
                  B(1, KBLO) = KBNEW
                  B(2, KBNEW)= KBLO
                  B(4, KBLO) = NT
                  NB=NB-1
                  B(1, KBHI) = AVAIL
                  AVAIL     = KBHI
C
                  call DTGADJ(NT,1,X,Y,NP,TRIANG,MT,B,MB, KBNEW, OK)
                  OK1 = OK1 .and. OK
c                     Done connecting point JNEW to JLO
c
                  KBHI = KBLO
                  KBLO = B(2, KBHI)
                  go to 130

            end if
         end if
C
C          LOOP FORWARD ALONG BOUNDARY FROM KBNEW TO FIND OTHER
C          POINTS TO WHICH POINT JNEW CAN BE CONNECTED.
C
         KBLO = B(1, KBNEW)
         KBHI = B(1, KBLO)
  150    if (KBHI.ne.0) then
            JHI= B(3, KBHI)
            JLO= B(3, KBLO)
            CROSS = -(Y(JLO)-Y(JHI)) * (X(JNEW)-X(JLO))
     *         +(X(JLO)-X(JHI)) * (Y(JNEW)-Y(JLO))
            if (CROSS .gt. 0) then
               if (NT .ge. NTLIM) then
C                       ERROR..  INSUFFICIENT SPACE IN TRIANG().
                  MODE = 5
                    go to 999
               endif
C          Connect point JNEW to JHI -- Update TRIANG()
                  NT=NT+1
                  call DTGUPD(B(3, KBLO),B(4, KBLO),TRI,TRIANG,MT,NT)
                  call DTGUPD(B(3, KBNEW),B(4, KBNEW),TRI,TRIANG,MT,NT)
                  call DTGSET(NT, B(4, KBLO), B(4, KBNEW), 0, JHI, JLO,
     1               JNEW, TRIANG, MT)
C
C                                 UPDATE B()
                  B(1, KBNEW) = KBHI
                  B(2, KBHI)  = KBNEW
                  B(4, KBNEW) = NT
                  NB=NB-1
                  B(1, KBLO)  = AVAIL
                  AVAIL      = KBLO
C
                  call DTGADJ(NT,1,X,Y,NP,TRIANG,MT,B,MB, KBNEW, OK)
                  OK1 = OK1 .and. OK
C                                 Done connecting point JNEW to JHI
c
                  KBLO = KBHI
                  KBHI = B(1, KBLO)
                  go to 150

            end if
         end if
  200 continue
      MBUSED = MBUSED + 1
C
  999 continue
c
c        As built, the linked boundary list runs (counterclockwise) from
c        b(*, 1) to b(*, 4), with b(*, 1) and b(*, 4) representing the
c        same point.  The back pointer in b(*, 1) and the forward
c        pointer in b(*, 4) are both zero, indicating ends of the list.
c        Change this so the linked list is fully circular, still
c        including b(*, 1), but detaching b(*, 4).  Let kbz be the value
c        of the back pointer in b(*, 4).  Set the back pointer in b(*,
c        1) to be kbz, and the forward pointer in b(*, kbz) to be 1.
c
      kbz = b(2, 4)
      b(2, 1) = kbz
      b(1, kbz) = 1

      if(MODE .eq. 0 .and. .not. OK1) MODE = 1
      INFO(1) = MODE
      INFO(2) = NB
      INFO(3) = MBUSED
      end
c     ==================================================================

      subroutine DTGUPD( JFIX, KB4, TRI, TRIANG, MT, NT)
c               Updates TRIANG, using JFIX and KB4
      integer JFIX, KB4, MT, NT, TRI(6), TRIANG(MT)
      integer I
c
      call DTGGET(KB4,TRI,TRIANG)
      do 10 I = 1, 3
         if (TRI(I+3) .eq. JFIX) then
            TRI(I) = NT
            go to 20
         end if
   10 continue
   20 call DTGPUT(KB4, TRI, TRIANG, MT)
      return
      end
c     ==================================================================
      DOUBLE PRECISION FUNCTION DTGANG  (XX,YY,AREF)
C>> 1995-09-26 CLL Editing for inclusion into MATH77.
C
c     Compute the pseudoangle  DTGANG  between the reference direction
c     whose pseudoangle is  AREF  and the vector  (XX,YY).
c     Pseudoangle is measured counterclockwise with period 8.
c     Require  0. .le. AREF .lt. 8. on input.
c     On return DTGANG will be in the interval 0. .le. DTGANG .lt. 8.
C
c     C.L.LAWSON, JPL, 1976 NOV 8
C
c     ------------------------------------------------------------------
      double precision xx,yy,aref,a
      double precision eight, four, two, zero
      parameter(eight = 8.0d0, four = 4.0d0, two = 2.0d0, zero = 0.0d0)
c     ------------------------------------------------------------------
      if (abs(xx) .gt. abs(yy)) then
         a = yy/xx
         if (xx .lt. zero) then
            a=a+four
         else
            if (yy .lt. zero) then
               a=a+ eight
            end if
         end if
      else
         if (yy .eq. zero) then
            a = zero
         else
            a = two - xx/yy
            if (yy .lt. zero) a = a+four
         end if
      end if
c     .                            A IS NOW THE PSEUDOANGLE OF (XX,YY)
c     .                            RELATIVE TO (1.,0.) AND LYING IN
c     .                            THE INTERVAL   0. .le. A .lt. 8.
      dtgang = a - aref
      if (dtgang .lt. zero) dtgang = dtgang + eight
      return
      end
c     ==================================================================
      SUBROUTINE DTGADJ(TIN,NIN,X,Y,NP,TRIANG,MT,B,MB,KBNEW, OK)
C>> 1995-09-26 CLL Editing for inclusion into MATH77.
C
C     C.L.LAWSON, JPL, 1976 NOV 30, Changed 1977 APR 7
C     Changed 1977 June 16 to get rid of internal arrays
c     used to stack triangles to be tested.
c
c     Given triangle TIN and the neighboring triangle that is in
c     position NIN relative to TIN,  Test for possible exchange of
c     their common edge.
c     Use the criterion of maximizing the smallest angle, which is the
c     same as the criterion of the empty circumcircle.
c     If the swap is made then certain other edges must be tested.
c     This subr does the tests, swaps edges as needed, and does
c     tests and swaps on other edges as needed.
c     The counter COUNT is used to guard against infinite looping
c     which should never happen anyway.  If count exceeds 100 the
c     subr sets OK = .false. indicating that the triangulation
c     produced may not be optimal.
c     ------------------------------------------------------------------
C                       Subroutine Arguments
c
c  TIN [integer, in]
c  NIN [integer, in]
c  X(),Y() [floating, in]
c  TRIANG() [integer, inout]
c  MT [integer, in]
c  B(1:4, ) [integer, inout]
c  MB [integer, in]
c  KBNEW [integer, in]  Index into the boundary array, B(,).
c  OK [logical, out]
c     ------------------------------------------------------------------
      integer mb
      INTEGER  ADD1(3)
      INTEGER  B(4, MB), COUNT
      integer i, k, kbkwd,KBNEW, kfwd
      integer more, mt
      integer nabor1, nabor2, nin, np
      INTEGER  QV1,QV2,QV3,QV4,QV1S(3),QV2S(3),QV3S(3),SUB1(3)
      INTEGER  t1, T11,T12,t2, T21,T22, tin
      INTEGER  TNEW, TOLD, TRI1(6),TRI2(6), TRIANG(MT)
      double precision a123, a234, a341, a412
      double precision diag13, diag24
      double precision dx12, dx23, dx34, dx41, dy12, dy23, dy34, dy41
      double precision s12, s23, s34, s41
      double precision tau13, tau24
      DOUBLE PRECISION X(NP),Y(NP)
      LOGICAL  OK, SWAP
      DATA QV1S(1),QV1S(2),QV1S(3)/ 5, 6, 4/
      DATA QV2S(1),QV2S(2),QV2S(3)/ 6, 4, 5/
      DATA QV3S(1),QV3S(2),QV3S(3)/ 4, 5, 6/
      DATA ADD1(1),ADD1(2),ADD1(3)/ 2, 3, 1/
      DATA SUB1(1),SUB1(2),SUB1(3)/ 3, 1, 2/
c     ------------------------------------------------------------------
      T1=TIN
      NABOR1 = NIN
      COUNT = 1
      MORE = 1
      OK = .TRUE.
   20 if ((MORE .gt. 0) .and. OK) then

         call DTGGET(T1,TRI1,TRIANG)
         T2 = TRI1(NABOR1)
         if ( T2 .ne. 0) then
C T2=0 MEANS SIDE NABOR OF TRIANGLE T1 IS A BOUNDARY EDGE SO THERE IS NO
C TRIANGLE T2.  HERE WE HAVE  T2 .ne. 0
            call DTGGET(T2,TRI2,TRIANG)
C  SET QV1,...,QV4 AS INDICES OF THE 4 VERTICES OF THE QUADRILATERAL
C  FORMED BY THE 2 TRIANGLES T1 AND T2
            I=QV1S(NABOR1)
            QV1=TRI1(I)
            I=QV2S(NABOR1)
            QV2=TRI1(I)
            I=QV3S(NABOR1)
            QV3=TRI1(I)
            do 40 NABOR2 = 1, 3
               if (TRI2(NABOR2) .eq. T1) then
                  I=QV2S(NABOR2)
                  QV4=TRI2(I)
                  go to 50
               end if
   40       continue
   50       continue

C The max-min angle test below is made using the max-min angle criterion
c as described in JPL internal TECH.  MEMO.  914-299, FEB., 1972 by
c C. L.  LAWSON.
C
            DX23 = X(QV3)-X(QV2)
            DY23 = Y(QV3)-Y(QV2)
            DX34 = X(QV4)-X(QV3)
            DY34 = Y(QV4)-Y(QV3)
            A234 = DX23*DY34 - DY23*DX34
            if (A234 .le. 0.) then
               SWAP = .FALSE.
            else
               DX12 = X(QV2)-X(QV1)
               DY12 = Y(QV2)-Y(QV1)
               DX41 = X(QV1)-X(QV4)
               DY41 = Y(QV1)-Y(QV4)
               A412 = DX41*DY12 - DY41*DX12
               if (A412 .le. 0.) then
                  SWAP = .FALSE.
               else
                  A341 = DX34*DY41 - DY34*DX41
                  A123 = A234 + A412 - A341
                  S12 = DX12**2 + DY12**2
                  S23 = DX23**2 + DY23**2
                  S34 = DX34**2 + DY34**2
                  S41 = DX41**2 + DY41**2
                  DIAG13 = (X(QV3)-X(QV1))**2 + (Y(QV3)-Y(QV1))**2
                  DIAG24 = (X(QV4)-X(QV2))**2 + (Y(QV4)-Y(QV2))**2
                  TAU13 = min( A123**2 / max(S12,S23),
     *               A341**2 / max(S34,S41)) / DIAG13
                  TAU24 = min( A234**2 / max(S23,S34),
     *               A412**2 / max(S41,S12)) / DIAG24
                  SWAP = TAU24 .gt. TAU13
               end if
            end if
c End of max-min angle test
         else
            SWAP = .FALSE.
         end if
         if ( SWAP) then
c  Remove the edge connecting QV1 and QV3.  Insert new edge connecting
C  QV2 and QV4.  Redefine triangle T1 to have vertices QV1, QV2, and
C  QV4.  Redefine triangle T2 to have vertices QV2, QV3, AND QV4.
            I = ADD1(NABOR1)
            T11 = TRI1(I)
            I = ADD1(I)
            T12 = TRI1(I)
C
            I   = ADD1(NABOR2)
            T21 = TRI2(I)
            I   = ADD1(I)
            T22 = TRI2(I)
C
            call DTGSET(T1,T11,T2,T22,QV1,QV2,QV4,TRIANG,MT)
            call DTGSET(T2,T12,T21,T1,QV2,QV3,QV4,TRIANG,MT)
C
            TNEW=T2
            if ( T12 .ne. 0) then
               call DTGGET( T12,TRI2, TRIANG)
               do 70 I = 1, 3
                  if (TRI2(I) .eq. T1) then
                     TRI2(I) = TNEW
                     go to 80
                  end if
   70          continue
   80          call DTGPUT( T12,TRI2, TRIANG, MT)
            else
C The following test of (QV2 .eq.  B(3, KBNEW)) will always be true if
C the outer algorithm calling this subroutine is the one originally
C planned, always adding new points from the exterior, not from the
C interior.  The false branch is provided to handle the more general
c case.
               if (QV2 .eq. B(3, KBNEW) ) then
                  B(4, KBNEW) = TNEW
               else
C                 WRITE(6,1001) QV2,B(3, KBNEW)
c Search alternately fwd and bkwd from KBNEW in the boundary list B()
c for the point QV2.  When found set the associated triangle to TNEW.
                  KFWD = KBNEW
                  KBKWD= B(2, KFWD )
   90             continue
                     if (KFWD .ne. 0) then
                        if (B(3, KFWD) .eq. QV2) then
                           B(4, KFWD) = TNEW
                           go to 100
                        end if
                        KFWD = B(1, KFWD)
                     end if
C
                     if (KBKWD .ne. 0) then
                        if (B(3, KBKWD) .eq. QV2) then
                           B(4, KBKWD) = TNEW
                           go to 100
                        end if
                        KBKWD = B(2, KBKWD)
                     end if
                  go to 90
               end if
            end if
C
  100       TNEW=T1
            if ( T22 .ne. 0) then
               call DTGGET( T22,TRI2, TRIANG)
               DO 110 I = 1, 3
                  if (TRI2(I) .eq. T2) then
                     TRI2(I) = TNEW
                     go to 120
                  end if
  110          continue
  120          call DTGPUT( T22,TRI2, TRIANG, MT)
            else
C              WRITE(6,1002) QV4
c Search alternately fwd and bkwd from KBNEW in the boundary list B()
c for the point QV2.  When found set the associated triangle to TNEW.
               KFWD = KBNEW
               KBKWD= B(2, KFWD )
  140          continue
                  if (KFWD .ne. 0) then
                     if (B(3, KFWD) .eq. QV4) then
                        B(4, KFWD) = TNEW
                        go to 150
                     end if
                     KFWD = B(1, KFWD)
                  end if
C
                  if (KBKWD .ne. 0) then
                     if (B(3, KBKWD) .eq. QV4) then
                        B(4, KBKWD) = TNEW
                        go to 150
                     end if
                     KBKWD = B(2, KBKWD)
                  end if
               go to 140
            end if
  150       NABOR1 = 3
            MORE = MORE+1
            COUNT = COUNT+1
            OK = COUNT .le. 100
         else
            MORE = MORE-1
            if (MORE .gt. 0) then
               TOLD = T1
               I = SUB1(NABOR1)
               T1=TRI1(I)
               call DTGGET(T1,TRI1,TRIANG)
               do 160 K = 1, 3
                  if (TRI1(K) .eq. TOLD) then
                     NABOR1 = SUB1(K)
                     go to 20
                  end if
  160          continue
            end if
         end if
         go to 20
      end if
      return
C
C1001 FORMAT(33H0ADJUST..  USING FIX-B WITH QV2 =,I5,5X,
C    * 12HB(3, KBNEW) =,I5)
C1002 FORMAT(33H0ADJUST..  USING FIX-B WITH QV4 =,I5)
      end
