ymx = yimx

    if ( y2 >= y1 ) then
      xmn = min ( x1, x2 )
    end if

    if ( y2 <= y1 ) then
      xmx = max ( x1, x2 )
    end if

    if ( x2 <= x1 ) then
      ymn = min ( y1, y2 )
    end if

    if ( x2 >= x1 ) then
      ymx = max ( y1, y2 )
    end if

    insd = 0

    do ixi = 1, nxi

      if ( xi(ixi) < xmn .or. xi(ixi) > xmx ) then
        if ( insd == 0 ) then
          go to 42
        end if
        iximx = ixi-1
        go to 43
      end if

      if ( insd /= 1 ) then
        insd = 1
        iximn = ixi
      end if

42    continue

    end do

    if ( insd == 0 ) go to 58

    iximx = nxi

43  continue

    do iyi = 1, nyi

      yii = yi(iyi)
      if(yii<ymn.or.yii>ymx)        go to 57

      do ixi = iximn,iximx

        xii = xi(ixi)
        l = 0
        if(vpdt(x1,y1,x2,y2,xii,yii))     46,45,56
   45       l = 1
   46       if(spdt(x2,y2,x1,y1,xii,yii))     56,47,48
   47       l = 1
   48       if(spdt(x1,y1,x2,y2,xii,yii))     56,49,50
   49       l = 1
   50       izi = nxi*(iyi-1)+ixi
        if(l==1)    go to 51
        ngp0 = ngp0+1
        jigp0 = jigp0+1
        igp(jigp0) = izi
        go to 56
 
51      continue
 
        do jigp1i = jigp1,nxinyi
          if(izi==igp(jigp1i))     go to 56
        end do
 
53      continue

        ngp1 = ngp1+1
        jigp1 = jigp1-1
        igp(jigp1) = izi

56      continue

      end do

57    continue

    end do

58  continue

    jngp0 = jngp0+1
    ngp(jngp0) = ngp0
    jngp1 = jngp1-1
    ngp(jngp1) = ngp1
!
!  In semi-infinite triangular area.
!
60  continue

    ngp0 = 0
    ngp1 = 0
    ilp1 = mod(il0,nl0)+1
    ilp1t3 = ilp1*3
    ip3 = ipl(ilp1t3-1)
    x3 = xd(ip3)
    y3 = yd(ip3)
    xmn = ximn
    xmx = ximx
    ymn = yimn
    ymx = yimx
    if(y3>=y2.and.y2>=y1)   xmn = x2
    if(y3<=y2.and.y2<=y1)   xmx = x2
    if(x3<=x2.and.x2<=x1)   ymn = y2
    if(x3>=x2.and.x2>=x1)   ymx = y2
    insd = 0

    do ixi = 1, nxi

      if ( xi(ixi) < xmn .or. xi(ixi) > xmx ) then
        if(insd==0)   go to 62
        iximx = ixi-1
        go to 63
      end if

      if ( insd /= 1 ) then
        insd = 1
        iximn = ixi
      end if

62    continue

    end do

    if(insd==0)     go to 78

    iximx = nxi
 
63  continue

    do iyi = 1, nyi
 
      yii = yi(iyi)
      if(yii<ymn.or.yii>ymx)        go to 77
 
      do ixi = iximn, iximx
 
        xii = xi(ixi)
        l = 0
        if ( spdt(x1,y1,x2,y2,xii,yii) )     66,65,76
   65       l = 1
   66       if(spdt(x3,y3,x2,y2,xii,yii))     70,67,76
   67       l = 1
   70       izi = nxi*(iyi-1)+ixi
 
        if ( l /= 1 ) then
          ngp0 = ngp0+1
          jigp0 = jigp0+1
          igp(jigp0) = izi
          go to 76
        end if
 
        do jigp1i = jigp1, nxinyi
          if(izi==igp(jigp1i)) go to 76
        end do
 
        ngp1 = ngp1+1
        jigp1 = jigp1-1
        igp(jigp1) = izi
 
76      continue

      end do
 
77    continue

    end do
 
78  continue

    jngp0 = jngp0+1
    ngp(jngp0) = ngp0
    jngp1 = jngp1-1
    ngp(jngp1) = ngp1
 
  end do
 
  return
end
subroutine idlctn ( ndp, xd, yd, nt, ipt, nl, ipl, xii, yii, iti, iwk, wk )
!
!*******************************************************************************
!
!! IDLCTN finds the triangle that contains a point.
!
!
!  Discusstion:
!
!    IDLCTN determines what triangle a given point (XII, YII) belongs to.  
!    When the given point does not lie inside the data area, IDLCTN 
!    determines the border line segment when the point lies in an outside
!    rectangular area, and two border line segments when the point
!    lies in an outside triangular area.
!
!  Parameters:
!
!    Input, integer NDP, the number of data points.
!
!    Input, real XD(NDP), YD(NDP), the X and Y coordinates of the data.
!
!    Input, integer NT, the number of triangles.
!
!    Input, integer IPT(3*NT), the point numbers of the vertexes of 
!    the triangles,
!
!    Input, integer NL, the number of border line segments.
!
!    Input, integer IPL(3*NL), the point numbers of the end points of 
!    the border line segments and their respective triangle numbers.
!
!    Input, real XII, YII, the coordinates of the point to be located.
!
!    Output, integer ITI, the triangle number, when the point is inside the
!    data area, or two border line segment numbers, il1 and il2,
!    coded to il1*(nt+nl)+il2, when the point is outside the data area.
!
!    Workspace, integer IWK(18*NDP).
!
!    Workspace, real WK(8*NDP).
!
  implicit none
!
  integer ndp
  integer nl
  integer nt
!
  integer i1
  integer i2
  integer i3
  integer idp
  integer idsc(9)
  integer il1
  integer il1t3
  integer il2
  integer ip1
  integer ip2
  integer ip3
  integer ipl(3*nl)
  integer ipt(3*nt)
  integer isc
  integer it0
  integer it0t3
  integer iti
  integer itipv
  integer itsc
  integer iwk(18*ndp)
  integer jiwk
  integer jwk
  integer nl0
  integer nt0
  integer ntl
  integer ntsc
  integer ntsci
  real spdt
  real u1
  real u2
  real u3
  real v1
  real v2
  real v3
  real vpdt
  real wk(8*ndp)
  real x0
  real x1
  real x2
  real x3
  real xd(ndp)
  real xii
  real xmn
  real xmx
  real xs1
  real xs2
  real y0
  real y1
  real y2
  real y3
  real yd(ndp)
  real yii
  real ymn
  real ymx
  real ys1
  real ys2
!
  save /idlc/
!
  common /idlc/ itipv,xs1,xs2,ys1,ys2,ntsc(9)
!
!  Statement functions
!
  spdt(u1,v1,u2,v2,u3,v3) = (u1-u2)*(u3-u2)+(v1-v2)*(v3-v2)

  vpdt(u1,v1,u2,v2,u3,v3) = (u1-u3)*(v2-v3)-(v1-v3)*(u2-u3)
!
!  Preliminary processing
!
  nt0 = nt
  nl0 = nl
  ntl = nt0+nl0
  x0 = xii
  y0 = yii
!
!  Processing for a new set of data points
!
  if ( itipv/=0)      go to 30
!
!  Divide the x-y plane into nine rectangular sections.
!
  xmn = xd(1)
  xmx = xd(1)
  ymn = yd(1)
  ymx = yd(1)
  do idp = 2, ndp
    xmn = min ( xd(idp), xmn )
    xmx = max ( xd(idp), xmx )
    ymn = min ( yd(idp), ymn )
    ymx = max ( yd(idp), ymx )
  end do
 
  xs1 = ( xmn + xmn + xmx ) / 3.0E+00
  xs2 = ( xmn + xmx + xmx ) / 3.0E+00
  ys1 = ( ymn + ymn + ymx ) / 3.0E+00
  ys2 = ( ymn + ymx + ymx ) / 3.0E+00
!
!  Determine and store in the iwk array, triangle numbers of
!  the triangles associated with each of the nine sections.
!
  ntsc(1:9) = 0
  idsc(1:9) = 0
 
  it0t3 = 0
  jwk = 0

  do it0 = 1, nt0

    it0t3 = it0t3+3
    i1 = ipt(it0t3-2)
    i2 = ipt(it0t3-1)
    i3 = ipt(it0t3)
    xmn = min ( xd(i1), xd(i2), xd(i3) )
    xmx = max ( xd(i1), xd(i2), xd(i3) )
    ymn = min ( yd(i1), yd(i2), yd(i3) )
    ymx = max ( yd(i1), yd(i2), yd(i3) )

    if ( ymn <= ys1 ) then
      if(xmn<=xs1)                   idsc(1) = 1
      if(xmx>=xs1.and.xmn<=xs2)    idsc(2) = 1
      if(xmx>=xs2)                   idsc(3) = 1
    end if

    if ( ymx >= ys1 .and. ymn <= ys2 ) then
      if(xmn<=xs1)                   idsc(4) = 1
      if(xmx>=xs1.and.xmn<=xs2)    idsc(5) = 1
      if(xmx>=xs2)                   idsc(6) = 1
    end if

    if(ymx<ys2)                   go to 25
    if(xmn<=xs1)                   idsc(7) = 1
    if(xmx>=xs1.and.xmn<=xs2)    idsc(8) = 1
    if(xmx>=xs2)                   idsc(9) = 1

25  continue

    do isc = 1, 9
      if ( idsc(isc) /= 0 ) then
        jiwk = 9*ntsc(isc)+isc
        iwk(jiwk) = it0
        ntsc(isc) = ntsc(isc)+1
        idsc(isc) = 0
      end if
    end do
!
!  Store in the wk array the minimum and maximum of the X and
!  Y coordinate values for each of the triangle.
!
    jwk = jwk+4
    wk(jwk-3) = xmn
    wk(jwk-2) = xmx
    wk(jwk-1) = ymn
    wk(jwk)   = ymx

  end do

  go to 60
!
!  Check if in the same triangle as previous.
!
30 continue

  it0 = itipv

  if(it0>nt0)      go to 40

  it0t3 = it0*3
  ip1 = ipt(it0t3-2)
  x1 = xd(ip1)
  y1 = yd(ip1)
  ip2 = ipt(it0t3-1)
  x2 = xd(ip2)
  y2 = yd(ip2)
  if(vpdt(x1,y1,x2,y2,x0,y0) < 0.0E+00 )      go to 60
  ip3 = ipt(it0t3)
  x3 = xd(ip3)
  y3 = yd(ip3)
  if(vpdt(x2,y2,x3,y3,x0,y0) < 0.0E+00 )      go to 60
  if(vpdt(x3,y3,x1,y1,x0,y0) < 0.0E+00 )      go to 60

  iti = it0
  itipv = it0

  return
!
!  Check if on the same border line segment.
!
40 continue

  il1 = it0 / ntl
  il2 = it0-il1*ntl
  il1t3 = il1*3
  ip1 = ipl(il1t3-2)
  x1 = xd(ip1)
  y1 = yd(ip1)
  ip2 = ipl(il1t3-1)
  x2 = xd(ip2)
  y2 = yd(ip2)
  if(il2/=il1)      go to 50
  if(spdt(x1,y1,x2,y2,x0,y0) < 0.0E+00 )      go to 60
  if(spdt(x2,y2,x1,y1,x0,y0) < 0.0E+00 )      go to 60
  if(vpdt(x1,y1,x2,y2,x0,y0) > 0.0E+00 )      go to 60
  
  iti = it0
  itipv = it0

  return
!
!  Check if between the same two border line segments.
!
50 continue

  if(spdt(x1,y1,x2,y2,x0,y0) > 0.0E+00 )      go to 60

  ip3 = ipl(3*il2-1)
  x3 = xd(ip3)
  y3 = yd(ip3)

  if ( spdt(x3,y3,x2,y2,x0,y0) <= 0.0E+00 )  then
    iti = it0
    itipv = it0
    return
  end if
!
!  Locate inside the data area.
!  Determine the section in which the point in question lies.
!
60 continue
 
  isc = 1

  if ( x0 >= xs1 ) then
    isc = isc+1
  end if

  if ( x0 >= xs2 ) then
    isc = isc+1
  end if

  if ( y0 >= ys1 ) then
    isc = isc+3
  end if

  if ( y0 >= ys2 ) then
    isc = isc+3
  end if
!
!  Search through the triangles associated with the section.
!
  ntsci = ntsc(isc)
  if(ntsci<=0)      go to 70
  jiwk = -9+isc

  do itsc = 1, ntsci

    jiwk = jiwk+9
    it0 = iwk(jiwk)
    jwk = it0*4
    if(x0<wk(jwk-3))    go to 61
    if(x0>wk(jwk-2))    go to 61
    if(y0<wk(jwk-1))    go to 61
    if(y0>wk(jwk))      go to 61
    it0t3 = it0*3
    ip1 = ipt(it0t3-2)
    x1 = xd(ip1)
    y1 = yd(ip1)
    ip2 = ipt(it0t3-1)
    x2 = xd(ip2)
    y2 = yd(ip2)
    if(vpdt(x1,y1,x2,y2,x0,y0)<0.0E+00 )    go to 61
    ip3 = ipt(it0t3)
    x3 = xd(ip3)
    y3 = yd(ip3)

    if ( vpdt(x2,y2,x3,y3,x0,y0) >= 0.0E+00 ) then

      if ( vpdt(x3,y3,x1,y1,x0,y0) >= 0.0E+00 ) then
        iti = it0
        itipv = it0
        return
      end if

    end if

61  continue

  end do
!
!  Locate outside the data area.
!
70 continue

  do il1 = 1, nl0

    il1t3 = il1*3
    ip1 = ipl(il1t3-2)
    x1 = xd(ip1)
    y1 = yd(ip1)
    ip2 = ipl(il1t3-1)
    x2 = xd(ip2)
    y2 = yd(ip2)
    if(spdt(x2,y2,x1,y1,x0,y0)<0.0E+00 )    go to 72
    if(spdt(x1,y1,x2,y2,x0,y0)<0.0E+00 )    go to 71
    if(vpdt(x1,y1,x2,y2,x0,y0)>0.0E+00 )    go to 72
    il2 = il1
    go to 75

   71   continue

    il2 = mod(il1,nl0)+1
    ip3 = ipl(3*il2-1)
    x3 = xd(ip3)
    y3 = yd(ip3)
    if(spdt(x3,y3,x2,y2,x0,y0)<=0.0E+00 )    go to 75

   72   continue

  end do
 
  it0 = 1
  iti = it0
  itipv = it0

  return

   75 continue

  it0 = il1*ntl+il2
  iti = it0
  itipv = it0

  return
end
subroutine idpdrv ( ndp, xd, yd, zd, nt, ipt, pd, wk )
!
!*******************************************************************************
!
!! IDPDRV estimates first and second partial derivatives at data points.
!
!
!  Parameters:
!
!    Input, integer NDP, the number of data points.
!
!    Input, real XD(NDP), YD(NDP), the X and Y coordinates of the data.
!
!    Input, real ZD(NDP), the data values.
!
!    Input, integer NT, the number of triangles.
!
!    Input, integer IPT(3*NT), the point numbers of the vertexes of the 
!    triangles.
!
!    Output, real PD(5*NDP), the estimated zx, zy, zxx, zxy, and zyy values