p21 = 0.0E+00
  p23 = -zuu(2)+zuu(1)
  p22 = -1.5E+00*p23

  itpv = it0
!
!  Convert XII and YII to UV system.
!
50 continue

  dx = xii-x0
  dy = yii-y0
  u = ap*dx+bp*dy
  v = cp*dx+dp*dy
!
!  Evaluate the polynomial.
!
  p0 = p00+v*(p01+v*(p02+v*(p03+v*(p04+v*p05))))
  p1 = p10+v*(p11+v*(p12+v*p13))
  p2 = p20+v*(p21+v*(p22+v*p23))
  zii = p0+u*(p1+u*p2)

  return
!
!  Calculation of ZII by extrapolation in the triangle.
!  Check if the necessary coefficients have been calculated.
!
60 continue

  if ( it0 /= itpv ) then
!
!  Load coordinate and partial derivative values at the vertex of the triangle.
!
    jipl = 3*il2-2
    idp = ipl(jipl)
    x0 = xd(idp)
    y0 = yd(idp)
    z0 = zd(idp)
    jpdd = 5*(idp-1)
 
    do kpd = 1, 5
      jpdd = jpdd+1
      pd(kpd) = pdd(jpdd)
    end do
!
!  Calculate the coefficients of the polynomial.
!
    p00 = z0
    p10 = pd(1)
    p01 = pd(2)
    p20 = 0.5E+00*pd(3)
    p11 = pd(4)
    p02 = 0.5E+00*pd(5)
    itpv = it0

  end if
!
!  Convert XII and YII to UV system.
!
  u = xii-x0
  v = yii-y0
!
!  Evaluate the polynomial.
!
  p0 = p00+v*(p01+v*p02)
  p1 = p10+v*p11
  zii = p0+u*(p1+u*p20)
 
  return
end
subroutine idsfft ( md, ndp, xd, yd, zd, nxi, nyi, nzi, xi, yi, zi, iwk, wk )
!
!*******************************************************************************
!
!! IDSFFT fits a smooth surface Z(X,Y) given irregular (X,Y,Z) data.
!
!
!  Discussion:
!
!    IDSFFT performs smooth surface fitting when the projections of the
!    data points in the (X,Y) plane are irregularly distributed.
!
!  Special conditions:
!
!    Inadequate work space IWK and WK may may cause incorrect results.
!
!    The data points must be distinct and their projections in the XY
!    plane must not be collinear, otherwise an error return occurs.
!
!  Parameters:
!
!    Input, integer MD, mode of computation (must be 1, 2, or 3,
!    else an error return will occur).
!
!    1, if this is the first call to this routine, or if the value of 
!    NDP has been changed from the previous call, or if the contents of 
!    the XD or YD arrays have been changed from the previous call.
!
!    2, if the values of NDP and the XD, YD arrays are unchanged from 
!    the previous call, but new values for XI, YI are being used.  If 
!    MD = 2 and NDP has been changed since the previous call to IDSFFT, 
!    an error return occurs.
!
!    3, if the values of NDP, NXI, NYI, XD, YD, XI, YI are unchanged 
!    from the previous call, i.e. if the only change on input to idsfft 
!    is in the ZD array.  If MD = 3 and NDP, nxi or nyi has been changed 
!    since the previous call to idsfft, an error return occurs.
!
!    Between the call with MD = 2 or MD = 3 and the preceding call, the 
!    iwk and wk work arrays should not be disturbed.
!
!    Input, integer NDP, the number of data points.  NDP must be at least 4.
!
!    Input, real XD(NDP), YD(NDP), the X and Y coordinates of the data.
!
!    Input, real ZD(NDP), the data values.
!
!    Input, integer NXI, NYI, the number of output grid points in the
!    X and Y directions.  NXI and NYI must each be at least 1.
!
!    Input, integer NZI, the first dimension of ZI.  NZI must be at
!    least NXI.
!
!    Input, real XI(NXI), YI(NYI), the X and Y coordinates of the grid
!    points.
!
!    Workspace, integer IWK(31*NDP+NXI*NYI).
!
!    Workspace, real WK(6*NDP).
!
!    Output, real ZI(NZI,NYI), contains the interpolated Z values at the
!    grid points.
!
  implicit none
!
  integer ndp
  integer nxi
  integer nyi
  integer nzi
!
  real ap
  real bp
  real cp
  real dp
  integer il1
  integer il2
  integer iti
  integer itpv
  integer iwk(31*ndp + nxi*nyi)
  integer ixi
  integer iyi
  integer izi
  integer jig0mn
  integer jig0mx
  integer jig1mn
  integer jig1mx
  integer jigp
  integer jngp
  integer jwigp
  integer jwigp0
  integer jwipl
  integer jwipt
  integer jwiwl
  integer jwiwp
  integer jwngp
  integer jwngp0
  integer jwwpd
  integer md
  integer ngp0
  integer ngp1
  integer nl
  integer nngp
  integer nt
  real p00
  real p01
  real p02
  real p03
  real p04
  real p05
  real p10
  real p11
  real p12
  real p13
  real p14
  real p20
  real p21
  real p22
  real p23
  real p30
  real p31
  real p32
  real p40
  real p41
  real p50
  real wk(6*ndp)
  real x0
  real xd(ndp)
  real xi(nxi)
  real y0
  real yd(ndp)
  real yi(nyi)
  real zd(ndp)
  real zi(nzi,nyi)
!
  save /idpt/
!
  common /idpt/ itpv,x0,y0,ap,bp,cp,dp, &
                p00,p10,p20,p30,p40,p50,p01,p11,p21,p31,p41, &
                p02,p12,p22,p32,p03,p13,p23,p04,p14,p05
!
!  Error check.
!
  if ( md < 1 .or. md > 3 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IDSFFT - Fatal error!'
    write(*,*)'  Input parameter MD out of range.'
    stop
  end if
 
  if ( ndp < 4 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IDSFFT - Fatal error!'
    write ( *, '(a)' ) '  Input parameter NDP out of range.'
    stop
  end if
 
  if ( nxi < 1 .or. nyi < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IDSFFT - Fatal error!'
    write ( *, '(a)' ) '  Input parameter NXI or NYI out of range.'
    stop
  end if
 
  if ( nxi > nzi ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IDSFFT - Fatal error!'
    write ( *, '(a)' ) '  Input parameter NZI is less than NXI.'
    stop
  end if
 
  if ( md <= 1 ) then

    iwk(1) = ndp

  else
 
    if ( ndp /= iwk(1) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'IDSFFT - Fatal error!'
      write ( *, '(a)' ) '  MD = 2 or 3 but ndp was changed since last call.'
      stop
    end if

  end if

  if ( md <= 2 ) then

    iwk(3) = nxi
    iwk(4) = nyi
    
  else
 
    if ( nxi /= iwk(3) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'IDSFFT - Fatal error!'
      write ( *, '(a)' ) 'MD = 3 but nxi was changed since last call.'
      stop
    end if
 
    if ( nyi /= iwk(4) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'IDSFFT - Fatal error!'
      write ( *, '(a)' ) '  MD = 3 but nyi was changed since last call.'
      stop
    end if

  end if
!
!  Allocation of storage areas in the IWK array.
!
  jwipt = 16
  jwiwl = 6*ndp+1
  jwngp0 = jwiwl-1
  jwipl = 24*ndp+1
  jwiwp = 30*ndp+1
  jwigp0 = 31*ndp
  jwwpd = 5*ndp+1
!
!  Triangulate the XY plane.
!
  if ( md == 1 ) then
 
    call idtang ( ndp, xd, yd, nt, iwk(jwipt), nl, iwk(jwipl), &
      iwk(jwiwl), iwk(jwiwp), wk )

    iwk(5) = nt
    iwk(6) = nl
 
    if ( nt == 0 ) then
      return
    end if
 
  else

    nt = iwk(5)
    nl = iwk(6)

  end if
!
!  Sort output grid points in ascending order of the triangle
!  number and the border line segment number.
! 
  if ( md <= 2 ) then
 
    call idgrid ( xd, yd, nt, iwk(jwipt), nl, iwk(jwipl), nxi, &
      nyi, xi, yi, iwk(jwngp0+1), iwk(jwigp0+1) )
 
  end if
!
!  Estimate partial derivatives at all data points.
!
  call idpdrv ( ndp, xd, yd, zd, nt, iwk(jwipt), wk, wk(jwwpd) )
!
!  Interpolate the ZI values.
!
  itpv = 0
  jig0mx = 0
  jig1mn = nxi*nyi+1
  nngp = nt+2*nl
 
  do jngp = 1, nngp

    iti = jngp

    if ( jngp > nt ) then
      il1 = (jngp-nt+1)/2
      il2 = (jngp-nt+2)/2
      if(il2>nl) then
        il2 = 1
      end if
      iti = il1*(nt+nl)+il2
    end if

    jwngp = jwngp0+jngp
    ngp0 = iwk(jwngp)

    if ( ngp0 /= 0 ) then

      jig0mn = jig0mx+1
      jig0mx = jig0mx+ngp0
 
      do jigp = jig0mn, jig0mx

        jwigp = jwigp0+jigp
        izi = iwk(jwigp)
        iyi = (izi-1)/nxi+1
        ixi = izi-nxi*(iyi-1)

        call idptip ( ndp, xd, yd, zd, nt, iwk(jwipt), nl, iwk(jwipl), &
          wk, iti, xi(ixi), yi(iyi), zi(ixi,iyi) )

      end do
 
    end if

    jwngp = jwngp0+2*nngp+1-jngp
    ngp1 = iwk(jwngp)

    if ( ngp1 /= 0 ) then

      jig1mx = jig1mn-1
      jig1mn = jig1mn-ngp1
 
      do jigp = jig1mn, jig1mx

        jwigp = jwigp0+jigp
        izi = iwk(jwigp)
        iyi = (izi-1)/nxi+1
        ixi = izi-nxi*(iyi-1)

        call idptip ( ndp, xd, yd, zd, nt, iwk(jwipt), nl, iwk(jwipl), &
          wk, iti, xi(ixi), yi(iyi), zi(ixi,iyi) )

      end do

    end if
 
  end do
 
  return
end
subroutine idtang ( ndp, xd, yd, nt, ipt, nl, ipl, iwl, iwp, wk )
!
!*******************************************************************************
!
!! IDTANG performs triangulation.
!
!
!  Discussion:
!
!    The routine divides the XY plane into a number of triangles according to
!    given data points in the plane, determines line segments that form
!    the border of data area, and determines the triangle numbers
!    corresponding to the border line segments.
!
!    At completion, point numbers of the vertexes of each triangle
!    are listed counter-clockwise.  Point numbers of the end points
!    of each border line segment are listed counter-clockwise,
!    listing order of the line segments being counter-clockwise.
!
!  Parameters:
!
!    Input, integer NDP, the number of data points.
!
!    Input, real XD(NDP), YD(NDP), the X and Y coordinates of the data.
!
!    Output, integer NT, the number of triangles,
!
!    Output, integer IPT(6*NDP-15), where the point numbers of the 
!    vertexes of the IT-th triangle are to be stored as entries
!    3*IT-2, 3*IT-1, and 3*IT, for IT = 1 to NT.
!
!    Output, integer NL, the number of border line segments.
!
!    Output, integer IPL(6*NDP), where the point numbers of the end 
!    points of the (il)th border line segment and its respective triangle
!    number are to be stored as the (3*il-2)nd, (3*il-1)st, and (3*il)th
!    elements, il = 1,2,..., nl.
!
!    Workspace, integer IWL(18*NDP),
!
!    Workspace, integer IWP(NDP),
!
!    Workspace, real WK(NDP).
!
  implicit none
!
  integer ndp
!
  real dsqf
  real dsqi
  real dsqmn
  real, parameter :: epsln = 1.0E-06
  integer idxchg
  integer il
  integer ilf
  integer iliv
  integer ilt3
  integer ilvs
  integer ip
  integer ip1
  integer ip1p1
  integer ip2
  integer ip3
  integer ipl(6*ndp)
  integer ipl1
  integer ipl2
  integer iplj1
  integer iplj2
  integer ipmn1
  integer ipmn2
  integer ipt(6*ndp-15)
  integer ipt1
  integer ipt2
  integer ipt3
  integer ipti
  integer ipti1
  integer ipti2
  integer irep
  integer it
  integer it1t3
  integer it2t3
  integer itf(2)
  integer its
  integer itt3
  integer itt3r
  integer iwl(18*ndp)
  integer iwp(ndp)
  integer ixvs
  integer ixvspv
  integer jl1
  integer jl2
  integer jlt3
  integer jp
  integer jp1
  integer jp2
  integer jpc
  integer jpmn
  integer jpmx
  integer jwl
  integer jwl1
  integer jwl1mn
  integer nl
  integer nl0
  integer nlf
  integer nlfc
  integer nlft2
  integer nln
  integer nlnt3
  integer nlsh
  integer nlsht3
  integer nlt3
  integer, parameter :: nrep = 100
  integer nt
  integer nt0
  integer ntf
  integer ntt3
  integer ntt3p3
  real sp
  real spdt
  real u1
  real u2
  real u3
  real v1
  real v2
  real v3
  real vp
  real vpdt
  real wk(ndp)
  real x1
  real x2
  real x3
  real xd(ndp)
  real xdmp
  real y1
  real y2
  real y3
  real yd(ndp)
  real ydmp
!
!  Statement functions
!
  dsqf(u1,v1,u2,v2) = (u2-u1)**2+(v2-v1)**2
  spdt(u1,v1,u2,v2,u3,v3) = (u2-u1)*(u3-u1)+(v2-v1)*(v3-v1)
  vpdt(u1,v1,u2,v2,u3,v3) = (v3-v1)*(u2-u1)-(u3-u1)*(v2-v1)
!
!  Preliminary processing
!
  if ( ndp < 4 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IDTANG - Fatal error!'
    write ( *, '(a)' ) '  Input parameter NDP out of range.'
    stop
  end if
!
!  Determine IPMN1 and IPMN2, the closest pair of data points.
!
  dsqmn = dsqf(xd(1),yd(1),xd(2),yd(2))
  ipmn1 = 1
  ipmn2 = 2
 
  do ip1 = 1, ndp-1
 
    x1 = xd(ip1)
    y1 = yd(ip1)
    ip1p1 = ip1+1
 
    do ip2 = ip1p1, ndp
 
      dsqi = dsqf(x1,y1,xd(ip2),yd(ip2))
 
      if ( dsqi == 0.0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'IDTANG - Fatal error!'
        write ( *, '(a)' ) '  Two of the input data points are identical.'
        stop
      end if
 
      if(dsqi<dsqmn) then
        dsqmn = dsqi
        ipmn1 = ip1
        ipmn2 = ip2
      end if
 
    end do
 
  end do
!
!  Compute the midpoint of the closest two data points.
!
  xdmp = (xd(ipmn1)+xd(ipmn2)) / 2.0E+00
  ydmp = (yd(ipmn1)+yd(ipmn2)) / 2.0E+00
!
!  Sort the other (NDP-2) data points in ascending order of
!  distance from the midpoint and store the sorted data point
!  numbers in the IWP array.
!
  jp1 = 2