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     istat = error message.
  c     status= Starlink error status
  
  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 this subroutine calls the idcldp, idgrid, idpdrv, idptip, and
  c idtang subroutines.