module dierckx

  use const_and_precisions, only : wp_
  implicit none

contains

  subroutine bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk, &
     iwrk,kwrk,ier)
!  subroutine bispev evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
!  ,my a bivariate spline s(x,y) of degrees kx and ky, given in the
!  b-spline representation.
!
!  calling sequence:
! call bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk,
!    * iwrk,kwrk,ier)
!
!  input parameters:
!   tx    : real array, length nx, which contains the position of the
!       knots in the x-direction.
!   nx    : integer, giving the total number of knots in the x-direction
!   ty    : real array, length ny, which contains the position of the
!       knots in the y-direction.
!   ny    : integer, giving the total number of knots in the y-direction
!   c     : real array, length (nx-kx-1)*(ny-ky-1), which contains the
!       b-spline coefficients.
!   kx,ky : integer values, giving the degrees of the spline.
!   x     : real array of dimension (mx).
!       before entry x(i) must be set to the x co-ordinate of the
!       i-th grid point along the x-axis.
!       tx(kx+1)<=x(i-1)<=x(i)<=tx(nx-kx), i=2,...,mx.
!   mx    : on entry mx must specify the number of grid points along
!       the x-axis. mx >=1.
!   y     : real array of dimension (my).
!       before entry y(j) must be set to the y co-ordinate of the
!       j-th grid point along the y-axis.
!       ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my.
!   my    : on entry my must specify the number of grid points along
!       the y-axis. my >=1.
!   wrk   : real array of dimension lwrk. used as workspace.
!   lwrk  : integer, specifying the dimension of wrk.
!       lwrk >= mx*(kx+1)+my*(ky+1)
!   iwrk  : integer array of dimension kwrk. used as workspace.
!   kwrk  : integer, specifying the dimension of iwrk. kwrk >= mx+my.
!
!  output parameters:
!   z     : real array of dimension (mx*my).
!       on succesful exit z(my*(i-1)+j) contains the value of s(x,y)
!       at the point (x(i),y(j)),i=1,...,mx;j=1,...,my.
!   ier   : integer error flag
!    ier=0 : normal return
!    ier=10: invalid input data (see restrictions)
!
!  restrictions:
!   mx >=1, my >=1, lwrk>=mx*(kx+1)+my*(ky+1), kwrk>=mx+my
!   tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx
!   ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my
!
!  other subroutines required:
!    fpbisp,fpbspl
!
!  references :
!    de boor c : on calculating with b-splines, j. approximation theory
!            6 (1972) 50-62.
!    cox m.g.  : the numerical evaluation of b-splines, j. inst. maths
!            applics 10 (1972) 134-149.
!    dierckx p. : curve and surface fitting with splines, monographs on
!             numerical analysis, oxford university press, 1993.
!
!  author :
!    p.dierckx
!    dept. computer science, k.u.leuven
!    celestijnenlaan 200a, b-3001 heverlee, belgium.
!    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
!  latest update : march 1987
!
    implicit none
! arguments
    integer, intent(in) :: nx, ny, kx, ky, mx, my, lwrk, kwrk
    integer, intent(out) :: ier
    integer, intent(inout) :: iwrk(kwrk)
    real(wp_), intent(in) :: tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(mx), y(my)
    real(wp_), intent(out) :: z(mx*my)
    real(wp_), intent(inout) :: wrk(lwrk)
! local variables
    integer :: i, iw, lwest
!  ..
!  before starting computations a data check is made. if the input data
!  are invalid control is immediately repassed to the calling program.
    ier = 10
    lwest = (kx+1)*mx+(ky+1)*my
    if(lwrk<lwest .or. kwrk<(mx+my) .or. mx<1 .or. my<1) return
    do i=2,mx
      if(x(i)<x(i-1)) return
    end do
    do i=2,my
      if(y(i)<y(i-1)) return
    end do
    ier = 0
    iw = mx*(kx+1)+1
    call fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk(1),wrk(iw), &
     iwrk(1),iwrk(mx+1))
  end subroutine bispev

  subroutine surfit(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest, &
      nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier)
! given the set of data points (x(i),y(i),z(i)) and the set of positive
! numbers w(i),i=1,...,m, subroutine surfit determines a smooth bivar-
! iate spline approximation s(x,y) of degrees kx and ky on the rect-
! angle xb <= x <= xe, yb <= y <= ye.
! if iopt = -1 surfit calculates the weighted least-squares spline
! according to a given set of knots.
! if iopt >= 0 the total numbers nx and ny of these knots and their
! position tx(j),j=1,...,nx and ty(j),j=1,...,ny are chosen automatic-
! ally by the routine. the smoothness of s(x,y) is then achieved by
! minimalizing the discontinuity jumps in the derivatives of s(x,y)
! across the boundaries of the subpanels (tx(i),tx(i+1))*(ty(j),ty(j+1).
! the amounth of smoothness is determined by the condition that f(p) =
! sum ((w(i)*(z(i)-s(x(i),y(i))))**2) be <= s, with s a given non-neg-
! ative constant, called the smoothing factor.
! the fit is given in the b-spline representation (b-spline coefficients
! c((ny-ky-1)*(i-1)+j),i=1,...,nx-kx-1;j=1,...,ny-ky-1) and can be eval-
! uated by means of subroutine bispev.
!
! calling sequence:
! call surfit(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,
!    *  nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ier)
!
! parameters:
!  iopt  : integer flag. on entry iopt must specify whether a weighted
!      least-squares spline (iopt=-1) or a smoothing spline (iopt=0
!      or 1) must be determined.
!      if iopt=0 the routine will start with an initial set of knots
!      tx(i)=xb,tx(i+kx+1)=xe,i=1,...,kx+1;ty(i)=yb,ty(i+ky+1)=ye,i=
!      1,...,ky+1. if iopt=1 the routine will continue with the set
!      of knots found at the last call of the routine.
!      attention: a call with iopt=1 must always be immediately pre-
!                 ceded by another call with iopt=1 or iopt=0.
!      unchanged on exit.
!  m     : integer. on entry m must specify the number of data points.
!      m >= (kx+1)*(ky+1). unchanged on exit.
!  x     : real array of dimension at least (m).
!  y     : real array of dimension at least (m).
!  z     : real array of dimension at least (m).
!      before entry, x(i),y(i),z(i) must be set to the co-ordinates
!      of the i-th data point, for i=1,...,m. the order of the data
!      points is immaterial. unchanged on exit.
!  w     : real array of dimension at least (m). before entry, w(i) must
!      be set to the i-th value in the set of weights. the w(i) must
!      be strictly positive. unchanged on exit.
!  xb,xe : real values. on entry xb,xe,yb and ye must specify the bound-
!  yb,ye   aries of the rectangular approximation domain.
!      xb<=x(i)<=xe,yb<=y(i)<=ye,i=1,...,m. unchanged on exit.
!  kx,ky : integer values. on entry kx and ky must specify the degrees
!      of the spline. 1<=kx,ky<=5. it is recommended to use bicubic
!      (kx=ky=3) splines. unchanged on exit.
!  s     : real. on entry (in case iopt>=0) s must specify the smoothing
!      factor. s >=0. unchanged on exit.
!      for advice on the choice of s see further comments
!  nxest : integer. unchanged on exit.
!  nyest : integer. unchanged on exit.
!      on entry, nxest and nyest must specify an upper bound for the
!      number of knots required in the x- and y-directions respect.
!      these numbers will also determine the storage space needed by
!      the routine. nxest >= 2*(kx+1), nyest >= 2*(ky+1).
!      in most practical situation nxest = kx+1+sqrt(m/2), nyest =
!      ky+1+sqrt(m/2) will be sufficient. see also further comments.
!  nmax  : integer. on entry nmax must specify the actual dimension of
!      the arrays tx and ty. nmax >= nxest, nmax >=nyest.
!      unchanged on exit.
!  eps   : real.
!      on entry, eps must specify a threshold for determining the
!      effective rank of an over-determined linear system of equat-
!      ions. 0 < eps < 1.  if the number of decimal digits in the
!      computer representation of a real number is q, then 10**(-q)
!      is a suitable value for eps in most practical applications.
!      unchanged on exit.
!  nx    : integer.
!      unless ier=10 (in case iopt >=0), nx will contain the total
!      number of knots with respect to the x-variable, of the spline
!      approximation returned. if the computation mode iopt=1 is
!      used, the value of nx should be left unchanged between sub-
!      sequent calls.
!      in case iopt=-1, the value of nx should be specified on entry
!  tx    : real array of dimension nmax.
!      on succesful exit, this array will contain the knots of the
!      spline with respect to the x-variable, i.e. the position of
!      the interior knots tx(kx+2),...,tx(nx-kx-1) as well as the
!      position of the additional knots tx(1)=...=tx(kx+1)=xb and
!      tx(nx-kx)=...=tx(nx)=xe needed for the b-spline representat.
!      if the computation mode iopt=1 is used, the values of tx(1),
!      ...,tx(nx) should be left unchanged between subsequent calls.
!      if the computation mode iopt=-1 is used, the values tx(kx+2),
!      ...tx(nx-kx-1) must be supplied by the user, before entry.
!      see also the restrictions (ier=10).
!  ny    : integer.
!      unless ier=10 (in case iopt >=0), ny will contain the total
!      number of knots with respect to the y-variable, of the spline
!      approximation returned. if the computation mode iopt=1 is
!      used, the value of ny should be left unchanged between sub-
!      sequent calls.
!      in case iopt=-1, the value of ny should be specified on entry
!  ty    : real array of dimension nmax.
!      on succesful exit, this array will contain the knots of the
!      spline with respect to the y-variable, i.e. the position of
!      the interior knots ty(ky+2),...,ty(ny-ky-1) as well as the
!      position of the additional knots ty(1)=...=ty(ky+1)=yb and
!      ty(ny-ky)=...=ty(ny)=ye needed for the b-spline representat.
!      if the computation mode iopt=1 is used, the values of ty(1),
!      ...,ty(ny) should be left unchanged between subsequent calls.
!      if the computation mode iopt=-1 is used, the values ty(ky+2),
!      ...ty(ny-ky-1) must be supplied by the user, before entry.
!      see also the restrictions (ier=10).
!  c     : real array of dimension at least (nxest-kx-1)*(nyest-ky-1).
!      on succesful exit, c contains the coefficients of the spline
!      approximation s(x,y)
!  fp    : real. unless ier=10, fp contains the weighted sum of
!      squared residuals of the spline approximation returned.
!  wrk1  : real array of dimension (lwrk1). used as workspace.
!      if the computation mode iopt=1 is used the value of wrk1(1)
!      should be left unchanged between subsequent calls.
!      on exit wrk1(2),wrk1(3),...,wrk1(1+(nx-kx-1)*(ny-ky-1)) will
!      contain the values d(i)/max(d(i)),i=1,...,(nx-kx-1)*(ny-ky-1)
!      with d(i) the i-th diagonal element of the reduced triangular
!      matrix for calculating the b-spline coefficients. it includes
!      those elements whose square is less than eps,which are treat-
!      ed as 0 in the case of presumed rank deficiency (ier<-2).
!  lwrk1 : integer. on entry lwrk1 must specify the actual dimension of
!      the array wrk1 as declared in the calling (sub)program.
!      lwrk1 must not be too small. let
!        u = nxest-kx-1, v = nyest-ky-1, km = max(kx,ky)+1,
!        ne = max(nxest,nyest), bx = kx*v+ky+1, by = ky*u+kx+1,
!        if(bx<=by) b1 = bx, b2 = b1+v-ky
!        if(bx>by) b1 = by, b2 = b1+u-kx  then
!      lwrk1 >= u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1
!  wrk2  : real array of dimension (lwrk2). used as workspace, but
!      only in the case a rank deficient system is encountered.
!  lwrk2 : integer. on entry lwrk2 must specify the actual dimension of
!      the array wrk2 as declared in the calling (sub)program.
!      lwrk2 > 0 . a save upper boundfor lwrk2 = u*v*(b2+1)+b2
!      where u,v and b2 are as above. if there are enough data
!      points, scattered uniformly over the approximation domain
!      and if the smoothing factor s is not too small, there is a
!      good chance that this extra workspace is not needed. a lot
!      of memory might therefore be saved by setting lwrk2=1.
!      (see also ier > 10)
!  iwrk  : integer array of dimension (kwrk). used as workspace.
!  kwrk  : integer. on entry kwrk must specify the actual dimension of
!      the array iwrk as declared in the calling (sub)program.
!      kwrk >= m+(nxest-2*kx-1)*(nyest-2*ky-1).
!  ier   : integer. unless the routine detects an error, ier contains a
!      non-positive value on exit, i.e.
!   ier=0  : normal return. the spline returned has a residual sum of
!        squares fp such that abs(fp-s)/s <= tol with tol a relat-
!        ive tolerance set to 0.001 by the program.
!   ier=-1 : normal return. the spline returned is an interpolating
!        spline (fp=0).
!   ier=-2 : normal return. the spline returned is the weighted least-
!        squares polynomial of degrees kx and ky. in this extreme
!        case fp gives the upper bound for the smoothing factor s.
!   ier<-2 : warning. the coefficients of the spline returned have been
!        computed as the minimal norm least-squares solution of a
!        (numerically) rank deficient system. (-ier) gives the rank.
!        especially if the rank deficiency which can be computed as
!        (nx-kx-1)*(ny-ky-1)+ier, is large the results may be inac-
!        curate. they could also seriously depend on the value of
!        eps.
!   ier=1  : error. the required storage space exceeds the available
!        storage space, as specified by the parameters nxest and
!        nyest.
!        probably causes : nxest or nyest too small. if these param-
!        eters are already large, it may also indicate that s is
!        too small
!        the approximation returned is the weighted least-squares
!        spline according to the current set of knots.
!        the parameter fp gives the corresponding weighted sum of
!        squared residuals (fp>s).
!   ier=2  : error. a theoretically impossible result was found during
!        the iteration proces for finding a smoothing spline with
!        fp = s. probably causes : s too small or badly chosen eps.
!        there is an approximation returned but the corresponding
!        weighted sum of squared residuals does not satisfy the
!        condition abs(fp-s)/s < tol.
!   ier=3  : error. the maximal number of iterations maxit (set to 20
!        by the program) allowed for finding a smoothing spline
!        with fp=s has been reached. probably causes : s too small
!        there is an approximation returned but the corresponding
!        weighted sum of squared residuals does not satisfy the
!        condition abs(fp-s)/s < tol.
!   ier=4  : error. no more knots can be added because the number of
!        b-spline coefficients (nx-kx-1)*(ny-ky-1) already exceeds
!        the number of data points m.
!        probably causes : either s or m too small.
!        the approximation returned is the weighted least-squares
!        spline according to the current set of knots.
!        the parameter fp gives the corresponding weighted sum of
!        squared residuals (fp>s).
!   ier=5  : error. no more knots can be added because the additional
!        knot would (quasi) coincide with an old one.
!        probably causes : s too small or too large a weight to an
!        inaccurate data point.
!        the approximation returned is the weighted least-squares
!        spline according to the current set of knots.
!        the parameter fp gives the corresponding weighted sum of
!        squared residuals (fp>s).
!   ier=10 : error. on entry, the input data are controlled on validity
!        the following restrictions must be satisfied.
!        -1<=iopt<=1, 1<=kx,ky<=5, m>=(kx+1)*(ky+1), nxest>=2*kx+2,
!        nyest>=2*ky+2, 0<eps<1, nmax>=nxest, nmax>=nyest,
!        xb<=x(i)<=xe, yb<=y(i)<=ye, w(i)>0, i=1,...,m
!        lwrk1 >= u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1
!        kwrk >= m+(nxest-2*kx-1)*(nyest-2*ky-1)
!        if iopt=-1: 2*kx+2<=nx<=nxest
!                    xb<tx(kx+2)<tx(kx+3)<...<tx(nx-kx-1)<xe
!                    2*ky+2<=ny<=nyest
!                    yb<ty(ky+2)<ty(ky+3)<...<ty(ny-ky-1)<ye
!        if iopt>=0: s>=0
!        if one of these conditions is found to be violated,control
!        is immediately repassed to the calling program. in that
!        case there is no approximation returned.
!   ier>10 : error. lwrk2 is too small, i.e. there is not enough work-
!        space for computing the minimal least-squares solution of
!        a rank deficient system of linear equations. ier gives the
!        requested value for lwrk2. there is no approximation re-
!        turned but, having saved the information contained in nx,
!        ny,tx,ty,wrk1, and having adjusted the value of lwrk2 and
!        the dimension of the array wrk2 accordingly, the user can
!        continue at the point the program was left, by calling
!        surfit with iopt=1.
!
! further comments:
!  by means of the parameter s, the user can control the tradeoff
!   between closeness of fit and smoothness of fit of the approximation.
!   if s is too large, the spline will be too smooth and signal will be
!   lost ; if s is too small the spline will pick up too much noise. in
!   the extreme cases the program will return an interpolating spline if
!   s=0 and the weighted least-squares polynomial (degrees kx,ky)if s is
!   very large. between these extremes, a properly chosen s will result
!   in a good compromise between closeness of fit and smoothness of fit.
!   to decide whether an approximation, corresponding to a certain s is
!   satisfactory the user is highly recommended to inspect the fits
!   graphically.
!   recommended values for s depend on the weights w(i). if these are
!   taken as 1/d(i) with d(i) an estimate of the standard deviation of
!   z(i), a good s-value should be found in the range (m-sqrt(2*m),m+
!   sqrt(2*m)). if nothing is known about the statistical error in z(i)
!   each w(i) can be set equal to one and s determined by trial and
!   error, taking account of the comments above. the best is then to
!   start with a very large value of s ( to determine the least-squares
!   polynomial and the corresponding upper bound fp0 for s) and then to
!   progressively decrease the value of s ( say by a factor 10 in the
!   beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
!   approximation shows more detail) to obtain closer fits.
!   to choose s very small is strongly discouraged. this considerably
!   increases computation time and memory requirements. it may also
!   cause rank-deficiency (ier<-2) and endager numerical stability.
!   to economize the search for a good s-value the program provides with
!   different modes of computation. at the first call of the routine, or
!   whenever he wants to restart with the initial set of knots the user
!   must set iopt=0.
!   if iopt=1 the program will continue with the set of knots found at
!   the last call of the routine. this will save a lot of computation
!   time if surfit is called repeatedly for different values of s.
!   the number of knots of the spline returned and their location will
!   depend on the value of s and on the complexity of the shape of the
!   function underlying the data. if the computation mode iopt=1
!   is used, the knots returned may also depend on the s-values at
!   previous calls (if these were smaller). therefore, if after a number
!   of trials with different s-values and iopt=1, the user can finally
!   accept a fit as satisfactory, it may be worthwhile for him to call
!   surfit once more with the selected value for s but now with iopt=0.
!   indeed, surfit may then return an approximation of the same quality
!   of fit but with fewer knots and therefore better if data reduction
!   is also an important objective for the user.
!   the number of knots may also depend on the upper bounds nxest and
!   nyest. indeed, if at a certain stage in surfit the number of knots
!   in one direction (say nx) has reached the value of its upper bound
!   (nxest), then from that moment on all subsequent knots are added
!   in the other (y) direction. this may indicate that the value of
!   nxest is too small. on the other hand, it gives the user the option
!   of limiting the number of knots the routine locates in any direction
!   for example, by setting nxest=2*kx+2 (the lowest allowable value for
!   nxest), the user can indicate that he wants an approximation which
!   is a simple polynomial of degree kx in the variable x.
!
!  other subroutines required:
!    fpback,fpbspl,fpsurf,fpdisc,fpgivs,fprank,fprati,fprota,fporde
!
!  references:
!   dierckx p. : an algorithm for surface fitting with spline functions
!            ima j. numer. anal. 1 (1981) 267-283.
!   dierckx p. : an algorithm for surface fitting with spline functions
!            report tw50, dept. computer science,k.u.leuven, 1980.
!   dierckx p. : curve and surface fitting with splines, monographs on
!            numerical analysis, oxford university press, 1993.
!
!  author:
!    p.dierckx
!    dept. computer science, k.u. leuven
!    celestijnenlaan 200a, b-3001 heverlee, belgium.
!    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
!  creation date : may 1979
!  latest update : march 1987
!
!  ..
    implicit none
!  ..scalar arguments..
    real(wp_) xb,xe,yb,ye,s,eps,fp
    integer iopt,m,kx,ky,nxest,nyest,nmax,nx,ny,lwrk1,lwrk2,kwrk,ier
!  ..array arguments..
    real(wp_) x(m),y(m),z(m),w(m),tx(nmax),ty(nmax), &
     c((nxest-kx-1)*(nyest-ky-1)),wrk1(lwrk1),wrk2(lwrk2)
    integer iwrk(kwrk)
!  ..local scalars..
    real(wp_) tol
    integer i,ib1,ib3,jb1,ki,kmax,km1,km2,kn,kwest,kx1,ky1,la,lbx, &
     lby,lco,lf,lff,lfp,lh,lq,lsx,lsy,lwest,maxit,ncest,nest,nek, &
     nminx,nminy,nmx,nmy,nreg,nrint,nxk,nyk
!  ..function references..
    integer max0
!  ..subroutine references..
!    fpsurf
!  ..
!  we set up the parameters tol and maxit.
    maxit = 20
    tol = 0.1e-02
!  before starting computations a data check is made. if the input data
!  are invalid,control is immediately repassed to the calling program.
    ier = 10
    if(eps<=0. .or. eps>=1.) return
    if(kx<=0 .or. kx>5) return
    kx1 = kx+1
    if(ky<=0 .or. ky>5) return
    ky1 = ky+1
    kmax = max0(kx,ky)
    km1 = kmax+1
    km2 = km1+1
    if(iopt<(-1) .or. iopt>1) return
    if(m<(kx1*ky1)) return
    nminx = 2*kx1
    if(nxest<nminx .or. nxest>nmax) return
    nminy = 2*ky1
    if(nyest<nminy .or. nyest>nmax) return
    nest = max0(nxest,nyest)
    nxk = nxest-kx1
    nyk = nyest-ky1
    ncest = nxk*nyk
    nmx = nxest-nminx+1
    nmy = nyest-nminy+1
    nrint = nmx+nmy
    nreg = nmx*nmy
    ib1 = kx*nyk+ky1
    jb1 = ky*nxk+kx1
    ib3 = kx1*nyk+1
    if(ib1>jb1) then
      ib1 = jb1
      ib3 = ky1*nxk+1
    end if
    lwest = ncest*(2+ib1+ib3)+2*(nrint+nest*km2+m*km1)+ib3
    kwest = m+nreg
    if(lwrk1<lwest .or. kwrk<kwest) return
    if(xb>=xe .or. yb>=ye) return
    do i=1,m
      if(w(i)<=0.) return
      if(x(i)<xb .or. x(i)>xe) return
      if(y(i)<yb .or. y(i)>ye) return
    end do
    if(iopt<0) then
      if(nx<nminx .or. nx>nxest) return
      nxk = nx-kx1
      tx(kx1) = xb
      tx(nxk+1) = xe
      do i=kx1,nxk
        if(tx(i+1)<=tx(i)) return
      end do
      if(ny<nminy .or. ny>nyest) return
      nyk = ny-ky1
      ty(ky1) = yb
      ty(nyk+1) = ye
      do i=ky1,nyk
        if(ty(i+1)<=ty(i)) return
      end do
    else
      if(s<0.) return
    end if
    ier = 0
!  we partition the working space and determine the spline approximation
    kn = 1
    ki = kn+m
    lq = 2
    la = lq+ncest*ib3
    lf = la+ncest*ib1
    lff = lf+ncest
    lfp = lff+ncest
    lco = lfp+nrint
    lh = lco+nrint
    lbx = lh+ib3
    nek = nest*km2
    lby = lbx+nek
    lsx = lby+nek
    lsy = lsx+m*km1
    call fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest, &
     eps,tol,maxit,nest,km1,km2,ib1,ib3,ncest,nrint,nreg,nx,tx, &
     ny,ty,c,fp,wrk1(1),wrk1(lfp),wrk1(lco),wrk1(lf),wrk1(lff), &
     wrk1(la),wrk1(lq),wrk1(lbx),wrk1(lby),wrk1(lsx),wrk1(lsy), &
     wrk1(lh),iwrk(ki),iwrk(kn),wrk2,lwrk2,ier)
  end subroutine surfit

  subroutine fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kxx,kyy,s,nxest, &
     nyest,eta,tol,maxit,nmax,km1,km2,ib1,ib3,nc,intest,nrest, &
     nx0,tx,ny0,ty,c,fp,fp0,fpint,coord,f,ff,a,q,bx,by,spx,spy,h, &
     idx,nummer,wrk,lwrk,ier)
!  ..
    implicit none
!  ..scalar arguments..
    real(wp_) xb,xe,yb,ye,s,eta,tol,fp,fp0
    integer iopt,m,kxx,kyy,nxest,nyest,maxit,nmax,km1,km2,ib1,ib3, &
     nc,intest,nrest,nx0,ny0,lwrk,ier
!  ..array arguments..
    real(wp_) x(m),y(m),z(m),w(m),tx(nmax),ty(nmax),c(nc),fpint(intest), &
     coord(intest),f(nc),ff(nc),a(nc,ib1),q(nc,ib3),bx(nmax,km2), &
     by(nmax,km2),spx(m,km1),spy(m,km1),h(ib3),wrk(lwrk)
    integer idx(nrest),nummer(m)
!  ..local scalars..
    real(wp_) acc,arg,cs,dmax,fac1,fac2,fpmax,fpms,f1,f2,f3,hxi,p,pinv, &
     piv,p1,p2,p3,sigma,sn,sq,store,wi,x0,x1,y0,y1,zi,eps, &
     rn,one,con1,con9,con4,ten
    integer i,iband,iband1,iband3,iband4,ibb,ichang,ich1,ich3,ii, &
     in,irot,iter,i1,i2,i3,j,jrot,jxy,j1,kx,kx1,kx2,ky,ky1,ky2,l, &
     la,lf,lh,lwest,lx,ly,l1,l2,n,ncof,nk1x,nk1y,nminx,nminy,nreg, &
     nrint,num,num1,nx,nxe,nxx,ny,nye,nyy,n1,rank
!  ..local arrays..
    real(wp_) hx(6),hy(6)
!  ..function references..
!  real(8) abs,sqrt
!  integer min0
!  ..subroutine references..
!    fpback,fpbspl,fpgivs,fpdisc,fporde,fprank,fprota
!  ..
!  set constants
    one = 0.1e+01
    con1 = 0.1e0
    con9 = 0.9e0
    con4 = 0.4e-01
    ten = 0.1e+02
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! part 1: determination of the number of knots and their position.     c
! ****************************************************************     c
! given a set of knots we compute the least-squares spline sinf(x,y),  c
! and the corresponding weighted sum of squared residuals fp=f(p=inf). c
! if iopt=-1  sinf(x,y) is the requested approximation.                c
! if iopt=0 or iopt=1 we check whether we can accept the knots:        c
!   if fp <=s we will continue with the current set of knots.          c
!   if fp > s we will increase the number of knots and compute the     c
!  corresponding least-squares spline until finally  fp<=s.        c
! the initial choice of knots depends on the value of s and iopt.      c
!   if iopt=0 we first compute the least-squares polynomial of degree  c
! kx in x and ky in y; nx=nminx=2*kx+2 and ny=nminy=2*ky+2.        c
! fp0=f(0) denotes the corresponding weighted sum of squared       c
! residuals                                                        c
!   if iopt=1 we start with the knots found at the last call of the    c
! routine, except for the case that s>=fp0; then we can compute    c
! the least-squares polynomial directly.                           c
! eventually the independent variables x and y (and the corresponding  c
! parameters) will be switched if this can reduce the bandwidth of the c
! system to be solved.                                                 c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!  ichang denotes whether(1) or not(-1) the directions have been inter-
!  changed.
    ichang = -1
    x0 = xb
    x1 = xe
    y0 = yb
    y1 = ye
    kx = kxx
    ky = kyy
    kx1 = kx+1
    ky1 = ky+1
    nxe = nxest
    nye = nyest
    eps = sqrt(eta)
    if(iopt<0) go to 20
!  calculation of acc, the absolute tolerance for the root of f(p)=s.
    acc = tol*s
    if(iopt==0) go to 10
    if(fp0>s) go to 20
!  initialization for the least-squares polynomial.
    10 continue
    nminx = 2*kx1
    nminy = 2*ky1
    nx = nminx
    ny = nminy
    ier = -2
    go to 30
    20 continue
    nx = nx0
    ny = ny0
!  main loop for the different sets of knots. m is a save upper bound
!  for the number of trials.
    30 continue
    do iter=1,m
!  find the position of the additional knots which are needed for the
!  b-spline representation of s(x,y).
      l = nx
      do i=1,kx1
        tx(i) = x0
        tx(l) = x1
        l = l-1
      end do
      l = ny
      do i=1,ky1
        ty(i) = y0
        ty(l) = y1
        l = l-1
      end do
!  find nrint, the total number of knot intervals and nreg, the number
!  of panels in which the approximation domain is subdivided by the
!  intersection of knots.
      nxx = nx-2*kx1+1
      nyy = ny-2*ky1+1
      nrint = nxx+nyy
      nreg = nxx*nyy
!  find the bandwidth of the observation matrix a.
!  if necessary, interchange the variables x and y, in order to obtain
!  a minimal bandwidth.
      iband1 = kx*(ny-ky1)+ky
      l = ky*(nx-kx1)+kx
      if(iband1>l) then
        iband1 = l
        ichang = -ichang
        do i=1,m
          store = x(i)
          x(i) = y(i)
          y(i) = store
        end do
        store = x0
        x0 = y0
        y0 = store
        store = x1
        x1 = y1
        y1 = store
        n = min0(nx,ny)
        do i=1,n
          store = tx(i)
          tx(i) = ty(i)
          ty(i) = store
        end do
        n1 = n+1
        if(nx<ny) then
          do i=n1,ny
            tx(i) = ty(i)
          end do
        else if (nx>ny) then
          do i=n1,nx
            ty(i) = tx(i)
          end do
        end if
        l = nx
        nx = ny
        ny = l
        l = nxe
        nxe = nye
        nye = l
        l = nxx
        nxx = nyy
        nyy = l
        l = kx
        kx = ky
        ky = l
        kx1 = kx+1
        ky1 = ky+1
      end if
      iband = iband1+1
!  arrange the data points according to the panel they belong to.
      call fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,idx,nreg)
!  find ncof, the number of b-spline coefficients.
      nk1x = nx-kx1
      nk1y = ny-ky1
      ncof = nk1x*nk1y
!  initialize the observation matrix a.
      do i=1,ncof
        f(i) = 0.
        do j=1,iband
          a(i,j) = 0.
        end do
      end do
!  initialize the sum of squared residuals.
      fp = 0.
!  fetch the data points in the new order. main loop for the
!  different panels.
      do num=1,nreg
!  fix certain constants for the current panel; jrot records the column
!  number of the first non-zero element in a row of the observation
!  matrix according to a data point of the panel.
        num1 = num-1
        lx = num1/nyy
        l1 = lx+kx1
        ly = num1-lx*nyy
        l2 = ly+ky1
        jrot = lx*nk1y+ly
!  test whether there are still data points in the panel.
        in = idx(num)
        do
          if(in==0) exit
!  fetch a new data point.
          wi = w(in)
          zi = z(in)*wi
!  evaluate for the x-direction, the (kx+1) non-zero b-splines at x(in).
          call fpbspl(tx,nx,kx,x(in),l1,hx)
!  evaluate for the y-direction, the (ky+1) non-zero b-splines at y(in).
          call fpbspl(ty,ny,ky,y(in),l2,hy)
!  store the value of these b-splines in spx and spy respectively.
          do i=1,kx1
            spx(in,i) = hx(i)
          end do
          do i=1,ky1
            spy(in,i) = hy(i)
          end do
!  initialize the new row of observation matrix.
          do i=1,iband
            h(i) = 0.
          end do
!  calculate the non-zero elements of the new row by making the cross
!  products of the non-zero b-splines in x- and y-direction.
          i1 = 0
          do i=1,kx1
            hxi = hx(i)
            j1 = i1
            do j=1,ky1
              j1 = j1+1
              h(j1) = hxi*hy(j)*wi
            end do
            i1 = i1+nk1y
          end do
!  rotate the row into triangle by givens transformations .
          irot = jrot
          do i=1,iband
            irot = irot+1
            piv = h(i)
            if(piv==0.) cycle
!  calculate the parameters of the givens transformation.
            call fpgivs(piv,a(irot,1),cs,sn)
!  apply that transformation to the right hand side.
            call fprota(cs,sn,zi,f(irot))
            if(i==iband) exit
!  apply that transformation to the left hand side.
            i2 = 1
            i3 = i+1
            do j=i3,iband
              i2 = i2+1
              call fprota(cs,sn,h(j),a(irot,i2))
            end do
          end do
!  add the contribution of the row to the sum of squares of residual
!  right hand sides.
          fp = fp+zi**2
!  find the number of the next data point in the panel.
          in = nummer(in)
        end do
      end do
!  find dmax, the maximum value for the diagonal elements in the reduced
!  triangle.
      dmax = 0.
      do i=1,ncof
        if(a(i,1)<=dmax) cycle
        dmax = a(i,1)
      end do
!  check whether the observation matrix is rank deficient.
      sigma = eps*dmax
      do i=1,ncof
        if(a(i,1)<=sigma) go to 280
      end do
!  backward substitution in case of full rank.
      call fpback(a,f,ncof,iband,c,nc)
      rank = ncof
      do i=1,ncof
        q(i,1) = a(i,1)/dmax
      end do
      go to 300
!  in case of rank deficiency, find the minimum norm solution.
!  check whether there is sufficient working space
      280 continue
      lwest = ncof*iband+ncof+iband
      if(lwrk<lwest) go to 780
      do i=1,ncof
        ff(i) = f(i)
        do j=1,iband
          q(i,j) = a(i,j)
        end do
      end do
      lf =1
      lh = lf+ncof
      la = lh+iband
      call fprank(q,ff,ncof,iband,nc,sigma,c,sq,rank,wrk(la), &
        wrk(lf),wrk(lh))
      do i=1,ncof
        q(i,1) = q(i,1)/dmax
      end do
!  add to the sum of squared residuals, the contribution of reducing
!  the rank.
      fp = fp+sq
      300 continue
      if(ier==(-2)) fp0 = fp
!  test whether the least-squares spline is an acceptable solution.
      if(iopt<0) go to 820
      fpms = fp-s
      if(abs(fpms)<=acc) then
        if(fp>0.0_wp_) then
          go to 820
        else
          go to 815
        end if
      end if
!  test whether we can accept the choice of knots.
      if(fpms<0.) exit
!  test whether we cannot further increase the number of knots.
      if(ncof>m) go to 790
      ier = 0
!  search where to add a new knot.
!  find for each interval the sum of squared residuals fpint for the
!  data points having the coordinate belonging to that knot interval.
!  calculate also coord which is the same sum, weighted by the position
!  of the data points considered.
      do i=1,nrint
        fpint(i) = 0.
        coord(i) = 0.
      end do
      do num=1,nreg
        num1 = num-1
        lx = num1/nyy
        l1 = lx+1
        ly = num1-lx*nyy
        l2 = ly+1+nxx
        jrot = lx*nk1y+ly
        in = idx(num)
        do
          if(in==0) exit
          store = 0.
          i1 = jrot
          do i=1,kx1
            hxi = spx(in,i)
            j1 = i1
            do j=1,ky1
              j1 = j1+1
              store = store+hxi*spy(in,j)*c(j1)
            end do
            i1 = i1+nk1y
          end do
          store = (w(in)*(z(in)-store))**2
          fpint(l1) = fpint(l1)+store
          coord(l1) = coord(l1)+store*x(in)
          fpint(l2) = fpint(l2)+store
          coord(l2) = coord(l2)+store*y(in)
          in = nummer(in)
        end do
      end do
!  find the interval for which fpint is maximal on the condition that
!  there still can be added a knot.
      do
        l = 0
        fpmax = 0.
        l1 = 1
        l2 = nrint
        if(nx==nxe) l1 = nxx+1
        if(ny==nye) l2 = nxx
        if(l1>l2) go to 810
        do i=l1,l2
          if(fpmax>=fpint(i)) cycle
          l = i
          fpmax = fpint(i)
        end do
!  test whether we cannot further increase the number of knots.
        if(l==0) go to 785
!  calculate the position of the new knot.
        arg = coord(l)/fpint(l)
!  test in what direction the new knot is going to be added.
        if(l<=nxx) then
!  addition in the x-direction.
          jxy = l+kx1
          fpint(l) = 0.
          fac1 = tx(jxy)-arg
          fac2 = arg-tx(jxy-1)
          if(fac1>(ten*fac2) .or. fac2>(ten*fac1)) cycle
          j = nx
          do i=jxy,nx
            tx(j+1) = tx(j)
            j = j-1
          end do
          tx(jxy) = arg
          nx = nx+1
        else
!  addition in the y-direction.
          jxy = l+ky1-nxx
          fpint(l) = 0.
          fac1 = ty(jxy)-arg
          fac2 = arg-ty(jxy-1)
          if(fac1>(ten*fac2) .or. fac2>(ten*fac1)) cycle
          j = ny
          do i=jxy,ny
            ty(j+1) = ty(j)
            j = j-1
          end do
          ty(jxy) = arg
          ny = ny+1
        end if
        exit
      end do
!  restart the computations with the new set of knots.
    end do
!  test whether the least-squares polynomial is a solution of our
!  approximation problem.
    if(ier==(-2)) go to 830
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! part 2: determination of the smoothing spline sp(x,y)                c
! *****************************************************                c
! we have determined the number of knots and their position. we now    c
! compute the b-spline coefficients of the smoothing spline sp(x,y).   c
! the observation matrix a is extended by the rows of a matrix,        c
! expressing that sp(x,y) must be a polynomial of degree kx in x and   c
! ky in y. the corresponding weights of these additional rows are set  c
! to 1./p.  iteratively we than have to determine the value of p       c
! such that f(p)=sum((w(i)*(z(i)-sp(x(i),y(i))))**2) be = s.           c
! we already know that the least-squares polynomial corresponds to     c
! p=0  and that the least-squares spline corresponds to p=infinity.    c
! the iteration process which is proposed here makes use of rational   c
! interpolation. since f(p) is a convex and strictly decreasing        c
! function of p, it can be approximated by a rational function r(p)=   c
! (u*p+v)/(p+w). three values of p(p1,p2,p3) with corresponding values c
! of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used to calculate the c
! new value of p such that r(p)=s. convergence is guaranteed by taking c
! f1 > 0 and f3 < 0.                                                   c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    kx2 = kx1+1
!  test whether there are interior knots in the x-direction.
!  and
!  evaluate the discotinuity jumps of the kx-th order derivative of
!  the b-splines at the knots tx(l),l=kx+2,...,nx-kx-1.
    if(nk1x/=kx1) call fpdisc(tx,nx,kx2,bx,nmax) 
    ky2 = ky1 + 1
!  test whether there are interior knots in the y-direction.
!  and
!  evaluate the discontinuity jumps of the ky-th order derivative of
!  the b-splines at the knots ty(l),l=ky+2,...,ny-ky-1.
    if(nk1y/=ky1) call fpdisc(ty,ny,ky2,by,nmax)
!  initial value for p.
    p1 = 0.
    f1 = fp0-s
    p3 = -one
    f3 = fpms
    p = 0.
    do i=1,ncof
      p = p+a(i,1)
    end do
    rn = ncof
    p = rn/p
!  find the bandwidth of the extended observation matrix.
    iband3 = kx1*nk1y
    iband4 = iband3 +1
    ich1 = 0
    ich3 = 0
!  iteration process to find the root of f(p)=s.
    do 770 iter=1,maxit
      pinv = one/p
!  store the triangularized observation matrix into q.
      do i=1,ncof
        ff(i) = f(i)
        do j=1,iband
          q(i,j) = a(i,j)
        end do
        ibb = iband+1
        do j=ibb,iband4
          q(i,j) = 0.
        end do
      end do
      if(nk1y/=ky1) then
!  extend the observation matrix with the rows of a matrix, expressing
!  that for x=cst. sp(x,y) must be a polynomial in y of degree ky.
        do i=ky2,nk1y
          ii = i-ky1
          do j=1,nk1x
!  initialize the new row.
            do l=1,iband
              h(l) = 0.
            end do
!  fill in the non-zero elements of the row. jrot records the column
!  number of the first non-zero element in the row.
            do l=1,ky2
              h(l) = by(ii,l)*pinv
            end do
            zi = 0.
            jrot = (j-1)*nk1y+ii
!  rotate the new row into triangle by givens transformations without
!  square roots.
            do irot=jrot,ncof
              piv = h(1)
              i2 = min0(iband1,ncof-irot)
              if(piv/=0.) then
!  calculate the parameters of the givens transformation.
                call fpgivs(piv,q(irot,1),cs,sn)
!  apply that givens transformation to the right hand side.
                call fprota(cs,sn,zi,ff(irot))
                if(i2==0) exit
!  apply that givens transformation to the left hand side.
                do l=1,i2
                  l1 = l+1
                  call fprota(cs,sn,h(l1),q(irot,l1))
                end do
              else
                if(i2<=0) exit
              end if
              do l=1,i2
                h(l) = h(l+1)
              end do
              h(i2+1) = 0.
            end do
          end do
        end do
      end if
      if(nk1x/=kx1) then
!  extend the observation matrix with the rows of a matrix expressing
!  that for y=cst. sp(x,y) must be a polynomial in x of degree kx.
        do i=kx2,nk1x
          ii = i-kx1
          do j=1,nk1y
!  initialize the new row
            do l=1,iband4
              h(l) = 0.
            end do
!  fill in the non-zero elements of the row. jrot records the column
!  number of the first non-zero element in the row.
            j1 = 1
            do l=1,kx2
              h(j1) = bx(ii,l)*pinv
              j1 = j1+nk1y
            end do
            zi = 0.
            jrot = (i-kx2)*nk1y+j
!  rotate the new row into triangle by givens transformations .
            do irot=jrot,ncof
              piv = h(1)
              i2 = min0(iband3,ncof-irot)
              if(piv/=0.) then
!  calculate the parameters of the givens transformation.
                call fpgivs(piv,q(irot,1),cs,sn)
!  apply that givens transformation to the right hand side.
                call fprota(cs,sn,zi,ff(irot))
                if(i2==0) exit
!  apply that givens transformation to the left hand side.
                do l=1,i2
                  l1 = l+1
                  call fprota(cs,sn,h(l1),q(irot,l1))
                end do
              else
                if(i2<=0) exit
              end if
              do l=1,i2
                h(l) = h(l+1)
              end do
              h(i2+1) = 0.
            end do
          end do
        end do
      end if
!  find dmax, the maximum value for the diagonal elements in the
!  reduced triangle.
      dmax = 0.
      do i=1,ncof
        if(q(i,1)<=dmax) cycle
        dmax = q(i,1)
      end do
!  check whether the matrix is rank deficient.
      sigma = eps*dmax
      do i=1,ncof
        if(q(i,1)<=sigma) go to 670
      end do
!  backward substitution in case of full rank.
      call fpback(q,ff,ncof,iband4,c,nc)
      rank = ncof
      go to 675
!  in case of rank deficiency, find the minimum norm solution.
      670 continue
      lwest = ncof*iband4+ncof+iband4
      if(lwrk<lwest) go to 780
      lf = 1
      lh = lf+ncof
      la = lh+iband4
      call fprank(q,ff,ncof,iband4,nc,sigma,c,sq,rank,wrk(la), &
       wrk(lf),wrk(lh))
      675 continue
      do i=1,ncof
        q(i,1) = q(i,1)/dmax
      end do
!  compute f(p).
      fp = 0.
      do num = 1,nreg
        num1 = num-1
        lx = num1/nyy
        ly = num1-lx*nyy
        jrot = lx*nk1y+ly
        in = idx(num)
        do
          if(in==0) exit
          store = 0.
          i1 = jrot
          do i=1,kx1
            hxi = spx(in,i)
            j1 = i1
            do j=1,ky1
             j1 = j1+1
              store = store+hxi*spy(in,j)*c(j1)
            end do
            i1 = i1+nk1y
          end do
          fp = fp+(w(in)*(z(in)-store))**2
          in = nummer(in)
        end do
      end do
!  test whether the approximation sp(x,y) is an acceptable solution.
      fpms = fp-s
      if(abs(fpms)<=acc) go to 820
!  test whether the maximum allowable number of iterations has been
!  reached.
      if(iter==maxit) go to 795
!  carry out one more step of the iteration process.
      p2 = p
      f2 = fpms
      if(ich3/=0) go to 740
      if((f2-f3)>acc) go to 730
!  our initial choice of p is too large.
      p3 = p2
      f3 = f2
      p = p*con4
      if(p<=p1) p = p1*con9 + p2*con1
      go to 770
   730    if(f2<0.) ich3 = 1
   740    if(ich1/=0) go to 760
      if((f1-f2)>acc) go to 750
!  our initial choice of p is too small
      p1 = p2
      f1 = f2
      p = p/con4
      if(p3<0.) go to 770
      if(p>=p3) p = p2*con1 + p3*con9
      go to 770
   750    if(f2>0.) ich1 = 1
!  test whether the iteration process proceeds as theoretically
!  expected.
   760    if(f2>=f1 .or. f2<=f3) go to 800
!  find the new value of p.
      p = fprati(p1,f1,p2,f2,p3,f3)
   770  continue
!  error codes and messages.
   780  ier = lwest
    go to 830
   785  ier = 5
    go to 830
   790  ier = 4
    go to 830
   795  ier = 3
    go to 830
   800  ier = 2
    go to 830
   810  ier = 1
    go to 830
   815  ier = -1
    fp = 0.
   820  if(ncof/=rank) ier = -rank
!  test whether x and y are in the original order.
   830  if(ichang<0) go to 930
!  if not, interchange x and y once more.
    l1 = 1
    do i=1,nk1x
      l2 = i
      do j=1,nk1y
        f(l2) = c(l1)
        l1 = l1+1
        l2 = l2+nk1x
      end do
    end do
    do i=1,ncof
      c(i) = f(i)
    end do
    do i=1,m
      store = x(i)
      x(i) = y(i)
      y(i) = store
    end do
    n = min0(nx,ny)
    do i=1,n
      store = tx(i)
      tx(i) = ty(i)
      ty(i) = store
    end do
    n1 = n+1
    if(nx<ny) then
      do i=n1,ny
        tx(i) = ty(i)
      end do
    else if(nx>ny) then
      do i=n1,nx
        ty(i) = tx(i)
      end do
    end if
    l = nx
    nx = ny
    ny = l
    930 continue
    if(iopt>=0) then
      nx0 = nx
      ny0 = ny
    end if
  end subroutine fpsurf

  subroutine fpback(a,z,n,k,c,nest)
!  subroutine fpback calculates the solution of the system of
!  equations a*c = z with a a n x n upper triangular matrix
!  of bandwidth k.
!  ..
    implicit none
! arguments
    integer, intent(in) :: n, k, nest
    real(wp_), intent(in) :: a(nest,k), z(n)
    real(wp_), intent(inout) :: c(n)
! local variables
    real(wp_) :: store
    integer :: i, i1, j, k1, l, m
!  ..
    k1 = k-1
    c(n) = z(n)/a(n,1)
    i = n-1
    if(i==0) return
    do j=2,n
      store = z(i)
      i1 = k1
      if(j<=k1) i1 = j-1
      m = i
      do l=1,i1
        m = m+1
        store = store-c(m)*a(i,l+1)
      end do
      c(i) = store/a(i,1)
      i = i-1
    end do
  end subroutine fpback

  subroutine fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wx,wy,lx,ly)
    implicit none
! arguments
    integer, intent(in) :: nx, ny, kx, ky, mx, my
    integer, intent(out) :: lx(mx), ly(my)
    real(wp_), intent(in) :: tx(nx), ty(ny), c((nx-kx-1)*(ny-ky-1)), &
        x(mx), y(my)
    real(wp_), intent(out) :: z(mx*my)
    real(wp_), intent(out) :: wx(mx,kx+1), wy(my,ky+1)
! local variables
    integer :: kx1, ky1, l, l1, l2, m, nkx1, nky1, i, j, i1, j1
    real(wp_) :: arg, sp, tb, te, h(6)
!  ..subroutine references..
!    fpbspl
!  ..
    kx1 = kx+1
    nkx1 = nx-kx1
    tb = tx(kx1)
    te = tx(nkx1+1)
    l = kx1
    l1 = l+1
    do i=1,mx
      arg = x(i)
      if(arg<tb) arg = tb
      if(arg>te) arg = te
      do
        if(arg<tx(l1) .or. l==nkx1) exit
        l = l1
        l1 = l+1
      end do
      call fpbspl(tx,nx,kx,arg,l,h)
      lx(i) = l-kx1
      do j=1,kx1
        wx(i,j) = h(j)
      end do
    end do
    ky1 = ky+1
    nky1 = ny-ky1
    tb = ty(ky1)
    te = ty(nky1+1)
    l = ky1
    l1 = l+1
    do i=1,my
      arg = y(i)
      if(arg<tb) arg = tb
      if(arg>te) arg = te
      do
        if(arg<ty(l1) .or. l==nky1) exit
        l = l1
        l1 = l+1
      end do
      call fpbspl(ty,ny,ky,arg,l,h)
      ly(i) = l-ky1
      do j=1,ky1
        wy(i,j) = h(j)
      end do
    end do
    m = 0
    do i=1,mx
      l = lx(i)*nky1
      do i1=1,kx1
        h(i1) = wx(i,i1)
      end do
      do j=1,my
        l1 = l+ly(j)
        sp = 0.0_wp_
        do i1=1,kx1
          l2 = l1
          do j1=1,ky1
            l2 = l2+1
            sp = sp+c(l2)*h(i1)*wy(j,j1)
          end do
          l1 = l1+nky1
        end do
        m = m+1
        z(m) = sp
      end do
    end do
  end subroutine fpbisp

  subroutine fpbspl(t,n,k,x,l,h)
!  subroutine fpbspl evaluates the (k+1) non-zero b-splines of
!  degree k at t(l) <= x < t(l+1) using the stable recurrence
!  relation of de boor and cox.
!  ..
    implicit none
! arguments
    integer, intent(in) :: n, k, l
    real(wp_), intent(in) :: x, t(n)
    real(wp_), intent(out) :: h(6)
! local variables
    integer :: i, j, li, lj
    real(wp_) :: f, hh(5)
! parameters
    real(wp_), parameter :: one = 0.1e1_wp_
!  ..
    h(1) = one
    do j=1,k
      do i=1,j
        hh(i) = h(i)
      end do
      h(1) = 0.0_wp_
      do i=1,j
        li = l+i
        lj = li-j
        f = hh(i)/(t(li)-t(lj))
        h(i) = h(i)+f*(t(li)-x)
        h(i+1) = f*(x-t(lj))
      end do
    end do
  end subroutine fpbspl

  subroutine fpchec(x,m,t,n,k,ier)
!  subroutine fpchec verifies the number and the position of the knots
!  t(j),j=1,2,...,n of a spline of degree k, in relation to the number
!  and the position of the data points x(i),i=1,2,...,m. if all of the
!  following conditions are fulfilled, the error parameter ier is set
!  to zero. if one of the conditions is violated ier is set to ten.
!  1) k+1 <= n-k-1 <= m
!  2) t(1) <= t(2) <= ... <= t(k+1)
!     t(n-k) <= t(n-k+1) <= ... <= t(n)
!  3) t(k+1) < t(k+2) < ... < t(n-k)
!  4) t(k+1) <= x(i) <= t(n-k)
!  5) the conditions specified by schoenberg and whitney must hold
!     for at least one subset of data points, i.e. there must be a
!     subset of data points y(j) such that
!         t(j) < y(j) < t(j+k+1), j=1,2,...,n-k-1
!  ..
    implicit none
! arguments
    integer, intent(in) :: m, n, k
    real(wp_), intent(in) :: x(m), t(n)
    integer, intent(out) :: ier
! local variables
    integer :: i, j, k1, k2, l, nk1, nk2, nk3
    real(wp_) :: tj, tl
!  ..
    k1 = k+1
    k2 = k1+1
    nk1 = n-k1
    nk2 = nk1+1
    ier = 10
!  check condition no 1
    if(nk1<k1 .or. nk1>m) return
!  check condition no 2
    j = n
    do i=1,k
      if(t(i)>t(i+1) .or. t(j)<t(j-1)) return
      j = j-1
    end do
!  check condition no 3
    do i=k2,nk2
      if(t(i)<=t(i-1)) return
    end do
!  check condition no 4
    if(x(1)<t(k1) .or. x(m)>t(nk2)) return
!  check condition no 5
    if(x(1)>=t(k2) .or. x(m)<=t(nk1)) return
    i = 1
    l = k2
    nk3 = nk1-1
    do j=2,nk3
      tj = t(j)
      l = l+1
      tl = t(l)
      do
        i = i+1
        if(i>=m) return
        if(x(i)>tj) exit
      end do
      if(x(i)>=tl) return
    end do
    ier = 0
  end subroutine fpchec

  subroutine fpdisc(t,n,k2,b,nest)
!  subroutine fpdisc calculates the discontinuity jumps of the kth
!  derivative of the b-splines of degree k at the knots t(k+2)..t(n-k-1)
    implicit none
! arguments
    integer, intent(in) :: n, k2, nest
    real(wp_), intent(in) :: t(n)
    real(wp_), intent(out) :: b(nest,k2)
! local variables
    real(wp_) :: an, fac, prod
    integer :: i, ik, j, jk, k, k1, l, lj, lk, lmk, lp, nk1, nrint
    real(wp_), dimension(12) :: h
!  ..
    k1 = k2-1
    k = k1-1
    nk1 = n-k1
    nrint = nk1-k
    an = nrint
    fac = an/(t(nk1+1)-t(k1))
    do l=k2,nk1
      lmk = l-k1
      do j=1,k1
        ik = j+k1
        lj = l+j
        lk = lj-k2
        h(j) = t(l)-t(lk)
        h(ik) = t(l)-t(lj)
      end do
      lp = lmk
      do j=1,k2
        jk = j
        prod = h(j)
        do i=1,k
          jk = jk+1
          prod = prod*h(jk)*fac
        end do
        lk = lp+k1
        b(lmk,j) = (t(lk)-t(lp))/prod
        lp = lp+1
      end do
    end do
  end subroutine fpdisc

  subroutine fprank(a,f,n,m,na,tol,c,sq,rank,aa,ff,h)
!  subroutine fprank finds the minimum norm solution of a least-
!  squares problem in case of rank deficiency.
!
!  input parameters:
!    a : array, which contains the non-zero elements of the observation
!    matrix after triangularization by givens transformations.
!    f : array, which contains the transformed right hand side.
!    n : integer,wich contains the dimension of a.
!    m : integer, which denotes the bandwidth of a.
!  tol : real value, giving a threshold to determine the rank of a.
!
!  output parameters:
!    c : array, which contains the minimum norm solution.
!   sq : real value, giving the contribution of reducing the rank
!    to the sum of squared residuals.
! rank : integer, which contains the rank of matrix a.
!
    implicit none
!  ..scalar arguments..
    integer n,m,na,rank
    real(wp_) tol,sq
!  ..array arguments..
    real(wp_) a(na,m),f(n),c(n),aa(n,m),ff(n),h(m)
!  ..local scalars..
    integer i,ii,ij,i1,i2,j,jj,j1,j2,j3,k,kk,m1,nl
    real(wp_) cs,fac,piv,sn,yi
    real(wp_) store,stor1,stor2,stor3
!  ..function references..
    integer min0
!  ..subroutine references..
!    fpgivs,fprota
!  ..
    m1 = m-1
!  the rank deficiency nl is considered to be the number of sufficient
!  small diagonal elements of a.
    nl = 0
    sq = 0.
    do i=1,n
      if(a(i,1)>tol) cycle
!  if a sufficient small diagonal element is found, we put it to
!  zero. the remainder of the row corresponding to that zero diagonal
!  element is then rotated into triangle by givens rotations .
!  the rank deficiency is increased by one.
      nl = nl+1
      if(i==n) cycle
      yi = f(i)
      do j=1,m1
        h(j) = a(i,j+1)
      end do
      h(m) = 0.
      i1 = i+1
      do ii=i1,n
        i2 = min0(n-ii,m1)
        piv = h(1)
        if(piv/=0.) then
          call fpgivs(piv,a(ii,1),cs,sn)
          call fprota(cs,sn,yi,f(ii))
          if(i2==0) exit
          do j=1,i2
            j1 = j+1
            call fprota(cs,sn,h(j1),a(ii,j1))
            h(j) = h(j1)
          end do
        else
          if(i2==0) exit
          do j=1,i2
            h(j) = h(j+1)
          end do
        end if
        h(i2+1) = 0.
      end do
!  add to the sum of squared residuals the contribution of deleting
!  the row with small diagonal element.
      sq = sq+yi**2
    end do
!  rank denotes the rank of a.
    rank = n-nl
!  let b denote the (rank*n) upper trapezoidal matrix which can be
!  obtained from the (n*n) upper triangular matrix a by deleting
!  the rows and interchanging the columns corresponding to a zero
!  diagonal element. if this matrix is factorized using givens
!  transformations as  b = (r) (u)  where
!    r is a (rank*rank) upper triangular matrix,
!    u is a (rank*n) orthonormal matrix
!  then the minimal least-squares solution c is given by c = b' v,
!  where v is the solution of the system  (r) (r)' v = g  and
!  g denotes the vector obtained from the old right hand side f, by
!  removing the elements corresponding to a zero diagonal element of a.
!  initialization.
    do i=1,rank
      do j=1,m
        aa(i,j) = 0.
      end do
    end do
!  form in aa the upper triangular matrix obtained from a by
!  removing rows and columns with zero diagonal elements. form in ff
!  the new right hand side by removing the elements of the old right
!  hand side corresponding to a deleted row.
    ii = 0
    do i=1,n
      if(a(i,1)<=tol) cycle
      ii = ii+1
      ff(ii) = f(i)
      aa(ii,1) = a(i,1)
      jj = ii
      kk = 1
      j = i
      j1 = min0(j-1,m1)
      do k=1,j1
        j = j-1
        if(a(j,1)<=tol) cycle
        kk = kk+1
        jj = jj-1
        aa(jj,kk) = a(j,k+1)
      end do
    end do
!  form successively in h the columns of a with a zero diagonal element.
    ii = 0
    do i=1,n
      ii = ii+1
      if(a(i,1)>tol) cycle
      ii = ii-1
      if(ii==0) cycle
      jj = 1
      j = i
      j1 = min0(j-1,m1)
      do k=1,j1
        j = j-1
        if(a(j,1)<=tol) cycle
        h(jj) = a(j,k+1)
        jj = jj+1
      end do
      do kk=jj,m
        h(kk) = 0.
      end do
!  rotate this column into aa by givens transformations.
      jj = ii
      do i1=1,ii
        j1 = min0(jj-1,m1)
        piv = h(1)
        if(piv==0.) then
          if(j1==0) exit
          do j2=1,j1
            j3 = j2+1
            h(j2) = h(j3)
          end do
        else
          call fpgivs(piv,aa(jj,1),cs,sn)
          if(j1==0) exit
          kk = jj
          do j2=1,j1
            j3 = j2+1
            kk = kk-1
            call fprota(cs,sn,h(j3),aa(kk,j3))
            h(j2) = h(j3)
          end do
        end if
        jj = jj-1
        h(j3) = 0.
      end do
    end do
!  solve the system (aa) (f1) = ff
    ff(rank) = ff(rank)/aa(rank,1)
    i = rank-1
    do j=2,rank
      store = ff(i)
      i1 = min0(j-1,m1)
      k = i
      do ii=1,i1
        k = k+1
        stor1 = ff(k)
        stor2 = aa(i,ii+1)
        store = store-stor1*stor2
      end do
      stor1 = aa(i,1)
      ff(i) = store/stor1
      i = i-1
    end do
!  solve the system  (aa)' (f2) = f1
    ff(1) = ff(1)/aa(1,1)
    do j=2,rank
      store = ff(j)
      i1 = min0(j-1,m1)
      k = j
      do ii=1,i1
        k = k-1
        stor1 = ff(k)
        stor2 = aa(k,ii+1)
        store = store-stor1*stor2
      end do
      stor1 = aa(j,1)
      ff(j) = store/stor1
    end do
!  premultiply f2 by the transpoze of a.
    k = 0
    do i=1,n
      store = 0.
      if(a(i,1)>tol) k = k+1
      j1 = min0(i,m)
      kk = k
      ij = i+1
      do j=1,j1
        ij = ij-1
        if(a(ij,1)<=tol) cycle
        stor1 = a(ij,j)
        stor2 = ff(kk)
        store = store+stor1*stor2
        kk = kk-1
      end do
      c(i) = store
    end do
!  add to the sum of squared residuals the contribution of putting
!  to zero the small diagonal elements of matrix (a).
    stor3 = 0.
    do i=1,n
      if(a(i,1)>tol) cycle
      store = f(i)
      i1 = min0(n-i,m1)
      do j=1,i1
        ij = i+j
        stor1 = c(ij)
        stor2 = a(i,j+1)
        store = store-stor1*stor2
      end do
      fac = a(i,1)*c(i)
      stor1 = a(i,1)
      stor2 = c(i)
      stor1 = stor1*stor2
      stor3 = stor3+stor1*(stor1-store-store)
    end do
    fac = stor3
    sq = sq+fac
  end subroutine fprank

  subroutine fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,idx,nreg)
!  subroutine fporde sorts the data points (x(i),y(i)),i=1,2,...,m
!  according to the panel tx(l)<=x<tx(l+1),ty(k)<=y<ty(k+1), they belong
!  to. for each panel a stack is constructed  containing the numbers
!  of data points lying inside; index(j),j=1,2,...,nreg points to the
!  first data point in the jth panel while nummer(i),i=1,2,...,m gives
!  the number of the next data point in the panel.
!  ..
    implicit none
!  ..scalar arguments..
    integer m,kx,ky,nx,ny,nreg
!  ..array arguments..
    real(wp_) x(m),y(m),tx(nx),ty(ny)
    integer nummer(m),idx(nreg)
!  ..local scalars..
    real(wp_) xi,yi
    integer i,im,k,kx1,ky1,k1,l,l1,nk1x,nk1y,num,nyy
!  ..
    kx1 = kx+1
    ky1 = ky+1
    nk1x = nx-kx1
    nk1y = ny-ky1
    nyy = nk1y-ky
    do i=1,nreg
      idx(i) = 0
    end do
    do im=1,m
      xi = x(im)
      yi = y(im)
      l = kx1
      l1 = l+1
      do
        if(xi<tx(l1) .or. l==nk1x) exit
        l = l1
        l1 = l+1
      end do
      k = ky1
      k1 = k+1
      do
        if(yi<ty(k1) .or. k==nk1y) exit
        k = k1
        k1 = k+1
      end do
      num = (l-kx1)*nyy+k-ky
      nummer(im) = idx(num)
      idx(num) = im
    end do
  end subroutine fporde

  subroutine fpgivs(piv,ww,cs,sn)
!  subroutine fpgivs calculates the parameters of a givens
!  transformation .
!  ..
    implicit none
! arguments
    real(wp_), intent(in)    :: piv
    real(wp_), intent(out)   :: cs, sn
    real(wp_), intent(inout) :: ww
! local variables
    real(wp_) :: dd, store
! parameters
    real(wp_), parameter :: one = 0.1e+01_wp_
!  ..
    store = abs(piv)
    if(store>=ww) dd = store*sqrt(one+(ww/piv)**2)
    if(store<ww) dd = ww*sqrt(one+(piv/ww)**2)
    cs = ww/dd
    sn = piv/dd
    ww = dd
  end subroutine fpgivs

  subroutine fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx, &
     ty,ny,p,c,nc,fp,fpx,fpy,mm,mynx,kx1,kx2,ky1,ky2,spx,spy,right,q, &
     ax,ay,bx,by,nrx,nry)
!  ..
    implicit none
! arguments
    integer, intent(in) :: kx, kx1, kx2, ky, ky1, ky2, mm, mx, my, mz, &
                           mynx, nx, ny, nc
    integer, intent(inout) :: ifsx, ifsy, ifbx, ifby, nrx(mx), nry(my)
    real(wp_), intent(in) :: p, x(mx), y(my), z(mz), tx(nx), ty(ny)
    real(wp_), intent(out) :: fp, c(nc), q(mynx), ax(nx,kx2), ay(ny,ky2), fpx(nx), fpy(ny)
    real(wp_), intent(inout) :: spx(mx,kx1), spy(my,ky1), bx(nx,kx2), by(ny,ky2)
    real(wp_), intent(out) :: right(mm)
! local variables
    real(wp_) :: arg, cs, fac, pinv, piv, sn, term
    integer :: i, ibandx, ibandy, ic, iq, irot, it, iz,i1, i2, i3,  &
     j, k, k1, k2, l, l1, l2, ncof, nk1x, nk1y, nrold, nroldx, nroldy,    &
     number, numx, numx1, numy, numy1, n1
    real(wp_), dimension(7) :: h
! parameters
    real(wp_), parameter :: one = 1.0_wp_, half = 0.5_wp_
!  ..subroutine references..
!    fpback,fpbspl,fpgivs,fpdisc,fprota
!  ..
!  the b-spline coefficients of the smoothing spline are calculated as
!  the least-squares solution of the over-determined linear system of
!  equations  (ay) c (ax)' = q       where
!
!           |   (spx)    |            |   (spy)    |
!    (ax) = | ---------- |     (ay) = | ---------- |
!           | (1/p) (bx) |            | (1/p) (by) |
!
!                            | z  ' 0 |
!                        q = | ------ |
!                            | 0  ' 0 |
!
!  with c      : the (ny-ky-1) x (nx-kx-1) matrix which contains the
!            b-spline coefficients.
!   z      : the my x mx matrix which contains the function values.
!   spx,spy: the mx x (nx-kx-1) and  my x (ny-ky-1) observation
!            matrices according to the least-squares problems in
!            the x- and y-direction.
!   bx,by  : the (nx-2*kx-1) x (nx-kx-1) and (ny-2*ky-1) x (ny-ky-1)
!            matrices which contain the discontinuity jumps of the
!            derivatives of the b-splines in the x- and y-direction.
    nk1x = nx-kx1
    nk1y = ny-ky1
    if(p>0.0_wp_) pinv = one/p
!  it depends on the value of the flags ifsx,ifsy,ifbx and ifby and on
!  the value of p whether the matrices (spx),(spy),(bx) and (by) still
!  must be determined.
    if(ifsx==0) then
!  calculate the non-zero elements of the matrix (spx) which is the
!  observation matrix according to the least-squares spline approximat-
!  ion problem in the x-direction.
      l = kx1
      l1 = kx2
      number = 0
      do it=1,mx
        arg = x(it)
        do
          if(arg<tx(l1) .or. l==nk1x) exit
          l = l1
          l1 = l+1
          number = number+1
        end do
        call fpbspl(tx,nx,kx,arg,l,h)
        do i=1,kx1
          spx(it,i) = h(i)
        end do
        nrx(it) = number
      end do
      ifsx = 1
    end if
    if(ifsy==0) then
!  calculate the non-zero elements of the matrix (spy) which is the
!  observation matrix according to the least-squares spline approximat-
!  ion problem in the y-direction.
      l = ky1
      l1 = ky2
      number = 0
      do it=1,my
        arg = y(it)
        do
          if(arg<ty(l1) .or. l==nk1y) exit
          l = l1
          l1 = l+1
          number = number+1
        end do
        call fpbspl(ty,ny,ky,arg,l,h)
        do i=1,ky1
          spy(it,i) = h(i)
        end do
        nry(it) = number
      end do
      ifsy = 1
    end if
    if(p>0.0_wp_) then
!  calculate the non-zero elements of the matrix (bx).
      if(ifbx==0 .and. nx/=2*kx1) then
        call fpdisc(tx,nx,kx2,bx,nx)
        ifbx = 1
      end if
!  calculate the non-zero elements of the matrix (by).
      if(ifby==0 .and. ny/=2*ky1) then
        call fpdisc(ty,ny,ky2,by,ny)
        ifby = 1
      end if
    end if
!  reduce the matrix (ax) to upper triangular form (rx) using givens
!  rotations. apply the same transformations to the rows of matrix q
!  to obtain the my x (nx-kx-1) matrix g.
!  store matrix (rx) into (ax) and g into q.
    l = my*nk1x
!  initialization.
    do i=1,l
      q(i) = 0.0_wp_
    end do
    do i=1,nk1x
      do j=1,kx2
        ax(i,j) = 0.0_wp_
      end do
    end do
    l = 0
    nrold = 0
!  ibandx denotes the bandwidth of the matrices (ax) and (rx).
    ibandx = kx1
    do it=1,mx
      number = nrx(it)
      do
        if(nrold/=number) then
          if(p<=0.0_wp_) then
            nrold = nrold+1
            cycle
          end if
          ibandx = kx2
!  fetch a new row of matrix (bx).
          n1 = nrold+1
          do j=1,kx2
            h(j) = bx(n1,j)*pinv
          end do
!  find the appropriate column of q.
          do j=1,my
            right(j) = 0.0_wp_
          end do
          irot = nrold
        else
!  fetch a new row of matrix (spx).
          h(ibandx) = 0.0_wp_
          do j=1,kx1
            h(j) = spx(it,j)
          end do
!  find the appropriate column of q.
          do j=1,my
            l = l+1
            right(j) = z(l)
          end do
          irot = number
        end if
!  rotate the new row of matrix (ax) into triangle.
        do i=1,ibandx
          irot = irot+1
          piv = h(i)
          if(piv==0.0_wp_) cycle
!  calculate the parameters of the givens transformation.
          call fpgivs(piv,ax(irot,1),cs,sn)
!  apply that transformation to the rows of matrix q.
          iq = (irot-1)*my
          do j=1,my
            iq = iq+1
            call fprota(cs,sn,right(j),q(iq))
          end do
!  apply that transformation to the columns of (ax).
          if(i==ibandx) exit
          i2 = 1
          i3 = i+1
          do j=i3,ibandx
            i2 = i2+1
            call fprota(cs,sn,h(j),ax(irot,i2))
          end do
        end do
        if(nrold==number) exit
        nrold = nrold+1
      end do
    end do
!  reduce the matrix (ay) to upper triangular form (ry) using givens
!  rotations. apply the same transformations to the columns of matrix g
!  to obtain the (ny-ky-1) x (nx-kx-1) matrix h.
!  store matrix (ry) into (ay) and h into c.
    ncof = nk1x*nk1y
!  initialization.
    do i=1,ncof
      c(i) = 0.0_wp_
    end do
    do i=1,nk1y
      do j=1,ky2
        ay(i,j) = 0.0_wp_
      end do
    end do
    nrold = 0
!  ibandy denotes the bandwidth of the matrices (ay) and (ry).
    ibandy = ky1
    do it=1,my
      number = nry(it)
      do
        if(nrold/=number) then
          if(p<=0.0_wp_) then
            nrold = nrold+1
            cycle
          end if
          ibandy = ky2
!  fetch a new row of matrix (by).
          n1 = nrold+1
          do j=1,ky2
            h(j) = by(n1,j)*pinv
          end do
!  find the appropiate row of g.
          do j=1,nk1x
            right(j) = 0.0_wp_
          end do
          irot = nrold
!  fetch a new row of matrix (spy)
        else
          h(ibandy) = 0.0_wp_
          do j=1,ky1
            h(j) = spy(it,j)
          end do
!  find the appropiate row of g.
          l = it
          do j=1,nk1x
            right(j) = q(l)
            l = l+my
          end do
          irot = number
        end if
!  rotate the new row of matrix (ay) into triangle.
        do i=1,ibandy
          irot = irot+1
          piv = h(i)
          if(piv==0.0_wp_) cycle
!  calculate the parameters of the givens transformation.
          call fpgivs(piv,ay(irot,1),cs,sn)
!  apply that transformation to the colums of matrix g.
          ic = irot
          do j=1,nk1x
            call fprota(cs,sn,right(j),c(ic))
            ic = ic+nk1y
          end do
!  apply that transformation to the columns of matrix (ay).
          if(i==ibandy) exit
          i2 = 1
          i3 = i+1
          do j=i3,ibandy
            i2 = i2+1
            call fprota(cs,sn,h(j),ay(irot,i2))
          end do
        end do
        if(nrold==number) exit
        nrold = nrold+1
      end do
    end do
!  backward substitution to obtain the b-spline coefficients as the
!  solution of the linear system    (ry) c (rx)' = h.
!  first step: solve the system  (ry) (c1) = h.
    k = 1
    do i=1,nk1x
      call fpback(ay,c(k),nk1y,ibandy,c(k),ny)
      k = k+nk1y
    end do
!  second step: solve the system  c (rx)' = (c1).
    k = 0
    do j=1,nk1y
      k = k+1
      l = k
      do i=1,nk1x
        right(i) = c(l)
        l = l+nk1y
      end do
      call fpback(ax,right,nk1x,ibandx,right,nx)
      l = k
      do i=1,nk1x
        c(l) = right(i)
        l = l+nk1y
      end do
    end do
!  calculate the quantities
!    res(i,j) = (z(i,j) - s(x(i),y(j)))**2 , i=1,2,..,mx;j=1,2,..,my
!    fp = sumi=1,mx(sumj=1,my(res(i,j)))
!    fpx(r) = sum''i(sumj=1,my(res(i,j))) , r=1,2,...,nx-2*kx-1
!              tx(r+kx) <= x(i) <= tx(r+kx+1)
!    fpy(r) = sumi=1,mx(sum''j(res(i,j))) , r=1,2,...,ny-2*ky-1
!              ty(r+ky) <= y(j) <= ty(r+ky+1)
    fp = 0.0_wp_
    do i=1,nx
      fpx(i) = 0.0_wp_
    end do
    do i=1,ny
      fpy(i) = 0.0_wp_
    end do
    nk1y = ny-ky1
    iz = 0
    nroldx = 0
!  main loop for the different grid points.
    do i1=1,mx
      numx = nrx(i1)
      numx1 = numx+1
      nroldy = 0
      do i2=1,my
        numy = nry(i2)
        numy1 = numy+1
        iz = iz+1
!  evaluate s(x,y) at the current grid point by making the sum of the
!  cross products of the non-zero b-splines at (x,y), multiplied with
!  the appropiate b-spline coefficients.
        term = 0.0_wp_
        k1 = numx*nk1y+numy
        do l1=1,kx1
          k2 = k1
          fac = spx(i1,l1)
          do l2=1,ky1
            k2 = k2+1
            term = term+fac*spy(i2,l2)*c(k2)
          end do
          k1 = k1+nk1y
        end do
!  calculate the squared residual at the current grid point.
        term = (z(iz)-term)**2
!  adjust the different parameters.
        fp = fp+term
        fpx(numx1) = fpx(numx1)+term
        fpy(numy1) = fpy(numy1)+term
        fac = term*half
        if(numy/=nroldy) then
          fpy(numy1) = fpy(numy1)-fac
          fpy(numy) = fpy(numy)+fac
        end if
        nroldy = numy
        if(numx/=nroldx) then
          fpx(numx1) = fpx(numx1)-fac
          fpx(numx) = fpx(numx)+fac
        end if
      end do
      nroldx = numx
    end do
  end subroutine fpgrre

  subroutine fpknot(x,m,t,n,fpint,nrdata,nrint,nest,istart)
!  subroutine fpknot locates an additional knot for a spline of degree
!  k and adjusts the corresponding parameters,i.e.
!    t     : the position of the knots.
!    n     : the number of knots.
!    nrint : the number of knotintervals.
!    fpint : the sum of squares of residual right hand sides
!        for each knot interval.
!    nrdata: the number of data points inside each knot interval.
!  istart indicates that the smallest data point at which the new knot
!  may be added is x(istart+1)
!  ..
    implicit none
! arguments
    integer, intent(in) :: m, nest, istart
    integer, intent(inout) :: n, nrint, nrdata(nest)
    real(wp_), intent(in) :: x(m)
    real(wp_), intent(inout) :: t(nest), fpint(nest)
! local variables
    real(wp_) :: an, am, fpmax
    integer :: ihalf, j, jbegin, jj, jk, jpoint, k, maxbeg, maxpt, next, nrx, number
!  ..
    k = (n-nrint-1)/2
!  search for knot interval t(number+k) <= x <= t(number+k+1) where
!  fpint(number) is maximal on the condition that nrdata(number)
!  not equals zero.
    fpmax = 0.0_wp_
    jbegin = istart
    do j=1,nrint
      jpoint = nrdata(j)
      if(fpmax<fpint(j) .and. jpoint/=0) then
        fpmax = fpint(j)
        number = j
        maxpt = jpoint
        maxbeg = jbegin
      end if
      jbegin = jbegin+jpoint+1
    end do
!  let coincide the new knot t(number+k+1) with a data point x(nrx)
!  inside the old knot interval t(number+k) <= x <= t(number+k+1).
    ihalf = maxpt/2+1
    nrx = maxbeg+ihalf
    next = number+1
    if(next<=nrint) then
!  adjust the different parameters.
      do j=next,nrint
        jj = next+nrint-j
        fpint(jj+1) = fpint(jj)
        nrdata(jj+1) = nrdata(jj)
        jk = jj+k
        t(jk+1) = t(jk)
      end do
    end if
    nrdata(number) = ihalf-1
    nrdata(next) = maxpt-ihalf
    am = maxpt
    an = nrdata(number)
    fpint(number) = fpmax*an/am
    an = nrdata(next)
    fpint(next) = fpmax*an/am
    jk = next+k
    t(jk) = x(nrx)
    n = n+1
    nrint = nrint+1
  end subroutine fpknot

  subroutine fpregr(iopt,x,mx,y,my,z,mz,xb,xe,yb,ye,kx,ky,s, &
     nxest,nyest,tol,maxit,nc,nx,tx,ny,ty,c,fp,fp0,fpold,reducx, &
     reducy,fpintx,fpinty,lastdi,nplusx,nplusy,nrx,nry,nrdatx,nrdaty, &
     wrk,lwrk,ier)
!  ..
    implicit none
! arguments
    integer, intent(in) :: kx,ky,mx,my,mz,nxest,nyest,maxit,nc,lwrk,iopt
    integer, intent(inout) :: nx,ny,lastdi,nplusx,nplusy, &
     nrdatx(nxest),nrdaty(nyest),nrx(mx),nry(my)
    integer, intent(out) :: ier
    real(wp_), intent(in) :: xb,xe,yb,ye,s,tol,x(mx),y(my),z(mz)
    real(wp_), intent(out) :: fp,c(nc)
    real(wp_), intent(inout) :: fp0,fpold,reducx,reducy,tx(nxest), &
     ty(nyest),fpintx(nxest),fpinty(nyest),wrk(lwrk)
! local variables
    real(wp_) :: acc,fpms,f1,f2,f3,p,p1,p2,p3,rn
    integer :: i,ich1,ich3,ifbx,ifby,ifsx,ifsy,iter,j,kx1,kx2,ky1,ky2,k3, &
     l,lax,lay,lbx,lby,lq,lri,lsx,lsy,mk1,mm,mpm,mynx,nmaxx,nmaxy,nminx,nminy, &
     nplx,nply,npl1,nrintx,nrinty,nxe,nxk,nye
! parameters
    real(wp_), parameter :: one=1.0_wp_,half=0.5_wp_,con1=0.1_wp_,con9=0.9_wp_,con4=0.4e-01_wp_

!  ..subroutine references..
!    fpgrre,fpknot
!  ..
!  we partition the working space.
    kx1 = kx+1
    ky1 = ky+1
    kx2 = kx1+1
    ky2 = ky1+1
    lsx = 1
    lsy = lsx+mx*kx1
    lri = lsy+my*ky1
    mm = max(nxest,my)
    lq = lri+mm
    mynx = nxest*my
    lax = lq+mynx
    nxk = nxest*kx2
    lbx = lax+nxk
    lay = lbx+nxk
    lby = lay+nyest*ky2
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! part 1: determination of the number of knots and their position.     c
! ****************************************************************     c
!  given a set of knots we compute the least-squares spline sinf(x,y), c
!  and the corresponding sum of squared residuals fp=f(p=inf).         c
!  if iopt=-1  sinf(x,y) is the requested approximation.               c
!  if iopt=0 or iopt=1 we check whether we can accept the knots:       c
!    if fp <=s we will continue with the current set of knots.         c
!    if fp > s we will increase the number of knots and compute the    c
!   corresponding least-squares spline until finally fp<=s.        c
!    the initial choice of knots depends on the value of s and iopt.   c
!    if s=0 we have spline interpolation; in that case the number of   c
!    knots equals nmaxx = mx+kx+1  and  nmaxy = my+ky+1.               c
!    if s>0 and                                                        c
! *iopt=0 we first compute the least-squares polynomial of degree  c
!  kx in x and ky in y; nx=nminx=2*kx+2 and ny=nymin=2*ky+2.       c
! *iopt=1 we start with the knots found at the last call of the    c
!  routine, except for the case that s > fp0; then we can compute  c
!  the least-squares polynomial directly.                          c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!  determine the number of knots for polynomial approximation.
    nminx = 2*kx1
    nminy = 2*ky1
    if(iopt>=0) then
!  acc denotes the absolute tolerance for the root of f(p)=s.
      acc = tol*s
!  find nmaxx and nmaxy which denote the number of knots in x- and y-
!  direction in case of spline interpolation.
      nmaxx = mx+kx1
      nmaxy = my+ky1
!  find nxe and nye which denote the maximum number of knots
!  allowed in each direction
      nxe = min(nmaxx,nxest)
      nye = min(nmaxy,nyest)
      if(s<=0.0_wp_) then
!  if s = 0, s(x,y) is an interpolating spline.
        nx = nmaxx
        ny = nmaxy
!  test whether the required storage space exceeds the available one.
        if(ny>nyest .or. nx>nxest) then
          ier = 1
          return
        end if
!  find the position of the interior knots in case of interpolation.
!  the knots in the x-direction.
        mk1 = mx-kx1
        if(mk1/=0) then
          k3 = kx/2
          i = kx1+1
          j = k3+2
          if(k3*2/=kx) then
            do l=1,mk1
              tx(i) = x(j)
              i = i+1
              j = j+1
            end do
          else
            do l=1,mk1
              tx(i) = (x(j)+x(j-1))*half
              i = i+1
              j = j+1
            end do
          end if
        end if
!  the knots in the y-direction.
        mk1 = my-ky1
        if(mk1/=0) then
          k3 = ky/2
          i = ky1+1
          j = k3+2
          if(k3*2/=ky) then
            do l=1,mk1
              ty(i) = y(j)
              i = i+1
              j = j+1
            end do
          else
            do l=1,mk1
              ty(i) = (y(j)+y(j-1))*half
              i = i+1
              j = j+1
            end do
          end if
        end if
      else
!  if s > 0 our initial choice of knots depends on the value of iopt.
        if(iopt/=0 .and. fp0>s) then
!  if iopt=1 and fp0 > s we start computing the least- squares spline
!  according to the set of knots found at the last call of the routine.
!  we determine the number of grid coordinates x(i) inside each knot
!  interval (tx(l),tx(l+1)).
          l = kx2
          j = 1
          nrdatx(1) = 0
          mpm = mx-1
          do i=2,mpm
            nrdatx(j) = nrdatx(j)+1
            if(x(i)<tx(l)) cycle
            nrdatx(j) = nrdatx(j)-1
            l = l+1
            j = j+1
            nrdatx(j) = 0
          end do
!  we determine the number of grid coordinates y(i) inside each knot
!  interval (ty(l),ty(l+1)).
          l = ky2
          j = 1
          nrdaty(1) = 0
          mpm = my-1
          do i=2,mpm
            nrdaty(j) = nrdaty(j)+1
            if(y(i)<ty(l)) cycle
            nrdaty(j) = nrdaty(j)-1
            l = l+1
            j = j+1
            nrdaty(j) = 0
          end do
        else
!  if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares
!  polynomial of degree kx in x and ky in y (which is a spline without
!  interior knots).
          nx = nminx
          ny = nminy
          nrdatx(1) = mx-2
          nrdaty(1) = my-2
          lastdi = 0
          nplusx = 0
          nplusy = 0
          fp0 = 0.0_wp_
          fpold = 0.0_wp_
          reducx = 0.0_wp_
          reducy = 0.0_wp_
        end if
      end if
    end if
    mpm = mx+my
    ifsx = 0
    ifsy = 0
    ifbx = 0
    ifby = 0
    p = -one
!  main loop for the different sets of knots.mpm=mx+my is a save upper
!  bound for the number of trials.
    do iter=1,mpm
      if(nx==nminx .and. ny==nminy) ier = -2
!  find nrintx (nrinty) which is the number of knot intervals in the
!  x-direction (y-direction).
      nrintx = nx-nminx+1
      nrinty = ny-nminy+1
!  find ncof, the number of b-spline coefficients for the current set
!  of knots.
!    nk1x = nx-kx1
!    nk1y = ny-ky1
!    ncof = nk1x*nk1y
!  find the position of the additional knots which are needed for the
!  b-spline representation of s(x,y).
      i = nx
      do j=1,kx1
        tx(j) = xb
        tx(i) = xe
        i = i-1
      end do
      i = ny
      do j=1,ky1
        ty(j) = yb
        ty(i) = ye
        i = i-1
      end do
!  find the least-squares spline sinf(x,y) and calculate for each knot
!  interval tx(j+kx)<=x<=tx(j+kx+1) (ty(j+ky)<=y<=ty(j+ky+1)) the sum
!  of squared residuals fpintx(j),j=1,2,...,nx-2*kx-1 (fpinty(j),j=1,2,
!  ...,ny-2*ky-1) for the data points having their absciss (ordinate)-
!  value belonging to that interval.
!  fp gives the total sum of squared residuals.
      call fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx,ty, &
      ny,p,c,nc,fp,fpintx,fpinty,mm,mynx,kx1,kx2,ky1,ky2,wrk(lsx), &
      wrk(lsy),wrk(lri),wrk(lq),wrk(lax),wrk(lay),wrk(lbx),wrk(lby), &
      nrx,nry)
      if(ier==(-2)) fp0 = fp
!  test whether the least-squares spline is an acceptable solution.
      if(iopt<0) return
      fpms = fp-s
      if(abs(fpms) < acc) return
!  if f(p=inf) < s, we accept the choice of knots.
      if(fpms<0.0_wp_) exit
!  if nx=nmaxx and ny=nmaxy, sinf(x,y) is an interpolating spline.
      if(nx==nmaxx .and. ny==nmaxy) then
        ier = -1
        fp = 0.0_wp_
        return
      end if
!  increase the number of knots.
!  if nx=nxe and ny=nye we cannot further increase the number of knots
!  because of the storage capacity limitation.
      if(nx==nxe .and. ny==nye) then
        ier = 1
        return
      end if
      ier = 0
!  adjust the parameter reducx or reducy according to the direction
!  in which the last added knots were located.
      if(lastdi<0) then
        reducx = fpold-fp
      else if(lastdi>0) then
        reducy = fpold-fp
      end if
!  store the sum of squared residuals for the current set of knots.
      fpold = fp
!  find nplx, the number of knots we should add in the x-direction.
      nplx = 1
      if(nx/=nminx) then
        npl1 = nplusx*2
        rn = nplusx
        if(reducx>acc) npl1 = int(rn*fpms/reducx)
        nplx = min(nplusx*2,max(npl1,nplusx/2,1))
      end if
!  find nply, the number of knots we should add in the y-direction.
      nply = 1
      if(ny/=nminy) then
        npl1 = nplusy*2
        rn = nplusy
        if(reducy>acc) npl1 = int(rn*fpms/reducy)
        nply = min0(nplusy*2,max0(npl1,nplusy/2,1))
      end if
      if (ny==nye .or. (nx/=nxe .and. ny/=nye .and.&
          (nplx<nply .or. (nplx==nply .and. lastdi>=0)))) then
!  addition in the x-direction.
        lastdi = -1
        nplusx = nplx
        ifsx = 0
        do l=1,nplusx
!  add a new knot in the x-direction
          call fpknot(x,mx,tx,nx,fpintx,nrdatx,nrintx,nxest,1)
!  test whether we cannot further increase the number of knots in the
!  x-direction.
          if (nx==nxe) exit
        end do
      else
!  addition in the y-direction.
        lastdi = 1
        nplusy = nply
        ifsy = 0
        do l=1,nplusy
!  add a new knot in the y-direction.
          call fpknot(y,my,ty,ny,fpinty,nrdaty,nrinty,nyest,1)
!  test whether we cannot further increase the number of knots in the
!  y-direction.
          if (ny==nye) exit
        end do
      end if
!  restart the computations with the new set of knots.
    end do
!  test whether the least-squares polynomial is a solution of our
!  approximation problem.
    if(ier==(-2)) return
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! part 2: determination of the smoothing spline sp(x,y)                c
! *****************************************************                c
!  we have determined the number of knots and their position. we now   c
!  compute the b-spline coefficients of the smoothing spline sp(x,y).  c
!  this smoothing spline varies with the parameter p in such a way thatc
!    f(p) = sumi=1,mx(sumj=1,my((z(i,j)-sp(x(i),y(j)))**2)             c
!  is a continuous, strictly decreasing function of p. moreover the    c
!  least-squares polynomial corresponds to p=0 and the least-squares   c
!  spline to p=infinity. iteratively we then have to determine the     c
!  positive value of p such that f(p)=s. the process which is proposed c
!  here makes use of rational interpolation. f(p) is approximated by a c
!  rational function r(p)=(u*p+v)/(p+w); three values of p (p1,p2,p3)  c
!  with corresponding values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s)c
!  are used to calculate the new value of p such that r(p)=s.          c
!  convergence is guaranteed by taking f1 > 0 and f3 < 0.              c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!  initial value for p.
    p1 = 0.0_wp_
    f1 = fp0-s
    p3 = -one
    f3 = fpms
    p = one
    ich1 = 0
    ich3 = 0
!  iteration process to find the root of f(p)=s.
    do iter = 1,maxit
!  find the smoothing spline sp(x,y) and the corresponding sum of
!  squared residuals fp.
      call fpgrre(ifsx,ifsy,ifbx,ifby,x,mx,y,my,z,mz,kx,ky,tx,nx,ty, &
      ny,p,c,nc,fp,fpintx,fpinty,mm,mynx,kx1,kx2,ky1,ky2,wrk(lsx), &
      wrk(lsy),wrk(lri),wrk(lq),wrk(lax),wrk(lay),wrk(lbx),wrk(lby), &
      nrx,nry)
!  test whether the approximation sp(x,y) is an acceptable solution.
      fpms = fp-s
      if(abs(fpms)<acc) return
!  test whether the maximum allowable number of iterations has been
!  reached.
      if(iter==maxit) then
        ier = 3
        return
      end if
!  carry out one more step of the iteration process.
      p2 = p
      f2 = fpms
      if(ich3==0) then
        if((f2-f3)<=acc) then
!  our initial choice of p is too large.
          p3 = p2
          f3 = f2
          p = p*con4
          if(p<=p1) p = p1*con9 + p2*con1
          cycle
        end if
        if(f2<0.0_wp_) ich3 = 1
      end if
      if(ich1==0) then
        if((f1-f2)<=acc) then
!  our initial choice of p is too small
          p1 = p2
          f1 = f2
          p = p/con4
          if(p3>=0.0_wp_ .and. p>=p3) p = p2*con1 + p3*con9
          cycle
        end if
!  test whether the iteration process proceeds as theoretically
!  expected.
        if(f2>0.0_wp_) ich1 = 1
      end if
      if(f2>=f1 .or. f2<=f3) then
        ier = 2
        return
      end if
!  find the new value of p.
      p = fprati(p1,f1,p2,f2,p3,f3)
    end do
  end subroutine fpregr

  subroutine fprota(cs,sn,a,b)
!  subroutine fprota applies a givens rotation to a and b.
!  ..
    implicit none
! arguments
    real(wp_), intent(in)    :: cs, sn
    real(wp_), intent(inout) :: a, b
! local variables
    real(wp_) :: stor1,stor2
!  ..
    stor1 = a
    stor2 = b
    b = cs*stor2+sn*stor1
    a = cs*stor1-sn*stor2
  end subroutine fprota

  function fprati(p1,f1,p2,f2,p3,f3)
!  given three points (p1,f1),(p2,f2) and (p3,f3), function fprati
!  gives the value of p such that the rational interpolating function
!  of the form r(p) = (u*p+v)/(p+w) equals zero at p.
!  ..
    implicit none
    real(wp_) :: fprati
! arguments
    real(wp_), intent(in) :: p2, f2
    real(wp_), intent(inout) :: p1, f1, p3, f3
! local variables
    real(wp_) :: h1, h2, h3, p
!  ..
    if(p3<=0.0_wp_) then
!  value of p in case p3 = infinity.
      p = (p1*(f1-f3)*f2-p2*(f2-f3)*f1)/((f1-f2)*f3)
    else
!  value of p in case p3 ^= infinity.
      h1 = f1*(f2-f3)
      h2 = f2*(f3-f1)
      h3 = f3*(f1-f2)
      p = -(p1*p2*h3+p2*p3*h1+p3*p1*h2)/(p1*h1+p2*h2+p3*h3)
!  adjust the value of p1,f1,p3 and f3 such that f1 > 0 and f3 < 0.
    end if
    if(f2>=0.0_wp_) then
      p1 = p2
      f1 = f2
    else
      p3 = p2
      f3 = f2
    end if
    fprati = p
  end function fprati

  subroutine regrid(iopt,mx,x,my,y,z,xb,xe,yb,ye,kx,ky,s, &
     nxest,nyest,nx,tx,ny,ty,c,fp,wrk,lwrk,iwrk,kwrk,ier)
! given the set of values z(i,j) on the rectangular grid (x(i),y(j)),
! i=1,...,mx;j=1,...,my, subroutine regrid determines a smooth bivar-
! iate spline approximation s(x,y) of degrees kx and ky on the rect-
! angle xb <= x <= xe, yb <= y <= ye.
! if iopt = -1 regrid calculates the least-squares spline according
! to a given set of knots.
! if iopt >= 0 the total numbers nx and ny of these knots and their
! position tx(j),j=1,...,nx and ty(j),j=1,...,ny are chosen automatic-
! ally by the routine. the smoothness of s(x,y) is then achieved by
! minimalizing the discontinuity jumps in the derivatives of s(x,y)
! across the boundaries of the subpanels (tx(i),tx(i+1))*(ty(j),ty(j+1).
! the amounth of smoothness is determined by the condition that f(p) =
! sum ((z(i,j)-s(x(i),y(j))))**2) be <= s, with s a given non-negative
! constant, called the smoothing factor.
! the fit is given in the b-spline representation (b-spline coefficients
! c((ny-ky-1)*(i-1)+j),i=1,...,nx-kx-1;j=1,...,ny-ky-1) and can be eval-
! uated by means of subroutine bispev.
!
! calling sequence:
! call regrid(iopt,mx,x,my,y,z,xb,xe,yb,ye,kx,ky,s,nxest,nyest,
!    *  nx,tx,ny,ty,c,fp,wrk,lwrk,iwrk,kwrk,ier)
!
! parameters:
!  iopt  : integer flag. on entry iopt must specify whether a least-
!      squares spline (iopt=-1) or a smoothing spline (iopt=0 or 1)
!      must be determined.
!      if iopt=0 the routine will start with an initial set of knots
!      tx(i)=xb,tx(i+kx+1)=xe,i=1,...,kx+1;ty(i)=yb,ty(i+ky+1)=ye,i=
!      1,...,ky+1. if iopt=1 the routine will continue with the set
!      of knots found at the last call of the routine.
!      attention: a call with iopt=1 must always be immediately pre-
!                 ceded by another call with iopt=1 or iopt=0 and
!                 s/=0.
!      unchanged on exit.
!  mx    : integer. on entry mx must specify the number of grid points
!      along the x-axis. mx > kx . unchanged on exit.
!  x     : real array of dimension at least (mx). before entry, x(i)
!      must be set to the x-co-ordinate of the i-th grid point
!      along the x-axis, for i=1,2,...,mx. these values must be
!      supplied in strictly ascending order. unchanged on exit.
!  my    : integer. on entry my must specify the number of grid points
!      along the y-axis. my > ky . unchanged on exit.
!  y     : real array of dimension at least (my). before entry, y(j)
!      must be set to the y-co-ordinate of the j-th grid point
!      along the y-axis, for j=1,2,...,my. these values must be
!      supplied in strictly ascending order. unchanged on exit.
!  z     : real array of dimension at least (mx*my).
!      before entry, z(my*(i-1)+j) must be set to the data value at
!      the grid point (x(i),y(j)) for i=1,...,mx and j=1,...,my.
!      unchanged on exit.
!  xb,xe : real values. on entry xb,xe,yb and ye must specify the bound-
!  yb,ye   aries of the rectangular approximation domain.
!      xb<=x(i)<=xe,i=1,...,mx; yb<=y(j)<=ye,j=1,...,my.
!      unchanged on exit.
!  kx,ky : integer values. on entry kx and ky must specify the degrees
!      of the spline. 1<=kx,ky<=5. it is recommended to use bicubic
!      (kx=ky=3) splines. unchanged on exit.
!  s     : real. on entry (in case iopt>=0) s must specify the smoothing
!      factor. s >=0. unchanged on exit.
!      for advice on the choice of s see further comments
!  nxest : integer. unchanged on exit.
!  nyest : integer. unchanged on exit.
!      on entry, nxest and nyest must specify an upper bound for the
!      number of knots required in the x- and y-directions respect.
!      these numbers will also determine the storage space needed by
!      the routine. nxest >= 2*(kx+1), nyest >= 2*(ky+1).
!      in most practical situation nxest = mx/2, nyest=my/2, will
!      be sufficient. always large enough are nxest=mx+kx+1, nyest=
!      my+ky+1, the number of knots needed for interpolation (s=0).
!      see also further comments.
!  nx    : integer.
!      unless ier=10 (in case iopt >=0), nx will contain the total
!      number of knots with respect to the x-variable, of the spline
!      approximation returned. if the computation mode iopt=1 is
!      used, the value of nx should be left unchanged between sub-
!      sequent calls.
!      in case iopt=-1, the value of nx should be specified on entry
!  tx    : real array of dimension nmax.
!      on succesful exit, this array will contain the knots of the
!      spline with respect to the x-variable, i.e. the position of
!      the interior knots tx(kx+2),...,tx(nx-kx-1) as well as the
!      position of the additional knots tx(1)=...=tx(kx+1)=xb and
!      tx(nx-kx)=...=tx(nx)=xe needed for the b-spline representat.
!      if the computation mode iopt=1 is used, the values of tx(1),
!      ...,tx(nx) should be left unchanged between subsequent calls.
!      if the computation mode iopt=-1 is used, the values tx(kx+2),
!      ...tx(nx-kx-1) must be supplied by the user, before entry.
!      see also the restrictions (ier=10).
!  ny    : integer.
!      unless ier=10 (in case iopt >=0), ny will contain the total
!      number of knots with respect to the y-variable, of the spline
!      approximation returned. if the computation mode iopt=1 is
!      used, the value of ny should be left unchanged between sub-
!      sequent calls.
!      in case iopt=-1, the value of ny should be specified on entry
!  ty    : real array of dimension nmax.
!      on succesful exit, this array will contain the knots of the
!      spline with respect to the y-variable, i.e. the position of
!      the interior knots ty(ky+2),...,ty(ny-ky-1) as well as the
!      position of the additional knots ty(1)=...=ty(ky+1)=yb and
!      ty(ny-ky)=...=ty(ny)=ye needed for the b-spline representat.
!      if the computation mode iopt=1 is used, the values of ty(1),
!      ...,ty(ny) should be left unchanged between subsequent calls.
!      if the computation mode iopt=-1 is used, the values ty(ky+2),
!      ...ty(ny-ky-1) must be supplied by the user, before entry.
!      see also the restrictions (ier=10).
!  c     : real array of dimension at least (nxest-kx-1)*(nyest-ky-1).
!      on succesful exit, c contains the coefficients of the spline
!      approximation s(x,y)
!  fp    : real. unless ier=10, fp contains the sum of squared
!      residuals of the spline approximation returned.
!  wrk   : real array of dimension (lwrk). used as workspace.
!      if the computation mode iopt=1 is used the values of wrk(1),
!      ...,wrk(4) should be left unchanged between subsequent calls.
!  lwrk  : integer. on entry lwrk must specify the actual dimension of
!      the array wrk as declared in the calling (sub)program.
!      lwrk must not be too small.
!       lwrk >= 4+nxest*(my+2*kx+5)+nyest*(2*ky+5)+mx*(kx+1)+
!        my*(ky+1) +u
!       where u is the larger of my and nxest.
!  iwrk  : integer array of dimension (kwrk). used as workspace.
!      if the computation mode iopt=1 is used the values of iwrk(1),
!      ...,iwrk(3) should be left unchanged between subsequent calls
!  kwrk  : integer. on entry kwrk must specify the actual dimension of
!      the array iwrk as declared in the calling (sub)program.
!      kwrk >= 3+mx+my+nxest+nyest.
!  ier   : integer. unless the routine detects an error, ier contains a
!      non-positive value on exit, i.e.
!   ier=0  : normal return. the spline returned has a residual sum of
!        squares fp such that abs(fp-s)/s <= tol with tol a relat-
!        ive tolerance set to 0.001 by the program.
!   ier=-1 : normal return. the spline returned is an interpolating
!        spline (fp=0).
!   ier=-2 : normal return. the spline returned is the least-squares
!        polynomial of degrees kx and ky. in this extreme case fp
!        gives the upper bound for the smoothing factor s.
!   ier=1  : error. the required storage space exceeds the available
!        storage space, as specified by the parameters nxest and
!        nyest.
!        probably causes : nxest or nyest too small. if these param-
!        eters are already large, it may also indicate that s is
!        too small
!        the approximation returned is the least-squares spline
!        according to the current set of knots. the parameter fp
!        gives the corresponding sum of squared residuals (fp>s).
!   ier=2  : error. a theoretically impossible result was found during
!        the iteration proces for finding a smoothing spline with
!        fp = s. probably causes : s too small.
!        there is an approximation returned but the corresponding
!        sum of squared residuals does not satisfy the condition
!        abs(fp-s)/s < tol.
!   ier=3  : error. the maximal number of iterations maxit (set to 20
!        by the program) allowed for finding a smoothing spline
!        with fp=s has been reached. probably causes : s too small
!        there is an approximation returned but the corresponding
!        sum of squared residuals does not satisfy the condition
!        abs(fp-s)/s < tol.
!   ier=10 : error. on entry, the input data are controlled on validity
!        the following restrictions must be satisfied.
!        -1<=iopt<=1, 1<=kx,ky<=5, mx>kx, my>ky, nxest>=2*kx+2,
!        nyest>=2*ky+2, kwrk>=3+mx+my+nxest+nyest,
!        lwrk >= 4+nxest*(my+2*kx+5)+nyest*(2*ky+5)+mx*(kx+1)+
!         my*(ky+1) +max(my,nxest),
!        xb<=x(i-1)<x(i)<=xe,i=2,..,mx,yb<=y(j-1)<y(j)<=ye,j=2,..,my
!        if iopt=-1: 2*kx+2<=nx<=min(nxest,mx+kx+1)
!                    xb<tx(kx+2)<tx(kx+3)<...<tx(nx-kx-1)<xe
!                    2*ky+2<=ny<=min(nyest,my+ky+1)
!                    yb<ty(ky+2)<ty(ky+3)<...<ty(ny-ky-1)<ye
!                the schoenberg-whitney conditions, i.e. there must
!                be subset of grid co-ordinates xx(p) and yy(q) such
!                that   tx(p) < xx(p) < tx(p+kx+1) ,p=1,...,nx-kx-1
!                       ty(q) < yy(q) < ty(q+ky+1) ,q=1,...,ny-ky-1
!        if iopt>=0: s>=0
!                    if s=0 : nxest>=mx+kx+1, nyest>=my+ky+1
!        if one of these conditions is found to be violated,control
!        is immediately repassed to the calling program. in that
!        case there is no approximation returned.
!
! further comments:
!   regrid does not allow individual weighting of the data-values.
!   so, if these were determined to widely different accuracies, then
!   perhaps the general data set routine surfit should rather be used
!   in spite of efficiency.
!   by means of the parameter s, the user can control the tradeoff
!   between closeness of fit and smoothness of fit of the approximation.
!   if s is too large, the spline will be too smooth and signal will be
!   lost ; if s is too small the spline will pick up too much noise. in
!   the extreme cases the program will return an interpolating spline if
!   s=0 and the least-squares polynomial (degrees kx,ky) if s is
!   very large. between these extremes, a properly chosen s will result
!   in a good compromise between closeness of fit and smoothness of fit.
!   to decide whether an approximation, corresponding to a certain s is
!   satisfactory the user is highly recommended to inspect the fits
!   graphically.
!   recommended values for s depend on the accuracy of the data values.
!   if the user has an idea of the statistical errors on the data, he
!   can also find a proper estimate for s. for, by assuming that, if he
!   specifies the right s, regrid will return a spline s(x,y) which
!   exactly reproduces the function underlying the data he can evaluate
!   the sum((z(i,j)-s(x(i),y(j)))**2) to find a good estimate for this s
!   for example, if he knows that the statistical errors on his z(i,j)-
!   values is not greater than 0.1, he may expect that a good s should
!   have a value not larger than mx*my*(0.1)**2.
!   if nothing is known about the statistical error in z(i,j), s must
!   be determined by trial and error, taking account of the comments
!   above. the best is then to start with a very large value of s (to
!   determine the least-squares polynomial and the corresponding upper
!   bound fp0 for s) and then to progressively decrease the value of s
!   ( say by a factor 10 in the beginning, i.e. s=fp0/10,fp0/100,...
!   and more carefully as the approximation shows more detail) to
!   obtain closer fits.
!   to economize the search for a good s-value the program provides with
!   different modes of computation. at the first call of the routine, or
!   whenever he wants to restart with the initial set of knots the user
!   must set iopt=0.
!   if iopt=1 the program will continue with the set of knots found at
!   the last call of the routine. this will save a lot of computation
!   time if regrid is called repeatedly for different values of s.
!   the number of knots of the spline returned and their location will
!   depend on the value of s and on the complexity of the shape of the
!   function underlying the data. if the computation mode iopt=1
!   is used, the knots returned may also depend on the s-values at
!   previous calls (if these were smaller). therefore, if after a number
!   of trials with different s-values and iopt=1, the user can finally
!   accept a fit as satisfactory, it may be worthwhile for him to call
!   regrid once more with the selected value for s but now with iopt=0.
!   indeed, regrid may then return an approximation of the same quality
!   of fit but with fewer knots and therefore better if data reduction
!   is also an important objective for the user.
!   the number of knots may also depend on the upper bounds nxest and
!   nyest. indeed, if at a certain stage in regrid the number of knots
!   in one direction (say nx) has reached the value of its upper bound
!   (nxest), then from that moment on all subsequent knots are added
!   in the other (y) direction. this may indicate that the value of
!   nxest is too small. on the other hand, it gives the user the option
!   of limiting the number of knots the routine locates in any direction
!   for example, by setting nxest=2*kx+2 (the lowest allowable value for
!   nxest), the user can indicate that he wants an approximation which
!   is a simple polynomial of degree kx in the variable x.
!
!  other subroutines required:
!    fpback,fpbspl,fpregr,fpdisc,fpgivs,fpgrre,fprati,fprota,fpchec,
!    fpknot
!
!  references:
!   dierckx p. : a fast algorithm for smoothing data on a rectangular
!            grid while using spline functions, siam j.numer.anal.
!            19 (1982) 1286-1304.
!   dierckx p. : a fast algorithm for smoothing data on a rectangular
!            grid while using spline functions, report tw53, dept.
!            computer science,k.u.leuven, 1980.
!   dierckx p. : curve and surface fitting with splines, monographs on
!            numerical analysis, oxford university press, 1993.
!
!  author:
!    p.dierckx
!    dept. computer science, k.u. leuven
!    celestijnenlaan 200a, b-3001 heverlee, belgium.
!    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
!  creation date : may 1979
!  latest update : march 1989
!
!  ..
    implicit none
! arguments
    integer, intent(in) :: iopt, mx, my, kx, ky, nxest, nyest, lwrk, kwrk
    integer, intent(out) :: ier
    integer, intent(inout) :: nx, ny, iwrk(kwrk)
    real(wp_), intent(in) :: xb, xe, yb, ye, x(mx), y(my), z(mx*my), s
    real(wp_), intent(out) :: fp, c((nxest-kx-1)*(nyest-ky-1))
    real(wp_), intent(inout) :: tx(nxest), ty(nyest), wrk(lwrk)
! local variables
    integer :: i, j, jwrk, kndx, kndy, knrx, knry, kwest, kx1, kx2, &
               ky1, ky2, lfpx, lfpy, lwest, lww, nc, nminx, nminy, mz
! parameters
    integer, parameter :: maxit = 20
    real(wp_), parameter :: tol = 0.1e-02_wp_
!  ..subroutine references..
!    fpregr,fpchec
!  ..
!  before starting computations a data check is made. if the input data
!  are invalid, control is immediately repassed to the calling program.
    ier = 10
    if(kx<=0 .or. kx>5) return
    kx1 = kx+1
    kx2 = kx1+1
    if(ky<=0 .or. ky>5) return
    ky1 = ky+1
    ky2 = ky1+1
    if(iopt<(-1) .or. iopt>1) return
    nminx = 2*kx1
    if(mx<kx1 .or. nxest<nminx) return
    nminy = 2*ky1
    if(my<ky1 .or. nyest<nminy) return
    mz = mx*my
    nc = (nxest-kx1)*(nyest-ky1)
    lwest = 4+nxest*(my+2*kx2+1)+nyest*(2*ky2+1)+mx*kx1+ &
     my*ky1+max(nxest,my)
    kwest = 3+mx+my+nxest+nyest
    if(lwrk<lwest .or. kwrk<kwest) return
    if(xb>x(1) .or. xe<x(mx)) return
    do i=2,mx
      if(x(i-1)>=x(i)) return
    end do
    if(yb>y(1) .or. ye<y(my)) return
    do i=2,my
      if(y(i-1)>=y(i)) return
    end do
    if(iopt<0) then
      if(nx<nminx .or. nx>nxest) return
      j = nx
      do i=1,kx1
        tx(i) = xb
        tx(j) = xe
        j = j-1
      end do
      call fpchec(x,mx,tx,nx,kx,ier)
      if(ier/=0) return
      if(ny<nminy .or. ny>nyest) return
      j = ny
      do i=1,ky1
        ty(i) = yb
        ty(j) = ye
        j = j-1
      end do
      call fpchec(y,my,ty,ny,ky,ier)
      if(ier/=0) return
    else
      if(s<0.0_wp_) return
      if(s==0.0_wp_ .and. (nxest<(mx+kx1) .or. nyest<(my+ky1)) ) &
        return
      ier = 0
    end if
!  we partition the working space and determine the spline approximation
    lfpx = 5
    lfpy = lfpx+nxest
    lww = lfpy+nyest
    jwrk = lwrk-4-nxest-nyest
    knrx = 4
    knry = knrx+mx
    kndx = knry+my
    kndy = kndx+nxest
    call fpregr(iopt,x,mx,y,my,z,mz,xb,xe,yb,ye,kx,ky,s,nxest,nyest, &
     tol,maxit,nc,nx,tx,ny,ty,c,fp,wrk(1),wrk(2),wrk(3),wrk(4), &
     wrk(lfpx),wrk(lfpy),iwrk(1),iwrk(2),iwrk(3),iwrk(knrx), &
     iwrk(knry),iwrk(kndx),iwrk(kndy),wrk(lww),jwrk,ier)
  end subroutine regrid

  subroutine parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z, &
     wrk,lwrk,iwrk,kwrk,ier)
!  subroutine parder evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
!  ,my the partial derivative ( order nux,nuy) of a bivariate spline
!  s(x,y) of degrees kx and ky, given in the b-spline representation.
!
!  calling sequence:
! call parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,wrk,lwrk,
!    * iwrk,kwrk,ier)
!
!  input parameters:
!   tx    : real array, length nx, which contains the position of the
!       knots in the x-direction.
!   nx    : integer, giving the total number of knots in the x-direction
!   ty    : real array, length ny, which contains the position of the
!       knots in the y-direction.
!   ny    : integer, giving the total number of knots in the y-direction
!   c     : real array, length (nx-kx-1)*(ny-ky-1), which contains the
!       b-spline coefficients.
!   kx,ky : integer values, giving the degrees of the spline.
!   nux   : integer values, specifying the order of the partial
!   nuy     derivative. 0<=nux<kx, 0<=nuy<ky.
!   x     : real array of dimension (mx).
!       before entry x(i) must be set to the x co-ordinate of the
!       i-th grid point along the x-axis.
!       tx(kx+1)<=x(i-1)<=x(i)<=tx(nx-kx), i=2,...,mx.
!   mx    : on entry mx must specify the number of grid points along
!       the x-axis. mx >=1.
!   y     : real array of dimension (my).
!       before entry y(j) must be set to the y co-ordinate of the
!       j-th grid point along the y-axis.
!       ty(ky+1)<=y(j-1)<=y(j)<=ty(ny-ky), j=2,...,my.
!   my    : on entry my must specify the number of grid points along
!       the y-axis. my >=1.
!   wrk   : real array of dimension lwrk. used as workspace.
!   lwrk  : integer, specifying the dimension of wrk.
!       lwrk >= mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1)
!   iwrk  : integer array of dimension kwrk. used as workspace.
!   kwrk  : integer, specifying the dimension of iwrk. kwrk >= mx+my.
!
!  output parameters:
!   z     : real array of dimension (mx*my).
!       on succesful exit z(my*(i-1)+j) contains the value of the
!       specified partial derivative of s(x,y) at the point
!       (x(i),y(j)),i=1,...,mx;j=1,...,my.
!   ier   : integer error flag
!    ier=0 : normal return
!    ier=10: invalid input data (see restrictions)
!
!  restrictions:
!   mx >=1, my >=1, 0 <= nux < kx, 0 <= nuy < ky, kwrk>=mx+my
!   lwrk>=mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1),
!   tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx
!   ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my
!
!  other subroutines required:
!    fpbisp,fpbspl
!
!  references :
!    de boor c : on calculating with b-splines, j. approximation theory
!            6 (1972) 50-62.
!   dierckx p. : curve and surface fitting with splines, monographs on
!            numerical analysis, oxford university press, 1993.
!
!  author :
!    p.dierckx
!    dept. computer science, k.u.leuven
!    celestijnenlaan 200a, b-3001 heverlee, belgium.
!    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
!  latest update : march 1989
!
    implicit none
! arguments
    integer, intent(in) :: nx, ny, kx, ky, nux, nuy, mx, my, lwrk, kwrk
    integer, intent(out) :: ier
    integer, intent(inout) :: iwrk(kwrk)
    real(wp_), intent(in) :: tx(nx), ty(ny), c((nx-kx-1)*(ny-ky-1)), &
                           x(mx), y(my)
    real(wp_), intent(out) :: z(mx*my)
    real(wp_), intent(inout) :: wrk(lwrk)
! local variables
    integer :: i, iwx, iwy, j, kkx, kky, kx1, ky1, lx, ly, lwest,   &
               l1, l2, m, m0, m1, nc, nkx1, nky1, nxx, nyy
    real(wp_) :: ak, fac
!  ..
!  before starting computations a data check is made. if the input data
!  are invalid control is immediately repassed to the calling program.
    ier = 10
    kx1 = kx+1
    ky1 = ky+1
    nkx1 = nx-kx1
    nky1 = ny-ky1
    nc = nkx1*nky1
    if(nux<0 .or. nux>=kx) return
    if(nuy<0 .or. nuy>=ky) return
    lwest = nc +(kx1-nux)*mx+(ky1-nuy)*my
    if(lwrk<lwest) return
    if(kwrk<(mx+my)) return
    if(mx<1) return
    do i=2,mx
      if(x(i)<x(i-1)) return
    end do
    if(my<1) return
    do i=2,my
      if(y(i)<y(i-1)) return
    end do
    ier = 0
    nxx = nkx1
    nyy = nky1
    kkx = kx
    kky = ky
!  the partial derivative of order (nux,nuy) of a bivariate spline of
!  degrees kx,ky is a bivariate spline of degrees kx-nux,ky-nuy.
!  we calculate the b-spline coefficients of this spline
    do i=1,nc
      wrk(i) = c(i)
    end do
    if(nux/=0) then
      lx = 1
      do j=1,nux
        ak = kkx
        nxx = nxx-1
        l1 = lx
        m0 = 1
        do i=1,nxx
          l1 = l1+1
          l2 = l1+kkx
          fac = tx(l2)-tx(l1)
          if(fac>0.0_wp_) then
            do m=1,nyy
              m1 = m0+nyy
              wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
              m0  = m0+1
            end do
          end if
        end do
        lx = lx+1
        kkx = kkx-1
      end do
    end if
    if(nuy/=0) then
      ly = 1
      do j=1,nuy
        ak = kky
        nyy = nyy-1
        l1 = ly
        do i=1,nyy
          l1 = l1+1
          l2 = l1+kky
          fac = ty(l2)-ty(l1)
          if(fac>0.0_wp_) then
            m0 = i
            do m=1,nxx
              m1 = m0+1
              wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
              m0  = m0+nky1
            end do
          end if
        end do
        ly = ly+1
        kky = kky-1
      end do
      m0 = nyy
      m1 = nky1
      do m=2,nxx
        do i=1,nyy
          m0 = m0+1
          m1 = m1+1
          wrk(m0) = wrk(m1)
        end do
        m1 = m1+nuy
      end do
    end if
!  we partition the working space and evaluate the partial derivative
    iwx = 1+nxx*nyy
    iwy = iwx+mx*(kx1-nux)
    call fpbisp(tx(nux+1),nx-2*nux,ty(nuy+1),ny-2*nuy,wrk,kkx,kky, &
     x,mx,y,my,z,wrk(iwx),wrk(iwy),iwrk(1),iwrk(mx+1))
  end subroutine parder

  subroutine coeff_parder(tx,nx,ty,ny,c,kx,ky,nux,nuy, &
     wrk,lwrk,ier)
!  subroutine parder evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
!  ,my the partial derivative ( order nux,nuy) of a bivariate spline
!  s(x,y) of degrees kx and ky, given in the b-spline representation.
!
!  calling sequence:
! call parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,mx,y,my,z,wrk,lwrk,
!    * iwrk,kwrk,ier)
!
!  input parameters:
!   tx    : real array, length nx, which contains the position of the
!       knots in the x-direction.
!   nx    : integer, giving the total number of knots in the x-direction
!   ty    : real array, length ny, which contains the position of the
!       knots in the y-direction.
!   ny    : integer, giving the total number of knots in the y-direction
!   c     : real array, length (nx-kx-1)*(ny-ky-1), which contains the
!       b-spline coefficients.
!   kx,ky : integer values, giving the degrees of the spline.
!   nux   : integer values, specifying the order of the partial
!   nuy     derivative. 0<=nux<kx, 0<=nuy<ky.
!   wrk   : real array of dimension lwrk. used as workspace.
!   lwrk  : integer, specifying the dimension of wrk.
!       lwrk >= mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1)
!
!  output parameters:
!   ier   : integer error flag
!    ier=0 : normal return
!    ier=10: invalid input data (see restrictions)
!
!  restrictions:
!   mx >=1, my >=1, 0 <= nux < kx, 0 <= nuy < ky
!   lwrk>=mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1),
!
!  other subroutines required:
!    fpbisp,fpbspl
!
!  references :
!    de boor c : on calculating with b-splines, j. approximation theory
!            6 (1972) 50-62.
!   dierckx p. : curve and surface fitting with splines, monographs on
!            numerical analysis, oxford university press, 1993.
!
!  author :
!    p.dierckx
!    dept. computer science, k.u.leuven
!    celestijnenlaan 200a, b-3001 heverlee, belgium.
!    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
!  latest update : march 1989
!
    implicit none
! arguments
    integer, intent(in) :: nx, ny, kx, ky, nux, nuy, lwrk
    integer, intent(out) :: ier
    real(wp_), intent(in) :: tx(nx), ty(ny), c((nx-kx-1)*(ny-ky-1))
    real(wp_), intent(inout) :: wrk(lwrk)
! local variables
    integer :: mx, my, i, j, kkx, kky, kx1, ky1, lx, ly, lwest,   &
               l1, l2, m, m0, m1, nc, nkx1, nky1, nxx, nyy
    real(wp_) :: ak, fac
!  ..
!  before starting computations a data check is made. if the input data
!  are invalid control is immediately repassed to the calling program.
    mx=1
    my=1
    ier = 10
    kx1 = kx+1
    ky1 = ky+1
    nkx1 = nx-kx1
    nky1 = ny-ky1
    nc = nkx1*nky1
    if(nux<0 .or. nux>=kx) return
    if(nuy<0 .or. nuy>=ky) return
    lwest = nc +(kx1-nux)*mx+(ky1-nuy)*my
    if(lwrk<lwest) return
    ier = 0
    nxx = nkx1
    nyy = nky1
    kkx = kx
    kky = ky
!  the partial derivative of order (nux,nuy) of a bivariate spline of
!  degrees kx,ky is a bivariate spline of degrees kx-nux,ky-nuy.
!  we calculate the b-spline coefficients of this spline
    do i=1,nc
      wrk(i) = c(i)
    end do
    if(nux/=0) then
      lx = 1
      do j=1,nux
        ak = kkx
        nxx = nxx-1
        l1 = lx
        m0 = 1
        do i=1,nxx
          l1 = l1+1
          l2 = l1+kkx
          fac = tx(l2)-tx(l1)
          if(fac>0.0_wp_) then
            do m=1,nyy
              m1 = m0+nyy
              wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
              m0  = m0+1
            end do
          end if
        end do
        lx = lx+1
        kkx = kkx-1
      end do
    end if
    if(nuy/=0) then
      ly = 1
      do j=1,nuy
        ak = kky
        nyy = nyy-1
        l1 = ly
        do i=1,nyy
          l1 = l1+1
          l2 = l1+kky
          fac = ty(l2)-ty(l1)
          if(fac>0.0_wp_) then
            m0 = i
            do m=1,nxx
              m1 = m0+1
              wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
              m0  = m0+nky1
            end do
          end if
        end do
        ly = ly+1
        kky = kky-1
      end do
      m0 = nyy
      m1 = nky1
      do m=2,nxx
        do i=1,nyy
          m0 = m0+1
          m1 = m1+1
          wrk(m0) = wrk(m1)
        end do
        m1 = m1+nuy
      end do
    end if
  end subroutine coeff_parder

  subroutine curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp, &
     wrk,lwrk,iwrk,ier)
!  given the set of data points (x(i),y(i)) and the set of positive
!  numbers w(i),i=1,2,...,m,subroutine curfit determines a smooth spline
!  approximation of degree k on the interval xb <= x <= xe.
!  if iopt=-1 curfit calculates the weighted least-squares spline
!  according to a given set of knots.
!  if iopt>=0 the number of knots of the spline s(x) and the position
!  t(j),j=1,2,...,n is chosen automatically by the routine. the smooth-
!  ness of s(x) is then achieved by minimalizing the discontinuity
!  jumps of the k-th derivative of s(x) at the knots t(j),j=k+2,k+3,...,
!  n-k-1. the amount of smoothness is determined by the condition that
!  f(p)=sum((w(i)*(y(i)-s(x(i))))**2) be <= s, with s a given non-
!  negative constant, called the smoothing factor.
!  the fit s(x) is given in the b-spline representation (b-spline coef-
!  ficients c(j),j=1,2,...,n-k-1) and can be evaluated by means of
!  subroutine splev.
!
!  calling sequence:
! call curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,wrk,
!    * lwrk,iwrk,ier)
!
!  parameters:
!   iopt  : integer flag. on entry iopt must specify whether a weighted
!       least-squares spline (iopt=-1) or a smoothing spline (iopt=
!       0 or 1) must be determined. if iopt=0 the routine will start
!       with an initial set of knots t(i)=xb, t(i+k+1)=xe, i=1,2,...
!       k+1. if iopt=1 the routine will continue with the knots
!       found at the last call of the routine.
!       attention: a call with iopt=1 must always be immediately
!       preceded by another call with iopt=1 or iopt=0.
!       unchanged on exit.
!   m     : integer. on entry m must specify the number of data points.
!       m > k. unchanged on exit.
!   x     : real array of dimension at least (m). before entry, x(i)
!       must be set to the i-th value of the independent variable x,
!       for i=1,2,...,m. these values must be supplied in strictly
!       ascending order. unchanged on exit.
!   y     : real array of dimension at least (m). before entry, y(i)
!       must be set to the i-th value of the dependent variable y,
!       for i=1,2,...,m. unchanged on exit.
!   w     : real array of dimension at least (m). before entry, w(i)
!       must be set to the i-th value in the set of weights. the
!       w(i) must be strictly positive. unchanged on exit.
!       see also further comments.
!   xb,xe : real values. on entry xb and xe must specify the boundaries
!       of the approximation interval. xb<=x(1), xe>=x(m).
!       unchanged on exit.
!   k     : integer. on entry k must specify the degree of the spline.
!       1<=k<=5. it is recommended to use cubic splines (k=3).
!       the user is strongly dissuaded from choosing k even,together
!       with a small s-value. unchanged on exit.
!   s     : real.on entry (in case iopt>=0) s must specify the smoothing
!       factor. s >=0. unchanged on exit.
!       for advice on the choice of s see further comments.
!   nest  : integer. on entry nest must contain an over-estimate of the
!       total number of knots of the spline returned, to indicate
!       the storage space available to the routine. nest >=2*k+2.
!       in most practical situation nest=m/2 will be sufficient.
!       always large enough is  nest=m+k+1, the number of knots
!       needed for interpolation (s=0). unchanged on exit.
!   n     : integer.
!       unless ier =10 (in case iopt >=0), n will contain the
!       total number of knots of the spline approximation returned.
!       if the computation mode iopt=1 is used this value of n
!       should be left unchanged between subsequent calls.
!       in case iopt=-1, the value of n must be specified on entry.
!   t     : real array of dimension at least (nest).
!       on succesful exit, this array will contain the knots of the
!       spline,i.e. the position of the interior knots t(k+2),t(k+3)
!       ...,t(n-k-1) as well as the position of the additional knots
!       t(1)=t(2)=...=t(k+1)=xb and t(n-k)=...=t(n)=xe needed for
!       the b-spline representation.
!       if the computation mode iopt=1 is used, the values of t(1),
!       t(2),...,t(n) should be left unchanged between subsequent
!       calls. if the computation mode iopt=-1 is used, the values
!       t(k+2),...,t(n-k-1) must be supplied by the user, before
!       entry. see also the restrictions (ier=10).
!   c     : real array of dimension at least (nest).
!       on succesful exit, this array will contain the coefficients
!       c(1),c(2),..,c(n-k-1) in the b-spline representation of s(x)
!   fp    : real. unless ier=10, fp contains the weighted sum of
!       squared residuals of the spline approximation returned.
!   wrk   : real array of dimension at least (m*(k+1)+nest*(7+3*k)).
!       used as working space. if the computation mode iopt=1 is
!       used, the values wrk(1),...,wrk(n) should be left unchanged
!       between subsequent calls.
!   lwrk  : integer. on entry,lwrk must specify the actual dimension of
!       the array wrk as declared in the calling (sub)program.lwrk
!       must not be too small (see wrk). unchanged on exit.
!   iwrk  : integer array of dimension at least (nest).
!       used as working space. if the computation mode iopt=1 is
!       used,the values iwrk(1),...,iwrk(n) should be left unchanged
!       between subsequent calls.
!   ier   : integer. unless the routine detects an error, ier contains a
!       non-positive value on exit, i.e.
!    ier=0  : normal return. the spline returned has a residual sum of
!         squares fp such that abs(fp-s)/s <= tol with tol a relat-
!         ive tolerance set to 0.001 by the program.
!    ier=-1 : normal return. the spline returned is an interpolating
!         spline (fp=0).
!    ier=-2 : normal return. the spline returned is the weighted least-
!         squares polynomial of degree k. in this extreme case fp
!         gives the upper bound fp0 for the smoothing factor s.
!    ier=1  : error. the required storage space exceeds the available
!         storage space, as specified by the parameter nest.
!         probably causes : nest too small. if nest is already
!         large (say nest > m/2), it may also indicate that s is
!         too small
!         the approximation returned is the weighted least-squares
!         spline according to the knots t(1),t(2),...,t(n). (n=nest)
!         the parameter fp gives the corresponding weighted sum of
!         squared residuals (fp>s).
!    ier=2  : error. a theoretically impossible result was found during
!         the iteration proces for finding a smoothing spline with
!         fp = s. probably causes : s too small.
!         there is an approximation returned but the corresponding
!         weighted sum of squared residuals does not satisfy the
!         condition abs(fp-s)/s < tol.
!    ier=3  : error. the maximal number of iterations maxit (set to 20
!         by the program) allowed for finding a smoothing spline
!         with fp=s has been reached. probably causes : s too small
!         there is an approximation returned but the corresponding
!         weighted sum of squared residuals does not satisfy the
!         condition abs(fp-s)/s < tol.
!    ier=10 : error. on entry, the input data are controlled on validity
!         the following restrictions must be satisfied.
!         -1<=iopt<=1, 1<=k<=5, m>k, nest>2*k+2, w(i)>0,i=1,2,...,m
!         xb<=x(1)<x(2)<...<x(m)<=xe, lwrk>=(k+1)*m+nest*(7+3*k)
!         if iopt=-1: 2*k+2<=n<=min(nest,m+k+1)
!                     xb<t(k+2)<t(k+3)<...<t(n-k-1)<xe
!                   the schoenberg-whitney conditions, i.e. there
!                   must be a subset of data points xx(j) such that
!                     t(j) < xx(j) < t(j+k+1), j=1,2,...,n-k-1
!         if iopt>=0: s>=0
!                     if s=0 : nest >= m+k+1
!         if one of these conditions is found to be violated,control
!         is immediately repassed to the calling program. in that
!         case there is no approximation returned.
!
!  further comments:
!   by means of the parameter s, the user can control the tradeoff
!   between closeness of fit and smoothness of fit of the approximation.
!   if s is too large, the spline will be too smooth and signal will be
!   lost ; if s is too small the spline will pick up too much noise. in
!   the extreme cases the program will return an interpolating spline if
!   s=0 and the weighted least-squares polynomial of degree k if s is
!   very large. between these extremes, a properly chosen s will result
!   in a good compromise between closeness of fit and smoothness of fit.
!   to decide whether an approximation, corresponding to a certain s is
!   satisfactory the user is highly recommended to inspect the fits
!   graphically.
!   recommended values for s depend on the weights w(i). if these are
!   taken as 1/d(i) with d(i) an estimate of the standard deviation of
!   y(i), a good s-value should be found in the range (m-sqrt(2*m),m+
!   sqrt(2*m)). if nothing is known about the statistical error in y(i)
!   each w(i) can be set equal to one and s determined by trial and
!   error, taking account of the comments above. the best is then to
!   start with a very large value of s ( to determine the least-squares
!   polynomial and the corresponding upper bound fp0 for s) and then to
!   progressively decrease the value of s ( say by a factor 10 in the
!   beginning, i.e. s=fp0/10, fp0/100,...and more carefully as the
!   approximation shows more detail) to obtain closer fits.
!   to economize the search for a good s-value the program provides with
!   different modes of computation. at the first call of the routine, or
!   whenever he wants to restart with the initial set of knots the user
!   must set iopt=0.
!   if iopt=1 the program will continue with the set of knots found at
!   the last call of the routine. this will save a lot of computation
!   time if curfit is called repeatedly for different values of s.
!   the number of knots of the spline returned and their location will
!   depend on the value of s and on the complexity of the shape of the
!   function underlying the data. but, if the computation mode iopt=1
!   is used, the knots returned may also depend on the s-values at
!   previous calls (if these were smaller). therefore, if after a number
!   of trials with different s-values and iopt=1, the user can finally
!   accept a fit as satisfactory, it may be worthwhile for him to call
!   curfit once more with the selected value for s but now with iopt=0.
!   indeed, curfit may then return an approximation of the same quality
!   of fit but with fewer knots and therefore better if data reduction
!   is also an important objective for the user.
!
!  other subroutines required:
!    fpback,fpbspl,fpchec,fpcurf,fpdisc,fpgivs,fpknot,fprati,fprota
!
!  references:
!   dierckx p. : an algorithm for smoothing, differentiation and integ-
!            ration of experimental data using spline functions,
!            j.comp.appl.maths 1 (1975) 165-184.
!   dierckx p. : a fast algorithm for smoothing data on a rectangular
!            grid while using spline functions, siam j.numer.anal.
!            19 (1982) 1286-1304.
!   dierckx p. : an improved algorithm for curve fitting with spline
!            functions, report tw54, dept. computer science,k.u.
!            leuven, 1981.
!   dierckx p. : curve and surface fitting with splines, monographs on
!            numerical analysis, oxford university press, 1993.
!
!  author:
!    p.dierckx
!    dept. computer science, k.u. leuven
!    celestijnenlaan 200a, b-3001 heverlee, belgium.
!    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
!  creation date : may 1979
!  latest update : march 1987
!
!  ..
    implicit none
! arguments
    integer, intent(in)  :: iopt, m, k, nest, lwrk
    integer, intent(out) :: ier
    integer, intent(inout) :: n, iwrk(nest)
    real(wp_), intent(in)    :: xb, xe, s, x(m), y(m), w(m)
    real(wp_), intent(out)   :: fp, c(nest)
    real(wp_), intent(inout) :: t(nest), wrk(lwrk)
! local variables
    integer :: i, ia, ib, ifp, ig, iq, iz, j, k1, k2, lwest, nmin
! parameters
    integer, parameter :: maxit = 20
    real(wp_), parameter :: tol = 0.1e-02_wp_
!  ..
!  before starting computations a data check is made. if the input data
!  are invalid, control is immediately repassed to the calling program.
    ier = 10
    if(k<=0 .or. k>5) return
    k1 = k+1
    k2 = k1+1
    if(iopt<(-1) .or. iopt>1) return
    nmin = 2*k1
    if(m<k1 .or. nest<nmin) return
    lwest = m*k1+nest*(7+3*k)
    if(lwrk<lwest) return
    if(xb>x(1) .or. xe<x(m) .or. w(1)<=0.0_wp_) return
    do i=2,m
       if(x(i-1)>=x(i) .or. w(i)<=0.0_wp_) return
    end do
    if(iopt<0) then
      if(n<nmin .or. n>nest) return
      j = n
      do i=1,k1
         t(i) = xb
         t(j) = xe
         j = j-1
      end do
      call fpchec(x,m,t,n,k,ier)
      if(ier/=0) return
    else
      if(s<0.0_wp_) return
      if(s==0.0_wp_ .and. nest<(m+k1)) return
      ier = 0
    end if
! we partition the working space and determine the spline approximation.
    ifp = 1
    iz = ifp+nest
    ia = iz+nest
    ib = ia+nest*k1
    ig = ib+nest*k2
    iq = ig+nest*k2
    call fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp, &
     wrk(ifp),wrk(iz),wrk(ia),wrk(ib),wrk(ig),wrk(iq),iwrk,ier)
  end subroutine curfit

  subroutine fpcurf(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2, &
     n,t,c,fp,fpint,z,a,b,g,q,nrdata,ier)
!  ..
    implicit none
! arguments
    integer, intent(in) :: iopt, m, k, nest, maxit, k1, k2
    integer, intent(out) :: ier
    integer, intent(inout) :: n, nrdata(nest)
    real(wp_), intent(in) :: xb, xe, s, tol, x(m), y(m), w(m)
    real(wp_), intent(out) :: fp, c(nest)
    real(wp_), intent(inout) :: t(nest), fpint(nest), z(nest), &
                      a(nest,k1), b(nest,k2), g(nest,k2), q(m,k1)
! local variables
    real(wp_) :: acc, cs, fpart, fpms, fpold, fp0, f1, f2, f3, p, pinv, &
               piv, p1, p2, p3, rn, sn, store, term, wi, xi, yi, h(7)
    integer :: i, ich1, ich3, it, iter, i1, i2, i3, j, k3, l, l0, mk1, new, &
               nk1, nmax, nmin, nplus, npl1, nrint, n8
    logical :: rstart
! parameters
    real(wp_), parameter :: one=1.0_wp_,con1=0.1_wp_,con9=0.9_wp_,con4=0.4e-01_wp_,half=0.5_wp_
!  ..function references
!  real(8) abs
!  integer max0,min0
!  ..subroutine references..
!    fpback,fpbspl,fpgivs,fpdisc,fpknot,fprota
!  ..
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!  part 1: determination of the number of knots and their position     c
!  **************************************************************      c
!  given a set of knots we compute the least-squares spline sinf(x),   c
!  and the corresponding sum of squared residuals fp=f(p=inf).         c
!  if iopt=-1 sinf(x) is the requested approximation.                  c
!  if iopt=0 or iopt=1 we check whether we can accept the knots:       c
!    if fp <=s we will continue with the current set of knots.         c
!    if fp > s we will increase the number of knots and compute the    c
!   corresponding least-squares spline until finally fp<=s.        c
!    the initial choice of knots depends on the value of s and iopt.   c
!    if s=0 we have spline interpolation; in that case the number of   c
!    knots equals nmax = m+k+1.                                        c
!    if s > 0 and                                                      c
!  iopt=0 we first compute the least-squares polynomial of         c
!  degree k; n = nmin = 2*k+2                                      c
!  iopt=1 we start with the set of knots found at the last         c
!  call of the routine, except for the case that s > fp0; then     c
!  we compute directly the least-squares polynomial of degree k.   c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!  determine nmin, the number of knots for polynomial approximation.
    nmin = 2*k1
    if(iopt>=0) then
!  calculation of acc, the absolute tolerance for the root of f(p)=s.
      acc = tol*s
!  determine nmax, the number of knots for spline interpolation.
      nmax = m+k1
      if(s<=0.0_wp_) then
!  if s=0, s(x) is an interpolating spline.
!  test whether the required storage space exceeds the available one.
        n = nmax
        if(nmax>nest) then
          ier = 1
          return
        end if
!  find the position of the interior knots in case of interpolation.
        mk1 = m-k1
        if(mk1/=0) then
          k3 = k/2
          i = k2
          j = k3+2
          if(k3*2/=k) then
            do l=1,mk1
              t(i) = x(j)
              i = i+1
              j = j+1
            end do
          else
            do l=1,mk1
              t(i) = (x(j)+x(j-1))*half
              i = i+1
              j = j+1
            end do
          end if
        end if
      else
!  if s>0 our initial choice of knots depends on the value of iopt.
!  if iopt=0 or iopt=1 and s>=fp0, we start computing the least-squares
!  polynomial of degree k which is a spline without interior knots.
!  if iopt=1 and fp0>s we start computing the least squares spline
!  according to the set of knots found at the last call of the routine.
        if (iopt==0 .or. n==nmin) then
          n = nmin
          fpold = 0.0_wp_
          nplus = 0
          nrdata(1) = m-2
        else
          fp0 = fpint(n)
          if (s>=fp0) then
            n = nmin
            fpold = 0.0_wp_
            nplus = 0
            nrdata(1) = m-2
          else
            fpold = fpint(n-1)
            nplus = nrdata(n)
          end if
        end if
      end if
    end if
    do
      rstart=.false.
!  main loop for the different sets of knots. m is a save upper bound
!  for the number of trials.
      do iter = 1,m
        if(n==nmin) ier = -2
!  find nrint, tne number of knot intervals.
        nrint = n-nmin+1
!  find the position of the additional knots which are needed for
!  the b-spline representation of s(x).
        nk1 = n-k1
        i = n
        do j=1,k1
          t(j) = xb
          t(i) = xe
          i = i-1
        end do
!  compute the b-spline coefficients of the least-squares spline
!  sinf(x). the observation matrix a is built up row by row and
!  reduced to upper triangular form by givens transformations.
!  at the same time fp=f(p=inf) is computed.
        fp = 0.0_wp_
!  initialize the observation matrix a.
        do i=1,nk1
          z(i) = 0.0_wp_
          do j=1,k1
            a(i,j) = 0.0_wp_
          end do
        end do
        l = k1
        do it=1,m
!  fetch the current data point x(it),y(it).
          xi = x(it)
          wi = w(it)
          yi = y(it)*wi
!  search for knot interval t(l) <= xi < t(l+1).
          do
            if(xi<t(l+1) .or. l==nk1) exit
            l = l+1
          end do
!  evaluate the (k+1) non-zero b-splines at xi and store them in q.
          call fpbspl(t,n,k,xi,l,h)
          do i=1,k1
            q(it,i) = h(i)
            h(i) = h(i)*wi
          end do
!  rotate the new row of the observation matrix into triangle.
          j = l-k1
          do i=1,k1
            j = j+1
            piv = h(i)
            if(piv==0.0_wp_) cycle
!  calculate the parameters of the givens transformation.
            call fpgivs(piv,a(j,1),cs,sn)
!  transformations to right hand side.
            call fprota(cs,sn,yi,z(j))
            if(i==k1) exit
            i2 = 1
            i3 = i+1
            do i1 = i3,k1
              i2 = i2+1
!  transformations to left hand side.
              call fprota(cs,sn,h(i1),a(j,i2))
            end do
          end do
!  add contribution of this row to the sum of squares of residual
!  right hand sides.
          fp = fp+yi**2
        end do
        if(ier==(-2)) fp0 = fp
        fpint(n) = fp0
        fpint(n-1) = fpold
        nrdata(n) = nplus
!  backward substitution to obtain the b-spline coefficients.
        call fpback(a,z,nk1,k1,c,nest)
!  test whether the approximation sinf(x) is an acceptable solution.
        if(iopt<0) return
        fpms = fp-s
        if(abs(fpms)<acc) return
!  if f(p=inf) < s accept the choice of knots.
        if(fpms<0.0_wp_) exit
!  if n = nmax, sinf(x) is an interpolating spline.
        if(n==nmax) then
          ier = -1
          return
        end if
!  increase the number of knots.
!  if n=nest we cannot increase the number of knots because of
!  the storage capacity limitation.
        if(n==nest) then
          ier = 1
          return
        end if
!  determine the number of knots nplus we are going to add.
        if(ier/=0) then
          nplus = 1
          ier = 0
        else
          npl1 = nplus*2
          rn = nplus
          if(fpold-fp>acc) npl1 = int(rn*fpms/(fpold-fp))
          nplus = min0(nplus*2,max0(npl1,nplus/2,1))
        end if
        fpold = fp
!  compute the sum((w(i)*(y(i)-s(x(i))))**2) for each knot interval
!  t(j+k) <= x(i) <= t(j+k+1) and store it in fpint(j),j=1,2,...nrint.
        fpart = 0.0_wp_
        i = 1
        l = k2
        new = 0
        do it=1,m
          if(x(it)>=t(l) .and. l<=nk1) then
            new = 1
            l = l+1
          end if
          term = 0.0_wp_
          l0 = l-k2
          do j=1,k1
            l0 = l0+1
            term = term+c(l0)*q(it,j)
          end do
          term = (w(it)*(term-y(it)))**2
          fpart = fpart+term
          if(new==0) cycle
          store = term*half
          fpint(i) = fpart-store
          i = i+1
          fpart = store
          new = 0
        end do
        fpint(nrint) = fpart
        do l=1,nplus
!  add a new knot.
          call fpknot(x,m,t,n,fpint,nrdata,nrint,nest,1)
!  if n=nmax we locate the knots as for interpolation.
          if(n==nmax) then
            rstart=.true.
            exit
          end if
!  test whether we cannot further increase the number of knots.
          if(n==nest) exit
        end do
        if(rstart) exit
!  restart the computations with the new set of knots.
      end do
      if(rstart) then
!  find the position of the interior knots in case of interpolation.
        mk1 = m-k1
        if (mk1/=0) then
          k3 = k/2
          i = k2
          j = k3+2
          if(k3*2/=k) then
            do l=1,mk1
              t(i) = x(j)
              i = i+1
              j = j+1
            end do
          else
            do l=1,mk1
              t(i) = (x(j)+x(j-1))*half
              i = i+1
              j = j+1
            end do
          end if
        end if
      else
        exit
      end if
    end do
!  test whether the least-squares kth degree polynomial is a solution
!  of our approximation problem.
    if(ier==(-2)) return
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!  part 2: determination of the smoothing spline sp(x).                c
!  ***************************************************                 c
!  we have determined the number of knots and their position.          c
!  we now compute the b-spline coefficients of the smoothing spline    c
!  sp(x). the observation matrix a is extended by the rows of matrix   c
!  b expressing that the kth derivative discontinuities of sp(x) at    c
!  the interior knots t(k+2),...t(n-k-1) must be zero. the corres-     c
!  ponding weights of these additional rows are set to 1/p.            c
!  iteratively we then have to determine the value of p such that      c
!  f(p)=sum((w(i)*(y(i)-sp(x(i))))**2) be = s. we already know that    c
!  the least-squares kth degree polynomial corresponds to p=0, and     c
!  that the least-squares spline corresponds to p=infinity. the        c
!  iteration process which is proposed here, makes use of rational     c
!  interpolation. since f(p) is a convex and strictly decreasing       c
!  function of p, it can be approximated by a rational function        c
!  r(p) = (u*p+v)/(p+w). three values of p(p1,p2,p3) with correspond-  c
!  ing values of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used      c
!  to calculate the new value of p such that r(p)=s. convergence is    c
!  guaranteed by taking f1>0 and f3<0.                                 c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!  evaluate the discontinuity jump of the kth derivative of the
!  b-splines at the knots t(l),l=k+2,...n-k-1 and store in b.
    call fpdisc(t,n,k2,b,nest)
!  initial value for p.
    p1 = 0.0_wp_
    f1 = fp0-s
    p3 = -one
    f3 = fpms
    p = 0.0_wp_
    do i=1,nk1
      p = p+a(i,1)
    end do
    rn = nk1
    p = rn/p
    ich1 = 0
    ich3 = 0
    n8 = n-nmin
!  iteration process to find the root of f(p) = s.
    do iter=1,maxit
!  the rows of matrix b with weight 1/p are rotated into the
!  triangularised observation matrix a which is stored in g.
      pinv = one/p
      do i=1,nk1
        c(i) = z(i)
        g(i,k2) = 0.0_wp_
        do j=1,k1
          g(i,j) = a(i,j)
        end do
      end do
      do it=1,n8
!  the row of matrix b is rotated into triangle by givens transformation
        do i=1,k2
          h(i) = b(it,i)*pinv
        end do
        yi = 0.0_wp_
        do j=it,nk1
          piv = h(1)
!  calculate the parameters of the givens transformation.
          call fpgivs(piv,g(j,1),cs,sn)
!  transformations to right hand side.
          call fprota(cs,sn,yi,c(j))
          if(j==nk1) exit
          i2 = k1
          if(j>n8) i2 = nk1-j
          do i=1,i2
!  transformations to left hand side.
            i1 = i+1
            call fprota(cs,sn,h(i1),g(j,i1))
            h(i) = h(i1)
          end do
          h(i2+1) = 0.0_wp_
        end do
      end do
!  backward substitution to obtain the b-spline coefficients.
      call fpback(g,c,nk1,k2,c,nest)
!  computation of f(p).
      fp = 0.0_wp_
      l = k2
      do it=1,m
        if(x(it)>=t(l) .and. l<=nk1) l = l+1
        l0 = l-k2
        term = 0.0_wp_
        do j=1,k1
          l0 = l0+1
          term = term+c(l0)*q(it,j)
        end do
        fp = fp+(w(it)*(term-y(it)))**2
      end do
!  test whether the approximation sp(x) is an acceptable solution.
      fpms = fp-s
      if(abs(fpms)<acc) return
!  test whether the maximal number of iterations is reached.
      if(iter==maxit) then
        ier = 3
        return
      end if
!  carry out one more step of the iteration process.
      p2 = p
      f2 = fpms
      if(ich3==0) then
        if((f2-f3)<=acc) then
!  our initial choice of p is too large.
          p3 = p2
          f3 = f2
          p = p*con4
          if(p<=p1) p=p1*con9 + p2*con1
          cycle
        end if
        if(f2<0.0_wp_) ich3=1
      end if
      if(ich1==0) then
        if((f1-f2)<=acc) then
!  our initial choice of p is too small
          p1 = p2
          f1 = f2
          p = p/con4
          if(p3>=0.0_wp_ .and. p>=p3) p = p2*con1 + p3*con9
          cycle
        end if
        if(f2>0.0_wp_) ich1=1
      end if
!  test whether the iteration process proceeds as theoretically
!  expected.
      if(f2>=f1 .or. f2<=f3) then
        ier = 2
        return
      end if
!  find the new value for p.
      p = fprati(p1,f1,p2,f2,p3,f3)
    end do
  end subroutine fpcurf

  subroutine splder(t,n,c,k,nu,x,y,m,wrk,ier)
!  subroutine splder evaluates in a number of points x(i),i=1,2,...,m
!  the derivative of order nu of a spline s(x) of degree k,given in
!  its b-spline representation.
!
!  calling sequence:
! call splder(t,n,c,k,nu,x,y,m,wrk,ier)
!
!  input parameters:
!    t    : array,length n, which contains the position of the knots.
!    n    : integer, giving the total number of knots of s(x).
!    c    : array,length n, which contains the b-spline coefficients.
!    k    : integer, giving the degree of s(x).
!    nu   : integer, specifying the order of the derivative. 0<=nu<=k
!    x    : array,length m, which contains the points where the deriv-
!       ative of s(x) must be evaluated.
!    m    : integer, giving the number of points where the derivative
!       of s(x) must be evaluated
!    wrk  : real array of dimension n. used as working space.
!
!  output parameters:
!    y    : array,length m, giving the value of the derivative of s(x)
!       at the different points.
!    ier  : error flag
!  ier = 0 : normal return
!  ier =10 : invalid input data (see restrictions)
!
!  restrictions:
!    0 <= nu <= k
!    m >= 1
!    t(k+1) <= x(i) <= x(i+1) <= t(n-k) , i=1,2,...,m-1.
!
!  other subroutines required: fpbspl
!
!  references :
!    de boor c : on calculating with b-splines, j. approximation theory
!            6 (1972) 50-62.
!    cox m.g.  : the numerical evaluation of b-splines, j. inst. maths
!            applics 10 (1972) 134-149.
!   dierckx p. : curve and surface fitting with splines, monographs on
!            numerical analysis, oxford university press, 1993.
!
!  author :
!    p.dierckx
!    dept. computer science, k.u.leuven
!    celestijnenlaan 200a, b-3001 heverlee, belgium.
!    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
!  latest update : march 1987
!
    implicit none
! arguments
    integer, intent(in) :: n, k, nu, m
    integer, intent(out) :: ier
    real(wp_), intent(in) :: t(n), c(n), x(m)
    real(wp_), intent(out) :: y(m), wrk(n)
! local variables
    integer :: i, j, kk, k1, k2, l, ll, l1, l2, nk1, nk2
    real(wp_) :: ak, arg, fac, sp, tb, te, h(6)
!  before starting computations a data check is made. if the input data
!  are invalid control is immediately repassed to the calling program.
    ier = 10
    if(nu<0 .or. nu>k) return
    if(m<1) return
    do i=2,m
      if(x(i)<x(i-1)) return
    end do
    ier = 0
!  fetch tb and te, the boundaries of the approximation interval.
    k1 = k+1
    nk1 = n-k1
    tb = t(k1)
    te = t(nk1+1)
!  the derivative of order nu of a spline of degree k is a spline of
!  degree k-nu,the b-spline coefficients wrk(i) of which can be found
!  using the recurrence scheme of de boor.
    l = 1
    kk = k
!  nn = n
    do i=1,nk1
      wrk(i) = c(i)
    end do
    if(nu/=0) then
      nk2 = nk1
      do j=1,nu
        ak = kk
        nk2 = nk2-1
        l1 = l
        do i=1,nk2
          l1 = l1+1
          l2 = l1+kk
          fac = t(l2)-t(l1)
          if(fac>0.0_wp_) wrk(i) = ak*(wrk(i+1)-wrk(i))/fac
        end do
        l = l+1
        kk = kk-1
      end do
      if(kk==0) then
!  if nu=k the derivative is a piecewise constant function
        j = 1
        do i=1,m
          arg = x(i)
          do
            if(arg<t(l+1) .or. l==nk1) exit
            l = l+1
            j = j+1
          end do
          y(i) = wrk(j)
        end do
        return
      end if
    end if
    l = k1
    l1 = l+1
    k2 = k1-nu
!  main loop for the different points.
    do i=1,m
!  fetch a new x-value arg.
      arg = x(i)
      if(arg<tb) arg = tb
      if(arg>te) arg = te
!  search for knot interval t(l) <= arg < t(l+1)
      do
        if(arg<t(l1) .or. l==nk1) exit
        l = l1
        l1 = l+1
      end do
!  evaluate the non-zero b-splines of degree k-nu at arg.
      call fpbspl(t,n,kk,arg,l,h)
!  find the value of the derivative at x=arg.
      sp = 0.0_wp_
      ll = l-k1
      do j=1,k2
        ll = ll+1
        sp = sp+wrk(ll)*h(j)
      end do
      y(i) = sp
    end do
  end subroutine splder

  subroutine splev(t,n,c,k,x,y,m,ier)
!  subroutine splev evaluates in a number of points x(i),i=1,2,...,m
!  a spline s(x) of degree k, given in its b-spline representation.
!
!  calling sequence:
! call splev(t,n,c,k,x,y,m,ier)
!
!  input parameters:
!    t    : array,length n, which contains the position of the knots.
!    n    : integer, giving the total number of knots of s(x).
!    c    : array,length n, which contains the b-spline coefficients.
!    k    : integer, giving the degree of s(x).
!    x    : array,length m, which contains the points where s(x) must
!       be evaluated.
!    m    : integer, giving the number of points where s(x) must be
!       evaluated.
!
!  output parameter:
!    y    : array,length m, giving the value of s(x) at the different
!       points.
!    ier  : error flag
!  ier = 0 : normal return
!  ier =10 : invalid input data (see restrictions)
!
!  restrictions:
!    m >= 1
!    t(k+1) <= x(i) <= x(i+1) <= t(n-k) , i=1,2,...,m-1.
!
!  other subroutines required: fpbspl.
!
!  references :
!    de boor c  : on calculating with b-splines, j. approximation theory
!             6 (1972) 50-62.
!    cox m.g.   : the numerical evaluation of b-splines, j. inst. maths
!             applics 10 (1972) 134-149.
!    dierckx p. : curve and surface fitting with splines, monographs on
!             numerical analysis, oxford university press, 1993.
!
!  author :
!    p.dierckx
!    dept. computer science, k.u.leuven
!    celestijnenlaan 200a, b-3001 heverlee, belgium.
!    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
!  latest update : march 1987
!
    implicit none
! arguments
    integer, intent(in) :: n, k, m
    integer, intent(out) :: ier
    real(wp_), intent(in) :: t(n), c(n), x(m)
    real(wp_), intent(out) :: y(m)
! local variables
    integer       :: i, j, k1, l, ll, l1, nk1
    real(wp_) :: arg, sp, tb, te, h(6)
!  ..
!  before starting computations a data check is made. if the input data
!  are invalid control is immediately repassed to the calling program.
    ier = 10
    if(m<1) return
    do i=2,m
      if(x(i)<x(i-1)) return
    end do
    ier = 0
!  fetch tb and te, the boundaries of the approximation interval.
    k1 = k+1
    nk1 = n-k1
    tb = t(k1)
    te = t(nk1+1)
    l = k1
    l1 = l+1
!  main loop for the different points.
    do i=1,m
!  fetch a new x-value arg.
      arg = x(i)
      if(arg<tb) arg = tb
      if(arg>te) arg = te
!  search for knot interval t(l) <= arg < t(l+1)
      do
        if(arg<t(l1) .or. l==nk1) exit
        l = l1
        l1 = l+1
      end do
!  evaluate the non-zero b-splines at arg.
      call fpbspl(t,n,k,arg,l,h)
!  find the value of s(x) at x=arg.
      sp = 0.0_wp_
      ll = l-k1
      do j=1,k1
        ll = ll+1
        sp = sp+c(ll)*h(j)
      end do
      y(i) = sp
    end do
  end subroutine splev

  subroutine sproota(val,t,n,c,zero,mest,m,ier)
!  subroutine sproot finds the zeros of a cubic spline s(x),which is
!  given in its normalized b-spline representation.
!
!  calling sequence:
! call sproot(t,n,c,zero,mest,m,ier)
!
!  input parameters:
!    t    : real array,length n, containing the knots of s(x).
!    n    : integer, containing the number of knots.  n>=8
!    c    : real array,length n, containing the b-spline coefficients.
!    mest : integer, specifying the dimension of array zero.
!
!  output parameters:
!    zero : real array,lenth mest, containing the zeros of s(x).
!    m    : integer,giving the number of zeros.
!    ier  : error flag:
!  ier = 0: normal return.
!  ier = 1: the number of zeros exceeds mest.
!  ier =10: invalid input data (see restrictions).
!
!  other subroutines required: fpcuro
!
!  restrictions:
!    1) n>= 8.
!    2) t(4) < t(5) < ... < t(n-4) < t(n-3).
!   t(1) <= t(2) <= t(3) <= t(4)
!   t(n-3) <= t(n-2) <= t(n-1) <= t(n)
!
!  author :
!    p.dierckx
!    dept. computer science, k.u.leuven
!    celestijnenlaan 200a, b-3001 heverlee, belgium.
!    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
!  latest update : march 1987
!
! ..
    implicit none
! ..scalar arguments..
    integer, intent(in) :: n,mest
    integer, intent(out) :: m,ier
    real(wp_), intent(in) :: val
!  ..array arguments..
    real(wp_), intent(in) :: t(n),c(n)
    real(wp_), intent(out) :: zero(mest)
!  ..local scalars..
    integer :: i,j,j1,l,n4
    real(wp_) :: ah,a0,a1,a2,a3,bh,b0,b1,c1,c2,c3,c4,c5,d4,d5,h1,h2, &
      t1,t2,t3,t4,t5,zz
    logical :: z0,z1,z2,z3,z4,nz0,nz1,nz2,nz3,nz4
!  ..local array..
    real(wp_) :: y(3)
!  ..
!  set some constants
    real(wp_), parameter :: two = 0.2e+01_wp_, three = 0.3e+01_wp_
!  before starting computations a data check is made. if the input data
!  are invalid, control is immediately repassed to the calling program.
    n4 = n-4
    ier = 10
    if(n<8) return
    j = n
    do i=1,3
      if(t(i)>t(i+1)) return
      if(t(j)<t(j-1)) return
      j = j-1
    end do
    do i=4,n4
      if(t(i)>=t(i+1)) return
    end do
!  the problem considered reduces to finding the zeros of the cubic
!  polynomials pl(x) which define the cubic spline in each knot
!  interval t(l)<=x<=t(l+1). a zero of pl(x) is also a zero of s(x) on
!  the condition that it belongs to the knot interval.
!  the cubic polynomial pl(x) is determined by computing s(t(l)),
!  s'(t(l)),s(t(l+1)) and s'(t(l+1)). in fact we only have to compute
!  s(t(l+1)) and s'(t(l+1)); because of the continuity conditions of
!  splines and their derivatives, the value of s(t(l)) and s'(t(l))
!  is already known from the foregoing knot interval.
    ier = 0
!  evaluate some constants for the first knot interval
    h1 = t(4)-t(3)
    h2 = t(5)-t(4)
    t1 = t(4)-t(2)
    t2 = t(5)-t(3)
    t3 = t(6)-t(4)
    t4 = t(5)-t(2)
    t5 = t(6)-t(3)
!  calculate a0 = s(t(4)) and ah = s'(t(4)).
    c1 = c(1)
    c2 = c(2)
    c3 = c(3)
    c4 = (c2-c1)/t4
    c5 = (c3-c2)/t5
    d4 = (h2*c1+t1*c2)/t4
    d5 = (t3*c2+h1*c3)/t5
    a0 = (h2*d4+h1*d5)/t2 - val
    ah = three*(h2*c4+h1*c5)/t2
    z1 = .true.
    if(ah<0.0_wp_) z1 = .false.
    nz1 = .not.z1
    m = 0
!  main loop for the different knot intervals.
    do l=4,n4
!  evaluate some constants for the knot interval t(l) <= x <= t(l+1).
      h1 = h2
      h2 = t(l+2)-t(l+1)
      t1 = t2
      t2 = t3
      t3 = t(l+3)-t(l+1)
      t4 = t5
      t5 = t(l+3)-t(l)
!  find a0 = s(t(l)), ah = s'(t(l)), b0 = s(t(l+1)) and bh = s'(t(l+1)).
      c1 = c2
      c2 = c3
      c3 = c(l)
      c4 = c5
      c5 = (c3-c2)/t5
      d4 = (h2*c1+t1*c2)/t4
      d5 = (h1*c3+t3*c2)/t5
      b0 = (h2*d4+h1*d5)/t2 - val
      bh = three*(h2*c4+h1*c5)/t2
!  calculate the coefficients a0,a1,a2 and a3 of the cubic polynomial
!  pl(x) = ql(y) = a0+a1*y+a2*y**2+a3*y**3 ; y = (x-t(l))/(t(l+1)-t(l)).
      a1 = ah*h1
      b1 = bh*h1
      a2 = three*(b0-a0)-b1-two*a1
      a3 = two*(a0-b0)+b1+a1
!  test whether or not pl(x) could have a zero in the range
!  t(l) <= x <= t(l+1).
      z3 = .true.
      if(b1<0.0_wp_) z3 = .false.
      nz3 = .not.z3
      if(a0*b0>0.0_wp_) then
        z0 = .true.
        if(a0<0.0_wp_) z0 = .false.
        nz0 = .not.z0
        z2 = .true.
        if(a2<0.0_wp_) z2 = .false.
        nz2 = .not.z2
        z4 = .true.
        if(3.0_wp_*a3+a2<0.0_wp_) z4 = .false.
        nz4 = .not.z4
      else
        z0 = .true.
        nz0 = .not.z0
        z2 = .true.
        nz2 = .not.z2
        z4 = .true.
        nz4 = .not.z4
      end if
      if(( z0.and.(nz1.and.( z3.or. z2.and.nz4).or.nz2.and. z3.and. z4) &
      .or.nz0.and.( z1.and.(nz3.or.nz2.and. z4).or. z2.and.nz3.and.nz4) &
         ) .or. (a0*b0<=0.0_wp_) ) then
!  find the zeros of ql(y).
        call fpcuro(a3,a2,a1,a0,y,j)
!  find which zeros of pl(x) are zeros of s(x).
        do i=1,j
          if(y(i)<0.0_wp_ .or. y(i)>1.0_wp_) cycle
!  test whether the number of zeros of s(x) exceeds mest.
          if(m>=mest) then
            ier = 1
            return
          end if
          m = m+1
          zero(m) = t(l)+h1*y(i)
        end do
      end if
      a0 = b0
      ah = bh
      z1 = z3
      nz1 = nz3
    end do
!  the zeros of s(x) are arranged in increasing order.
    do i=2,m
      j = i
      do
        j1 = j-1
        if(j1==0) exit
        if(zero(j)>=zero(j1)) exit
        zz = zero(j)
        zero(j) = zero(j1)
        zero(j1) = zz
        j = j1
      end do
    end do
    j = m
    m = 1
    do i=2,j
      if(zero(i)==zero(m)) cycle
      m = m+1
      zero(m) = zero(i)
    end do
  end subroutine sproota

  subroutine profil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
!  if iopt=0 subroutine profil calculates the b-spline coefficients of
!  the univariate spline f(y) = s(u,y) with s(x,y) a bivariate spline of
!  degrees kx and ky, given in the b-spline representation.
!  if iopt = 1 it calculates the b-spline coefficients of the univariate
!  spline g(x) = s(x,u)
!
!  calling sequence:
! call profil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
!
!  input parameters:
!   iopt  : integer flag, specifying whether the profile f(y) (iopt=0)
!       or the profile g(x) (iopt=1) must be determined.
!   tx    : real array, length nx, which contains the position of the
!       knots in the x-direction.
!   nx    : integer, giving the total number of knots in the x-direction
!   ty    : real array, length ny, which contains the position of the
!       knots in the y-direction.
!   ny    : integer, giving the total number of knots in the y-direction
!   c     : real array, length (nx-kx-1)*(ny-ky-1), which contains the
!       b-spline coefficients.
!   kx,ky : integer values, giving the degrees of the spline.
!   u     : real value, specifying the requested profile.
!       tx(kx+1)<=u<=tx(nx-kx), if iopt=0.
!       ty(ky+1)<=u<=ty(ny-ky), if iopt=1.
!   nu    : on entry nu must specify the dimension of the array cu.
!       nu >= ny if iopt=0, nu >= nx if iopt=1.
!
!  output parameters:
!   cu    : real array of dimension (nu).
!       on succesful exit this array contains the b-spline
!   ier   : integer error flag
!    ier=0 : normal return
!    ier=10: invalid input data (see restrictions)
!
!  restrictions:
!   if iopt=0 : tx(kx+1) <= u <= tx(nx-kx), nu >=ny.
!   if iopt=1 : ty(ky+1) <= u <= ty(ny-ky), nu >=nx.
!
!  other subroutines required:
!    fpbspl
!
!  author :
!    p.dierckx
!    dept. computer science, k.u.leuven
!    celestijnenlaan 200a, b-3001 heverlee, belgium.
!    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
!
!  latest update : march 1987
!
    implicit none
!  ..scalar arguments..
    integer,intent(in) :: iopt,nx,ny,kx,ky,nu
    integer,intent(out) :: ier
    real(wp_), intent(in) :: u
!  ..array arguments..
    real(wp_),intent(in) :: tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1))
    real(wp_),intent(out) :: cu(nu)
!  ..local scalars..
    integer :: i,j,kx1,ky1,l,l1,m,m0,nkx1,nky1
    real(wp_) :: summ
!  ..local array
    real(wp_) :: h(6)
!  ..
!  before starting computations a data check is made. if the input data
!  are invalid control is immediately repassed to the calling program.
    kx1 = kx+1
    ky1 = ky+1
    nkx1 = nx-kx1
    nky1 = ny-ky1
    ier = 10
    if(iopt==0) then
      if(nu<ny) return
      if(u<tx(kx1) .or. u>tx(nkx1+1)) return
!  the b-splinecoefficients of f(y) = s(u,y).
      ier = 0
      l = kx1
      l1 = l+1
      do
        if(u<tx(l1) .or. l==nkx1) exit
        l = l1
        l1 = l+1
      end do
      call fpbspl(tx,nx,kx,u,l,h)
      m0 = (l-kx1)*nky1+1
      do i=1,nky1
        m = m0
        summ = 0.0_wp_
        do j=1,kx1
          summ = summ+h(j)*c(m)
          m = m+nky1
        end do
        cu(i) = summ
        m0 = m0+1
      end do
    else
      if(nu<nx) return
      if(u<ty(ky1) .or. u>ty(nky1+1)) return
!  the b-splinecoefficients of g(x) = s(x,u).
      ier = 0
      l = ky1
      l1 = l+1
      do
        if(u<ty(l1) .or. l==nky1) exit
        l = l1
        l1 = l+1
      end do
      call fpbspl(ty,ny,ky,u,l,h)
      m0 = l-ky
      do i=1,nkx1
        m = m0
        summ = 0.0_wp_
        do j=1,ky1
          summ = summ+h(j)*c(m)
          m = m+1
        end do
        cu(i) = summ
        m0 = m0+nky1
      end do
    end if
  end subroutine profil

  subroutine fpcuro(a,b,c,d,x,n)
!  subroutine fpcuro finds the real zeros of a cubic polynomial
!  p(x) = a*x**3+b*x**2+c*x+d.
!
!  calling sequence:
! call fpcuro(a,b,c,d,x,n)
!
!  input parameters:
!    a,b,c,d: real values, containing the coefficients of p(x).
!
!  output parameters:
!    x      : real array,length 3, which contains the real zeros of p(x)
!    n      : integer, giving the number of real zeros of p(x).
!  ..
    implicit none
!  ..scalar arguments..
    real(wp_), intent(in) :: a,b,c,d
    integer, intent(out) :: n
!  ..array argument..
    real(wp_), intent(out) :: x(3)
!  ..local scalars..
    integer :: i
    real(wp_) :: a1,b1,c1,df,disc,d1,f,pi3,p3,q,r,step,u,u1,u2,y
!  ..function references..
!  real(8) abs,max,atan,atan2,cos,sign,sqrt
!  set constants
    real(wp_), parameter :: two =0.2e+01_wp_,three=0.3e+01_wp_,four=0.4e+01_wp_
    real(wp_), parameter :: ovfl=0.1e+05_wp_,half =0.5e+00_wp_,tent=0.1e+00_wp_
    real(wp_), parameter :: e3 = tent/0.3_wp_
    pi3 = atan(0.1e+01_wp_)/0.75_wp_
    a1 = abs(a)
    b1 = abs(b)
    c1 = abs(c)
    d1 = abs(d)
!  test whether p(x) is a third degree polynomial.
    if(max(b1,c1,d1)>=a1*ovfl) then
!  test whether p(x) is a second degree polynomial.
      if(max(c1,d1)>=b1*ovfl) then
!  test whether p(x) is a first degree polynomial.
        if(d1>=c1*ovfl) then
!  p(x) is a constant function.
          n = 0
          return
        end if
!  p(x) is a first degree polynomial.
        n = 1
        x(1) = -d/c
      else
!  p(x) is a second degree polynomial.
        disc = c*c-four*b*d
        n = 0
        if(disc<0.0_wp_) return
        n = 2
        u = sqrt(disc)
        b1 = b+b
        x(1) = (-c+u)/b1
        x(2) = (-c-u)/b1
      end if
    else
!  p(x) is a third degree polynomial.
      b1 = b/a*e3
      c1 = c/a
      d1 = d/a
      q = c1*e3-b1*b1
      r = b1*b1*b1+(d1-b1*c1)*half
      disc = q*q*q+r*r
      if(disc<=0.0_wp_) then
        u = sqrt(abs(q))
        if(r<0.0_wp_) u = -u
        p3 = atan2(sqrt(-disc),abs(r))*e3
        u2 = u+u
        n = 3
        x(1) = -u2*cos(p3)-b1
        x(2) = u2*cos(pi3-p3)-b1
        x(3) = u2*cos(pi3+p3)-b1
      else
        u = sqrt(disc)
        u1 = -r+u
        u2 = -r-u
        n = 1
        x(1) = sign(abs(u1)**e3,u1)+sign(abs(u2)**e3,u2)-b1
      end if
    end if
!  apply a newton iteration to improve the accuracy of the roots.
    do i=1,n
      y = x(i)
      f = ((a*y+b)*y+c)*y+d
      df = (three*a*y+two*b)*y+c
      step = 0.0_wp_
      if(abs(f)<abs(df)*tent) step = f/df
      x(i) = y-step
    end do
  end subroutine fpcuro

end module dierckx