do ip1 = 1, ndp
    if ( ip1 /= ipmn1 .and. ip1 /= ipmn2 ) then
      jp1 = jp1+1
      iwp(jp1) = ip1
      wk(jp1) = dsqf(xdmp,ydmp,xd(ip1),yd(ip1))
    end if
  end do
 
  do jp1 = 3, ndp-1
 
    dsqmn = wk(jp1)
    jpmn = jp1
 
    do jp2 = jp1, ndp
      if(wk(jp2)<dsqmn) then
        dsqmn = wk(jp2)
        jpmn = jp2
      end if
    end do
 
    its = iwp(jp1)
    iwp(jp1) = iwp(jpmn)
    iwp(jpmn) = its
    wk(jpmn) = wk(jp1)
 
  end do
!
!  If necessary, modify the ordering in such a way that the
!  first three data points are not collinear.
!
  x1 = xd(ipmn1)
  y1 = yd(ipmn1)
  x2 = xd(ipmn2)
  y2 = yd(ipmn2)
 
  do jp = 3, ndp
    ip = iwp(jp)
    sp = spdt(xd(ip),yd(ip),x1,y1,x2,y2)
    vp = vpdt(xd(ip),yd(ip),x1,y1,x2,y2)
    if ( abs(vp) > ( abs(sp) * epsln ) )   go to 37
  end do
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'IDTANG - Fatal error!'
  write ( *, '(a)' ) '  All collinear data points.'
  stop
 
   37 continue
 
  if ( jp /= 3 ) then
 
    jpmx = jp
 
    do jpc = 4, jpmx
      jp = jpmx+4-jpc
      iwp(jp) = iwp(jp-1)
    end do
 
    iwp(3) = ip
 
  end if
!
!  Form the first triangle.  
!
!  Store point numbers of the vertexes of the triangle in the IPT array, 
!  store point numbers of the border line segments and the triangle number in
!  the IPL array.
!
  ip1 = ipmn1
  ip2 = ipmn2
  ip3 = iwp(3)
 
  if ( vpdt(xd(ip1),yd(ip1),xd(ip2),yd(ip2),xd(ip3),yd(ip3)) < 0.0E+00 ) then
    ip1 = ipmn2
    ip2 = ipmn1
  end if
 
  nt0 = 1
  ntt3 = 3
  ipt(1) = ip1
  ipt(2) = ip2
  ipt(3) = ip3
  nl0 = 3
  nlt3 = 9
  ipl(1) = ip1
  ipl(2) = ip2
  ipl(3) = 1
  ipl(4) = ip2
  ipl(5) = ip3
  ipl(6) = 1
  ipl(7) = ip3
  ipl(8) = ip1
  ipl(9) = 1
!
!  Add the remaining data points, one by one.
!
  do jp1 = 4, ndp

    ip1 = iwp(jp1)
    x1 = xd(ip1)
    y1 = yd(ip1)
!
!  Determine the first invisible and visible border line segments, iliv and
!  ilvs.
!
    do il = 1, nl0

      ip2 = ipl(3*il-2)
      ip3 = ipl(3*il-1)
      x2 = xd(ip2)
      y2 = yd(ip2)
      x3 = xd(ip3)
      y3 = yd(ip3)
      sp = spdt(x1,y1,x2,y2,x3,y3)
      vp = vpdt(x1,y1,x2,y2,x3,y3)

      if ( il == 1 ) then
        ixvs = 0
        if(vp<=(abs(sp)*(-epsln)))   ixvs = 1
        iliv = 1
        ilvs = 1
        go to 53
      end if

      ixvspv = ixvs

      if ( vp <= (abs(sp)*(-epsln)) ) then
        ixvs = 1
        if(ixvspv==1)      go to 53
        ilvs = il
        if(iliv/=1)        go to 54
        go to 53
      end if

      ixvs = 0

      if ( ixvspv /= 0 ) then
        iliv = il
        if(ilvs/=1)        go to 54
      end if

53     continue

    end do

    if(iliv==1.and.ilvs==1)  ilvs = nl0

54   continue

    if(ilvs<iliv)  ilvs = ilvs+nl0
!
!  Shift (rotate) the IPL array to have the invisible border
!  line segments contained in the first part of the array.
!
55   continue
 
    if ( iliv /= 1 ) then
 
      nlsh = iliv-1
      nlsht3 = nlsh*3
 
      do jl1 = 1,nlsht3
        jl2 = jl1+nlt3
        ipl(jl2) = ipl(jl1)
      end do
 
      do jl1 = 1,nlt3
        jl2 = jl1+nlsht3
        ipl(jl1) = ipl(jl2)
      end do
 
      ilvs = ilvs-nlsh
 
    end if
!
!  Add triangles to the IPT array, 
!  update border line segments in the IPL array, 
!  set flags for the border line segments to be reexamined in the IWL array.
!
    jwl = 0
 
    do il = ilvs, nl0
 
      ilt3 = il*3
      ipl1 = ipl(ilt3-2)
      ipl2 = ipl(ilt3-1)
      it   = ipl(ilt3)
!
!  Add a triangle to the IPT array.
!
      nt0 = nt0+1
      ntt3 = ntt3+3
      ipt(ntt3-2) = ipl2
      ipt(ntt3-1) = ipl1
      ipt(ntt3)   = ip1
!
!  Update border line segments in the IPL array.
!
      if ( il == ilvs ) then
        ipl(ilt3-1) = ip1
        ipl(ilt3)   = nt0
      end if
 
      if ( il == nl0 ) then
        nln = ilvs+1
        nlnt3 = nln*3
        ipl(nlnt3-2) = ip1
        ipl(nlnt3-1) = ipl(1)
        ipl(nlnt3)   = nt0
      end if
!
!  Determine the vertex that does not lie on the border
!  line segments.
!
      itt3 = it*3
      ipti = ipt(itt3-2)
 
      if ( ipti == ipl1 .or. ipti == ipl2 ) then
        ipti = ipt(itt3-1)
        if ( ipti == ipl1 .or. ipti == ipl2 ) then
          ipti = ipt(itt3)
        end if
      end if
!
!  Check if the exchange is necessary.
!
      if ( idxchg(xd,yd,ip1,ipti,ipl1,ipl2) /= 0 ) then
!
!  Modify the IPT array.
!
        ipt(itt3-2) = ipti
        ipt(itt3-1) = ipl1
        ipt(itt3)   = ip1
        ipt(ntt3-1) = ipti
        if(il==ilvs)  ipl(ilt3) = it
        if(il==nl0.and.ipl(3)==it)      ipl(3) = nt0
!
!  Set flags in the IWL array.
!
        jwl = jwl+4
        iwl(jwl-3) = ipl1
        iwl(jwl-2) = ipti
        iwl(jwl-1) = ipti
        iwl(jwl)   = ipl2

      end if
 
    end do
 
    nl0 = nln
    nlt3 = nlnt3
    nlf = jwl/2
 
    if ( nlf == 0 ) then
      go to 79
    end if
!
!  Improve triangulation.
!
    ntt3p3 = ntt3+3

    do irep = 1, nrep

      do ilf = 1,nlf

        ipl1 = iwl(2*ilf-1)
        ipl2 = iwl(2*ilf)
!
!  Locate in the ipt array two triangles on both sides of
!  the flagged line segment.
!
        ntf = 0

        do itt3r = 3,ntt3,3
          itt3 = ntt3p3-itt3r
          ipt1 = ipt(itt3-2)
          ipt2 = ipt(itt3-1)
          ipt3 = ipt(itt3)
          if(ipl1/=ipt1.and.ipl1/=ipt2.and. ipl1/=ipt3)      go to 71
          if(ipl2/=ipt1.and.ipl2/=ipt2.and. ipl2/=ipt3)      go to 71
          ntf = ntf+1
          itf(ntf) = itt3/3
          if(ntf==2)     go to 72
71         continue
        end do

        if ( ntf < 2 )       go to 76
!
!  Determine the vertexes of the triangles that do not lie
!  on the line segment.
!
72       continue

        it1t3 = itf(1)*3
        ipti1 = ipt(it1t3-2)
        if(ipti1/=ipl1.and.ipti1/=ipl2)    go to 73
        ipti1 = ipt(it1t3-1)

        if ( ipti1 == ipl1 .or. ipti1 == ipl2 ) then
          ipti1 = ipt(it1t3)
        end if

73       continue

        it2t3 = itf(2)*3
        ipti2 = ipt(it2t3-2)
        if(ipti2/=ipl1.and.ipti2/=ipl2)    go to 74
        ipti2 = ipt(it2t3-1)
        if(ipti2/=ipl1.and.ipti2/=ipl2)    go to 74
        ipti2 = ipt(it2t3)
!
!  Check if the exchange is necessary.
!
74       continue

        if(idxchg(xd,yd,ipti1,ipti2,ipl1,ipl2)==0) then
          go to 76
        end if
!
!  Modify the IPT array.
!
        ipt(it1t3-2) = ipti1
        ipt(it1t3-1) = ipti2
        ipt(it1t3)   = ipl1
        ipt(it2t3-2) = ipti2
        ipt(it2t3-1) = ipti1
        ipt(it2t3)   = ipl2
!
!  Set new flags.
!
        jwl = jwl+8
        iwl(jwl-7) = ipl1
        iwl(jwl-6) = ipti1
        iwl(jwl-5) = ipti1
        iwl(jwl-4) = ipl2
        iwl(jwl-3) = ipl2
        iwl(jwl-2) = ipti2
        iwl(jwl-1) = ipti2
        iwl(jwl)   = ipl1
        do jlt3 = 3,nlt3,3
          iplj1 = ipl(jlt3-2)
          iplj2 = ipl(jlt3-1)

          if((iplj1==ipl1.and.iplj2==ipti2).or. &
             (iplj2==ipl1.and.iplj1==ipti2)) then
                               ipl(jlt3) = itf(1)
          end if

          if((iplj1==ipl2.and.iplj2==ipti1).or. &
             (iplj2==ipl2.and.iplj1==ipti1)) then
                              ipl(jlt3) = itf(2)
          end if

        end do

76       continue

      end do
 
      nlfc = nlf
      nlf = jwl/2
!
!  Reset the IWL array for the next round.
!
      if ( nlf == nlfc ) go to 79

      jwl1mn = 2*nlfc+1
      nlft2 = nlf*2
 
      do jwl1 = jwl1mn,nlft2
        jwl = jwl1+1-jwl1mn
        iwl(jwl) = iwl(jwl1)
      end do
 
      nlf = jwl/2

    end do

79  continue

  end do
!
!  Rearrange the IPT array so that the vertexes of each triangle
!  are listed counter-clockwise.
!
  do itt3 = 3, ntt3, 3
 
    ip1 = ipt(itt3-2)
    ip2 = ipt(itt3-1)
    ip3 = ipt(itt3)
 
    if(vpdt(xd(ip1),yd(ip1),xd(ip2),yd(ip2),xd(ip3),yd(ip3)) < 0.0E+00 ) then
      ipt(itt3-2) = ip2
      ipt(itt3-1) = ip1
    end if
 
  end do
 
  nt = nt0
  nl = nl0
 
  return
end
function idxchg ( x, y, i1, i2, i3, i4 )
!
!*******************************************************************************
!
!! IDXCHG determines whether two triangles should be exchanged.
!
!
!  Discussion:
!
!    The max-min-angle criterion of C L Lawson is used.
!
!  Parameters:
!
!    Input, real X(*), Y(*), the coordinates of the data points.
!
!    Input, integer I1, I2, I3, I4, are the point numbers of
!    four points P1, P2, P3, and P4 that form a quadrilateral,
!    with P3 and P4 connected diagonally.
!
!    Output, integer IDXCHG, reports whether the triangles should be
!    exchanged:
!    0, no exchange is necessary.
!    1, an exchange is necessary.
!
  implicit none
!
  real a1sq
  real a2sq
  real a3sq
  real a4sq
  real c1sq
  real c3sq
  real, parameter :: epsln = 1.0E-06
  integer i1
  integer i2
  integer i3
  integer i4
  integer idx
  integer idxchg
  real s1sq
  real s2sq
  real s3sq
  real s4sq
  real u1
  real u2
  real u3
  real u4
  real x(*)
  real x1
  real x2
  real x3
  real x4
  real y(*)
  real y1
  real y2
  real y3
  real y4
!
!  Preliminary processing
!
  x1 = x(i1)
  y1 = y(i1)
  x2 = x(i2)
  y2 = y(i2)
  x3 = x(i3)
  y3 = y(i3)
  x4 = x(i4)
  y4 = y(i4)

  idx = 0
 
  u3 = (y2-y3)*(x1-x3)-(x2-x3)*(y1-y3)
  u4 = (y1-y4)*(x2-x4)-(x1-x4)*(y2-y4)
 
  if ( u3 * u4 > 0.0E+00 ) then
 
    u1 = (y3-y1)*(x4-x1)-(x3-x1)*(y4-y1)
    u2 = (y4-y2)*(x3-x2)-(x4-x2)*(y3-y2)

    a1sq = (x1-x3)**2+(y1-y3)**2
    a4sq = (x4-x1)**2+(y4-y1)**2
    c1sq = (x3-x4)**2+(y3-y4)**2
    a2sq = (x2-x4)**2+(y2-y4)**2
    a3sq = (x3-x2)**2+(y3-y2)**2
    c3sq = (x2-x1)**2+(y2-y1)**2

    s1sq = u1*u1 / (c1sq*max(a1sq,a4sq))
    s2sq = u2*u2 / (c1sq*max(a2sq,a3sq))
    s3sq = u3*u3 / (c3sq*max(a3sq,a1sq))
    s4sq = u4*u4 / (c3sq*max(a4sq,a2sq))
 
    if ( min ( s3sq, s4sq ) - min ( s1sq, s2sq ) > epsln ) then
      idx = 1
    end if
 
  end if
 
  idxchg = idx
 
  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end