/****************************************************************************
                  INTERNATIONAL AVS CENTER
	(This disclaimer must remain at the top of all files)

WARRANTY DISCLAIMER

This module and the files associated with it are distributed free of charge.
It is placed in the public domain and permission is granted for anyone to use,
duplicate, modify, and redistribute it unless otherwise noted.  Some modules
may be copyrighted.  You agree to abide by the conditions also included in
the AVS Licensing Agreement, version 1.0, located in the main module
directory located at the International AVS Center ftp site and to include
the AVS Licensing Agreement when you distribute any files downloaded from 
that site.

The International AVS Center, MCNC, the AVS Consortium and the individual
submitting the module and files associated with said module provide absolutely
NO WARRANTY OF ANY KIND with respect to this software.  The entire risk as to
the quality and performance of this software is with the user.  IN NO EVENT
WILL The International AVS Center, MCNC, the AVS Consortium and the individual
submitting the module and files associated with said module BE LIABLE TO
ANYONE FOR ANY DAMAGES ARISING FROM THE USE OF THIS SOFTWARE, INCLUDING,
WITHOUT LIMITATION, DAMAGES RESULTING FROM LOST DATA OR LOST PROFITS, OR ANY
SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES.

This AVS module and associated files are public domain software unless
otherwise noted.  Permission is hereby granted to do whatever you like with
it, subject to the conditions that may exist in copyrighted materials. Should
you wish to make a contribution toward the improvement, modification, or
general performance of this module, please send us your comments:  why you
liked or disliked it, how you use it, and most important, how it helps your
work. We will receive your comments at avs@ncsc.org.

Please send AVS module bug reports to avs@ncsc.org.

******************************************************************************/
C
C NOTE: this software was obtained from netlib@ornl.gov.  All
C relevant disclaimers apply.
C 
C 22 may 1991 w.bethel  lawrence berkeley laboratory
C
C 31 jan 1992 w.bethel
C   fixed some problems relating to internal hard-coded limits which
C   turned out to be unreasonable.
C 
      subroutine  idbvip(md,ncp,ndp,xd,yd,zd,nip,xi,yi,zi,
     1                   iwk,wk)
c this subroutine performs bivariate interpolation when the pro-
c jections of the data points in the x-y plane are irregularly
c distributed in the plane.
c the input parameters are
c     md  = mode of computation (must be 1, 2, or 3),
c         = 1 for new ncp and/or new xd-yd,
c         = 2 for old ncp, old xd-yd, new xi-yi,
c         = 3 for old ncp, old xd-yd, old xi-yi,
c     ncp = number of additional data points used for esti-
c           mating partial derivatives at each data point
c           (must be 2 or greater, but smaller than ndp),
c     ndp = number of data points (must be 4 or greater),
c     xd  = array of dimension ndp containing the x
c           coordinates of the data points,
c     yd  = array of dimension ndp containing the y
c           coordinates of the data points,
c     zd  = array of dimension ndp containing the z
c           coordinates of the data points,
c     nip = number of output points at which interpolation
c           is to be performed (must be 1 or greater),
c     xi  = array of dimension nip containing the x
c           coordinates of the output points,
c     yi  = array of dimension nip containing the y
c           coordinates of the output points.
c the output parameter is
c     zi  = array of dimension nip where interpolated z
c           values are to be stored.
c the other parameters are
c     iwk = integer array of dimension
c              max0(31,27+ncp)*ndp+nip
c           used internally as a work area,
c     wk  = array of dimension 8*ndp used internally as a
c           work area.
c the very first call to this subroutine and the call with a new
c ncp value, a new ndp value, and/or new contents of the xd and
c yd arrays must be made with md=1.  the call with md=2 must be
c preceded by another call with the same ncp and ndp values and
c with the same contents of the xd and yd arrays.  the call with
c md=3 must be preceded by another call with the same ncp, ndp,
c and nip values and with the same contents of the xd, yd, xi,
c and yi arrays.  between the call with md=2 or md=3 and its
c preceding call, the iwk and wk arrays must not be disturbed.
c use of a value between 3 and 5 (inclusive) for ncp is recom-
c mended unless there are evidences that dictate otherwise.
c the lun constant in the data initialization statement is the
c logical unit number of the standard output unit and is,
c therefore, system dependent.
c this subroutine calls the idcldp, idlctn, idpdrv, idptip, and
c idtang subroutines.
c declaration statements
      dimension   xd(100),yd(100),zd(100),xi(1000),yi(1000),
     1            zi(1000),iwk(4100),wk(800)
      common/idlc/nit
      common/idpi/itpv
      data  lun/6/
c setting of some input parameters to local variables.
c (for md=1,2,3)
   10 md0=md
      ncp0=ncp
      ndp0=ndp
      nip0=nip
c error check.  (for md=1,2,3)
   20 if(md0.lt.1.or.md0.gt.3)           go to 90
      if(ncp0.lt.2.or.ncp0.ge.ndp0)      go to 90
      if(ndp0.lt.4)                      go to 90
      if(nip0.lt.1)                      go to 90
      if(md0.ge.2)        go to 21
      iwk(1)=ncp0
      iwk(2)=ndp0
      go to 22
   21 ncppv=iwk(1)
      ndppv=iwk(2)
      if(ncp0.ne.ncppv)   go to 90
      if(ndp0.ne.ndppv)   go to 90
   22 if(md0.ge.3)        go to 23
      iwk(3)=nip
      go to 30
   23 nippv=iwk(3)
      if(nip0.ne.nippv)   go to 90
c allocation of storage areas in the iwk array.  (for md=1,2,3)
   30 jwipt=16
      jwiwl=6*ndp0+1
      jwiwk=jwiwl
      jwipl=24*ndp0+1
      jwiwp=30*ndp0+1
      jwipc=27*ndp0+1
      jwit0=max0(31,27+ncp0)*ndp0
c triangulates the x-y plane.  (for md=1)
   40 if(md0.gt.1)   go to 50
      call idtang(ndp0,xd,yd,nt,iwk(jwipt),nl,iwk(jwipl),
     1            iwk(jwiwl),iwk(jwiwp),wk)
      iwk(5)=nt
      iwk(6)=nl
      if(nt.eq.0)    return
c determines ncp points closest to each data point.  (for md=1)
   50 if(md0.gt.1)   go to 60
      call idcldp(ndp0,xd,yd,ncp0,iwk(jwipc))
      if(iwk(jwipc).eq.0)      return
c locates all points at which interpolation is to be performed.
c (for md=1,2)
   60 if(md0.eq.3)   go to 70
      nit=0
      jwit=jwit0
      do 61  iip=1,nip0
        jwit=jwit+1
        call idlctn(ndp0,xd,yd,nt,iwk(jwipt),nl,iwk(jwipl),
     1            xi(iip),yi(iip),iwk(jwit),iwk(jwiwk),wk)
   61 continue
c estimates partial derivatives at all data points.
c (for md=1,2,3)
   70 call idpdrv(ndp0,xd,yd,zd,ncp0,iwk(jwipc),wk)
c interpolates the zi values.  (for md=1,2,3)
   80 itpv=0
      jwit=jwit0
      do 81  iip=1,nip0
        jwit=jwit+1
        call idptip(xd,yd,zd,nt,iwk(jwipt),nl,iwk(jwipl),wk,
     1              iwk(jwit),xi(iip),yi(iip),zi(iip))
   81 continue
      return
c error exit
   90 write (lun,2090) md0,ncp0,ndp0,nip0
      return
c format statement for error message
 2090 format(1x/41h ***   improper input parameter value(s)./
     1   7h   md =,i4,10x,5hncp =,i6,10x,5hndp =,i6,
     2   10x,5hnip =,i6/
     3   35h error detected in routine   idbvip/)
      end
C
      subroutine  idcldp(ndp,xd,yd,ncp,ipc)
c this subroutine selects several data points that are closest
c to each of the data point.
c the input parameters are
c     ndp = number of data points,
c     xd,yd = arrays of dimension ndp containing the x and y
c           coordinates of the data points,
c     ncp = number of data points closest to each data
c           points.
c the output parameter is
c     ipc = integer array of dimension ncp*ndp, where the
c           point numbers of ncp data points closest to
c           each of the ndp data points are to be stored.
c this subroutine arbitrarily sets a restriction that ncp must
c not exceed 25.
c the lun constant in the data initialization statement is the
c logical unit number of the standard output unit and is,
c therefore, system dependent.
c declaration statements
c      dimension   xd(100),yd(100),ipc(400)
C
C changed 31 jan 92 w.bethel to remove restriction of ncp max == 25
C     the arbitrary max size is now hard coded to 1,000,000
c
c      dimension   dsq0(25),ipc0(25)
c      data  ncpmx/25/, lun/6/
      dimension   xd(1),yd(1),ipc(ncp*ndp)
      dimension   dsq0(1000000),ipc0(100000)
      data  ncpmx/1000000/, lun/6/
c statement function
      dsqf(u1,v1,u2,v2)=(u2-u1)**2+(v2-v1)**2
c preliminary processing
   10 ndp0=ndp
      ncp0=ncp
      if(ndp0.lt.2)  go to 90
      if(ncp0.lt.1.or.ncp0.gt.ncpmx.or.ncp0.ge.ndp0)    go to 90
c calculation
   20 do 59  ip1=1,ndp0
c - selects ncp points.
        x1=xd(ip1)
        y1=yd(ip1)
        j1=0
        dsqmx=0.0
        do 22  ip2=1,ndp0
          if(ip2.eq.ip1)  go to 22
          dsqi=dsqf(x1,y1,xd(ip2),yd(ip2))
          j1=j1+1
          dsq0(j1)=dsqi
          ipc0(j1)=ip2
          if(dsqi.le.dsqmx)    go to 21
          dsqmx=dsqi
          jmx=j1
   21     if(j1.ge.ncp0)  go to 23
   22   continue
   23   ip2mn=ip2+1
        if(ip2mn.gt.ndp0)      go to 30
        do 25  ip2=ip2mn,ndp0
          if(ip2.eq.ip1)  go to 25
          dsqi=dsqf(x1,y1,xd(ip2),yd(ip2))
          if(dsqi.ge.dsqmx)    go to 25
          dsq0(jmx)=dsqi
          ipc0(jmx)=ip2
          dsqmx=0.0
          do 24  j1=1,ncp0
            if(dsq0(j1).le.dsqmx)   go to 24
            dsqmx=dsq0(j1)
            jmx=j1
   24     continue
   25   continue
c - checks if all the ncp+1 points are collinear.
   30   ip2=ipc0(1)
        dx12=xd(ip2)-x1
        dy12=yd(ip2)-y1
        do 31  j3=2,ncp0
          ip3=ipc0(j3)
          dx13=xd(ip3)-x1
          dy13=yd(ip3)-y1
          if((dy13*dx12-dx13*dy12).ne.0.0)    go to 50
   31   continue
c - searches for the closest noncollinear point.
   40   nclpt=0
        do 43  ip3=1,ndp0
          if(ip3.eq.ip1)       go to 43
          do 41  j4=1,ncp0
            if(ip3.eq.ipc0(j4))     go to 43
   41     continue
          dx13=xd(ip3)-x1
          dy13=yd(ip3)-y1
          if((dy13*dx12-dx13*dy12).eq.0.0)    go to 43
          dsqi=dsqf(x1,y1,xd(ip3),yd(ip3))
          if(nclpt.eq.0)       go to 42
          if(dsqi.ge.dsqmn)    go to 43
   42     nclpt=1
          dsqmn=dsqi
          ip3mn=ip3
   43   continue
        if(nclpt.eq.0)    go to 91
        dsqmx=dsqmn
        ipc0(jmx)=ip3mn
c - replaces the local array for the output array.
   50   j1=(ip1-1)*ncp0
        do 51  j2=1,ncp0
          j1=j1+1
          ipc(j1)=ipc0(j2)
   51   continue
   59 continue
      return
c error exit
   90 write (lun,2090)
      go to 92
   91 write (lun,2091)
   92 write (lun,2092)  ndp0,ncp0
      ipc(1)=0
      return
c format statements for error messages
 2090 format(1x/41h ***   improper input parameter value(s).)
 2091 format(1x/33h ***   all collinear data points.)
 2092 format(8h   ndp =,i5,5x,5hncp =,i5/
     1   35h error detected in routine   idcldp/)
      end
      subroutine idgrid(xd, yd, nt, ipt, nl, ipl, nxi, nyi, xi, yi,
     *  ngp, igp)
c this subroutine organizes grid points for surface fitting by
c sorting them in ascending order of triangle numbers and of the
c border line segment number.
c the input parameters are
c     xd,yd = arrays of dimension ndp containing the x and y
c           coordinates of the data points, where ndp is the
c           number of the data points,
c     nt  = number of triangles,
c     ipt = integer array of dimension 3*nt containing the
c           point numbers of the vertexes of the triangles,
c     nl  = number of border line segments,
c     ipl = integer array of dimension 3*nl containing the
c           point numbers of the end points of the border
c           line segments and their respective triangle
c           numbers,
c     nxi = number of grid points in the x coordinate,
c     nyi = number of grid points in the y coordinate,
c     xi,yi = arrays of dimension nxi and nyi containing
c           the x and y coordinates of the grid points,
c           respectively.
c the output parameters are
c     ngp = integer array of dimension 2*(nt+2*nl) where the
c           number of grid points that belong to each of the
c           triangles or of the border line segments are to
c           be stored,
c     igp = integer array of dimension nxi*nyi where the
c           grid point numbers are to be stored in ascending
c           order of the triangle number and the border line
c           segment number.
c declaration statements
      dimension xd(100), yd(100), ipt(585), ipl(300), xi(101),
     *  yi(101), ngp(800), igp(10201)
c statement functions
      side(u1,v1,u2,v2,u3,v3) = (u1-u3)*(v2-v3) - (v1-v3)*(u2-u3)
      spdt(u1,v1,u2,v2,u3,v3) = (u1-u2)*(u3-u2) + (v1-v2)*(v3-v2)
c preliminary processing
      nt0 = nt
      nl0 = nl
      nxi0 = nxi
      nyi0 = nyi
      nxinyi = nxi0*nyi0
      ximn = amin1(xi(1),xi(nxi0))
      ximx = amax1(xi(1),xi(nxi0))
      yimn = amin1(yi(1),yi(nyi0))
      yimx = amax1(yi(1),yi(nyi0))
c determines grid points inside the data area.
      jngp0 = 0
      jngp1 = 2*(nt0+2*nl0) + 1
      jigp0 = 0
      jigp1 = nxinyi + 1
      do 160 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 = amin1(x1,x2,x3)
        xmx = amax1(x1,x2,x3)
        ymn = amin1(y1,y2,y3)
        ymx = amax1(y1,y2,y3)
        insd = 0
        do 20 ixi=1,nxi0
          if (xi(ixi).ge.xmn .and. xi(ixi).le.xmx) go to 10
          if (insd.eq.0) go to 20
          iximx = ixi - 1
          go to 30
   10     if (insd.eq.1) go to 20
          insd = 1
          iximn = ixi
   20   continue
        if (insd.eq.0) go to 150
        iximx = nxi0
   30   do 140 iyi=1,nyi0
          yii = yi(iyi)
          if (yii.lt.ymn .or. yii.gt.ymx) go to 140
          do 130 ixi=iximn,iximx
            xii = xi(ixi)
            l = 0
            if (side(x1,y1,x2,y2,xii,yii)) 130, 40, 50
   40       l = 1
   50       if (side(x2,y2,x3,y3,xii,yii)) 130, 60, 70
   60       l = 1
   70       if (side(x3,y3,x1,y1,xii,yii)) 130, 80, 90
   80       l = 1
   90       izi = nxi0*(iyi-1) + ixi
            if (l.eq.1) go to 100
            ngp0 = ngp0 + 1
            jigp0 = jigp0 + 1
            igp(jigp0) = izi
            go to 130
  100       if (jigp1.gt.nxinyi) go to 120
            do 110 jigp1i=jigp1,nxinyi
              if (izi.eq.igp(jigp1i)) go to 130
  110       continue
  120       ngp1 = ngp1 + 1
            jigp1 = jigp1 - 1
            igp(jigp1) = izi
  130     continue
  140   continue
  150   jngp0 = jngp0 + 1
        ngp(jngp0) = ngp0
        jngp1 = jngp1 - 1
        ngp(jngp1) = ngp1
  160 continue
c determines grid points outside the data area.
c - in semi-infinite rectangular area.
      do 450 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
        ymx = yimx
        if (y2.ge.y1) xmn = amin1(x1,x2)
        if (y2.le.y1) xmx = amax1(x1,x2)
        if (x2.le.x1) ymn = amin1(y1,y2)
        if (x2.ge.x1) ymx = amax1(y1,y2)
        insd = 0
        do 180 ixi=1,nxi0
          if (xi(ixi).ge.xmn .and. xi(ixi).le.xmx) go to 170
          if (insd.eq.0) go to 180
          iximx = ixi - 1
          go to 190
  170     if (insd.eq.1) go to 180
          insd = 1
          iximn = ixi
  180   continue
        if (insd.eq.0) go to 310
        iximx = nxi0
  190   do 300 iyi=1,nyi0
          yii = yi(iyi)
          if (yii.lt.ymn .or. yii.gt.ymx) go to 300
          do 290 ixi=iximn,iximx
            xii = xi(ixi)
            l = 0
            if (side(x1,y1,x2,y2,xii,yii)) 210, 200, 290
  200       l = 1
  210       if (spdt(x2,y2,x1,y1,xii,yii)) 290, 220, 230
  220       l = 1
  230       if (spdt(x1,y1,x2,y2,xii,yii)) 290, 240, 250
  240       l = 1
  250       izi = nxi0*(iyi-1) + ixi
            if (l.eq.1) go to 260
            ngp0 = ngp0 + 1
            jigp0 = jigp0 + 1
            igp(jigp0) = izi
            go to 290
  260       if (jigp1.gt.nxinyi) go to 280
            do 270 jigp1i=jigp1,nxinyi
              if (izi.eq.igp(jigp1i)) go to 290
  270       continue
  280       ngp1 = ngp1 + 1
            jigp1 = jigp1 - 1
            igp(jigp1) = izi
  290     continue
  300   continue
  310   jngp0 = jngp0 + 1
        ngp(jngp0) = ngp0
        jngp1 = jngp1 - 1
        ngp(jngp1) = ngp1
c - in semi-infinite triangular area.
        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.ge.y2 .and. y2.ge.y1) xmn = x2
        if (y3.le.y2 .and. y2.le.y1) xmx = x2
        if (x3.le.x2 .and. x2.le.x1) ymn = y2
        if (x3.ge.x2 .and. x2.ge.x1) ymx = y2
        insd = 0
        do 330 ixi=1,nxi0
          if (xi(ixi).ge.xmn .and. xi(ixi).le.xmx) go to 320
          if (insd.eq.0) go to 330
          iximx = ixi - 1
          go to 340
  320     if (insd.eq.1) go to 330
          insd = 1
          iximn = ixi
  330   continue
        if (insd.eq.0) go to 440
        iximx = nxi0
  340   do 430 iyi=1,nyi0
          yii = yi(iyi)
          if (yii.lt.ymn .or. yii.gt.ymx) go to 430
          do 420 ixi=iximn,iximx
            xii = xi(ixi)
            l = 0
            if (spdt(x1,y1,x2,y2,xii,yii)) 360, 350, 420
  350       l = 1
  360       if (spdt(x3,y3,x2,y2,xii,yii)) 380, 370, 420
  370       l = 1
  380       izi = nxi0*(iyi-1) + ixi
            if (l.eq.1) go to 390
            ngp0 = ngp0 + 1
            jigp0 = jigp0 + 1
            igp(jigp0) = izi
            go to 420
  390       if (jigp1.gt.nxinyi) go to 410
            do 400 jigp1i=jigp1,nxinyi
              if (izi.eq.igp(jigp1i)) go to 420
  400       continue
  410       ngp1 = ngp1 + 1
            jigp1 = jigp1 - 1
            igp(jigp1) = izi
  420     continue
  430   continue
  440   jngp0 = jngp0 + 1
        ngp(jngp0) = ngp0
        jngp1 = jngp1 - 1
        ngp(jngp1) = ngp1
  450 continue
      return
      end
      subroutine idlctn(ndp, xd, yd, nt, ipt, nl, ipl, xii, yii, iti,
     *  iwk, wk)
c this subroutine locates a point, i.e., determines to what tri-
c angle a given point (xii,yii) belongs.  when the given point
c does not lie inside the data area, this subroutine determines
c the border line segment when the point lies in an outside
c rectangular area, and two border line segments when the point
c lies in an outside triangular area.
c the input parameters are
c     ndp = number of data points,
c     xd,yd = arrays of dimension ndp containing the x and y
c           coordinates of the data points,
c     nt  = number of triangles,
c     ipt = integer array of dimension 3*nt containing the
c           point numbers of the vertexes of the triangles,
c     nl  = number of border line segments,
c     ipl = integer array of dimension 3*nl containing the
c           point numbers of the end points of the border
c           line segments and their respective triangle
c           numbers,
c     xii,yii = x and y coordinates of the point to be
c           located.
c the output parameter is
c     iti = triangle number, when the point is inside the
c           data area, or
c           two border line segment numbers, il1 and il2,
c           coded to il1*(nt+nl)+il2, when the point is
c           outside the data area.
c the other parameters are
c     iwk = integer array of dimension 18*ndp used inter-
c           nally as a work area,
c     wk  = array of dimension 8*ndp used internally as a
c           work area.
c declaration statements
      dimension xd(100), yd(100), ipt(585), ipl(300), iwk(1800),
     *  wk(800)
      dimension ntsc(9), idsc(9)
      common /idlc/ nit
c statement functions
      side(u1,v1,u2,v2,u3,v3) = (u1-u3)*(v2-v3) - (v1-v3)*(u2-u3)
      spdt(u1,v1,u2,v2,u3,v3) = (u1-u2)*(u3-u2) + (v1-v2)*(v3-v2)
c preliminary processing
      ndp0 = ndp
      nt0 = nt
      nl0 = nl
      ntl = nt0 + nl0
      x0 = xii
      y0 = yii
c processing for a new set of data points
      if (nit.ne.0) go to 80
      nit = 1
c - divides the x-y plane into nine rectangular sections.
      xmn = xd(1)
      xmx = xmn
      ymn = yd(1)
      ymx = ymn
      do 10 idp=2,ndp0
        xi = xd(idp)
        yi = yd(idp)
        xmn = amin1(xi,xmn)
        xmx = amax1(xi,xmx)
        ymn = amin1(yi,ymn)
        ymx = amax1(yi,ymx)
   10 continue
      xs1 = (xmn+xmn+xmx)/3.0
      xs2 = (xmn+xmx+xmx)/3.0
      ys1 = (ymn+ymn+ymx)/3.0
      ys2 = (ymn+ymx+ymx)/3.0
c - determines and stores in the iwk array triangle numbers of
c - the triangles associated with each of the nine sections.
      do 20 isc=1,9
        ntsc(isc) = 0
        idsc(isc) = 0
   20 continue
      it0t3 = 0
      jwk = 0
      do 70 it0=1,nt0
        it0t3 = it0t3 + 3
        i1 = ipt(it0t3-2)
        i2 = ipt(it0t3-1)
        i3 = ipt(it0t3)
        xmn = amin1(xd(i1),xd(i2),xd(i3))
        xmx = amax1(xd(i1),xd(i2),xd(i3))
        ymn = amin1(yd(i1),yd(i2),yd(i3))
        ymx = amax1(yd(i1),yd(i2),yd(i3))
        if (ymn.gt.ys1) go to 30
        if (xmn.le.xs1) idsc(1) = 1
        if (xmx.ge.xs1 .and. xmn.le.xs2) idsc(2) = 1
        if (xmx.ge.xs2) idsc(3) = 1
   30   if (ymx.lt.ys1 .or. ymn.gt.ys2) go to 40
        if (xmn.le.xs1) idsc(4) = 1
        if (xmx.ge.xs1 .and. xmn.le.xs2) idsc(5) = 1
        if (xmx.ge.xs2) idsc(6) = 1
   40   if (ymx.lt.ys2) go to 50
        if (xmn.le.xs1) idsc(7) = 1
        if (xmx.ge.xs1 .and. xmn.le.xs2) idsc(8) = 1
        if (xmx.ge.xs2) idsc(9) = 1
   50   do 60 isc=1,9
          if (idsc(isc).eq.0) go to 60
          jiwk = 9*ntsc(isc) + isc
          iwk(jiwk) = it0
          ntsc(isc) = ntsc(isc) + 1
          idsc(isc) = 0
   60   continue
c - stores in the wk array the minimum and maximum of the x and
c - 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
   70 continue
      go to 110
c checks if in the same triangle as previous.
   80 it0 = itipv
      if (it0.gt.nt0) go to 90
      it0t3 = it0*3
      ip1 = ipt(it0t3-2)
      x1 = xd(ip1)
      y1 = yd(ip1)
      ip2 = ipt(it0t3-1)
      x2 = xd(ip2)
      y2 = yd(ip2)
      if (side(x1,y1,x2,y2,x0,y0).lt.0.0) go to 110
      ip3 = ipt(it0t3)
      x3 = xd(ip3)
      y3 = yd(ip3)
      if (side(x2,y2,x3,y3,x0,y0).lt.0.0) go to 110
      if (side(x3,y3,x1,y1,x0,y0).lt.0.0) go to 110
      go to 170
c checks if on the same border line segment.
   90 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.ne.il1) go to 100
      if (spdt(x1,y1,x2,y2,x0,y0).lt.0.0) go to 110
      if (spdt(x2,y2,x1,y1,x0,y0).lt.0.0) go to 110
      if (side(x1,y1,x2,y2,x0,y0).gt.0.0) go to 110
      go to 170
c checks if between the same two border line segments.
  100 if (spdt(x1,y1,x2,y2,x0,y0).gt.0.0) go to 110
      ip3 = ipl(3*il2-1)
      x3 = xd(ip3)
      y3 = yd(ip3)
      if (spdt(x3,y3,x2,y2,x0,y0).le.0.0) go to 170
c locates inside the data area.
c - determines the section in which the point in question lies.
  110 isc = 1
      if (x0.ge.xs1) isc = isc + 1
      if (x0.ge.xs2) isc = isc + 1
      if (y0.ge.ys1) isc = isc + 3
      if (y0.ge.ys2) isc = isc + 3
c - searches through the triangles associated with the section.
      ntsci = ntsc(isc)
      if (ntsci.le.0) go to 130
      jiwk = -9 + isc
      do 120 itsc=1,ntsci
        jiwk = jiwk + 9
        it0 = iwk(jiwk)
        jwk = it0*4
        if (x0.lt.wk(jwk-3)) go to 120
        if (x0.gt.wk(jwk-2)) go to 120
        if (y0.lt.wk(jwk-1)) go to 120
        if (y0.gt.wk(jwk)) go to 120
        it0t3 = it0*3
        ip1 = ipt(it0t3-2)
        x1 = xd(ip1)
        y1 = yd(ip1)
        ip2 = ipt(it0t3-1)
        x2 = xd(ip2)
        y2 = yd(ip2)
        if (side(x1,y1,x2,y2,x0,y0).lt.0.0) go to 120
        ip3 = ipt(it0t3)
        x3 = xd(ip3)
        y3 = yd(ip3)
        if (side(x2,y2,x3,y3,x0,y0).lt.0.0) go to 120
        if (side(x3,y3,x1,y1,x0,y0).lt.0.0) go to 120
        go to 170
  120 continue
c locates outside the data area.
  130 do 150 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).lt.0.0) go to 150
        if (spdt(x1,y1,x2,y2,x0,y0).lt.0.0) go to 140
        if (side(x1,y1,x2,y2,x0,y0).gt.0.0) go to 150
        il2 = il1
        go to 160
  140   il2 = mod(il1,nl0) + 1
        ip3 = ipl(3*il2-1)
        x3 = xd(ip3)
        y3 = yd(ip3)
        if (spdt(x3,y3,x2,y2,x0,y0).le.0.0) go to 160
  150 continue
      it0 = 1
      go to 170
  160 it0 = il1*ntl + il2
c normal exit
  170 iti = it0
      itipv = it0
      return
      end
      subroutine  idpdrv(ndp,xd,yd,zd,ncp,ipc,pd)
c this subroutine estimates partial derivatives of the first and
c second order at the data points.
c the input parameters are
c     ndp = number of data points,
c     xd,yd,zd = arrays of dimension ndp containing the x,
c           y, and z coordinates of the data points,
c     ncp = number of additional data points used for esti-
c           mating partial derivatives at each data point,
c     ipc = integer array of dimension ncp*ndp containing
c           the point numbers of ncp data points closest to
c           each of the ndp data points.
c the output parameter is
c     pd  = array of dimension 5*ndp, where the estimated
c           zx, zy, zxx, zxy, and zyy values at the data
c           points are to be stored.
c declaration statements
      dimension   xd(100),yd(100),zd(100),ipc(400),pd(500)
      real        nmx,nmy,nmz,nmxx,nmxy,nmyx,nmyy
c preliminary processing
   10 ndp0=ndp
      ncp0=ncp
      ncpm1=ncp0-1
c estimation of zx and zy
  20  do 24  ip0=1,ndp0
        x0=xd(ip0)
        y0=yd(ip0)
        z0=zd(ip0)
        nmx=0.0
        nmy=0.0
        nmz=0.0
        jipc0=ncp0*(ip0-1)
        do 23  ic1=1,ncpm1
          jipc=jipc0+ic1
          ipi=ipc(jipc)
          dx1=xd(ipi)-x0
          dy1=yd(ipi)-y0
          dz1=zd(ipi)-z0
          ic2mn=ic1+1
          do 22  ic2=ic2mn,ncp0
            jipc=jipc0+ic2
            ipi=ipc(jipc)
            dx2=xd(ipi)-x0
            dy2=yd(ipi)-y0
            dnmz=dx1*dy2-dy1*dx2
            if(dnmz.eq.0.0)    go to 22
            dz2=zd(ipi)-z0
            dnmx=dy1*dz2-dz1*dy2
            dnmy=dz1*dx2-dx1*dz2
            if(dnmz.ge.0.0)    go to 21
            dnmx=-dnmx
            dnmy=-dnmy
            dnmz=-dnmz
   21       nmx=nmx+dnmx
            nmy=nmy+dnmy
            nmz=nmz+dnmz
   22     continue
   23   continue
        jpd0=5*ip0
        pd(jpd0-4)=-nmx/nmz
        pd(jpd0-3)=-nmy/nmz
   24 continue
c estimation of zxx, zxy, and zyy
   30 do 34  ip0=1,ndp0
        jpd0=jpd0+5
        x0=xd(ip0)
        jpd0=5*ip0
        y0=yd(ip0)
        zx0=pd(jpd0-4)
        zy0=pd(jpd0-3)
        nmxx=0.0
        nmxy=0.0
        nmyx=0.0
        nmyy=0.0
        nmz =0.0
        jipc0=ncp0*(ip0-1)
        do 33  ic1=1,ncpm1
          jipc=jipc0+ic1
          ipi=ipc(jipc)
          dx1=xd(ipi)-x0
          dy1=yd(ipi)-y0
          jpd=5*ipi
          dzx1=pd(jpd-4)-zx0
          dzy1=pd(jpd-3)-zy0
          ic2mn=ic1+1
          do 32  ic2=ic2mn,ncp0
            jipc=jipc0+ic2
            ipi=ipc(jipc)
            dx2=xd(ipi)-x0
            dy2=yd(ipi)-y0
            dnmz =dx1*dy2 -dy1*dx2
            if(dnmz.eq.0.0)    go to 32
            jpd=5*ipi
            dzx2=pd(jpd-4)-zx0
            dzy2=pd(jpd-3)-zy0
            dnmxx=dy1*dzx2-dzx1*dy2
            dnmxy=dzx1*dx2-dx1*dzx2
            dnmyx=dy1*dzy2-dzy1*dy2
            dnmyy=dzy1*dx2-dx1*dzy2
            if(dnmz.ge.0.0)    go to 31
            dnmxx=-dnmxx
            dnmxy=-dnmxy
            dnmyx=-dnmyx
            dnmyy=-dnmyy
            dnmz =-dnmz
   31       nmxx=nmxx+dnmxx
            nmxy=nmxy+dnmxy
            nmyx=nmyx+dnmyx
            nmyy=nmyy+dnmyy
            nmz =nmz +dnmz
   32     continue
   33   continue
        pd(jpd0-2)=-nmxx/nmz
        pd(jpd0-1)=-(nmxy+nmyx)/(2.0*nmz)
        pd(jpd0)  =-nmyy/nmz
   34 continue
      return
      end
      subroutine  idptip(xd,yd,zd,nt,ipt,nl,ipl,pdd,iti,xii,yii,
     1                   zii)
c this subroutine performs punctual interpolation or extrapola-
c tion, i.e., determines the z value at a point.
c the input parameters are
c     xd,yd,zd = arrays of dimension ndp containing the x,
c           y, and z coordinates of the data points, where
c           ndp is the number of the data points,
c     nt  = number of triangles,
c     ipt = integer array of dimension 3*nt containing the
c           point numbers of the vertexes of the triangles,
c     nl  = number of border line segments,
c     ipl = integer array of dimension 3*nl containing the
c           point numbers of the end points of the border
c           line segments and their respective triangle
c           numbers,
c     pdd = array of dimension 5*ndp containing the partial
c           derivatives at the data points,
c     iti = triangle number of the triangle in which lies
c           the point for which interpolation is to be
c           performed,
c     xii,yii = x and y coordinates of the point for which
c           interpolation is to be performed.
c the output parameter is
c     zii = interpolated z value.
c declaration statements
      dimension   xd(100),yd(100),zd(100),ipt(585),ipl(300),
     1            pdd(500)
      common/idpi/itpv
      dimension   x(3),y(3),z(3),pd(15),
     1            zu(3),zv(3),zuu(3),zuv(3),zvv(3)
      real        lu,lv
      equivalence (p5,p50)
c preliminary processing
   10 it0=iti
      ntl=nt+nl
      if(it0.le.ntl)      go to 20
      il1=it0/ntl
      il2=it0-il1*ntl
      if(il1.eq.il2)      go to 40
      go to 60
c calculation of zii by interpolation.
c checks if the necessary coefficients have been calculated.
   20 if(it0.eq.itpv)     go to 30
c loads coordinate and partial derivative values at the
c vertexes.
   21 jipt=3*(it0-1)
      jpd=0
      do 23  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 22  kpd=1,5
          jpd=jpd+1
          jpdd=jpdd+1
          pd(jpd)=pdd(jpdd)
   22   continue
   23 continue
c determines the coefficients for the coordinate system
c transformation from the x-y system to the u-v system
c and vice versa.
   24 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
c converts the partial derivatives at the vertexes of the
c triangle for the u-v coordinate system.
   25 aa=a*a
      act2=2.0*a*c
      cc=c*c
      ab=a*b
      adbc=ad+bc
      cd=c*d
      bb=b*b
      bdt2=2.0*b*d
      dd=d*d
      do 26  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)
   26 continue
c calculates the coefficients of the polynomial.
   27 p00=z(1)
      p10=zu(1)
      p01=zv(1)
      p20=0.5*zuu(1)
      p11=zuv(1)
      p02=0.5*zvv(1)
      h1=z(2)-p00-p10-p20
      h2=zu(2)-p10-zuu(1)
      h3=zuu(2)-zuu(1)
      p30= 10.0*h1-4.0*h2+0.5*h3
      p40=-15.0*h1+7.0*h2    -h3
      p50=  6.0*h1-3.0*h2+0.5*h3
      h1=z(3)-p00-p01-p02
      h2=zv(3)-p01-zvv(1)
      h3=zvv(3)-zvv(1)
      p03= 10.0*h1-4.0*h2+0.5*h3
      p04=-15.0*h1+7.0*h2    -h3
      p05=  6.0*h1-3.0*h2+0.5*h3
      lu=sqrt(aa+cc)
      lv=sqrt(bb+dd)
      thxu=atan2(c,a)
      thuv=atan2(d,b)-thxu
      csuv=cos(thuv)
      p41=5.0*lv*csuv/lu*p50
      p14=5.0*lu*csuv/lv*p05
      h1=zv(2)-p01-p11-p41
      h2=zuv(2)-p11-4.0*p41
      p21= 3.0*h1-h2
      p31=-2.0*h1+h2
      h1=zu(3)-p10-p11-p14
      h2=zuv(3)-p11-4.0*p14
      p12= 3.0*h1-h2
      p13=-2.0*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.0*bc+2.0*ad)
      g2=cc*ac*(3.0*ad+2.0*bc)
      h1=-aa*aa*aa*(5.0*aa*bb*p50+(4.0*bc+ad)*p41)
     1   -cc*cc*cc*(5.0*cc*dd*p05+(4.0*ad+bc)*p14)
      h2=0.5*zvv(2)-p02-p12
      h3=0.5*zuu(3)-p20-p21
      p22=(g1*h2+g2*h3-h1)/(g1+g2)
      p32=h2-p22
      p23=h3-p22
      itpv=it0
c converts xii and yii to u-v system.
   30 dx=xii-x0
      dy=yii-y0
      u=ap*dx+bp*dy
      v=cp*dx+dp*dy
c evaluates the polynomial.
   31 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
      zii=p0+u*(p1+u*(p2+u*(p3+u*(p4+u*p5))))
      return
c calculation of zii by extrapolation in the rectangle.
c checks if the necessary coefficients have been calculated.
   40 if(it0.eq.itpv)     go to 50
c loads coordinate and partial derivative values at the end
c points of the border line segment.
   41 jipl=3*(il1-1)
      jpd=0
      do 43  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 42  kpd=1,5
          jpd=jpd+1
          jpdd=jpdd+1
          pd(jpd)=pdd(jpdd)
   42   continue
   43 continue
c determines the coefficients for the coordinate system
c transformation from the x-y system to the u-v system
c and vice versa.
   44 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
c converts the partial derivatives at the end points of the
c border line segment for the u-v coordinate system.
   45 aa=a*a
      act2=2.0*a*c
      cc=c*c
      ab=a*b
      adbc=ad+bc
      cd=c*d
      bb=b*b
      bdt2=2.0*b*d
      dd=d*d
      do 46  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)
   46 continue
c calculates the coefficients of the polynomial.
   47 p00=z(1)
      p10=zu(1)
      p01=zv(1)
      p20=0.5*zuu(1)
      p11=zuv(1)
      p02=0.5*zvv(1)
      h1=z(2)-p00-p01-p02
      h2=zv(2)-p01-zvv(1)
      h3=zvv(2)-zvv(1)
      p03= 10.0*h1-4.0*h2+0.5*h3
      p04=-15.0*h1+7.0*h2    -h3
      p05=  6.0*h1-3.0*h2+0.5*h3
      h1=zu(2)-p10-p11
      h2=zuv(2)-p11
      p12= 3.0*h1-h2
      p13=-2.0*h1+h2
      p21=0.0
      p23=-zuu(2)+zuu(1)
      p22=-1.5*p23
      itpv=it0
c converts xii and yii to u-v system.
   50 dx=xii-x0
      dy=yii-y0
      u=ap*dx+bp*dy
      v=cp*dx+dp*dy
c evaluates the polynomial.
   51 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
c calculation of zii by extrapolation in the triangle.
c checks if the necessary coefficients have been calculated.
   60 if(it0.eq.itpv)     go to 70
c loads coordinate and partial derivative values at the vertex
c of the triangle.
   61 jipl=3*il2-2
      idp=ipl(jipl)
      x(1)=xd(idp)
      y(1)=yd(idp)
      z(1)=zd(idp)
      jpdd=5*(idp-1)
      do 62  kpd=1,5
        jpdd=jpdd+1
        pd(kpd)=pdd(jpdd)
   62 continue
c calculates the coefficients of the polynomial.
   67 p00=z(1)
      p10=pd(1)
      p01=pd(2)
      p20=0.5*pd(3)
      p11=pd(4)
      p02=0.5*pd(5)
      itpv=it0
c converts xii and yii to u-v system.
   70 u=xii-x(1)
      v=yii-y(1)
c evaluates the polynomial.
   71 p0=p00+v*(p01+v*p02)
      p1=p10+v*p11
      zii=p0+u*(p1+u*p20)
      return
      end
      subroutine  idsfft(md,ncp,ndp,xd,yd,zd,nxi,nyi,xi,yi,zi,
     1                   iwk,wk)
c this subroutine performs smooth surface fitting when the pro-
c jections of the data points in the x-y plane are irregularly
c distributed in the plane.
c the input parameters are
c     md  = mode of computation (must be 1, 2, or 3),
c         = 1 for new ncp and/or new xd-yd,
c         = 2 for old ncp, old xd-yd, new xi-yi,
c         = 3 for old ncp, old xd-yd, old xi-yi,
c     ncp = number of additional data points used for esti-
c           mating partial derivatives at each data point
c           (must be 2 or greater, but smaller than ndp),
c     ndp = number of data points (must be 4 or greater),
c     xd  = array of dimension ndp containing the x
c           coordinates of the data points,
c     yd  = array of dimension ndp containing the y
c           coordinates of the data points,
c     zd  = array of dimension ndp containing the z
c           coordinates of the data points,
c     nxi = number of output grid points in the x coordinate
c           (must be 1 or greater),
c     nyi = number of output grid points in the y coordinate
c           (must be 1 or greater),
c     xi  = array of dimension nxi containing the x
c           coordinates of the output grid points,
c     yi  = array of dimension nyi containing the y
c           coordinates of the output grid points.
c the output parameter is
c     zi  = doubly-dimensioned array of dimension (nxi,nyi),
c           where the interpolated z values at the output
c           grid points are to be stored.
c the other parameters are
c     iwk = integer array of dimension
c              max0(31,27+ncp)*ndp+nxi*nyi
c           used internally as a work area,
c     wk  = array of dimension 5*ndp used internally as a
c           work area.
c the very first call to this subroutine and the call with a new
c ncp value, a new ndp value, and/or new contents of the xd and
c yd arrays must be made with md=1.  the call with md=2 must be
c preceded by another call with the same ncp and ndp values and
c with the same contents of the xd and yd arrays.  the call with
c md=3 must be preceded by another call with the same ncp, ndp,
c nxi, and nyi values and with the same contents of the xd, yd,
c xi, and yi arrays.  between the call with md=2 or md=3 and its
c preceding call, the iwk and wk arrays must not be disturbed.
c use of a value between 3 and 5 (inclusive) for ncp is recom-
c mended unless there are evidences that dictate otherwise.
c the lun constant in the data initialization statement is the
c logical unit number of the standard output unit and is,
c therefore, system dependent.
c this subroutine calls the idcldp, idgrid, idpdrv, idptip, and
c idtang subroutines.
c declaration statements
c      dimension   xd(100),yd(100),zd(100),xi(101),yi(101),
c     1            zi(10201),iwk(13301),wk(500)
      dimension   xd(1),yd(1),zd(1),xi(1),yi(1),
     1            zi(1),iwk(1),wk(1)
      common/idpi/itpv
      data  lun/6/
c setting of some input parameters to local variables.
c (for md=1,2,3)
   10 md0=md
      ncp0=ncp
      ndp0=ndp
      nxi0=nxi
      nyi0=nyi
c error check.  (for md=1,2,3)
   20 if(md0.lt.1.or.md0.gt.3)           go to 90
      if(ncp0.lt.2.or.ncp0.ge.ndp0)      go to 90
      if(ndp0.lt.4)                      go to 90
      if(nxi0.lt.1.or.nyi0.lt.1)         go to 90
      if(md0.ge.2)        go to 21
      iwk(1)=ncp0
      iwk(2)=ndp0
      go to 22
   21 ncppv=iwk(1)
      ndppv=iwk(2)
      if(ncp0.ne.ncppv)   go to 90
      if(ndp0.ne.ndppv)   go to 90
   22 if(md0.ge.3)        go to 23
      iwk(3)=nxi0
      iwk(4)=nyi0
      go to 30
   23 nxipv=iwk(3)
      nyipv=iwk(4)
      if(nxi0.ne.nxipv)   go to 90
      if(nyi0.ne.nyipv)   go to 90
c allocation of storage areas in the iwk array.  (for md=1,2,3)
   30 jwipt=16
      jwiwl=6*ndp0+1
      jwngp0=jwiwl-1
      jwipl=24*ndp0+1
      jwiwp=30*ndp0+1
      jwipc=27*ndp0+1
      jwigp0=max0(31,27+ncp0)*ndp0
c triangulates the x-y plane.  (for md=1)
   40 if(md0.gt.1)   go to 50
      call idtang(ndp0,xd,yd,nt,iwk(jwipt),nl,iwk(jwipl),
     1            iwk(jwiwl),iwk(jwiwp),wk)
      iwk(5)=nt
      iwk(6)=nl
      if(nt.eq.0)    return
c determines ncp points closest to each data point.  (for md=1)
   50 if(md0.gt.1)   go to 60
      call idcldp(ndp0,xd,yd,ncp0,iwk(jwipc))
      if(iwk(jwipc).eq.0)      return
c sorts output grid points in ascending order of the triangle
c number and the border line segment number.  (for md=1,2)
   60 if(md0.eq.3)   go to 70
      call idgrid(xd,yd,nt,iwk(jwipt),nl,iwk(jwipl),nxi0,nyi0,
     1            xi,yi,iwk(jwngp0+1),iwk(jwigp0+1))
c estimates partial derivatives at all data points.
c (for md=1,2,3)
   70 call idpdrv(ndp0,xd,yd,zd,ncp0,iwk(jwipc),wk)
c interpolates the zi values.  (for md=1,2,3)
   80 itpv=0
      jig0mx=0
      jig1mn=nxi0*nyi0+1
      nngp=nt+2*nl
      do 89  jngp=1,nngp
        iti=jngp
        if(jngp.le.nt)    go to 81
        il1=(jngp-nt+1)/2
        il2=(jngp-nt+2)/2
        if(il2.gt.nl)     il2=1
        iti=il1*(nt+nl)+il2
   81   jwngp=jwngp0+jngp
        ngp0=iwk(jwngp)
        if(ngp0.eq.0)     go to 86
        jig0mn=jig0mx+1
        jig0mx=jig0mx+ngp0
        do 82  jigp=jig0mn,jig0mx
          jwigp=jwigp0+jigp
          izi=iwk(jwigp)
          iyi=(izi-1)/nxi0+1
          ixi=izi-nxi0*(iyi-1)
          call idptip(xd,yd,zd,nt,iwk(jwipt),nl,iwk(jwipl),wk,
     1                iti,xi(ixi),yi(iyi),zi(izi))
   82   continue
   86   jwngp=jwngp0+2*nngp+1-jngp
        ngp1=iwk(jwngp)
        if(ngp1.eq.0)     go to 89
        jig1mx=jig1mn-1
        jig1mn=jig1mn-ngp1
        do 87  jigp=jig1mn,jig1mx
          jwigp=jwigp0+jigp
          izi=iwk(jwigp)
          iyi=(izi-1)/nxi0+1
          ixi=izi-nxi0*(iyi-1)
          call idptip(xd,yd,zd,nt,iwk(jwipt),nl,iwk(jwipl),wk,
     1                iti,xi(ixi),yi(iyi),zi(izi))
   87   continue
   89 continue
      return
c error exit
   90 write (lun,2090) md0,ncp0,ndp0,nxi0,nyi0
      return
c format statement for error message
 2090 format(1x/41h ***   improper input parameter value(s)./
     1   7h   md =,i4,10x,5hncp =,i6,10x,5hndp =,i6,
     2   10x,5hnxi =,i6,10x,5hnyi =,i6/
     3   35h error detected in routine   idsfft/)
      end
      subroutine  idtang(ndp,xd,yd,nt,ipt,nl,ipl,iwl,iwp,wk)
c this subroutine performs triangulation.  it divides the x-y
c plane into a number of triangles according to given data
c points in the plane, determines line segments that form the
c border of data area, and determines the triangle numbers
c corresponding to the border line segments.
c at completion, point numbers of the vertexes of each triangle
c are listed counter-clockwise.  point numbers of the end points
c of each border line segment are listed counter-clockwise,
c listing order of the line segments being counter-clockwise.
c the lun constant in the data initialization statement is the
c logical unit number of the standard output unit and is,
c therefore, system dependent.
c this subroutine calls the idxchg function.
c the input parameters are
c     ndp = number of data points,
c     xd  = array of dimension ndp containing the
c           x coordinates of the data points,
c     yd  = array of dimension ndp containing the
c           y coordinates of the data points.
c the output parameters are
c     nt  = number of triangles,
c     ipt = integer array of dimension 6*ndp-15, where the
c           point numbers of the vertexes of the (it)th
c           triangle are to be stored as the (3*it-2)nd,
c           (3*it-1)st, and (3*it)th elements,
c           it=1,2,...,nt,
c     nl  = number of border line segments,
c     ipl = integer array of dimension 6*ndp, where the
c           point numbers of the end points of the (il)th
c           border line segment and its respective triangle
c           number are to be stored as the (3*il-2)nd,
c           (3*il-1)st, and (3*il)th elements,
c           il=1,2,..., nl.
c the other parameters are
c     iwl = integer array of dimension 18*ndp used
c           internally as a work area,
c     iwp = integer array of dimension ndp used
c           internally as a work area,
c     wk  = array of dimension ndp used internally as a
c           work area.
c declaration statements
      dimension   xd(100),yd(100),ipt(585),ipl(600),
     1            iwl(1800),iwp(100),wk(100)
      dimension   itf(2)
      data  ratio/1.0e-6/, nrep/100/, lun/6/
c statement functions
      dsqf(u1,v1,u2,v2)=(u2-u1)**2+(v2-v1)**2
      side(u1,v1,u2,v2,u3,v3)=(v3-v1)*(u2-u1)-(u3-u1)*(v2-v1)
c preliminary processing
   10 ndp0=ndp
      ndpm1=ndp0-1
      if(ndp0.lt.4)       go to 90
c determines the closest pair of data points and their midpoint.
   20 dsqmn=dsqf(xd(1),yd(1),xd(2),yd(2))
      ipmn1=1
      ipmn2=2
      do 22  ip1=1,ndpm1
        x1=xd(ip1)
        y1=yd(ip1)
        ip1p1=ip1+1
        do 21  ip2=ip1p1,ndp0
          dsqi=dsqf(x1,y1,xd(ip2),yd(ip2))
          if(dsqi.eq.0.0)      go to 91
          if(dsqi.ge.dsqmn)    go to 21
          dsqmn=dsqi
          ipmn1=ip1
          ipmn2=ip2
   21   continue
   22 continue
      dsq12=dsqmn
      xdmp=(xd(ipmn1)+xd(ipmn2))/2.0
      ydmp=(yd(ipmn1)+yd(ipmn2))/2.0
c sorts the other (ndp-2) data points in ascending order of
c distance from the midpoint and stores the sorted data point
c numbers in the iwp array.
   30 jp1=2
      do 31  ip1=1,ndp0
        if(ip1.eq.ipmn1.or.ip1.eq.ipmn2)      go to 31
        jp1=jp1+1
        iwp(jp1)=ip1
        wk(jp1)=dsqf(xdmp,ydmp,xd(ip1),yd(ip1))
   31 continue
      do 33  jp1=3,ndpm1
        dsqmn=wk(jp1)
        jpmn=jp1
        do 32  jp2=jp1,ndp0
          if(wk(jp2).ge.dsqmn)      go to 32
          dsqmn=wk(jp2)
          jpmn=jp2
   32   continue
        its=iwp(jp1)
        iwp(jp1)=iwp(jpmn)
        iwp(jpmn)=its
        wk(jpmn)=wk(jp1)
   33 continue
c if necessary, modifies the ordering in such a way that the
c first three data points are not collinear.
   35 ar=dsq12*ratio
      x1=xd(ipmn1)
      y1=yd(ipmn1)
      dx21=xd(ipmn2)-x1
      dy21=yd(ipmn2)-y1
      do 36  jp=3,ndp0
        ip=iwp(jp)
        if(abs((yd(ip)-y1)*dx21-(xd(ip)-x1)*dy21).gt.ar)
     1               go to 37
   36 continue
      go to 92
   37 if(jp.eq.3)    go to 40
      jpmx=jp
      jp=jpmx+1
      do 38  jpc=4,jpmx
        jp=jp-1
        iwp(jp)=iwp(jp-1)
   38 continue
      iwp(3)=ip
c forms the first triangle.  stores point numbers of the ver-
c texes of the triangle in the ipt array, and stores point num-
c bers of the border line segments and the triangle number in
c the ipl array.
   40 ip1=ipmn1
      ip2=ipmn2
      ip3=iwp(3)
      if(side(xd(ip1),yd(ip1),xd(ip2),yd(ip2),xd(ip3),yd(ip3))
     1     .ge.0.0)       go to 41
      ip1=ipmn2
      ip2=ipmn1
   41 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
c adds the remaining (ndp-3) data points, one by one.
   50 do 79  jp1=4,ndp0
        ip1=iwp(jp1)
        x1=xd(ip1)
        y1=yd(ip1)
c - determines the visible border line segments.
        ip2=ipl(1)
        jpmn=1
        dxmn=xd(ip2)-x1
        dymn=yd(ip2)-y1
        dsqmn=dxmn**2+dymn**2
        armn=dsqmn*ratio
        jpmx=1
        dxmx=dxmn
        dymx=dymn
        dsqmx=dsqmn
        armx=armn
        do 52  jp2=2,nl0
          ip2=ipl(3*jp2-2)
          dx=xd(ip2)-x1
          dy=yd(ip2)-y1
          ar=dy*dxmn-dx*dymn
          if(ar.gt.armn)       go to 51
          dsqi=dx**2+dy**2
          if(ar.ge.(-armn).and.dsqi.ge.dsqmn)      go to 51
          jpmn=jp2
          dxmn=dx
          dymn=dy
          dsqmn=dsqi
          armn=dsqmn*ratio
   51     ar=dy*dxmx-dx*dymx
          if(ar.lt.(-armx))    go to 52
          dsqi=dx**2+dy**2
          if(ar.le.armx.and.dsqi.ge.dsqmx)    go to 52
          jpmx=jp2
          dxmx=dx
          dymx=dy
          dsqmx=dsqi
          armx=dsqmx*ratio
   52   continue
        if(jpmx.lt.jpmn)  jpmx=jpmx+nl0
        nsh=jpmn-1
        if(nsh.le.0)      go to 60
c - shifts (rotates) the ipl array to have the invisible border
c - line segments contained in the first part of the ipl array.
        nsht3=nsh*3
        do 53  jp2t3=3,nsht3,3
          jp3t3=jp2t3+nlt3
          ipl(jp3t3-2)=ipl(jp2t3-2)
          ipl(jp3t3-1)=ipl(jp2t3-1)
          ipl(jp3t3)  =ipl(jp2t3)
   53   continue
        do 54  jp2t3=3,nlt3,3
          jp3t3=jp2t3+nsht3
          ipl(jp2t3-2)=ipl(jp3t3-2)
          ipl(jp2t3-1)=ipl(jp3t3-1)
          ipl(jp2t3)  =ipl(jp3t3)
   54   continue
        jpmx=jpmx-nsh
c - adds triangles to the ipt array, updates border line
c - segments in the ipl array, and sets flags for the border
c - line segments to be reexamined in the iwl array.
   60   jwl=0
        do 64  jp2=jpmx,nl0
          jp2t3=jp2*3
          ipl1=ipl(jp2t3-2)
          ipl2=ipl(jp2t3-1)
          it  =ipl(jp2t3)
c - - adds a triangle to the ipt array.
          nt0=nt0+1
          ntt3=ntt3+3
          ipt(ntt3-2)=ipl2
          ipt(ntt3-1)=ipl1
          ipt(ntt3)  =ip1
c - - updates border line segments in the ipl array.
          if(jp2.ne.jpmx)      go to 61
          ipl(jp2t3-1)=ip1
          ipl(jp2t3)  =nt0
   61     if(jp2.ne.nl0)       go to 62
          nln=jpmx+1
          nlnt3=nln*3
          ipl(nlnt3-2)=ip1
          ipl(nlnt3-1)=ipl(1)
          ipl(nlnt3)  =nt0
c - - determines the vertex that does not lie on the border
c - - line segments.
   62     itt3=it*3
          ipti=ipt(itt3-2)
          if(ipti.ne.ipl1.and.ipti.ne.ipl2)   go to 63
          ipti=ipt(itt3-1)
          if(ipti.ne.ipl1.and.ipti.ne.ipl2)   go to 63
          ipti=ipt(itt3)
c - - checks if the exchange is necessary.
   63     if(idxchg(xd,yd,ip1,ipti,ipl1,ipl2).eq.0)     go to 64
c - - modifies the ipt array when necessary.
          ipt(itt3-2)=ipti
          ipt(itt3-1)=ipl1
          ipt(itt3)  =ip1
          ipt(ntt3-1)=ipti
          if(jp2.eq.jpmx)      ipl(jp2t3)=it
          if(jp2.eq.nl0.and.ipl(3).eq.it)     ipl(3)=nt0
c - - sets flags in the iwl array.
          jwl=jwl+4
          iwl(jwl-3)=ipl1
          iwl(jwl-2)=ipti
          iwl(jwl-1)=ipti
          iwl(jwl)  =ipl2
   64   continue
        nl0=nln
        nlt3=nlnt3
        nlf=jwl/2
        if(nlf.eq.0)      go to 79
c - improves triangulation.
   70   ntt3p3=ntt3+3
        do 78  irep=1,nrep
          do 76  ilf=1,nlf
            ilft2=ilf*2
            ipl1=iwl(ilft2-1)
            ipl2=iwl(ilft2)
c - - locates in the ipt array two triangles on both sides of
c - - the flagged line segment.
            ntf=0
            do 71  itt3r=3,ntt3,3
              itt3=ntt3p3-itt3r
              ipt1=ipt(itt3-2)
              ipt2=ipt(itt3-1)
              ipt3=ipt(itt3)
              if(ipl1.ne.ipt1.and.ipl1.ne.ipt2.and.
     1           ipl1.ne.ipt3)      go to 71
              if(ipl2.ne.ipt1.and.ipl2.ne.ipt2.and.
     1           ipl2.ne.ipt3)      go to 71
              ntf=ntf+1
              itf(ntf)=itt3/3
              if(ntf.eq.2)     go to 72
   71       continue
            if(ntf.lt.2)       go to 76
c - - determines the vertexes of the triangles that do not lie
c - - on the line segment.
   72       it1t3=itf(1)*3
            ipti1=ipt(it1t3-2)
            if(ipti1.ne.ipl1.and.ipti1.ne.ipl2)    go to 73
            ipti1=ipt(it1t3-1)
            if(ipti1.ne.ipl1.and.ipti1.ne.ipl2)    go to 73
            ipti1=ipt(it1t3)
   73       it2t3=itf(2)*3
            ipti2=ipt(it2t3-2)
            if(ipti2.ne.ipl1.and.ipti2.ne.ipl2)    go to 74
            ipti2=ipt(it2t3-1)
            if(ipti2.ne.ipl1.and.ipti2.ne.ipl2)    go to 74
            ipti2=ipt(it2t3)
c - - checks if the exchange is necessary.
   74       if(idxchg(xd,yd,ipti1,ipti2,ipl1,ipl2).eq.0)
     1         go to 76
c - - modifies the ipt array when necessary.
            ipt(it1t3-2)=ipti1
            ipt(it1t3-1)=ipti2
            ipt(it1t3)  =ipl1
            ipt(it2t3-2)=ipti2
            ipt(it2t3-1)=ipti1
            ipt(it2t3)  =ipl2
c - - sets 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 75  jlt3=3,nlt3,3
              iplj1=ipl(jlt3-2)
              iplj2=ipl(jlt3-1)
              if((iplj1.eq.ipl1.and.iplj2.eq.ipti2).or.
     1           (iplj2.eq.ipl1.and.iplj1.eq.ipti2))
     2                         ipl(jlt3)=itf(1)
              if((iplj1.eq.ipl2.and.iplj2.eq.ipti1).or.
     1           (iplj2.eq.ipl2.and.iplj1.eq.ipti1))
     2                         ipl(jlt3)=itf(2)
   75       continue
   76     continue
          nlfc=nlf
          nlf=jwl/2
          if(nlf.eq.nlfc)      go to 79
c - - resets the iwl array for the next round.
          jwl=0
          jwl1mn=(nlfc+1)*2
          nlft2=nlf*2
          do 77  jwl1=jwl1mn,nlft2,2
            jwl=jwl+2
            iwl(jwl-1)=iwl(jwl1-1)
            iwl(jwl)  =iwl(jwl1)
   77     continue
          nlf=jwl/2
   78   continue
   79 continue
c rearranges the ipt array so that the vertexes of each triangle
c are listed counter-clockwise.
   80 do 81  itt3=3,ntt3,3
        ip1=ipt(itt3-2)
        ip2=ipt(itt3-1)
        ip3=ipt(itt3)
        if(side(xd(ip1),yd(ip1),xd(ip2),yd(ip2),xd(ip3),yd(ip3))
     1       .ge.0.0)     go to 81
        ipt(itt3-2)=ip2
        ipt(itt3-1)=ip1
   81 continue
      nt=nt0
      nl=nl0
      return
c error exit
   90 write (lun,2090)  ndp0
      go to 93
   91 write (lun,2091)  ndp0,ip1,ip2,x1,y1
      go to 93
   92 write (lun,2092)  ndp0
   93 write (lun,2093)
      nt=0
      return
c format statements
 2090 format(1x/23h ***   ndp less than 4./8h   ndp =,i5)
 2091 format(1x/29h ***   identical data points./
     1   8h   ndp =,i5,5x,5hip1 =,i5,5x,5hip2 =,i5,
     2   5x,4hxd =,e12.4,5x,4hyd =,e12.4)
 2092 format(1x/33h ***   all collinear data points./
     1   8h   ndp =,i5)
 2093 format(35h error detected in routine   idtang/)
      end
      function  idxchg(x,y,i1,i2,i3,i4)
c this function determines whether or not the exchange of two
c triangles is necessary on the basis of max-min-angle criterion
c by c. l. lawson.
c the input parameters are
c     x,y = arrays containing the coordinates of the data
c           points,
c     i1,i2,i3,i4 = point numbers of four points p1, p2,
c           p3, and p4 that form a quadrilateral with p3
c           and p4 connected diagonally.
c this function returns an integer value 1 (one) when an ex-
c change is necessary, and 0 (zero) otherwise.
c declaration statements
      dimension   x(100),y(100)
      equivalence (c2sq,c1sq),(a3sq,b2sq),(b3sq,a1sq),
     1            (a4sq,b1sq),(b4sq,a2sq),(c4sq,c3sq)
c preliminary processing
   10 x1=x(i1)
      y1=y(i1)
      x2=x(i2)
      y2=y(i2)
      x3=x(i3)
      y3=y(i3)
      x4=x(i4)
      y4=y(i4)
c calculation
   20 idx=0
      u3=(y2-y3)*(x1-x3)-(x2-x3)*(y1-y3)
      u4=(y1-y4)*(x2-x4)-(x1-x4)*(y2-y4)
      if(u3*u4.le.0.0)    go to 30
      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
      b1sq=(x4-x1)**2+(y4-y1)**2
      c1sq=(x3-x4)**2+(y3-y4)**2
      a2sq=(x2-x4)**2+(y2-y4)**2
      b2sq=(x3-x2)**2+(y3-y2)**2
      c3sq=(x2-x1)**2+(y2-y1)**2
      s1sq=u1*u1/(c1sq*amax1(a1sq,b1sq))
      s2sq=u2*u2/(c2sq*amax1(a2sq,b2sq))
      s3sq=u3*u3/(c3sq*amax1(a3sq,b3sq))
      s4sq=u4*u4/(c4sq*amax1(a4sq,b4sq))
      if(amin1(s1sq,s2sq).lt.amin1(s3sq,s4sq))     idx=1
   30 idxchg=idx
      return
      end