subroutine idbvip ( md, ndp, xd, yd, zd, nip, xi, yi, zi, iwk, wk )
!
!*******************************************************************************
!
!! IDBVIP performs bivariate interpolation of irregular X, Y data.
!
!
!  Discussion:
!
!    The data points must be distinct and their projections in the
!    X-Y plane must not be collinear, otherwise an error return
!    occurs.
!
!    Inadequate work space IWK and WK may cause incorrect results.
!
!  Latest revision:
!
!    January, 1985
!
!  Purpose:
!
!    To provide bivariate interpolation and smooth surface fitting for 
!    values given at irregularly distributed points.
!
!    The resulting interpolating function and its first-order partial 
!    derivatives are continuous.
!
!    The method employed is local, i.e. a change in the data in one area 
!    of the plane does not affect the interpolating function except in
!    that local area.  This is advantageous over global interpolation 
!    methods.
!
!    Also, the method gives exact results when all points lie in a plane. 
!    This is advantageous over other methods such as two-dimensional
!    Fourier series interpolation.
!
!  Usage:
!
!    This package contains two user entries, IDBVIP and IDSFFT, both 
!    requiring input data to be given at points
!      ( X(I), Y(I) ), I = 1,...,N.
!
!    If the user desires the interpolated data to be output at grid 
!    points, i.e. at points
!      ( XI(I), YI(J) ), I = 1,...,NXI, J=1,...,NYI,
!    then IDSFFT should be used.  This is useful for generating an 
!    interpolating surface.
!
!    The other user entry point, IDBVIP, will produce interpolated 
!    values at scattered points
!      ( XI(I), YI(I) ), i = 1,...,NIP.  
!    This is useful for filling in missing data points on a grid.
!
!  History:
!
!    The original version of BIVAR was written by Hiroshi Akima in 
!    August 1975 and rewritten by him in late 1976.  It was incorporated 
!    into NCAR's public software libraries in January 1977.  In August 
!    1984 a new version of BIVAR, incorporating changes described in the 
!    Rocky Mountain Journal of Mathematics article cited below, was 
!    obtained from Dr Akima by Michael Pernice of NCAR's Scientific 
!    Computing Division, who evaluated it and made it available in February, 
!    1985.
!
!  Accuracy:
!
!    Accurate to machine precision on the input data points.  Accuracy at 
!    other points greatly depends on the input data.
!
!  References:
!
!    Hiroshi Akima,  
!    A Method of Bivariate Interpolation and Smooth Surface Fitting
!      for Values Given at Irregularly Distributed Points,
!    ACM Transactions on Mathematical Software, 
!    Volume 4, Number 2, June 1978.
!
!    Hiroshi Akima,  
!    On Estimating Partial Derivatives for Bivariate Interpolation
!      of Scattered Data,
!    Rocky Mountain Journal of Mathematics,
!    Volume 14, Number 1, Winter 1984.
!
!  Method:
!
!    The XY plane is divided into triangular cells, each cell having 
!    projections of three data points in the plane as its vertices, and
!    a bivariate quintic polynomial in X and Y is fitted to each 
!    triangular cell.
!
!    The coefficients in the fitted quintic polynomials are determined 
!    by continuity requirements and by estimates of partial derivatives 
!    at the vertices and along the edges of the triangles.  The method 
!    described in the rocky mountain journal reference guarantees that 
!    the generated surface depends continuously on the triangulation.
!
!    The resulting interpolating function is invariant under the following 
!    types of linear coordinate transformations:
!      1) a rotation of the XY coordinate system
!      2) linear scale transformation of the Z axis
!      3) tilting of the XY plane, i.e. new coordinates (u,v,w) given by
!           u = x
!           v = y
!           w = z + a*x + b*y
!         where a, b are arbitrary constants.
!
!    complete details of the method are given in the reference publications.
!
!  Parameters:
!
!    Input, integer MD, mode of computation.   MD must be 1,
!    2, or 3, else an error return occurs.
!
!    1: if this is the first call to this subroutine, 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 and 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 IDBVIP, an error return occurs.
!
!    3: if the values of NDP, NIP, XD, YD, XI, YI are unchanged from
!    the previous call, i.e. if the only change on input to idbvip is
!    in the ZD array.  If MD = 3 and NDP or NIP has been changed since
!    the previous call to IDBVIP, 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 (must be 4 or
!    greater, else an error return occurs).
!
!    Input, real XD(NDP), Y(NDP), the X and Y coordinates of the data points.
!
!    Input, real ZD(NDP), the data values at the data points.
!
!    Input, integer NIP, the number of output points at which
!    interpolation is to be performed (must be 1 or greater, else an
!    error return occurs).
!
!    Input, real XI(NIP), YI(NIP), the coordinates of the points at which
!    interpolation is to be performed.
!
!    Output, real ZI(NIP), the interpolated data values.
!
!    Workspace, integer IWK(31*NDP+NIP).
!
!    Workspace, real WK(8*NDP).
!
  implicit none
!
  integer ndp
  integer nip
!
  real ap
  real bp
  real cp
  real dp
  integer iip
  integer itipv
  integer itpv
  integer iwk(31*ndp + nip)
  integer jwipl
  integer jwipt
  integer jwit
  integer jwit0
  integer jwiwk
  integer jwiwl
  integer jwiwp
  integer jwwpd
  integer md
  integer nl
  integer nt
  integer ntsc
  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(8*ndp)
  real x0
  real xd(ndp)
  real xi(nip)
  real xs1
  real xs2
  real y0
  real yd(ndp)
  real yi(nip)
  real ys1
  real ys2
  real zd(ndp)
  real zi(nip)
!
  save /idlc/
  save /idpt/
!
  common /idlc/ itipv,xs1,xs2,ys1,ys2,ntsc(9)
  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)' ) 'IDBVIP - Fatal error!'
    write ( *, '(a)' ) '  Input parameter MD out of range.'
    stop
  end if
 
  if ( ndp < 4 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IDBVIP - Fatal error!'
    write ( *, '(a)' ) '  Input parameter NDP out of range.'
    stop
  end if
 
  if ( nip < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'IDBVIP - Fatal error!'
    write ( *, '(a)' ) '  Input parameter NIP out of range.'
    stop
  end if
 
  if ( md == 1 ) then
    iwk(1) = ndp
  else
    if ( ndp /= iwk(1) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'IDBVIP - 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) = nip
  else
    if ( nip < iwk(3) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'IDBVIP - Fatal error!'
      write ( *, '(a)' ) '  MD = 3 but NIP was changed since last call.'
      stop
    end if
  end if
!
!  Allocation of storage areas in the IWK array.
!
  jwipt = 16
  jwiwl = 6*ndp+1
  jwiwk = jwiwl
  jwipl = 24*ndp+1
  jwiwp = 30*ndp+1
  jwit0 = 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
!
!  Locate all points at which interpolation is to be performed.
!
  if ( md <= 2 ) then

    itipv = 0
    jwit = jwit0

    do iip = 1, nip

      jwit = jwit+1

      call idlctn ( ndp, xd, yd, nt, iwk(jwipt), nl, iwk(jwipl), &
        xi(iip), yi(iip), iwk(jwit), iwk(jwiwk), wk )

    end do

  end if
!
!  Estimate the partial derivatives at all data points.
!
  call idpdrv ( ndp, xd, yd, zd, nt, iwk(jwipt), wk, wk(jwwpd) )
!
!  Interpolate the ZI values.
!
  itpv = 0
  jwit = jwit0

  do iip = 1, nip

    jwit = jwit + 1

    call idptip ( ndp, xd, yd, zd, nt, iwk(jwipt), nl, iwk(jwipl), wk, &
      iwk(jwit), xi(iip), yi(iip), zi(iip) )

  end do
 
  return
end
subroutine idgrid ( xd, yd, nt, ipt, nl, ipl, nxi, nyi, xi, yi, ngp, igp )
!
!*******************************************************************************
!
!! IDGRID organizes grid points for surface fitting.
!
!
!  Discussion:
!
!    IDGRID sorts the points in ascending order of triangle numbers and 
!    of the border line segment number.
!
!  Parameters:
!
!    Input, real XD(NDP), YD(NDP), the X and Y coordinates of the data 
!    points.
!
!    Input, integer NT, the number of triangles.
!
!    Input, integer IPT(3*NT), the indices of the triangle vertexes.
!
!    Input, integer NL, the number of border line segments.
!
!    Input, integer IPL(3*NL), containing the point numbers of the end points 
!    of the border line segments and their respective triangle numbers,
!
!    Input, integer NXI, NYI, the number of grid points in the X and Y
!    coordinates.
!
!    Input, real XI(NXI), YI(NYI), the coordinates of the grid points.
!
!    Output, integer NGP(2*(NT+2*NL)) where the
!    number of grid points that belong to each of the
!    triangles or of the border line segments are to be stored.
!
!    Output, integer IGP(NXI*NYI), where the grid point numbers are to be 
!    stored in ascending order of the triangle number and the border line
!    segment number.
!
  implicit none
!
  integer nl
  integer nt
  integer nxi
  integer nyi
!
  integer igp(nxi*nyi)
  integer il0
  integer il0t3
  integer ilp1
  integer ilp1t3
  integer insd
  integer ip1
  integer ip2
  integer ip3
  integer ipl(3*nl)
  integer ipt(3*nt)
  integer it0
  integer it0t3
  integer ixi
  integer iximn
  integer iximx
  integer iyi
  integer izi
  integer jigp0
  integer jigp1
  integer jigp1i
  integer jngp0
  integer jngp1
  integer l
  integer ngp(2*(nt+2*nl))
  integer ngp0
  integer ngp1
  integer nl0
  integer nt0
  integer nxinyi
  real spdt
  real u1
  real u2
  real u3
  real v1
  real v2
  real v3
  real vpdt
  real x1
  real x2
  real x3
  real xd(*)
  real xi(nxi)
  real xii
  real ximn
  real ximx
  real xmn
  real xmx
  real y1
  real y2
  real y3
  real yd(*)
  real yi(nyi)
  real yii
  real yimn
  real yimx
  real ymn
  real ymx
!
!  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
  nxinyi = nxi * nyi
  ximn = min ( xi(1), xi(nxi) )
  ximx = max ( xi(1), xi(nxi) )
  yimn = min ( yi(1), yi(nyi) )
  yimx = max ( yi(1), yi(nyi) )
!
!  Determine grid points inside the data area.
!
  jngp0 = 0
  jngp1 = 2*(nt0+2*nl0)+1
  jigp0 = 0
  jigp1 = nxinyi + 1
 
  do it0 = 1, nt0

    ngp0 = 0
    ngp1 = 0
    it0t3 = it0*3
    ip1 = ipt(it0t3-2)
    ip2 = ipt(it0t3-1)
    ip3 = ipt(it0t3)
    x1 = xd(ip1)
    y1 = yd(ip1)
    x2 = xd(ip2)
    y2 = yd(ip2)
    x3 = xd(ip3)
    y3 = yd(ip3)
    xmn = min ( x1, x2, x3 )
    xmx = max ( x1, x2, x3 )
    ymn = min ( y1, y2, y3 )
    ymx = max ( y1, y2, y3 )
    insd = 0

    do ixi = 1, nxi

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

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

    end do
 
    if ( insd == 0 ) then
      go to 38
    end if

    iximx = nxi

23  continue

    do iyi = 1, nyi

      yii = yi(iyi)

      if ( yii < ymn .or. yii > ymx ) then
        go to 37
      end if

      do ixi = iximn, iximx

        xii = xi(ixi)
        l = 0
        if ( vpdt(x1,y1,x2,y2,xii,yii) ) 36,25,26

25      continue

        l = 1
26      continue
        if(vpdt(x2,y2,x3,y3,xii,yii))     36,27,28
27      continue
        l = 1
28      continue
        if(vpdt(x3,y3,x1,y1,xii,yii))     36,29,30
29      continue
        l = 1
30      continue
        izi = nxi*(iyi-1)+ixi

        if ( l == 1 ) go to 31

        ngp0 = ngp0+1
        jigp0 = jigp0+1
        igp(jigp0) = izi
        go to 36

31      continue
 
        do jigp1i = jigp1,nxinyi
          if ( izi == igp(jigp1i) ) then
            go to 36
          end if
        end do
 
        ngp1 = ngp1+1
        jigp1 = jigp1-1
        igp(jigp1) = izi

36      continue

      end do

37    continue

    end do

38  continue

    jngp0 = jngp0+1
    ngp(jngp0) = ngp0
    jngp1 = jngp1-1
    ngp(jngp1) = ngp1

  end do
!
!  Determine grid points outside the data area.
!  in semi-infinite rectangular area.
!
  do il0 = 1, nl0

    ngp0 = 0
    ngp1 = 0
    il0t3 = il0*3
    ip1 = ipl(il0t3-2)
    ip2 = ipl(il0t3-1)
    x1 = xd(ip1)
    y1 = yd(ip1)
    x2 = xd(ip2)
    y2 = yd(ip2)

    xmn = ximn
    xmx = ximx
    ymn = yimn