!    at the ith data point are to be stored as  the (5*i-4)th, (5*i-3)rd, 
!    (5*i-2)nd, (5*i-1)st and (5*i)th elements, respectively, where i = 
!    1, 2, ..., ndp.
!
!    Workspace, real WK(NDP).
!
  implicit none
!
  integer ndp
  integer nt
!
  real d12
  real d23
  real d31
  real dx1
  real dx2
  real dy1
  real dy2
  real dz1
  real dz2
  real dzx1
  real dzx2
  real dzy1
  real dzy2
  real, parameter :: epsln = 1.0E-06
  integer idp
  integer ipt(3*nt)
  integer ipti(3)
  integer it
  integer iv
  integer jpd
  integer jpd0
  integer jpdmx
  integer jpt
  integer jpt0
  integer nt0
  real pd(5*ndp)
  real vpx
  real vpxx
  real vpxy
  real vpy
  real vpyx
  real vpyy
  real vpz
  real vpzmn
  real w1(3)
  real w2(3)
  real wi
  real wk(ndp)
  real xd(ndp)
  real xv(3)
  real yd(ndp)
  real yv(3)
  real zd(ndp)
  real zv(3)
  real zxv(3)
  real zyv(3)
!
!  Preliminary processing.
!
  nt0 = nt
!
!  Clear the PD array.
!
  jpdmx = 5*ndp
 
  pd(1:jpdmx) = 0.0E+00
 
  wk(1:ndp) = 0.0E+00
!
!  Estimate ZX and ZY.
!
  do it = 1, nt0
 
    jpt0 = 3*(it-1)
 
    do iv = 1, 3
      jpt = jpt0+iv
      idp = ipt(jpt)
      ipti(iv) = idp
      xv(iv) = xd(idp)
      yv(iv) = yd(idp)
      zv(iv) = zd(idp)
    end do
 
    dx1 = xv(2)-xv(1)
    dy1 = yv(2)-yv(1)
    dz1 = zv(2)-zv(1)
    dx2 = xv(3)-xv(1)
    dy2 = yv(3)-yv(1)
    dz2 = zv(3)-zv(1)
    vpx = dy1*dz2-dz1*dy2
    vpy = dz1*dx2-dx1*dz2
    vpz = dx1*dy2-dy1*dx2
    vpzmn = abs(dx1*dx2+dy1*dy2)*epsln
 
    if ( abs(vpz) > vpzmn ) then
 
      d12 = sqrt((xv(2)-xv(1))**2+(yv(2)-yv(1))**2)
      d23 = sqrt((xv(3)-xv(2))**2+(yv(3)-yv(2))**2)
      d31 = sqrt((xv(1)-xv(3))**2+(yv(1)-yv(3))**2)
      w1(1) = 1.0E+00 / (d31*d12)
      w1(2) = 1.0E+00 / (d12*d23)
      w1(3) = 1.0E+00 / (d23*d31)
      w2(1) = vpz*w1(1)
      w2(2) = vpz*w1(2)
      w2(3) = vpz*w1(3)
 
      do iv = 1, 3
        idp = ipti(iv)
        jpd0 = 5*(idp-1)
        wi = (w1(iv)**2)*w2(iv)
        pd(jpd0+1) = pd(jpd0+1)+vpx*wi
        pd(jpd0+2) = pd(jpd0+2)+vpy*wi
        wk(idp) = wk(idp)+vpz*wi
      end do
 
    end if
 
  end do
 
  do idp = 1, ndp
    jpd0 = 5*(idp-1)
    pd(jpd0+1) = -pd(jpd0+1)/wk(idp)
    pd(jpd0+2) = -pd(jpd0+2)/wk(idp)
  end do
!
!  Estimate ZXX, ZXY, and ZYY.
!
  do it = 1, nt0
 
    jpt0 = 3*(it-1)
 
    do iv = 1, 3
      jpt = jpt0+iv
      idp = ipt(jpt)
      ipti(iv) = idp
      xv(iv) = xd(idp)
      yv(iv) = yd(idp)
      jpd0 = 5*(idp-1)
      zxv(iv) = pd(jpd0+1)
      zyv(iv) = pd(jpd0+2)
    end do
 
    dx1 = xv(2)-xv(1)
    dy1 = yv(2)-yv(1)
    dzx1 = zxv(2)-zxv(1)
    dzy1 = zyv(2)-zyv(1)
    dx2 = xv(3)-xv(1)
    dy2 = yv(3)-yv(1)
    dzx2 = zxv(3)-zxv(1)
    dzy2 = zyv(3)-zyv(1)
    vpxx = dy1*dzx2-dzx1*dy2
    vpxy = dzx1*dx2-dx1*dzx2
    vpyx = dy1*dzy2-dzy1*dy2
    vpyy = dzy1*dx2-dx1*dzy2
    vpz = dx1*dy2-dy1*dx2
    vpzmn = abs(dx1*dx2+dy1*dy2)*epsln
 
    if ( abs(vpz) > vpzmn ) then
 
      d12 = sqrt((xv(2)-xv(1))**2+(yv(2)-yv(1))**2)
      d23 = sqrt((xv(3)-xv(2))**2+(yv(3)-yv(2))**2)
      d31 = sqrt((xv(1)-xv(3))**2+(yv(1)-yv(3))**2)
      w1(1) = 1.0E+00 /(d31*d12)
      w1(2) = 1.0E+00 /(d12*d23)
      w1(3) = 1.0E+00 /(d23*d31)
      w2(1) = vpz*w1(1)
      w2(2) = vpz*w1(2)
      w2(3) = vpz*w1(3)
 
      do iv = 1, 3
        idp = ipti(iv)
        jpd0 = 5*(idp-1)
        wi = (w1(iv)**2)*w2(iv)
        pd(jpd0+3) = pd(jpd0+3)+vpxx*wi
        pd(jpd0+4) = pd(jpd0+4)+(vpxy+vpyx)*wi
        pd(jpd0+5) = pd(jpd0+5)+vpyy*wi
      end do
 
    end if
 
  end do
 
  do idp = 1, ndp
    jpd0 = 5*(idp-1)
    pd(jpd0+3) = -pd(jpd0+3)/wk(idp)
    pd(jpd0+4) = -pd(jpd0+4)/(2.0*wk(idp))
    pd(jpd0+5) = -pd(jpd0+5)/wk(idp)
  end do
 
  return
end
subroutine idptip ( ndp,xd, yd, zd, nt, ipt, nl, ipl, pdd, iti, xii, yii, zii )
!
!*******************************************************************************
!
!! IDPTIP performs interpolation, determining a value of Z given X and Y.
!
!
!  Modified:
!
!    19 February 2001
!
!  Parameters:
!
!    Input, integer NDP, the number of data values.
!
!    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, ipt = integer array of dimension 3*nt containing 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 PDD(5*NDP). the partial derivatives at the data points,
!
!    Input, integer ITI, triangle number of the triangle in which lies
!    the point for which interpolation is to be performed,
!
!    Input, real XII, YII, the X and Y coordinates of the point for which
!    interpolation is to be performed.
!
!    Output, real ZII, the interpolated Z value.
!
  implicit none
!
  integer ndp
  integer nl
  integer nt
!
  real a
  real aa
  real ab
  real ac
  real act2
  real ad
  real adbc
  real ap
  real b
  real bb
  real bc
  real bdt2
  real bp
  real c
  real cc
  real cd
  real cp
  real csuv
  real d
  real dd
  real dlt
  real dp
  real dx
  real dy
  real g1
  real g2
  real h1
  real h2
  real h3
  integer i
  integer idp
  integer il1
  integer il2
  integer ipl(3*nl)
  integer ipt(3*nt)
  integer it0
  integer iti
  integer itpv
  integer jipl
  integer jipt
  integer jpd
  integer jpdd
  integer kpd
  integer ntl
  real lu
  real lv
  real p0
  real p00
  real p01
  real p02
  real p03
  real p04
  real p05
  real p1
  real p10
  real p11
  real p12
  real p13
  real p14
  real p2
  real p20
  real p21
  real p22
  real p23
  real p3
  real p30
  real p31
  real p32
  real p4
  real p40
  real p41
  real p5
  real p50
  real pd(15)
  real pdd(5*ndp)
  real thsv
  real thus
  real thuv
  real thxu
  real u
  real v
  real x(3)
  real x0
  real xd(*)
  real xii
  real y(3)
  real y0
  real yd(*)
  real yii
  real z(3)
  real z0
  real zd(*)
  real zii
  real zu(3)
  real zuu(3)
  real zuv(3)
  real zv(3)
  real zvv(3)
!
  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
!
!  Preliminary processing
!
  it0 = iti
  ntl = nt+nl

  if ( it0 > ntl ) then
    il1 = it0/ntl
    il2 = it0-il1*ntl
    if(il1==il2)      go to 40
    go to 60
  end if
!
!  Calculation of ZII by interpolation.
!  Check if the necessary coefficients have been calculated.
!
  if ( it0 == itpv )     go to 30
!
!  Load coordinate and partial derivative values at the vertexes.
!
  jipt = 3*(it0-1)
  jpd = 0
 
  do i = 1, 3
 
    jipt = jipt+1
    idp = ipt(jipt)
    x(i) = xd(idp)
    y(i) = yd(idp)
    z(i) = zd(idp)
    jpdd = 5*(idp-1)
 
    do kpd = 1, 5
      jpd = jpd+1
      jpdd = jpdd+1
      pd(jpd) = pdd(jpdd)
    end do
 
  end do
!
!  Determine the coefficients for the coordinate system
!  transformation from the XY system to the UV system and vice versa.
!
  x0 = x(1)
  y0 = y(1)
  a = x(2)-x0
  b = x(3)-x0
  c = y(2)-y0
  d = y(3)-y0
  ad = a*d
  bc = b*c
  dlt = ad-bc
  ap =  d/dlt
  bp = -b/dlt
  cp = -c/dlt
  dp =  a/dlt
!
!  Convert the partial derivatives at the vertexes of the
!  triangle for the UV coordinate system.
!
  aa = a*a
  act2 = 2.0E+00 *a*c
  cc = c*c
  ab = a*b
  adbc = ad+bc
  cd = c*d
  bb = b*b
  bdt2 = 2.0E+00 *b*d
  dd = d*d
 
  do i = 1, 3
    jpd = 5*i
    zu(i) = a*pd(jpd-4)+c*pd(jpd-3)
    zv(i) = b*pd(jpd-4)+d*pd(jpd-3)
    zuu(i) = aa*pd(jpd-2)+act2*pd(jpd-1)+cc*pd(jpd)
    zuv(i) = ab*pd(jpd-2)+adbc*pd(jpd-1)+cd*pd(jpd)
    zvv(i) = bb*pd(jpd-2)+bdt2*pd(jpd-1)+dd*pd(jpd)
  end do
!
!  Calculate the coefficients of the polynomial.
!
  p00 = z(1)
  p10 = zu(1)
  p01 = zv(1)
  p20 = 0.5E+00 * zuu(1)
  p11 = zuv(1)
  p02 = 0.5E+00 * zvv(1)
  h1 = z(2)-p00-p10-p20
  h2 = zu(2)-p10-zuu(1)
  h3 = zuu(2)-zuu(1)
  p30 =  10.0E+00 * h1 - 4.0E+00 * h2 + 0.5E+00 * h3
  p40 = -15.0E+00 * h1 + 7.0E+00 * h2           - h3
  p50 =   6.0E+00 * h1 - 3.0E+00 * h2 + 0.5E+00 * h3
  h1 = z(3)-p00-p01-p02
  h2 = zv(3)-p01-zvv(1)
  h3 = zvv(3)-zvv(1)
  p03 =  10.0E+00 * h1 - 4.0E+00 * h2 + 0.5E+00 * h3
  p04 = -15.0E+00 * h1 + 7.0E+00 * h2    -h3
  p05 =   6.0E+00 * h1 - 3.0E+00 * h2 + 0.5E+00 * h3
  lu = sqrt(aa+cc)
  lv = sqrt(bb+dd)
  thxu = atan2(c,a)
  thuv = atan2(d,b)-thxu
  csuv = cos(thuv)
  p41 = 5.0E+00*lv*csuv/lu*p50
  p14 = 5.0E+00*lu*csuv/lv*p05
  h1 = zv(2)-p01-p11-p41
  h2 = zuv(2)-p11-4.0E+00 * p41
  p21 =  3.0E+00 * h1-h2
  p31 = -2.0E+00 * h1+h2
  h1 = zu(3)-p10-p11-p14
  h2 = zuv(3)-p11- 4.0E+00 * p14
  p12 =  3.0E+00 * h1-h2
  p13 = -2.0E+00 * h1+h2
  thus = atan2(d-c,b-a)-thxu
  thsv = thuv-thus
  aa =  sin(thsv)/lu
  bb = -cos(thsv)/lu
  cc =  sin(thus)/lv
  dd =  cos(thus)/lv
  ac = aa*cc
  ad = aa*dd
  bc = bb*cc
  g1 = aa * ac*(3.0E+00*bc+2.0E+00*ad)
  g2 = cc * ac*(3.0E+00*ad+2.0E+00*bc)
  h1 = -aa*aa*aa*(5.0E+00*aa*bb*p50+(4.0E+00*bc+ad)*p41) &
       -cc*cc*cc*(5.0E+00*cc*dd*p05+(4.0E+00*ad+bc)*p14)
  h2 = 0.5E+00 * zvv(2)-p02-p12
  h3 = 0.5E+00 * zuu(3)-p20-p21
  p22 = (g1*h2+g2*h3-h1)/(g1+g2)
  p32 = h2-p22
  p23 = h3-p22
  itpv = it0
!
!  Convert XII and YII to UV system.
!
30 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+v*p14)))
  p2 = p20+v*(p21+v*(p22+v*p23))
  p3 = p30+v*(p31+v*p32)
  p4 = p40+v*p41
  p5 = p50
  zii = p0+u*(p1+u*(p2+u*(p3+u*(p4+u*p5))))
  return
!
!  Calculation of ZII by extrapolation in the rectangle.
!  Check if the necessary coefficients have been calculated.
!
40 continue

  if(it0==itpv)     go to 50
!
!  Load coordinate and partial derivative values at the end
!  points of the border line segment.
!
  jipl = 3*(il1-1)
  jpd = 0
 
  do i = 1, 2
 
    jipl = jipl+1
    idp = ipl(jipl)
    x(i) = xd(idp)
    y(i) = yd(idp)
    z(i) = zd(idp)
    jpdd = 5*(idp-1)
 
    do kpd = 1, 5
      jpd = jpd+1
      jpdd = jpdd+1
      pd(jpd) = pdd(jpdd)
    end do
 
  end do
!
!  Determine the coefficients for the coordinate system
!  transformation from the XY system to the UV system
!  and vice versa.
!
  x0 = x(1)
  y0 = y(1)
  a = y(2)-y(1)
  b = x(2)-x(1)
  c = -b
  d = a
  ad = a*d
  bc = b*c
  dlt = ad-bc
  ap =  d/dlt
  bp = -b/dlt
  cp = -bp
  dp =  ap
!
!  Convert the partial derivatives at the end points of the
!  border line segment for the UV coordinate system.
!
  aa = a*a
  act2 = 2.0E+00 * a * c
  cc = c*c
  ab = a*b
  adbc = ad+bc
  cd = c*d
  bb = b*b
  bdt2 = 2.0E+00 * b * d
  dd = d*d
 
  do i = 1, 2
    jpd = 5*i
    zu(i) = a*pd(jpd-4)+c*pd(jpd-3)
    zv(i) = b*pd(jpd-4)+d*pd(jpd-3)
    zuu(i) = aa*pd(jpd-2)+act2*pd(jpd-1)+cc*pd(jpd)
    zuv(i) = ab*pd(jpd-2)+adbc*pd(jpd-1)+cd*pd(jpd)
    zvv(i) = bb*pd(jpd-2)+bdt2*pd(jpd-1)+dd*pd(jpd)
  end do
!
!  Calculate the coefficients of the polynomial.
!
  p00 = z(1)
  p10 = zu(1)
  p01 = zv(1)
  p20 = 0.5E+00 * zuu(1)
  p11 = zuv(1)
  p02 = 0.5E+00 * zvv(1)

  h1 = z(2)-p00-p01-p02
  h2 = zv(2)-p01-zvv(1)
  h3 = zvv(2)-zvv(1)

  p03 =  10.0E+00 * h1 - 4.0E+00*h2+0.5E+00*h3
  p04 = -15.0E+00 * h1 + 7.0E+00*h2    -h3
  p05 =   6.0E+00 * h1 - 3.0E+00*h2+0.5E+00*h3

  h1 = zu(2)-p10-p11
  h2 = zuv(2)-p11

  p12 =  3.0E+00*h1-h2
  p13 = -2.0E+00*h1+h2