tools_stripack.F90

Go to the documentation of this file.

!----------------------------------------------------------------------
! Module: tools_stripack
!> STRIPACK routines
! Source: https://people.sc.fsu.edu/~jburkardt/f_src/stripack/stripack.html
! Author: Robert Renka
! Original licensing: none
! Modified by Benjamin Menetrier for BUMP
! Licensing: this code is distributed under the CeCILL-C license
! Copyright © 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
!----------------------------------------------------------------------
module tools_stripack
 
use tools_kinds, only: kind_real
use tools_repro, only: eq,inf,sup,supeq,small
use type_mpl, only: mpl_type
 
implicit none
 
private
public :: addnod,bnodes,inside,trfind,trlist,trmesh
 
contains
 
subroutine addnod ( mpl, nst, k, x, y, z, list, lptr, lend, lnew, ier )
 
!*****************************************************************************80
!
! Subroutine: addnod
!> Add a node to a triangulation
!
!  Discussion:
!
!    This subroutine adds node K to a triangulation of the
!    convex hull of nodes 1, ..., K-1, producing a triangulation
!    of the convex hull of nodes 1, ..., K.
!
!    The algorithm consists of the following steps:  node K
!    is located relative to the triangulation (TRFIND), its
!    index is added to the data structure (INTADD or BDYADD),
!    and a sequence of swaps (SWPTST and SWAP) are applied to
!    the arcs opposite K so that all arcs incident on node K
!    and opposite node K are locally optimal (satisfy the circumcircle test).
!
!    Thus, if a Delaunay triangulation of nodes 1 through K-1 is input,
!    a Delaunay triangulation of nodes 1 through K will be output.
!
!  Modified:
!
!    15 May 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer NST, the index of a node at which TRFIND
!    begins its search.  Search time depends on the proximity of this node to
!    K.  If NST < 1, the search is begun at node K-1.
!
!    Input, integer K, the nodal index (index for X, Y, Z, and
!    LEND) of the new node to be added.  4 <= K.
!
!    Input, real ( kind_real ) X(K), Y(K), Z(K), the coordinates of the nodes.
!
!    Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(K),
!    LNEW.  On input, the data structure associated with the triangulation of
!    nodes 1 to K-1.  On output, the data has been updated to include node
!    K.  The array lengths are assumed to be large enough to add node K.
!    Refer to TRMESH.
!
!    Output, integer IER, error indicator:
!     0 if no errors were encountered.
!    -1 if K is outside its valid range on input.
!    -2 if all nodes (including K) are collinear (lie on a common geodesic).
!     L if nodes L and K coincide for some L < K.
!
!  Local parameters:
!
!    B1,B2,B3 = Unnormalized barycentric coordinates returned by TRFIND.
!    I1,I2,I3 = Vertex indexes of a triangle containing K
!    IN1 =      Vertex opposite K:  first neighbor of IO2
!               that precedes IO1.  IN1,IO1,IO2 are in
!               counterclockwise order.
!    IO1,IO2 =  Adjacent neighbors of K defining an arc to
!               be tested for a swap
!    IST =      Index of node at which TRFIND begins its search
!    KK =       Local copy of K
!    KM1 =      K-1
!    L =        Vertex index (I1, I2, or I3) returned in IER
!               if node K coincides with a vertex
!    LP =       LIST pointer
!    LPF =      LIST pointer to the first neighbor of K
!    LPO1 =     LIST pointer to IO1
!    LPO1S =    Saved value of LPO1
!    P =        Cartesian coordinates of node K
!
  implicit none
 
  type(mpl_type),intent(inout) :: mpl !< MPI data
 
  integer k
 
  real ( kind_real ) b1
  real ( kind_real ) b2
  real ( kind_real ) b3
  integer i1
  integer i2
  integer i3
  integer ier
  integer in1
  integer io1
  integer io2
  integer ist
  integer kk
  integer km1
  integer l
  integer lend(k)
  integer list(*)
  integer lnew
  integer lp
  integer lpf
  integer lpo1
  integer lpo1s
  integer lptr(*)
  integer nst
  real ( kind_real ) p(3)
  logical output
  real ( kind_real ) x(k)
  real ( kind_real ) y(k)
  real ( kind_real ) z(k)
  character(len=1024),parameter :: subr = 'addnod'
 
  kk = k
 
  if ( kk < 4 ) then
    ier = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ADDNOD - Fatal error!'
    write ( *, '(a)' ) '  K < 4.'
    call mpl%abort(subr,'stop in stripack')
  end if
!
!  Initialization:
!
  km1 = kk - 1
  ist = nst
  if ( ist < 1 ) then
    ist = km1
  end if
 
  p(1) = x(kk)
  p(2) = y(kk)
  p(3) = z(kk)
!
!  Find a triangle (I1,I2,I3) containing K or the rightmost
!  (I1) and leftmost (I2) visible boundary nodes as viewed
!  from node K.
!
 
  call trfind ( ist, p, km1, x, y, z, list, lptr, lend, b1, b2, b3, &
    i1, i2, i3 )
!
!  Test for collinear or duplicate nodes.
!
  if ( i1 == 0 ) call mpl%abort(subr,'colinear points')
  if ( i3 /= 0 ) then
 
    l = i1
 
    if ( eq(p(1),x(l)) .and. eq(p(2),y(l))  .and. eq(p(3),z(l)) ) then
      ier = l
      return
    end if
 
    l = i2
 
    if ( eq(p(1),x(l)) .and. eq(p(2),y(l))  .and. eq(p(3),z(l)) ) then
      ier = l
      return
    end if
 
    l = i3
    if ( eq(p(1),x(l)) .and. eq(p(2),y(l))  .and. eq(p(3),z(l)) ) then
      ier = l
      return
    end if
 
    call intadd ( kk, i1, i2, i3, list, lptr, lend, lnew )
 
  else
 
    if ( i1 /= i2 ) then
      call bdyadd ( kk, i1,i2, list, lptr, lend, lnew )
    else
      call covsph ( kk, i1, list, lptr, lend, lnew )
    end if
 
  end if
 
  ier = 0
!
!  Initialize variables for optimization of the triangulation.
!
  lp = lend(kk)
  lpf = lptr(lp)
  io2 = list(lpf)
  lpo1 = lptr(lpf)
  io1 = abs( list(lpo1) )
!
!  Begin loop: find the node opposite K.
!
  do
 
    call lstptr ( lend(io1), io2, list, lptr, lp )
 
    if ( 0 <= list(lp) ) then
 
      lp = lptr(lp)
      in1 = abs( list(lp) )
!
!  Swap test:  if a swap occurs, two new arcs are
!  opposite K and must be tested.
!
      lpo1s = lpo1
 
      call swptst ( in1, kk, io1, io2, x, y, z, output )
      if ( .not. output ) then
 
        if ( lpo1 == lpf .or. list(lpo1) < 0 ) then
          exit
        end if
 
        io2 = io1
        lpo1 = lptr(lpo1)
        io1 = abs( list(lpo1) )
        cycle
 
      end if
 
      call swap ( in1, kk, io1, io2, list, lptr, lend, lpo1 )
!
!  A swap is not possible because KK and IN1 are already
!  adjacent.  This error in SWPTST only occurs in the
!  neutral case and when there are nearly duplicate nodes.
!
      if ( lpo1 /= 0 ) then
        io1 = in1
        cycle
      end if
 
      lpo1 = lpo1s
 
    end if
!
!  No swap occurred.  Test for termination and reset IO2 and IO1.
!
    if ( lpo1 == lpf .or. list(lpo1) < 0 ) then
      exit
    end if
 
    io2 = io1
    lpo1 = lptr(lpo1)
    io1 = abs( list(lpo1) )
 
  end do
 
  return
end subroutine addnod
subroutine bdyadd ( kk, i1, i2, list, lptr, lend, lnew )
 
!*****************************************************************************80
!
! Subroutine: bdyadd
!> Add a boundary node to a triangulation
!
!  Discussion:
!
!    This subroutine adds a boundary node to a triangulation
!    of a set of KK-1 points on the unit sphere.  The data
!    structure is updated with the insertion of node KK, but no
!    optimization is performed.
!
!    This routine is identical to the similarly named routine
!    in TRIPACK.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer KK, the index of a node to be connected to
!    the sequence of all visible boundary nodes.  1 <= KK and
!    KK must not be equal to I1 or I2.
!
!    Input, integer I1, the first (rightmost as viewed from KK)
!    boundary node in the triangulation that is visible from
!    node KK (the line segment KK-I1 intersects no arcs.
!
!    Input, integer I2, the last (leftmost) boundary node that
!    is visible from node KK.  I1 and I2 may be determined by TRFIND.
!
!    Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N),
!    LNEW, the triangulation data structure created by TRMESH.
!    Nodes I1 and I2 must be included
!    in the triangulation.  On output, the data structure is updated with
!    the addition of node KK.  Node KK is connected to I1, I2, and
!    all boundary nodes in between.
!
!  Local parameters:
!
!    K =     Local copy of KK
!    LP =    LIST pointer
!    LSAV =  LIST pointer
!    N1,N2 = Local copies of I1 and I2, respectively
!    NEXT =  Boundary node visible from K
!    NSAV =  Boundary node visible from K
!
  implicit none
 
  integer i1
  integer i2
  integer k
  integer kk
  integer lend(*)
  integer list(*)
  integer lnew
  integer lp
  integer lptr(*)
  integer lsav
  integer n1
  integer n2
  integer next
  integer nsav
 
  k = kk
  n1 = i1
  n2 = i2
!
!  Add K as the last neighbor of N1.
!
  lp = lend(n1)
  lsav = lptr(lp)
  lptr(lp) = lnew
  list(lnew) = -k
  lptr(lnew) = lsav
  lend(n1) = lnew
  lnew = lnew + 1
  next = -list(lp)
  list(lp) = next
  nsav = next
!
!  Loop on the remaining boundary nodes between N1 and N2,
!  adding K as the first neighbor.
!
  do
 
    lp = lend(next)
    call insert ( k, lp, list, lptr, lnew )
 
    if ( next == n2 ) then
      exit
    end if
 
    next = -list(lp)
    list(lp) = next
 
  end do
!
!  Add the boundary nodes between N1 and N2 as neighbors of node K.
!
  lsav = lnew
  list(lnew) = n1
  lptr(lnew) = lnew + 1
  lnew = lnew + 1
  next = nsav
 
  do
 
    if ( next == n2 ) then
      exit
    end if
 
    list(lnew) = next
    lptr(lnew) = lnew + 1
    lnew = lnew + 1
    lp = lend(next)
    next = list(lp)
 
  end do
 
  list(lnew) = -n2
  lptr(lnew) = lsav
  lend(k) = lnew
  lnew = lnew + 1
 
  return
end subroutine bdyadd
subroutine bnodes ( n, list, lptr, lend, nodes, nb, na, nt )
 
!*****************************************************************************80
!
! Subroutine: bnodes
!> Return the boundary nodes of a triangulation
!
!  Discussion:
!
!    Given a triangulation of N nodes on the unit sphere created by TRMESH,
!    this subroutine returns an array containing the indexes (if any) of
!    the counterclockwise sequence of boundary nodes, that is, the nodes on
!    the boundary of the convex hull of the set of nodes.  The
!    boundary is empty if the nodes do not lie in a single
!    hemisphere.  The numbers of boundary nodes, arcs, and
!    triangles are also returned.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer N, the number of nodes in the triangulation.
!    3 <= N.
!
!    Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), the
!    data structure defining the triangulation, created by TRMESH.
!
!    Output, integer NODES(*), the ordered sequence of NB boundary
!    node indexes in the range 1 to N.  For safety, the dimension of NODES
!    should be N.
!
!    Output, integer NB, the number of boundary nodes.
!
!    Output, integer NA, NT, the number of arcs and triangles,
!    respectively, in the triangulation.
!
!  Local parameters:
!
!    K =   NODES index
!    LP =  LIST pointer
!    N0 =  Boundary node to be added to NODES
!    NN =  Local copy of N
!    NST = First element of nodes (arbitrarily chosen to be
!          the one with smallest index)
!
  implicit none
 
  integer n
 
  integer i
  integer k
  integer lend(n)
  integer list(6*(n-2))
  integer lp
  integer lptr(6*(n-2))
  integer n0
  integer na
  integer nb
  integer nn
  integer nodes(*)
  integer nst
  integer nt
 
  nn = n
!
!  Search for a boundary node.
!
  nst = 0
 
  do i = 1, nn
 
    lp = lend(i)
 
    if ( list(lp) < 0 ) then
      nst = i
      exit
    end if
 
  end do
!
!  The triangulation contains no boundary nodes.
!
  if ( nst == 0 ) then
    nb = 0
    na = 3 * ( nn - 2 )
    nt = 2 * ( nn - 2 )
    return
  end if
!
!  NST is the first boundary node encountered.
!
!  Initialize for traversal of the boundary.
!
  nodes(1) = nst
  k = 1
  n0 = nst
!
!  Traverse the boundary in counterclockwise order.
!
  do
 
    lp = lend(n0)
    lp = lptr(lp)
    n0 = list(lp)
 
    if ( n0 == nst ) then
      exit
    end if
 
    k = k + 1
    nodes(k) = n0
 
  end do
!
!  Store the counts.
!
  nb = k
  nt = 2 * n - nb - 2
  na = nt + n - 1
 
  return
end subroutine bnodes
subroutine covsph ( kk, n0, list, lptr, lend, lnew )
 
!*****************************************************************************80
!
! Subroutine: covsph
!> Connect an exterior node to boundary nodes, covering the sphere
!
!  Discussion:
!
!    This subroutine connects an exterior node KK to all
!    boundary nodes of a triangulation of KK-1 points on the
!    unit sphere, producing a triangulation that covers the
!    sphere.  The data structure is updated with the addition
!    of node KK, but no optimization is performed.  All
!    boundary nodes must be visible from node KK.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer KK = Index of the node to be connected to the
!    set of all boundary nodes.  4 <= KK.
!
!    Input, integer N0 = Index of a boundary node (in the range
!    1 to KK-1).  N0 may be determined by TRFIND.
!
!    Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N),
!    LNEW, the triangulation data structure created by TRMESH.  Node N0 must
!    be included in the triangulation.  On output, updated with the addition
!    of node KK as the last entry.  The updated triangulation contains no
!    boundary nodes.
!
!  Local parameters:
!
!    K =     Local copy of KK
!    LP =    LIST pointer
!    LSAV =  LIST pointer
!    NEXT =  Boundary node visible from K
!    NST =   Local copy of N0
!
  implicit none
 
  integer k
  integer kk
  integer lend(*)
  integer list(*)
  integer lnew
  integer lp
  integer lptr(*)
  integer lsav
  integer n0
  integer next
  integer nst
 
  k = kk
  nst = n0
!
!  Traverse the boundary in clockwise order, inserting K as
!  the first neighbor of each boundary node, and converting
!  the boundary node to an interior node.
!
  next = nst
 
  do
 
    lp = lend(next)
    call insert ( k, lp, list, lptr, lnew )
    next = -list(lp)
    list(lp) = next
 
    if ( next == nst ) then
      exit
    end if
 
  end do
!
!  Traverse the boundary again, adding each node to K's adjacency list.
!
  lsav = lnew
 
  do
 
    lp = lend(next)
    list(lnew) = next
    lptr(lnew) = lnew + 1
    lnew = lnew + 1
    next = list(lp)
 
    if ( next == nst ) then
      exit
    end if
 
  end do
 
  lptr(lnew-1) = lsav
  lend(k) = lnew - 1
 
  return
end subroutine covsph
subroutine det ( x1, y1, z1, x2, y2, z2, x0, y0, z0, output )
!*****************************************************************************80
!
! Subroutine: det
!> Compute 3D determinant
!
!  DET(X1,...,Z0) >= 0 if and only if (X0,Y0,Z0) is in the
!  (closed) left hemisphere defined by the plane containing (0,0,0),
!  (X1,Y1,Z1), and (X2,Y2,Z2), where left is defined relative to an
!  observer at (X1,Y1,Z1) facing (X2,Y2,Z2).
!
  implicit none
 
  real ( kind_real ) x1
  real ( kind_real ) y1
  real ( kind_real ) z1
  real ( kind_real ) x2
  real ( kind_real ) y2
  real ( kind_real ) z2
  real ( kind_real ) x0
  real ( kind_real ) y0
  real ( kind_real ) z0
  real ( kind_real ) t1
  real ( kind_real ) t2
  real ( kind_real ) t3
  real ( kind_real ) output
 
  t1 = x0*(y1*z2-y2*z1)
  t2 = y0*(x1*z2-x2*z1)
  t3 = z0*(x1*y2-x2*y1)
  output = t1 - t2 + t3
 
  ! Indistinguishability threshold for cross-plateform reproducibility
  if (small(output,t1).or.small(output,t2).or.small(output,t3)) output = 0.0
end subroutine det
subroutine insert ( k, lp, list, lptr, lnew )
 
!*****************************************************************************80
!
! Subroutine: insert
!> Insert K as a neighbor of N1
!
!  Discussion:
!
!    This subroutine inserts K as a neighbor of N1 following
!    N2, where LP is the LIST pointer of N2 as a neighbor of
!    N1.  Note that, if N2 is the last neighbor of N1, K will
!    become the first neighbor (even if N1 is a boundary node).
!
!    This routine is identical to the similarly named routine in TRIPACK.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer K, the index of the node to be inserted.
!
!    Input, integer LP, the LIST pointer of N2 as a neighbor of N1.
!
!    Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LNEW,
!    the data structure defining the triangulation, created by TRMESH.
!    On output, updated with the addition of node K.
!
  implicit none
 
  integer k
  integer list(*)
  integer lnew
  integer lp
  integer lptr(*)
  integer lsav
 
  lsav = lptr(lp)
  lptr(lp) = lnew
  list(lnew) = k
  lptr(lnew) = lsav
  lnew = lnew + 1
 
  return
end subroutine insert
function inside ( p, lv, xv, yv, zv, nv, listv, ier )
 
!*****************************************************************************80
!
! Function: inside
!> Determine if a point is inside a polygonal region
!
!  Discussion:
!
!    This function locates a point P relative to a polygonal
!    region R on the surface of the unit sphere, returning
!    INSIDE = TRUE if and only if P is contained in R.  R is
!    defined by a cyclically ordered sequence of vertices which
!    form a positively-oriented simple closed curve.  Adjacent
!    vertices need not be distinct but the curve must not be
!    self-intersecting.  Also, while polygon edges are by definition
!    restricted to a single hemisphere, R is not so
!    restricted.  Its interior is the region to the left as the
!    vertices are traversed in order.
!
!    The algorithm consists of selecting a point Q in R and
!    then finding all points at which the great circle defined
!    by P and Q intersects the boundary of R.  P lies inside R
!    if and only if there is an even number of intersection
!    points between Q and P.  Q is taken to be a point immediately
!    to the left of a directed boundary edge -- the first
!    one that results in no consistency-check failures.
!
!    If P is close to the polygon boundary, the problem is
!    ill-conditioned and the decision may be incorrect.  Also,
!    an incorrect decision may result from a poor choice of Q
!    (if, for example, a boundary edge lies on the great circle
!    defined by P and Q).  A more reliable result could be
!    obtained by a sequence of calls to INSIDE with the vertices
!    cyclically permuted before each call (to alter the
!    choice of Q).
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, real ( kind_real ) P(3), the coordinates of the point (unit vector)
!    to be located.
!
!    Input, integer LV, the length of arrays XV, YV, and ZV.
!
!    Input, real ( kind_real ) XV(LV), YV(LV), ZV(LV), the coordinates of unit
!    vectors (points on the unit sphere).
!
!    Input, integer NV, the number of vertices in the polygon.
!    3 <= NV <= LV.
!
!    Input, integer LISTV(NV), the indexes (for XV, YV, and ZV)
!    of a cyclically-ordered (and CCW-ordered) sequence of vertices that
!    define R.  The last vertex (indexed by LISTV(NV)) is followed by the
!    first (indexed by LISTV(1)).  LISTV entries must be in the range 1 to LV.
!
!    Output, logical INSIDE, TRUE if and only if P lies inside R unless
!    IER /= 0, in which case the value is not altered.
!
!    Output, integer IER, error indicator:
!    0, if no errors were encountered.
!    1, if LV or NV is outside its valid range.
!    2, if a LISTV entry is outside its valid range.
!    3, if the polygon boundary was found to be self-intersecting.  This
!      error will not necessarily be detected.
!    4, if every choice of Q (one for each boundary edge) led to failure of
!      some internal consistency check.  The most likely cause of this error
!      is invalid input:  P = (0,0,0), a null or self-intersecting polygon, etc.
!
!  Local parameters:
!
!    B =         Intersection point between the boundary and
!                the great circle defined by P and Q.
!
!    BP,BQ =     <B,P> and <B,Q>, respectively, maximized over
!                intersection points B that lie between P and
!                Q (on the shorter arc) -- used to find the
!                closest intersection points to P and Q
!    CN =        Q X P = normal to the plane of P and Q
!    D =         Dot product <B,P> or <B,Q>
!    EPS =       Parameter used to define Q as the point whose
!                orthogonal distance to (the midpoint of)
!                boundary edge V1->V2 is approximately EPS/
!                (2*Cos(A/2)), where <V1,V2> = Cos(A).
!    EVEN =      TRUE iff an even number of intersection points
!                lie between P and Q (on the shorter arc)
!    I1,I2 =     Indexes (LISTV elements) of a pair of adjacent
!                boundary vertices (endpoints of a boundary
!                edge)
!    IERR =      Error flag for calls to INTRSC (not tested)
!    IMX =       Local copy of LV and maximum value of I1 and I2
!    K =         DO-loop index and LISTV index
!    K0 =        LISTV index of the first endpoint of the
!                boundary edge used to compute Q
!    LFT1,LFT2 = Logical variables associated with I1 and I2 in
!                the boundary traversal:  TRUE iff the vertex
!                is strictly to the left of Q->P (<V,CN> > 0)
!    N =         Local copy of NV
!    NI =        Number of intersections (between the boundary
!                curve and the great circle P-Q) encountered
!    PINR =      TRUE iff P is to the left of the directed
!                boundary edge associated with the closest
!                intersection point to P that lies between P
!                and Q (a left-to-right intersection as
!                viewed from Q), or there is no intersection
!                between P and Q (on the shorter arc)
!    PN,QN =     P X CN and CN X Q, respectively:  used to
!                locate intersections B relative to arc Q->P
!    Q =         (V1 + V2 + EPS*VN/VNRM)/QNRM, where V1->V2 is
!                the boundary edge indexed by LISTV(K0) ->
!                LISTV(K0+1)
!    QINR =      TRUE iff Q is to the left of the directed
!                boundary edge associated with the closest
!                intersection point to Q that lies between P
!                and Q (a right-to-left intersection as
!                viewed from Q), or there is no intersection
!                between P and Q (on the shorter arc)
!    QNRM =      Euclidean norm of V1+V2+EPS*VN/VNRM used to
!                compute (normalize) Q
!    V1,V2 =     Vertices indexed by I1 and I2 in the boundary
!                traversal
!    VN =        V1 X V2, where V1->V2 is the boundary edge
!                indexed by LISTV(K0) -> LISTV(K0+1)
!    VNRM =      Euclidean norm of VN
!
  implicit none
 
  integer lv
  integer nv
 
  real ( kind_real ) b(3)
  real ( kind_real ) bp
  real ( kind_real ) bq
  real ( kind_real ) cn(3)
  real ( kind_real ) d
  real ( kind_real ), parameter :: eps = 0.001_kind_real
  logical even
  integer i1
  integer i2
  integer ier
  integer ierr
  integer imx
  logical inside
  integer k
  integer k0
  logical lft1
  logical lft2
  integer listv(nv)
  integer n
  integer ni
  real ( kind_real ) p(3)
  logical pinr
  real ( kind_real ) pn(3)
  real ( kind_real ) q(3)
  logical qinr
  real ( kind_real ) qn(3)
  real ( kind_real ) qnrm
  real ( kind_real ) v1(3)
  real ( kind_real ) v2(3)
  real ( kind_real ) vn(3)
  real ( kind_real ) vnrm
  real ( kind_real ) xv(lv)
  real ( kind_real ) yv(lv)
  real ( kind_real ) zv(lv)
!
!  Default value
!
  inside = .false.
!
!  Store local parameters.
!
  imx = lv
  n = nv
!
!  Test for error 1.
!
  if ( n < 3 .or. imx < n ) then
    ier = 1
    return
  end if
!
!  Initialize K0.
!
  k0 = 0
  i1 = listv(1)
 
  if ( i1 < 1 .or. imx < i1 ) then
    ier = 2
    return
  end if
!
!  Increment K0 and set Q to a point immediately to the left
!  of the midpoint of edge V1->V2 = LISTV(K0)->LISTV(K0+1):
!  Q = (V1 + V2 + EPS*VN/VNRM)/QNRM, where VN = V1 X V2.
!
1 continue
 
  k0 = k0 + 1
 
  if ( n < k0 ) then
    ier = 4
    return
  end if
 
  i1 = listv(k0)
 
  if ( k0 < n ) then
    i2 = listv(k0+1)
  else
    i2 = listv(1)
  end if
 
  if ( i2 < 1 .or. imx < i2 ) then
    ier = 2
    return
  end if
 
  vn(1) = yv(i1) * zv(i2) - zv(i1) * yv(i2)
  vn(2) = zv(i1) * xv(i2) - xv(i1) * zv(i2)
  vn(3) = xv(i1) * yv(i2) - yv(i1) * xv(i2)
  vnrm = sqrt( sum( vn(1:3)**2 ) )
 
  if ( eq(vnrm,0.0_kind_real) ) then
    go to 1
  end if
 
  q(1) = xv(i1) + xv(i2) + eps * vn(1) / vnrm
  q(2) = yv(i1) + yv(i2) + eps * vn(2) / vnrm
  q(3) = zv(i1) + zv(i2) + eps * vn(3) / vnrm
 
  qnrm = sqrt( sum( q(1:3)**2 ) )
 
  q(1) = q(1) / qnrm
  q(2) = q(2) / qnrm
  q(3) = q(3) / qnrm
!
!  Compute CN = Q X P, PN = P X CN, and QN = CN X Q.
!
  cn(1) = q(2) * p(3) - q(3) * p(2)
  cn(2) = q(3) * p(1) - q(1) * p(3)
  cn(3) = q(1) * p(2) - q(2) * p(1)
 
  if ( eq(cn(1),0.0_kind_real) .and. eq(cn(2),0.0_kind_real)  .and. eq(cn(2),0.0_kind_real) ) then
    go to 1
  end if
 
  pn(1) = p(2) * cn(3) - p(3) * cn(2)
  pn(2) = p(3) * cn(1) - p(1) * cn(3)
  pn(3) = p(1) * cn(2) - p(2) * cn(1)
  qn(1) = cn(2) * q(3) - cn(3) * q(2)
  qn(2) = cn(3) * q(1) - cn(1) * q(3)
  qn(3) = cn(1) * q(2) - cn(2) * q(1)
!
!  Initialize parameters for the boundary traversal.
!
  ni = 0
  even = .true.
  bp = -2.0_kind_real
  bq = -2.0_kind_real
  pinr = .true.
  qinr = .true.
  i2 = listv(n)
 
  if ( i2 < 1 .or. imx < i2 ) then
    ier = 2
    return
  end if
 
  lft2 = sup(cn(1) * xv(i2) + cn(2) * yv(i2) + cn(3) * zv(i2),0.0_kind_real)
!
!  Loop on boundary arcs I1->I2.
!
  do k = 1, n
 
    i1 = i2
    lft1 = lft2
    i2 = listv(k)
 
    if ( i2 < 1 .or. imx < i2 ) then
      ier = 2
      return
    end if
 
    lft2 = sup(cn(1) * xv(i2) + cn(2) * yv(i2) + cn(3) * zv(i2),0.0_kind_real)
 
    if ( lft1 .eqv. lft2 ) then
      cycle
    end if
!
!  I1 and I2 are on opposite sides of Q->P.  Compute the
!  point of intersection B.
!
    ni = ni + 1
    v1(1) = xv(i1)
    v1(2) = yv(i1)
    v1(3) = zv(i1)
    v2(1) = xv(i2)
    v2(2) = yv(i2)
    v2(3) = zv(i2)
 
    call intrsc ( v1, v2, cn, b, ierr )
!
!  B is between Q and P (on the shorter arc) iff
!  B Forward Q->P and B Forward P->Q       iff
!  <B,QN> > 0 and 0 < <B,PN>.
!
    if ( sup(dot_product( b(1:3), qn(1:3) ),0.0_kind_real) .and. &
         sup(dot_product( b(1:3), pn(1:3) ),0.0_kind_real) ) then
!
!  Update EVEN, BQ, QINR, BP, and PINR.
!
      even = .not. even
      d = dot_product( b(1:3), q(1:3) )
 
      if ( sup(d,bq) ) then
        bq = d
        qinr = lft2
      end if
 
      d = dot_product( b(1:3), p(1:3) )
 
      if ( sup(d,bq) ) then
        bp = d
        pinr = lft1
      end if
 
    end if
 
  end do
!
!  Test for consistency:  NI must be even and QINR must be TRUE.
!
  if ( ni /= 2 * ( ni / 2 ) .or. .not. qinr ) then
    go to 1
  end if
!
!  Test for error 3:  different values of PINR and EVEN.
!
  if ( pinr .neqv. even ) then
    ier = 3
    return
  end if
 
  ier = 0
  inside = even
 
  return
end function inside
subroutine intadd ( kk, i1, i2, i3, list, lptr, lend, lnew )
 
!*****************************************************************************80
!
! Subroutine: intadd
!> Add an interior node to a triangulation
!
!  Discussion:
!
!    This subroutine adds an interior node to a triangulation
!    of a set of points on the unit sphere.  The data structure
!    is updated with the insertion of node KK into the triangle
!    whose vertices are I1, I2, and I3.  No optimization of the
!    triangulation is performed.
!
!    This routine is identical to the similarly named routine in TRIPACK.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer KK, the index of the node to be inserted.
!    1 <= KK and KK must not be equal to I1, I2, or I3.
!
!    Input, integer I1, I2, I3, indexes of the
!    counterclockwise-ordered sequence of vertices of a triangle which contains
!    node KK.
!
!    Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), LNEW,
!    the data structure defining the triangulation, created by TRMESH.  Triangle
!    (I1,I2,I3) must be included in the triangulation.
!    On output, updated with the addition of node KK.  KK
!    will be connected to nodes I1, I2, and I3.
!
!  Local parameters:
!
!    K =        Local copy of KK
!    LP =       LIST pointer
!    N1,N2,N3 = Local copies of I1, I2, and I3
!
  implicit none
 
  integer i1
  integer i2
  integer i3
  integer k
  integer kk
  integer lend(*)
  integer list(*)
  integer lnew
  integer lp
  integer lptr(*)
  integer n1
  integer n2
  integer n3
 
  k = kk
!
!  Initialization.
!
  n1 = i1
  n2 = i2
  n3 = i3
!
!  Add K as a neighbor of I1, I2, and I3.
!
  call lstptr ( lend(n1), n2, list, lptr, lp )
  call insert ( k, lp, list, lptr, lnew )
 
  call lstptr ( lend(n2), n3, list, lptr, lp )
  call insert ( k, lp, list, lptr, lnew )
 
  call lstptr ( lend(n3), n1, list, lptr, lp )
  call insert ( k, lp, list, lptr, lnew )
!
!  Add I1, I2, and I3 as neighbors of K.
!
  list(lnew) = n1
  list(lnew+1) = n2
  list(lnew+2) = n3
  lptr(lnew) = lnew + 1
  lptr(lnew+1) = lnew + 2
  lptr(lnew+2) = lnew
  lend(k) = lnew + 2
  lnew = lnew + 3
 
  return
end subroutine intadd
subroutine intrsc ( p1, p2, cn, p, ier )
 
!*****************************************************************************80
!
! Subroutine: intrsc
!> Find the intersection of two great circles
!
!  Discussion:
!
!    Given a great circle C and points P1 and P2 defining an
!    arc A on the surface of the unit sphere, where A is the
!    shorter of the two portions of the great circle C12
!    associated with P1 and P2, this subroutine returns the point
!    of intersection P between C and C12 that is closer to A.
!    Thus, if P1 and P2 lie in opposite hemispheres defined by
!    C, P is the point of intersection of C with A.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, real ( kind_real ) P1(3), P2(3), the coordinates of unit vectors.
!
!    Input, real ( kind_real ) CN(3), the coordinates of a nonzero vector
!    which defines C as the intersection of the plane whose normal is CN
!    with the unit sphere.  Thus, if C is to be the great circle defined
!    by P and Q, CN should be P X Q.
!
!    Output, real ( kind_real ) P(3), point of intersection defined above
!    unless IER is not 0, in which case P is not altered.
!
!    Output, integer IER, error indicator.
!    0, if no errors were encountered.
!    1, if <CN,P1> = <CN,P2>.  This occurs iff P1 = P2 or CN = 0 or there are
!      two intersection points at the same distance from A.
!    2, if P2 = -P1 and the definition of A is therefore ambiguous.
!
!  Local parameters:
!
!    D1 =  <CN,P1>
!    D2 =  <CN,P2>
!    I =   DO-loop index
!    PP =  P1 + T*(P2-P1) = Parametric representation of the
!          line defined by P1 and P2
!    PPN = Norm of PP
!    T =   D1/(D1-D2) = Parameter value chosen so that PP lies
!          in the plane of C
!
  implicit none
 
  real ( kind_real ) cn(3)
  real ( kind_real ) d1
  real ( kind_real ) d2
  integer ier
  real ( kind_real ) p(3)
  real ( kind_real ) p1(3)
  real ( kind_real ) p2(3)
  real ( kind_real ) pp(3)
  real ( kind_real ) ppn
  real ( kind_real ) t
 
  d1 = dot_product( cn(1:3), p1(1:3) )
  d2 = dot_product( cn(1:3), p2(1:3) )
 
  if ( eq(d1,d2) ) then
    ier = 1
    return
  end if
!
!  Solve for T such that <PP,CN> = 0 and compute PP and PPN.
!
  t = d1 / ( d1 - d2 )
 
  pp(1:3) = p1(1:3) + t * ( p2(1:3) - p1(1:3) )
 
  ppn = dot_product( pp(1:3), pp(1:3) )
!
!  PPN = 0 iff PP = 0 iff P2 = -P1 (and T = .5).
!
  if ( eq(ppn,0.0_kind_real) ) then
    ier = 2
    return
  end if
 
  ppn = sqrt( ppn )
!
!  Compute P = PP/PPN.
!
  p(1:3) = pp(1:3) / ppn
 
  ier = 0
 
  return
end subroutine intrsc
subroutine jrand ( n, ix, iy, iz, output )
 
!*****************************************************************************80
!
! Subroutine: jrand
!> Return a random integer between 1 and N
!
!  Discussion:
!
!   This function returns a uniformly distributed pseudorandom integer
!   in the range 1 to N.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Brian Wichmann, David Hill,
!    An Efficient and Portable Pseudo-random Number Generator,
!    Applied Statistics,
!    Volume 31, Number 2, 1982, pages 188-190.
!
!  Parameters:
!
!    Input, integer N, the maximum value to be returned.
!
!    Input/output, integer IX, IY, IZ = seeds initialized to
!    values in the range 1 to 30,000 before the first call to JRAND, and
!    not altered between subsequent calls (unless a sequence of random
!    numbers is to be repeated by reinitializing the seeds).
!
!    Output, integer JRAND, a random integer in the range 1 to N.
!
!  Local parameters:
!
!    U = Pseudo-random number uniformly distributed in the interval (0,1).
!    X = Pseudo-random number in the range 0 to 3 whose fractional part is U.
!
  implicit none
 
  integer ix
  integer iy
  integer iz
  integer output
  integer n
  real ( kind_real ) u
  real ( kind_real ) x
 
  ix = mod( 171 * ix, 30269 )
  iy = mod( 172 * iy, 30307 )
  iz = mod( 170 * iz, 30323 )
 
  x = ( real( ix, kind_real ) / 30269.0_kind_real ) &
    + ( real( iy, kind_real ) / 30307.0_kind_real ) &
    + ( real( iz, kind_real ) / 30323.0_kind_real )
 
  u = x - int( x )
  output = int( real( n, kind_real ) * u ) + 1
 
  return
end subroutine jrand
subroutine left ( x1, y1, z1, x2, y2, z2, x0, y0, z0, output )
 
!*****************************************************************************80
!
! Subroutine: left
!> Determin whether a node is to the left of a plane through the origin
!
!  Discussion:
!
!    This function determines whether node N0 is in the
!    (closed) left hemisphere defined by the plane containing
!    N1, N2, and the origin, where left is defined relative to
!    an observer at N1 facing N2.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, real ( kind_real ) X1, Y1, Z1 = Coordinates of N1.
!
!    Input, real ( kind_real ) X2, Y2, Z2 = Coordinates of N2.
!
!    Input, real ( kind_real ) X0, Y0, Z0 = Coordinates of N0.
!
!    Output, logical LEFT = TRUE if and only if N0 is in the closed
!    left hemisphere.
!
  implicit none
 
  logical output
  real ( kind_real ) x0
  real ( kind_real ) x1
  real ( kind_real ) x2
  real ( kind_real ) y0
  real ( kind_real ) y1
  real ( kind_real ) y2
  real ( kind_real ) z0
  real ( kind_real ) z1
  real ( kind_real ) z2
  real ( kind_real ) zz
!
!  LEFT = TRUE iff <N0,N1 X N2> = det(N0,N1,N2) >= 0.
!
 
  call det ( x1, y1, z1, x2, y2, z2, x0, y0, z0, zz )
  output = sup(zz,0.0_kind_real)
 
  return
end subroutine left
subroutine lstptr ( lpl, nb, list, lptr, output )
 
!*****************************************************************************80
!
! Subroutine: lstptr
!> Return the index of NB in the adjacency list
!
!  Discussion:
!
!    This function returns the index (LIST pointer) of NB in
!    the adjacency list for N0, where LPL = LEND(N0).
!
!    This function is identical to the similarly named function in TRIPACK.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer LPL, is LEND(N0).
!
!    Input, integer NB, index of the node whose pointer is to
!    be returned.  NB must be connected to N0.
!
!    Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), the data
!    structure defining the triangulation, created by TRMESH.
!
!    Output, integer LSTPTR, pointer such that LIST(LSTPTR) = NB or
!    LIST(LSTPTR) = -NB, unless NB is not a neighbor of N0, in which
!    case LSTPTR = LPL.
!
!  Local parameters:
!
!    LP = LIST pointer
!    ND = Nodal index
!
  implicit none
 
  integer list(*)
  integer lp
  integer lpl
  integer lptr(*)
  integer output
  integer nb
  integer nd
 
  lp = lptr(lpl)
 
  do
 
    nd = list(lp)
 
    if ( nd == nb ) then
      exit
    end if
 
    lp = lptr(lp)
 
    if ( lp == lpl ) then
      exit
    end if
 
  end do
 
  output = lp
 
  return
end subroutine lstptr
subroutine swap ( in1, in2, io1, io2, list, lptr, lend, lp21 )
 
!*****************************************************************************80
!
! Subroutine: swap
!> Replace the diagonal arc of a quadrilateral with the other diagonal
!
!  Discussion:
!
!    Given a triangulation of a set of points on the unit
!    sphere, this subroutine replaces a diagonal arc in a
!    strictly convex quadrilateral (defined by a pair of adja-
!    cent triangles) with the other diagonal.  Equivalently, a
!    pair of adjacent triangles is replaced by another pair
!    having the same union.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer IN1, IN2, IO1, IO2, nodal indexes of the
!    vertices of the quadrilateral.  IO1-IO2 is replaced by IN1-IN2.
!    (IO1,IO2,IN1) and (IO2,IO1,IN2) must be triangles on input.
!
!    Input/output, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N),
!    the data structure defining the triangulation, created by TRMESH.
!    On output, updated with the swap; triangles (IO1,IO2,IN1) an (IO2,IO1,IN2)
!    are replaced by (IN1,IN2,IO2) and (IN2,IN1,IO1) unless LP21 = 0.
!
!    Output, integer LP21, index of IN1 as a neighbor of IN2 after
!    the swap is performed unless IN1 and IN2 are adjacent on input, in which
!    case LP21 = 0.
!
!  Local parameters:
!
!    LP, LPH, LPSAV = LIST pointers
!
  implicit none
 
  integer in1
  integer in2
  integer io1
  integer io2
  integer lend(*)
  integer list(*)
  integer lp
  integer lp21
  integer lph
  integer lpsav
  integer lptr(*)
!
!  Test for IN1 and IN2 adjacent.
!
  call lstptr ( lend(in1), in2, list, lptr, lp )
 
  if ( abs( list(lp) ) == in2 ) then
    lp21 = 0
    return
  end if
!
!  Delete IO2 as a neighbor of IO1.
!
  call lstptr ( lend(io1), in2, list, lptr, lp )
  lph = lptr(lp)
  lptr(lp) = lptr(lph)
!
!  If IO2 is the last neighbor of IO1, make IN2 the last neighbor.
!
  if ( lend(io1) == lph ) then
    lend(io1) = lp
  end if
!
!  Insert IN2 as a neighbor of IN1 following IO1 using the hole created above.
!
  call lstptr ( lend(in1), io1, list, lptr, lp )
  lpsav = lptr(lp)
  lptr(lp) = lph
  list(lph) = in2
  lptr(lph) = lpsav
!
!  Delete IO1 as a neighbor of IO2.
!
  call lstptr ( lend(io2), in1, list, lptr, lp )
  lph = lptr(lp)
  lptr(lp) = lptr(lph)
!
!  If IO1 is the last neighbor of IO2, make IN1 the last neighbor.
!
  if ( lend(io2) == lph ) then
    lend(io2) = lp
  end if
!
!  Insert IN1 as a neighbor of IN2 following IO2.
!
  call lstptr ( lend(in2), io2, list, lptr, lp )
  lpsav = lptr(lp)
  lptr(lp) = lph
  list(lph) = in1
  lptr(lph) = lpsav
  lp21 = lph
 
  return
end subroutine swap
subroutine swptst ( n1, n2, n3, n4, x, y, z, output )
 
!*****************************************************************************80
!
! Subroutine: swptst
!> Decide whether to replace a diagonal arc by the other
!
!  Discussion:
!
!    This function decides whether or not to replace a
!    diagonal arc in a quadrilateral with the other diagonal.
!    The decision will be to swap (SWPTST = TRUE) if and only
!    if N4 lies above the plane (in the half-space not contain-
!    ing the origin) defined by (N1,N2,N3), or equivalently, if
!    the projection of N4 onto this plane is interior to the
!    circumcircle of (N1,N2,N3).  The decision will be for no
!    swap if the quadrilateral is not strictly convex.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer N1, N2, N3, N4, indexes of the four nodes
!    defining the quadrilateral with N1 adjacent to N2, and (N1,N2,N3) in
!    counterclockwise order.  The arc connecting N1 to N2 should be replaced
!    by an arc connecting N3 to N4 if SWPTST = TRUE.  Refer to subroutine SWAP.
!
!    Input, real ( kind_real ) X(N), Y(N), Z(N), the coordinates of the nodes.
!
!    Output, logical SWPTST, TRUE if and only if the arc connecting N1
!    and N2 should be swapped for an arc connecting N3 and N4.
!
!  Local parameters:
!
!    DX1,DY1,DZ1 = Coordinates of N4->N1
!    DX2,DY2,DZ2 = Coordinates of N4->N2
!    DX3,DY3,DZ3 = Coordinates of N4->N3
!    X4,Y4,Z4 =    Coordinates of N4
!
  implicit none
 
  real ( kind_real ) dx1
  real ( kind_real ) dx2
  real ( kind_real ) dx3
  real ( kind_real ) dy1
  real ( kind_real ) dy2
  real ( kind_real ) dy3
  real ( kind_real ) dz1
  real ( kind_real ) dz2
  real ( kind_real ) dz3
  integer n1
  integer n2
  integer n3
  integer n4
  logical output
  real ( kind_real ) x(*)
  real ( kind_real ) x4
  real ( kind_real ) y(*)
  real ( kind_real ) y4
  real ( kind_real ) z(*)
  real ( kind_real ) z4
  real ( kind_real ) zz
 
  x4 = x(n4)
  y4 = y(n4)
  z4 = z(n4)
  dx1 = x(n1) - x4
  dx2 = x(n2) - x4
  dx3 = x(n3) - x4
  dy1 = y(n1) - y4
  dy2 = y(n2) - y4
  dy3 = y(n3) - y4
  dz1 = z(n1) - z4
  dz2 = z(n2) - z4
  dz3 = z(n3) - z4
!
!  N4 lies above the plane of (N1,N2,N3) iff N3 lies above
!  the plane of (N2,N1,N4) iff Det(N3-N4,N2-N4,N1-N4) =
!  (N3-N4,N2-N4 X N1-N4) > 0.
!
 
  call det ( dx2, dy2, dz2, dx1, dy1, dz1, dx3, dy3, dz3, zz )
  output = sup(zz,0.0_kind_real)
 
  return
end subroutine swptst
subroutine trfind ( nst, p, n, x, y, z, list, lptr, lend, b1, b2, b3, i1, &
  i2, i3 )
 
!*****************************************************************************80
!
! Subroutine: trfind
!> Locate a point relative to a triangulation
!
!  Discussion:
!
!    This subroutine locates a point P relative to a triangulation
!    created by TRMESH.  If P is contained in
!    a triangle, the three vertex indexes and barycentric
!    coordinates are returned.  Otherwise, the indexes of the
!    visible boundary nodes are returned.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer NST, index of a node at which TRFIND begins
!    its search.  Search time depends on the proximity of this node to P.
!
!    Input, real ( kind_real ) P(3), the x, y, and z coordinates (in that order)
!    of the point P to be located.
!
!    Input, integer N, the number of nodes in the triangulation.
!    3 <= N.
!
!    Input, real ( kind_real ) X(N), Y(N), Z(N), the coordinates of the
!    triangulation nodes (unit vectors).
!
!    Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), the
!    data structure defining the triangulation, created by TRMESH.
!
!    Output, real ( kind_real ) B1, B2, B3, the unnormalized barycentric
!    coordinates of the central projection of P onto the underlying planar
!    triangle if P is in the convex hull of the nodes.  These parameters
!    are not altered if I1 = 0.
!
!    Output, integer I1, I2, I3, the counterclockwise-ordered
!    vertex indexes of a triangle containing P if P is contained in a triangle.
!    If P is not in the convex hull of the nodes, I1 and I2 are the rightmost
!    and leftmost (boundary) nodes that are visible from P, and I3 = 0.  (If
!    all boundary nodes are visible from P, then I1 and I2 coincide.)
!    I1 = I2 = I3 = 0 if P and all of the nodes are coplanar (lie on a
!    common great circle.
!
!  Local parameters:
!
!    EPS =      Machine precision
!    IX,IY,IZ = Integer seeds for JRAND
!    LP =       LIST pointer
!    N0,N1,N2 = Nodes in counterclockwise order defining a
!               cone (with vertex N0) containing P, or end-
!               points of a boundary edge such that P Right
!               N1->N2
!    N1S,N2S =  Initially-determined values of N1 and N2
!    N3,N4 =    Nodes opposite N1->N2 and N2->N1, respectively
!    NEXT =     Candidate for I1 or I2 when P is exterior
!    NF,NL =    First and last neighbors of N0, or first
!               (rightmost) and last (leftmost) nodes
!               visible from P when P is exterior to the
!               triangulation
!    PTN1 =     Scalar product <P,N1>
!    PTN2 =     Scalar product <P,N2>
!    Q =        (N2 X N1) X N2  or  N1 X (N2 X N1) -- used in
!               the boundary traversal when P is exterior
!    S12 =      Scalar product <N1,N2>
!    TOL =      Tolerance (multiple of EPS) defining an upper
!               bound on the magnitude of a negative bary-
!               centric coordinate (B1 or B2) for P in a
!               triangle -- used to avoid an infinite number
!               of restarts with 0 <= B3 < EPS and B1 < 0 or
!               B2 < 0 but small in magnitude
!    XP,YP,ZP = Local variables containing P(1), P(2), and P(3)
!    X0,Y0,Z0 = Dummy arguments for DET
!    X1,Y1,Z1 = Dummy arguments for DET
!    X2,Y2,Z2 = Dummy arguments for DET
!
  implicit none
 
  integer n
 
  real ( kind_real ) b1
  real ( kind_real ) b2
  real ( kind_real ) b3
  real ( kind_real ) output
  real ( kind_real ) eps
  integer i1
  integer i2
  integer i3
  integer ( kind = 4 ), save :: ix = 1
  integer ( kind = 4 ), save :: iy = 2
  integer ( kind = 4 ), save :: iz = 3
  integer lend(n)
  integer list(6*(n-2))
  integer lp
  integer lptr(6*(n-2))
  integer n0
  integer n1
  integer n1s
  integer n2
  integer n2s
  integer n3
  integer n4
  integer next
  integer nf
  integer nl
  integer nst
  real ( kind_real ) p(3)
  real ( kind_real ) ptn1
  real ( kind_real ) ptn2
  real ( kind_real ) q(3)
  real ( kind_real ) s12
  real ( kind_real ) tol
  real ( kind_real ) x(n)
  real ( kind_real ) xp
  real ( kind_real ) y(n)
  real ( kind_real ) yp
  real ( kind_real ) z(n)
  real ( kind_real ) zp
!
!  Initialize variables.
!
  xp = p(1)
  yp = p(2)
  zp = p(3)
  n0 = nst
 
  if ( n0 < 1 .or. n < n0 ) then
    call jrand ( n, ix, iy, iz, n0 )
  end if
!
!  Compute the relative machine precision EPS and TOL.
!
  eps = epsilon( eps )
  tol = 100.0_kind_real * eps
!
!  Set NF and NL to the first and last neighbors of N0, and initialize N1 = NF.
!
2 continue
 
  lp = lend(n0)
  nl = list(lp)
  lp = lptr(lp)
  nf = list(lp)
  n1 = nf
!
!  Find a pair of adjacent neighbors N1,N2 of N0 that define
!  a wedge containing P:  P LEFT N0->N1 and P RIGHT N0->N2.
!
  if ( 0 < nl ) then
!
!  N0 is an interior node.  Find N1.
!
3   continue
 
    call det(x(n0),y(n0),z(n0),x(n1),y(n1),z(n1),xp,yp,zp,output)
    if ( inf(output,0.0_kind_real) ) then
      lp = lptr(lp)
      n1 = list(lp)
      if ( n1 == nl ) then
        go to 6
      end if
      go to 3
    end if
  else
!
!  N0 is a boundary node.  Test for P exterior.
!
    nl = -nl
!
!  Is P to the right of the boundary edge N0->NF?
!
    call det(x(n0),y(n0),z(n0),x(nf),y(nf),z(nf), xp,yp,zp,output)
    if ( inf(output,0.0_kind_real) ) then
      n1 = n0
      n2 = nf
      go to 9
    end if
!
!  Is P to the right of the boundary edge NL->N0?
!
    call det(x(nl),y(nl),z(nl),x(n0),y(n0),z(n0),xp,yp,zp,output)
    if ( inf(output,0.0_kind_real) ) then
      n1 = nl
      n2 = n0
      go to 9
    end if
 
  end if
!
!  P is to the left of arcs N0->N1 and NL->N0.  Set N2 to the
!  next neighbor of N0 (following N1).
!
4 continue
 
    lp = lptr(lp)
    n2 = abs( list(lp) )
 
    call det(x(n0),y(n0),z(n0),x(n2),y(n2),z(n2),xp,yp,zp,output)
    if ( inf(output,0.0_kind_real) ) then
      go to 7
    end if
 
    n1 = n2
 
    if ( n1 /= nl ) then
      go to 4
    end if
 
  call det ( x(n0), y(n0), z(n0), x(nf), y(nf), z(nf), xp, yp, zp, output )
  if ( inf(output,0.0_kind_real) ) then
    go to 6
  end if
!
!  P is left of or on arcs N0->NB for all neighbors NB
!  of N0.  Test for P = +/-N0.
!
  if ( inf(abs( x(n0 ) * xp + y(n0) * yp + z(n0) * zp),1.0-4.0*eps) ) then
!
!  All points are collinear iff P Left NB->N0 for all
!  neighbors NB of N0.  Search the neighbors of N0.
!  Note:  N1 = NL and LP points to NL.
!
    do
 
      call det(x(n1),y(n1),z(n1),x(n0),y(n0),z(n0),xp,yp,zp,output)
      if ( inf(output,0.0_kind_real) ) then
        exit
      end if
 
      lp = lptr(lp)
      n1 = abs( list(lp) )
 
      if ( n1 == nl ) then
        i1 = 0
        i2 = 0
        i3 = 0
        return
      end if
 
    end do
 
  end if
!
!  P is to the right of N1->N0, or P = +/-N0.  Set N0 to N1 and start over.
!
  n0 = n1
  go to 2
!
!  P is between arcs N0->N1 and N0->NF.
!
6 continue
 
  n2 = nf
!
!  P is contained in a wedge defined by geodesics N0-N1 and
!  N0-N2, where N1 is adjacent to N2.  Save N1 and N2 to
!  test for cycling.
!
7 continue
 
  n3 = n0
  n1s = n1
  n2s = n2
!
!  Top of edge-hopping loop:
!
8 continue
 
  call det ( x(n1),y(n1),z(n1),x(n2),y(n2),z(n2),xp,yp,zp,b3 )
 
  if ( inf(b3,0.0_kind_real) ) then
!
!  Set N4 to the first neighbor of N2 following N1 (the
!  node opposite N2->N1) unless N1->N2 is a boundary arc.
!
    call lstptr ( lend(n2), n1, list, lptr, lp )
 
    if ( list(lp) < 0 ) then
      go to 9
    end if
 
    lp = lptr(lp)
    n4 = abs( list(lp) )
!
!  Define a new arc N1->N2 which intersects the geodesic N0-P.
!
    call det ( x(n0),y(n0),z(n0),x(n4),y(n4),z(n4),xp,yp,zp,output )
    if ( inf(output,0.0_kind_real) ) then
      n3 = n2
      n2 = n4
      n1s = n1
      if ( n2 /= n2s .and. n2 /= n0 ) then
        go to 8
      end if
    else
      n3 = n1
      n1 = n4
      n2s = n2
      if ( n1 /= n1s .and. n1 /= n0 ) then
        go to 8
      end if
    end if
!
!  The starting node N0 or edge N1-N2 was encountered
!  again, implying a cycle (infinite loop).  Restart
!  with N0 randomly selected.
!
    call jrand ( n, ix, iy, iz, n0 )
    go to 2
 
  end if
!
!  P is in (N1,N2,N3) unless N0, N1, N2, and P are collinear
!  or P is close to -N0.
!
  if ( supeq(b3,eps) ) then
!
!  B3 /= 0.
!
    call det(x(n2),y(n2),z(n2),x(n3),y(n3),z(n3),xp,yp,zp,b1)
    call det(x(n3),y(n3),z(n3),x(n1),y(n1),z(n1),xp,yp,zp,b2)
!
!  Restart with N0 randomly selected.
!
    if ( inf(b1,-tol) .or. inf(b2,-tol) ) then
      call jrand ( n, ix, iy, iz, n0 )
      go to 2
    end if
 
  else
!
!  B3 = 0 and thus P lies on N1->N2. Compute
!  B1 = Det(P,N2 X N1,N2) and B2 = Det(P,N1,N2 X N1).
!
    b3 = 0.0_kind_real
    s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
    ptn1 = xp * x(n1) + yp * y(n1) + zp * z(n1)
    ptn2 = xp * x(n2) + yp * y(n2) + zp * z(n2)
    b1 = ptn1 - s12 * ptn2
    b2 = ptn2 - s12 * ptn1
!
!  Restart with N0 randomly selected.
!
    if ( inf(b1,-tol) .or. inf(b2,-tol) ) then
      call jrand ( n, ix, iy, iz, n0 )
      go to 2
    end if
 
  end if
!
!  P is in (N1,N2,N3).
!
  i1 = n1
  i2 = n2
  i3 = n3
  b1 = max( b1, 0.0_kind_real )
  b2 = max( b2, 0.0_kind_real )
  return
!
!  P Right N1->N2, where N1->N2 is a boundary edge.
!  Save N1 and N2, and set NL = 0 to indicate that
!  NL has not yet been found.
!
9 continue
 
  n1s = n1
  n2s = n2
  nl = 0
!
!  Counterclockwise Boundary Traversal:
!
10 continue
 
  lp = lend(n2)
  lp = lptr(lp)
  next = list(lp)
 
  call det(x(n2),y(n2),z(n2),x(next),y(next),z(next),xp,yp,zp,output)
  if ( supeq(output,0.0_kind_real) ) then
!
!  N2 is the rightmost visible node if P Forward N2->N1
!  or NEXT Forward N2->N1.  Set Q to (N2 X N1) X N2.
!
    s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
 
    q(1) = x(n1) - s12 * x(n2)
    q(2) = y(n1) - s12 * y(n2)
    q(3) = z(n1) - s12 * z(n2)
 
    if ( sup(xp * q(1) + yp * q(2) + zp * q(3),0.0_kind_real) ) then
      go to 11
    end if
 
    if ( sup(x(next) * q(1) + y(next) * q(2) + z(next) * q(3),0.0_kind_real) ) then
      go to 11
    end if
!
!  N1, N2, NEXT, and P are nearly collinear, and N2 is
!  the leftmost visible node.
!
    nl = n2
  end if
!
!  Bottom of counterclockwise loop:
!
  n1 = n2
  n2 = next
 
  if ( n2 /= n1s ) then
    go to 10
  end if
!
!  All boundary nodes are visible from P.
!
  i1 = n1s
  i2 = n1s
  i3 = 0
  return
!
!  N2 is the rightmost visible node.
!
11 continue
 
  nf = n2
 
  if ( nl == 0 ) then
!
!  Restore initial values of N1 and N2, and begin the search
!  for the leftmost visible node.
!
    n2 = n2s
    n1 = n1s
!
!  Clockwise Boundary Traversal:
!
12  continue
 
    lp = lend(n1)
    next = -list(lp)
 
    call det ( x(next), y(next), z(next), x(n1), y(n1), z(n1), xp, yp, zp, output )
    if ( supeq(output,0.0_kind_real)  ) then
!
!  N1 is the leftmost visible node if P or NEXT is
!  forward of N1->N2.  Compute Q = N1 X (N2 X N1).
!
      s12 = x(n1) * x(n2) + y(n1) * y(n2) + z(n1) * z(n2)
      q(1) = x(n2) - s12 * x(n1)
      q(2) = y(n2) - s12 * y(n1)
      q(3) = z(n2) - s12 * z(n1)
 
      if ( sup(xp * q(1) + yp * q(2) + zp * q(3),0.0_kind_real) ) then
        go to 13
      end if
 
      if ( sup(x(next) * q(1) + y(next) * q(2) + z(next) * q(3),0.0_kind_real) ) then
        go to 13
      end if
!
!  P, NEXT, N1, and N2 are nearly collinear and N1 is the rightmost
!  visible node.
!
      nf = n1
    end if
!
!  Bottom of clockwise loop:
!
    n2 = n1
    n1 = next
 
    if ( n1 /= n1s ) then
      go to 12
    end if
!
!  All boundary nodes are visible from P.
!
    i1 = n1
    i2 = n1
    i3 = 0
    return
!
!  N1 is the leftmost visible node.
!
13   continue
 
    nl = n1
 
  end if
!
!  NF and NL have been found.
!
  i1 = nf
  i2 = nl
  i3 = 0
 
  return
end subroutine trfind
subroutine trlist ( n, list, lptr, lend, nrow, nt, ltri, ier )
 
!*****************************************************************************80
!
! Subroutine: trlist
!> Convert a triangulation data structure to a triangle list
!
!  Discussion:
!
!    This subroutine converts a triangulation data structure
!    from the linked list created by TRMESH to a triangle list.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer N, the number of nodes in the triangulation.
!    3 <= N.
!
!    Input, integer LIST(6*(N-2)), LPTR(6*(N-2)), LEND(N), linked
!    list data structure defining the triangulation.  Refer to TRMESH.
!
!    Input, integer NROW, the number of rows (entries per triangle)
!    reserved for the triangle list LTRI.  The value must be 6 if only the
!    vertex indexes and neighboring triangle indexes are to be stored, or 9
!    if arc indexes are also to be assigned and stored.  Refer to LTRI.
!
!    Output, integer NT, the number of triangles in the
!    triangulation unless IER /=0, in which case NT = 0.  NT = 2N-NB-2 if
!    NB >= 3 or 2N-4 if NB = 0, where NB is the number of boundary nodes.
!
!    Output, integer LTRI(NROW,*).  The second dimension of LTRI
!    must be at least NT, where NT will be at most 2*N-4.  The J-th column
!    contains the vertex nodal indexes (first three rows), neighboring triangle
!    indexes (second three rows), and, if NROW = 9, arc indexes (last three
!    rows) associated with triangle J for J = 1,...,NT.  The vertices are
!    ordered counterclockwise with the first vertex taken to be the one
!    with smallest index.  Thus, LTRI(2,J) and LTRI(3,J) are larger than
!    LTRI(1,J) and index adjacent neighbors of node LTRI(1,J).  For
!    I = 1,2,3, LTRI(I+3,J) and LTRI(I+6,J) index the triangle and arc,
!    respectively, which are opposite (not shared by) node LTRI(I,J), with
!    LTRI(I+3,J) = 0 if LTRI(I+6,J) indexes a boundary arc.  Vertex indexes
!    range from 1 to N, triangle indexes from 0 to NT, and, if included,
!    arc indexes from 1 to NA, where NA = 3N-NB-3 if NB >= 3 or 3N-6 if
!    NB = 0.  The triangles are ordered on first (smallest) vertex indexes.
!
!    Output, integer IER, error indicator.
!    0, if no errors were encountered.
!    1, if N or NROW is outside its valid range on input.
!    2, if the triangulation data structure (LIST,LPTR,LEND) is invalid.
!      Note, however, that these arrays are not completely tested for validity.
!
!  Local parameters:
!
!    ARCS =     Logical variable with value TRUE iff are
!               indexes are to be stored
!    I,J =      LTRI row indexes (1 to 3) associated with
!               triangles KT and KN, respectively
!    I1,I2,I3 = Nodal indexes of triangle KN
!    ISV =      Variable used to permute indexes I1,I2,I3
!    KA =       Arc index and number of currently stored arcs
!    KN =       Index of the triangle that shares arc I1-I2 with KT
!    KT =       Triangle index and number of currently stored triangles
!    LP =       LIST pointer
!    LP2 =      Pointer to N2 as a neighbor of N1
!    LPL =      Pointer to the last neighbor of I1
!    LPLN1 =    Pointer to the last neighbor of N1
!    N1,N2,N3 = Nodal indexes of triangle KT
!    NM2 =      N-2
!
  implicit none
 
  integer n
  integer nrow
 
  logical arcs
  integer i
  integer i1
  integer i2
  integer i3
  integer ier
  integer isv
  integer j
  integer ka
  integer kn
  integer kt
  integer lend(n)
  integer list(6*(n-2))
  integer lp
  integer lp2
  integer lpl
  integer lpln1
  integer lptr(6*(n-2))
  integer ltri(nrow,*)
  integer n1
  integer n2
  integer n3
  integer nm2
  integer nt
!
!  Test for invalid input parameters.
!
  if ( n < 3 .or. ( nrow /= 6 .and. nrow /= 9 ) ) then
    nt = 0
    ier = 1
    return
  end if
!
!  Initialize parameters for loop on triangles KT = (N1,N2,
!  N3), where N1 < N2 and N1 < N3.
!
!  ARCS = TRUE iff arc indexes are to be stored.
!  KA,KT = Numbers of currently stored arcs and triangles.
!  NM2 = Upper bound on candidates for N1.
!
  arcs = nrow == 9
  ka = 0
  kt = 0
  nm2 = n-2
!
!  Loop on nodes N1.
!
  do n1 = 1, nm2
!
!  Loop on pairs of adjacent neighbors (N2,N3).  LPLN1 points
!  to the last neighbor of N1, and LP2 points to N2.
!
    lpln1 = lend(n1)
    lp2 = lpln1
 
1   continue
 
      lp2 = lptr(lp2)
      n2 = list(lp2)
      lp = lptr(lp2)
      n3 = abs( list(lp) )
 
      if ( n2 < n1 .or. n3 < n1 ) then
        go to 8
      end if
!
!  Add a new triangle KT = (N1,N2,N3).
!
      kt = kt + 1
      ltri(1,kt) = n1
      ltri(2,kt) = n2
      ltri(3,kt) = n3
!
!  Loop on triangle sides (I2,I1) with neighboring triangles
!  KN = (I1,I2,I3).
!
      do i = 1, 3
 
        if ( i == 1 ) then
          i1 = n3
          i2 = n2
        else if ( i == 2 ) then
          i1 = n1
          i2 = n3
        else
          i1 = n2
          i2 = n1
        end if
        j = 0
!
!  Set I3 to the neighbor of I1 that follows I2 unless
!  I2->I1 is a boundary arc.
!
        lpl = lend(i1)
        lp = lptr(lpl)
 
        do
 
          if ( list(lp) == i2 ) then
            go to 3
          end if
 
          lp = lptr(lp)
 
          if ( lp == lpl ) then
            exit
          end if
 
        end do
!
!  Invalid triangulation data structure:  I1 is a neighbor of
!  I2, but I2 is not a neighbor of I1.
!
        if ( abs( list(lp) ) /= i2 ) then
          nt = 0
          ier = 2
          return
        end if
!
!  I2 is the last neighbor of I1.  Bypass the search for a neighboring
!  triangle if I2->I1 is a boundary arc.
!
        kn = 0
 
        if ( list(lp) < 0 ) then
          go to 6
        end if
!
!  I2->I1 is not a boundary arc, and LP points to I2 as
!  a neighbor of I1.
!
3       continue
 
        lp = lptr(lp)
        i3 = abs( list(lp) )
!
!  Find J such that LTRI(J,KN) = I3 (not used if KT < KN),
!  and permute the vertex indexes of KN so that I1 is smallest.
!
        if ( i1 < i2 .and. i1 < i3 ) then
          j = 3
        else if ( i2 < i3 ) then
          j = 2
          isv = i1
          i1 = i2
          i2 = i3
          i3 = isv
        else
          j = 1
          isv = i1
          i1 = i3
          i3 = i2
          i2 = isv
        end if
!
!  Test for KT < KN (triangle index not yet assigned).
!
        if ( n1 < i1 ) then
          cycle
        end if
!
!  Find KN, if it exists, by searching the triangle list in
!  reverse order.
!
        do kn = kt-1, 1, -1
          if ( ltri(1,kn) == i1 .and. &
               ltri(2,kn) == i2 .and. &
               ltri(3,kn) == i3 ) then
            go to 5
          end if
        end do
 
        cycle
!
!  Store KT as a neighbor of KN.
!
5       continue
 
        ltri(j+3,kn) = kt
!
!  Store KN as a neighbor of KT, and add a new arc KA.
!
6       continue
 
        ltri(i+3,kt) = kn
 
        if ( arcs ) then
          ka = ka + 1
          ltri(i+6,kt) = ka
          if ( kn /= 0 ) then
            ltri(j+6,kn) = ka
          end if
        end if
 
        end do
!
!  Bottom of loop on triangles.
!
8     continue
 
      if ( lp2 /= lpln1 ) then
        go to 1
      end if
 
  end do
 
  nt = kt
  ier = 0
 
  return
end subroutine trlist
subroutine trmesh ( mpl, n, x, y, z, list, lptr, lend, lnew, near, next, dist, ier )
 
!*****************************************************************************80
!
! Subroutine: trmesh
!> Create a Delaunay triangulation on the unit sphere
!
!  Discussion:
!
!    This subroutine creates a Delaunay triangulation of a
!    set of N arbitrarily distributed points, referred to as
!    nodes, on the surface of the unit sphere.  The Delaunay
!    triangulation is defined as a set of (spherical) triangles
!    with the following five properties:
!
!     1)  The triangle vertices are nodes.
!     2)  No triangle contains a node other than its vertices.
!     3)  The interiors of the triangles are pairwise disjoint.
!     4)  The union of triangles is the convex hull of the set
!           of nodes (the smallest convex set that contains
!           the nodes).  If the nodes are not contained in a
!           single hemisphere, their convex hull is the
!           entire sphere and there are no boundary nodes.
!           Otherwise, there are at least three boundary nodes.
!     5)  The interior of the circumcircle of each triangle
!           contains no node.
!
!    The first four properties define a triangulation, and the
!    last property results in a triangulation which is as close
!    as possible to equiangular in a certain sense and which is
!    uniquely defined unless four or more nodes lie in a common
!    plane.  This property makes the triangulation well-suited
!    for solving closest-point problems and for triangle-based
!    interpolation.
!
!    Provided the nodes are randomly ordered, the algorithm
!    has expected time complexity O(N*log(N)) for most nodal
!    distributions.  Note, however, that the complexity may be
!    as high as O(N**2) if, for example, the nodes are ordered
!    on increasing latitude.
!
!    Spherical coordinates (latitude and longitude) may be
!    converted to Cartesian coordinates by Subroutine TRANS.
!
!    The following is a list of the software package modules
!    which a user may wish to call directly:
!
!    ADDNOD - Updates the triangulation by appending a new node.
!
!    BNODES - Returns an array containing the indexes of the
!             boundary nodes (if any) in counterclockwise
!             order.  Counts of boundary nodes, triangles,
!             and arcs are also returned.
!
!    EDGE   - Forces an arbitrary pair of nodes to be connected
!             by an arc in the triangulation.
!
!    INSIDE - Locates a point relative to a polygon on the
!             surface of the sphere.
!
!    INTRSC - Returns the point of intersection between a
!             pair of great circle arcs.
!
!    JRAND  - Generates a uniformly distributed pseudo-random integer.
!
!    LEFT   - Locates a point relative to a great circle.
!
!    TRANS  - Transforms spherical coordinates into Cartesian
!             coordinates on the unit sphere for input to
!             Subroutine TRMESH.
!
!    TRLIST - Converts the triangulation data structure to a
!             triangle list more suitable for use in a finite
!             element code.
!
!    TRMESH - Creates a Delaunay triangulation of a set of
!             nodes.
!
!  Modified:
!
!    16 June 2007
!
!  Author:
!
!    Robert Renka
!
!  Reference:
!
!    Robert Renka,
!    Algorithm 772: STRIPACK,
!    Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere,
!    ACM Transactions on Mathematical Software,
!    Volume 23, Number 3, September 1997, pages 416-434.
!
!  Parameters:
!
!    Input, integer N, the number of nodes in the triangulation.
!    3 <= N.
!
!    Input, real ( kind_real ) X(N), Y(N), Z(N), the coordinates of distinct
!    nodes.  (X(K),Y(K), Z(K)) is referred to as node K, and K is referred
!    to as a nodal index.  It is required that X(K)**2 + Y(K)**2 + Z(K)**2 = 1
!    for all K.  The first three nodes must not be collinear (lie on a
!    common great circle).
!
!    Output, integer LIST(6*(N-2)), nodal indexes which, along
!    with LPTR, LEND, and LNEW, define the triangulation as a set of N
!    adjacency lists; counterclockwise-ordered sequences of neighboring nodes
!    such that the first and last neighbors of a boundary node are boundary
!    nodes (the first neighbor of an interior node is arbitrary).  In order to
!    distinguish between interior and boundary nodes, the last neighbor of
!    each boundary node is represented by the negative of its index.
!
!    Output, integer LPTR(6*(N-2)), = Set of pointers (LIST
!    indexes) in one-to-one correspondence with the elements of LIST.
!    LIST(LPTR(I)) indexes the node which follows LIST(I) in cyclical
!    counterclockwise order (the first neighbor follows the last neighbor).
!
!    Output, integer LEND(N), pointers to adjacency lists.
!    LEND(K) points to the last neighbor of node K.  LIST(LEND(K)) < 0 if and
!    only if K is a boundary node.
!
!    Output, integer LNEW, pointer to the first empty location
!    in LIST and LPTR (list length plus one).  LIST, LPTR, LEND, and LNEW are
!    not altered if IER < 0, and are incomplete if 0 < IER.
!
!    Workspace, integer NEAR(N),
!    used to efficiently determine the nearest triangulation node to each
!    unprocessed node for use by ADDNOD.
!
!    Workspace, integer NEXT(N),
!    used to efficiently determine the nearest triangulation node to each
!    unprocessed node for use by ADDNOD.
!
!    Workspace, real ( kind_real ) DIST(N),
!    used to efficiently determine the nearest triangulation node to each
!    unprocessed node for use by ADDNOD.
!
!    Output, integer IER, error indicator:
!     0, if no errors were encountered.
!    -1, if N < 3 on input.
!    -2, if the first three nodes are collinear.
!     L, if nodes L and M coincide for some L < M.  The data structure
!      represents a triangulation of nodes 1 to M-1 in this case.
!
!  Local parameters:
!
!    D =        (Negative cosine of) distance from node K to node I
!    D1,D2,D3 = Distances from node K to nodes 1, 2, and 3, respectively
!    I,J =      Nodal indexes
!    I0 =       Index of the node preceding I in a sequence of
!               unprocessed nodes:  I = NEXT(I0)
!    K =        Index of node to be added and DO-loop index: 3 < K
!    LP =       LIST index (pointer) of a neighbor of K
!    LPL =      Pointer to the last neighbor of K
!    NEXTI =    NEXT(I)
!    NN =       Local copy of N
!
  implicit none
 
  type(mpl_type),intent(inout) :: mpl !< MPI data
 
  integer n
 
  real ( kind_real ) d
  real ( kind_real ) d1
  real ( kind_real ) d2
  real ( kind_real ) d3
  real ( kind_real ) dist(n)
  integer i
  integer i0
  integer ier
  integer j
  integer k
  logical output_1, output_2
  integer lend(n)
  integer list(6*(n-2))
  integer lnew
  integer lp
  integer lpl
  integer lptr(6*(n-2))
  integer near(n)
  integer next(n)
  integer nexti
  integer nn
  real ( kind_real ) x(n)
  real ( kind_real ) y(n)
  real ( kind_real ) z(n)
  character(len=1024),parameter :: subr = 'trmesh'
 
  nn = n
 
  if ( nn < 3 ) then
    ier = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TRMESH - Fatal error!'
    write ( *, '(a)' ) '  N < 3.'
    call mpl%abort(subr,'stop in stripack')
  end if
!
!  Initialize
!
  list = 0
  lptr = 0
  lend = 0
!
!  Store the first triangle in the linked list.
!
  call left (x(1),y(1),z(1),x(2),y(2),z(2),x(3),y(3),z(3),output_1)
  call left (x(2),y(2),z(2),x(1),y(1),z(1),x(3),y(3),z(3),output_2)
 
  if ( .not. output_1 ) then
!
!  The first triangle is (3,2,1) = (2,1,3) = (1,3,2).
!
    list(1) = 3
    lptr(1) = 2
    list(2) = -2
    lptr(2) = 1
    lend(1) = 2
 
    list(3) = 1
    lptr(3) = 4
    list(4) = -3
    lptr(4) = 3
    lend(2) = 4
 
    list(5) = 2
    lptr(5) = 6
    list(6) = -1
    lptr(6) = 5
    lend(3) = 6
 
  else if ( .not. output_2 ) then
!
!  The first triangle is (1,2,3):  3 Strictly Left 1->2,
!  i.e., node 3 lies in the left hemisphere defined by arc 1->2.
!
    list(1) = 2
    lptr(1) = 2
    list(2) = -3
    lptr(2) = 1
    lend(1) = 2
 
    list(3) = 3
    lptr(3) = 4
    list(4) = -1
    lptr(4) = 3
    lend(2) = 4
 
    list(5) = 1
    lptr(5) = 6
    list(6) = -2
    lptr(6) = 5
    lend(3) = 6
!
!  The first three nodes are collinear.
!
  else
 
    ier = -2
    return
 
  end if
!
!  Initialize LNEW and test for N = 3.
!
  lnew = 7
 
  if ( nn == 3 ) then
    ier = 0
    return
  end if
!
!  A nearest-node data structure (NEAR, NEXT, and DIST) is
!  used to obtain an expected-time (N*log(N)) incremental
!  algorithm by enabling constant search time for locating
!  each new node in the triangulation.
!
!  For each unprocessed node K, NEAR(K) is the index of the
!  triangulation node closest to K (used as the starting
!  point for the search in Subroutine TRFIND) and DIST(K)
!  is an increasing function of the arc length (angular
!  distance) between nodes K and NEAR(K):  -Cos(a) for arc
!  length a.
!
!  Since it is necessary to efficiently find the subset of
!  unprocessed nodes associated with each triangulation
!  node J (those that have J as their NEAR entries), the
!  subsets are stored in NEAR and NEXT as follows:  for
!  each node J in the triangulation, I = NEAR(J) is the
!  first unprocessed node in J's set (with I = 0 if the
!  set is empty), L = NEXT(I) (if 0 < I) is the second,
!  NEXT(L) (if 0 < L) is the third, etc.  The nodes in each
!  set are initially ordered by increasing indexes (which
!  maximizes efficiency) but that ordering is not main-
!  tained as the data structure is updated.
!
!  Initialize the data structure for the single triangle.
!
  near(1) = 0
  near(2) = 0
  near(3) = 0
 
  do k = nn, 4, -1
 
    d1 = -( x(k) * x(1) + y(k) * y(1) + z(k) * z(1) )
    d2 = -( x(k) * x(2) + y(k) * y(2) + z(k) * z(2) )
    d3 = -( x(k) * x(3) + y(k) * y(3) + z(k) * z(3) )
 
    if ( inf(d1,d2) .and. inf(d1,d3) ) then
      near(k) = 1
      dist(k) = d1
      next(k) = near(1)
      near(1) = k
    else if ( inf(d2,d1) .and. inf(d2,d3) ) then
      near(k) = 2
      dist(k) = d2
      next(k) = near(2)
      near(2) = k
    else
      near(k) = 3
      dist(k) = d3
      next(k) = near(3)
      near(3) = k
    end if
  end do
!
!  Add the remaining nodes.
!
  do k = 4, nn
 
    call addnod ( mpl, near(k), k, x, y, z, list, lptr, lend, lnew, ier )
 
    if ( ier /= 0 ) then
       write ( *, '(a)' ) 'TRMESH - ADDNOD - Fatal error!'
       write ( *, '(a,i8,a,i8)' ) '  Node ', k, ' is equal to already used node ', ier
       call mpl%abort(subr,'stop in stripack')
    end if
!
!  Remove K from the set of unprocessed nodes associated with NEAR(K).
!
    i = near(k)
    i0 = i
 
    if ( near(i) == k ) then
 
      near(i) = next(k)
 
    else
 
      i = near(i)
 
      do
 
        i0 = i
        i = next(i0)
 
        if ( i == k ) then
          exit
        end if
 
      end do
 
      next(i0) = next(k)
 
    end if
 
    near(k) = 0
!
!  Loop on neighbors J of node K.
!
    lpl = lend(k)
    lp = lpl
 
3     continue
 
    lp = lptr(lp)
    j = abs( list(lp) )
!
!  Loop on elements I in the sequence of unprocessed nodes
!  associated with J:  K is a candidate for replacing J
!  as the nearest triangulation node to I.  The next value
!  of I in the sequence, NEXT(I), must be saved before I
!  is moved because it is altered by adding I to K's set.
!
    i = near(j)
 
    do
 
      if ( i == 0 ) then
        exit
      end if
 
      nexti = next(i)
!
!  Test for the distance from I to K less than the distance
!  from I to J. Indistinguishability threshold for cross-plateform reproducibility
!
      d = - ( x(i) * x(k) + y(i) * y(k) + z(i) * z(k) )
      if ( inf(d,dist(i)) ) then
!
!  Replace J by K as the nearest triangulation node to I:
!  update NEAR(I) and DIST(I), and remove I from J's set
!  of unprocessed nodes and add it to K's set.
!
        near(i) = k
        dist(i) = d
 
        if ( i == near(j) ) then
          near(j) = nexti
        else
          next(i0) = nexti
        end if
 
        next(i) = near(k)
        near(k) = i
      else
        i0 = i
      end if
 
      i = nexti
 
    end do
!
!  Bottom of loop on neighbors J.
!
    if ( lp /= lpl ) then
      go to 3
    end if
  end do
 
  return
end subroutine trmesh
 
end module tools_stripack