1      subroutine spgrid(iopt,ider,mu,u,mv,v,r,r0,r1,s,nuest,nvest,
2     * nu,tu,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
3c  given the function values r(i,j) on the latitude-longitude grid
4c  (u(i),v(j)), i=1,...,mu ; j=1,...,mv , spgrid determines a smooth
5c  bicubic spline approximation on the rectangular domain 0<=u<=pi,
6c  vb<=v<=ve (vb = v(1), ve=vb+2*pi).
7c  this approximation s(u,v) will satisfy the properties
8c
9c    (1) s(0,v) = s(0,0) = dr(1)
10c
11c        d s(0,v)           d s(0,0)           d s(0,pi/2)
12c    (2) -------- = cos(v)* -------- + sin(v)* -----------
13c        d u                d u                d u
14c
15c                 = cos(v)*dr(2)+sin(v)*dr(3)
16c                                                     vb <= v <= ve
17c    (3) s(pi,v) = s(pi,0) = dr(4)
18c
19c        d s(pi,v)           d s(pi,0)           d s(pi,pi/2)
20c    (4) -------- = cos(v)*  --------- + sin(v)* ------------
21c        d u                 d u                 d u
22c
23c                 = cos(v)*dr(5)+sin(v)*dr(6)
24c
25c  and will be periodic in the variable v, i.e.
26c
27c         j           j
28c        d s(u,vb)   d s(u,ve)
29c    (5) --------- = ---------   0 <=u<= pi , j=0,1,2
30c           j           j
31c        d v         d v
32c
33c  the number of knots of s(u,v) and their position tu(i),i=1,2,...,nu;
34c  tv(j),j=1,2,...,nv, is chosen automatically by the routine. the
35c  smoothness of s(u,v) is achieved by minimalizing the discontinuity
36c  jumps of the derivatives of the spline at the knots. the amount of
37c  smoothness of s(u,v) is determined by the condition that
38c  fp=sumi=1,mu(sumj=1,mv((r(i,j)-s(u(i),v(j)))**2))+(r0-s(0,v))**2
39c  + (r1-s(pi,v))**2 <= s, with s a given non-negative constant.
40c  the fit s(u,v) is given in its b-spline representation and can be
41c  evaluated by means of routine bispev
42c
43c calling sequence:
44c     call spgrid(iopt,ider,mu,u,mv,v,r,r0,r1,s,nuest,nvest,nu,tu,
45c    *  ,nv,tv,c,fp,wrk,lwrk,iwrk,kwrk,ier)
46c
47c parameters:
48c  iopt  : integer array of dimension 3, specifying different options.
49c          unchanged on exit.
50c  iopt(1):on entry iopt(1) must specify whether a least-squares spline
51c          (iopt(1)=-1) or a smoothing spline (iopt(1)=0 or 1) must be
52c          determined.
53c          if iopt(1)=0 the routine will start with an initial set of
54c          knots tu(i)=0,tu(i+4)=pi,i=1,...,4;tv(i)=v(1)+(i-4)*2*pi,
55c          i=1,...,8.
56c          if iopt(1)=1 the routine will continue with the set of knots
57c          found at the last call of the routine.
58c          attention: a call with iopt(1)=1 must always be immediately
59c          preceded by another call with iopt(1) = 1 or iopt(1) = 0.
60c  iopt(2):on entry iopt(2) must specify the requested order of conti-
61c          nuity at the pole u=0.
62c          if iopt(2)=0 only condition (1) must be fulfilled and
63c          if iopt(2)=1 conditions (1)+(2) must be fulfilled.
64c  iopt(3):on entry iopt(3) must specify the requested order of conti-
65c          nuity at the pole u=pi.
66c          if iopt(3)=0 only condition (3) must be fulfilled and
67c          if iopt(3)=1 conditions (3)+(4) must be fulfilled.
68c  ider  : integer array of dimension 4, specifying different options.
69c          unchanged on exit.
70c  ider(1):on entry ider(1) must specify whether (ider(1)=0 or 1) or not
71c          (ider(1)=-1) there is a data value r0 at the pole u=0.
72c          if ider(1)=1, r0 will be considered to be the right function
73c          value, and it will be fitted exactly (s(0,v)=r0).
74c          if ider(1)=0, r0 will be considered to be a data value just
75c          like the other data values r(i,j).
76c  ider(2):on entry ider(2) must specify whether (ider(2)=1) or not
77c          (ider(2)=0) the approximation has vanishing derivatives
78c          dr(2) and dr(3) at the pole u=0  (in case iopt(2)=1)
79c  ider(3):on entry ider(3) must specify whether (ider(3)=0 or 1) or not
80c          (ider(3)=-1) there is a data value r1 at the pole u=pi.
81c          if ider(3)=1, r1 will be considered to be the right function
82c          value, and it will be fitted exactly (s(pi,v)=r1).
83c          if ider(3)=0, r1 will be considered to be a data value just
84c          like the other data values r(i,j).
85c  ider(4):on entry ider(4) must specify whether (ider(4)=1) or not
86c          (ider(4)=0) the approximation has vanishing derivatives
87c          dr(5) and dr(6) at the pole u=pi (in case iopt(3)=1)
88c  mu    : integer. on entry mu must specify the number of grid points
89c          along the u-axis. unchanged on exit.
90c          mu >= 1, mu >=mumin=4-i0-i1-ider(2)-ider(4) with
91c            i0=min(1,ider(1)+1), i1=min(1,ider(3)+1)
92c  u     : real array of dimension at least (mu). before entry, u(i)
93c          must be set to the u-co-ordinate of the i-th grid point
94c          along the u-axis, for i=1,2,...,mu. these values must be
95c          supplied in strictly ascending order. unchanged on exit.
96c          0 < u(i) < pi.
97c  mv    : integer. on entry mv must specify the number of grid points
98c          along the v-axis. mv > 3 . unchanged on exit.
99c  v     : real array of dimension at least (mv). before entry, v(j)
100c          must be set to the v-co-ordinate of the j-th grid point
101c          along the v-axis, for j=1,2,...,mv. these values must be
102c          supplied in strictly ascending order. unchanged on exit.
103c          -pi <= v(1) < pi , v(mv) < v(1)+2*pi.
104c  r     : real array of dimension at least (mu*mv).
105c          before entry, r(mv*(i-1)+j) must be set to the data value at
106c          the grid point (u(i),v(j)) for i=1,...,mu and j=1,...,mv.
107c          unchanged on exit.
108c  r0    : real value. on entry (if ider(1) >=0 ) r0 must specify the
109c          data value at the pole u=0. unchanged on exit.
110c  r1    : real value. on entry (if ider(1) >=0 ) r1 must specify the
111c          data value at the pole u=pi. unchanged on exit.
112c  s     : real. on entry (if iopt(1)>=0) s must specify the smoothing
113c          factor. s >=0. unchanged on exit.
114c          for advice on the choice of s see further comments
115c  nuest : integer. unchanged on exit.
116c  nvest : integer. unchanged on exit.
117c          on entry, nuest and nvest must specify an upper bound for the
118c          number of knots required in the u- and v-directions respect.
119c          these numbers will also determine the storage space needed by
120c          the routine. nuest >= 8, nvest >= 8.
121c          in most practical situation nuest = mu/2, nvest=mv/2, will
122c          be sufficient. always large enough are nuest=mu+6+iopt(2)+
123c          iopt(3), nvest = mv+7, the number of knots needed for
124c          interpolation (s=0). see also further comments.
125c  nu    : integer.
126c          unless ier=10 (in case iopt(1)>=0), nu will contain the total
127c          number of knots with respect to the u-variable, of the spline
128c          approximation returned. if the computation mode iopt(1)=1 is
129c          used, the value of nu should be left unchanged between sub-
130c          sequent calls. in case iopt(1)=-1, the value of nu should be
131c          specified on entry.
132c  tu    : real array of dimension at least (nuest).
133c          on successful exit, this array will contain the knots of the
134c          spline with respect to the u-variable, i.e. the position of
135c          the interior knots tu(5),...,tu(nu-4) as well as the position
136c          of the additional knots tu(1)=...=tu(4)=0 and tu(nu-3)=...=
137c          tu(nu)=pi needed for the b-spline representation.
138c          if the computation mode iopt(1)=1 is used,the values of tu(1)
139c          ...,tu(nu) should be left unchanged between subsequent calls.
140c          if the computation mode iopt(1)=-1 is used, the values tu(5),
141c          ...tu(nu-4) must be supplied by the user, before entry.
142c          see also the restrictions (ier=10).
143c  nv    : integer.
144c          unless ier=10 (in case iopt(1)>=0), nv will contain the total
145c          number of knots with respect to the v-variable, of the spline
146c          approximation returned. if the computation mode iopt(1)=1 is
147c          used, the value of nv should be left unchanged between sub-
148c          sequent calls. in case iopt(1) = -1, the value of nv should
149c          be specified on entry.
150c  tv    : real array of dimension at least (nvest).
151c          on successful exit, this array will contain the knots of the
152c          spline with respect to the v-variable, i.e. the position of
153c          the interior knots tv(5),...,tv(nv-4) as well as the position
154c          of the additional knots tv(1),...,tv(4) and tv(nv-3),...,
155c          tv(nv) needed for the b-spline representation.
156c          if the computation mode iopt(1)=1 is used,the values of tv(1)
157c          ...,tv(nv) should be left unchanged between subsequent calls.
158c          if the computation mode iopt(1)=-1 is used, the values tv(5),
159c          ...tv(nv-4) must be supplied by the user, before entry.
160c          see also the restrictions (ier=10).
161c  c     : real array of dimension at least (nuest-4)*(nvest-4).
162c          on successful exit, c contains the coefficients of the spline
163c          approximation s(u,v)
164c  fp    : real. unless ier=10, fp contains the sum of squared
165c          residuals of the spline approximation returned.
166c  wrk   : real array of dimension (lwrk). used as workspace.
167c          if the computation mode iopt(1)=1 is used the values of
168c          wrk(1),..,wrk(12) should be left unchanged between subsequent
169c          calls.
170c  lwrk  : integer. on entry lwrk must specify the actual dimension of
171c          the array wrk as declared in the calling (sub)program.
172c          lwrk must not be too small.
173c           lwrk >= 12+nuest*(mv+nvest+3)+nvest*24+4*mu+8*mv+q
174c           where q is the larger of (mv+nvest) and nuest.
175c  iwrk  : integer array of dimension (kwrk). used as workspace.
176c          if the computation mode iopt(1)=1 is used the values of
177c          iwrk(1),.,iwrk(5) should be left unchanged between subsequent
178c          calls.
179c  kwrk  : integer. on entry kwrk must specify the actual dimension of
180c          the array iwrk as declared in the calling (sub)program.
181c          kwrk >= 5+mu+mv+nuest+nvest.
182c  ier   : integer. unless the routine detects an error, ier contains a
183c          non-positive value on exit, i.e.
184c   ier=0  : normal return. the spline returned has a residual sum of
185c            squares fp such that abs(fp-s)/s <= tol with tol a relat-
186c            ive tolerance set to 0.001 by the program.
187c   ier=-1 : normal return. the spline returned is an interpolating
188c            spline (fp=0).
189c   ier=-2 : normal return. the spline returned is the least-squares
190c            constrained polynomial. in this extreme case fp gives the
191c            upper bound for the smoothing factor s.
192c   ier=1  : error. the required storage space exceeds the available
193c            storage space, as specified by the parameters nuest and
194c            nvest.
195c            probably causes : nuest or nvest too small. if these param-
196c            eters are already large, it may also indicate that s is
197c            too small
198c            the approximation returned is the least-squares spline
199c            according to the current set of knots. the parameter fp
200c            gives the corresponding sum of squared residuals (fp>s).
201c   ier=2  : error. a theoretically impossible result was found during
202c            the iteration process for finding a smoothing spline with
203c            fp = s. probably causes : s too small.
204c            there is an approximation returned but the corresponding
205c            sum of squared residuals does not satisfy the condition
206c            abs(fp-s)/s < tol.
207c   ier=3  : error. the maximal number of iterations maxit (set to 20
208c            by the program) allowed for finding a smoothing spline
209c            with fp=s has been reached. probably causes : s too small
210c            there is an approximation returned but the corresponding
211c            sum of squared residuals does not satisfy the condition
212c            abs(fp-s)/s < tol.
213c   ier=10 : error. on entry, the input data are controlled on validity
214c            the following restrictions must be satisfied.
215c            -1<=iopt(1)<=1, 0<=iopt(2)<=1, 0<=iopt(3)<=1,
216c            -1<=ider(1)<=1, 0<=ider(2)<=1, ider(2)=0 if iopt(2)=0.
217c            -1<=ider(3)<=1, 0<=ider(4)<=1, ider(4)=0 if iopt(3)=0.
218c            mu >= mumin (see above), mv >= 4, nuest >=8, nvest >= 8,
219c            kwrk>=5+mu+mv+nuest+nvest,
220c            lwrk >= 12+nuest*(mv+nvest+3)+nvest*24+4*mu+8*mv+
221c             max(nuest,mv+nvest)
222c            0< u(i-1)<u(i)< pi,i=2,..,mu,
223c            -pi<=v(1)< pi, v(1)<v(i-1)<v(i)<v(1)+2*pi, i=3,...,mv
224c            if iopt(1)=-1: 8<=nu<=min(nuest,mu+6+iopt(2)+iopt(3))
225c                           0<tu(5)<tu(6)<...<tu(nu-4)< pi
226c                           8<=nv<=min(nvest,mv+7)
227c                           v(1)<tv(5)<tv(6)<...<tv(nv-4)<v(1)+2*pi
228c                    the schoenberg-whitney conditions, i.e. there must
229c                    be subset of grid co-ordinates uu(p) and vv(q) such
230c                    that   tu(p) < uu(p) < tu(p+4) ,p=1,...,nu-4
231c                     (iopt(2)=1 and iopt(3)=1 also count for a uu-value
232c                           tv(q) < vv(q) < tv(q+4) ,q=1,...,nv-4
233c                     (vv(q) is either a value v(j) or v(j)+2*pi)
234c            if iopt(1)>=0: s>=0
235c                       if s=0: nuest>=mu+6+iopt(2)+iopt(3), nvest>=mv+7
236c            if one of these conditions is found to be violated,control
237c            is immediately repassed to the calling program. in that
238c            case there is no approximation returned.
239c
240c further comments:
241c   spgrid does not allow individual weighting of the data-values.
242c   so, if these were determined to widely different accuracies, then
243c   perhaps the general data set routine sphere should rather be used
244c   in spite of efficiency.
245c   by means of the parameter s, the user can control the tradeoff
246c   between closeness of fit and smoothness of fit of the approximation.
247c   if s is too large, the spline will be too smooth and signal will be
248c   lost ; if s is too small the spline will pick up too much noise. in
249c   the extreme cases the program will return an interpolating spline if
250c   s=0 and the constrained least-squares polynomial(degrees 3,0)if s is
251c   very large. between these extremes, a properly chosen s will result
252c   in a good compromise between closeness of fit and smoothness of fit.
253c   to decide whether an approximation, corresponding to a certain s is
254c   satisfactory the user is highly recommended to inspect the fits
255c   graphically.
256c   recommended values for s depend on the accuracy of the data values.
257c   if the user has an idea of the statistical errors on the data, he
258c   can also find a proper estimate for s. for, by assuming that, if he
259c   specifies the right s, spgrid will return a spline s(u,v) which
260c   exactly reproduces the function underlying the data he can evaluate
261c   the sum((r(i,j)-s(u(i),v(j)))**2) to find a good estimate for this s
262c   for example, if he knows that the statistical errors on his r(i,j)-
263c   values is not greater than 0.1, he may expect that a good s should
264c   have a value not larger than mu*mv*(0.1)**2.
265c   if nothing is known about the statistical error in r(i,j), s must
266c   be determined by trial and error, taking account of the comments
267c   above. the best is then to start with a very large value of s (to
268c   determine the least-squares polynomial and the corresponding upper
269c   bound fp0 for s) and then to progressively decrease the value of s
270c   ( say by a factor 10 in the beginning, i.e. s=fp0/10,fp0/100,...
271c   and more carefully as the approximation shows more detail) to
272c   obtain closer fits.
273c   to economize the search for a good s-value the program provides with
274c   different modes of computation. at the first call of the routine, or
275c   whenever he wants to restart with the initial set of knots the user
276c   must set iopt(1)=0.
277c   if iopt(1) = 1 the program will continue with the knots found at
278c   the last call of the routine. this will save a lot of computation
279c   time if spgrid is called repeatedly for different values of s.
280c   the number of knots of the spline returned and their location will
281c   depend on the value of s and on the complexity of the shape of the
282c   function underlying the data. if the computation mode iopt(1) = 1
283c   is used, the knots returned may also depend on the s-values at
284c   previous calls (if these were smaller). therefore, if after a number
285c   of trials with different s-values and iopt(1)=1,the user can finally
286c   accept a fit as satisfactory, it may be worthwhile for him to call
287c   spgrid once more with the chosen value for s but now with iopt(1)=0.
288c   indeed, spgrid may then return an approximation of the same quality
289c   of fit but with fewer knots and therefore better if data reduction
290c   is also an important objective for the user.
291c   the number of knots may also depend on the upper bounds nuest and
292c   nvest. indeed, if at a certain stage in spgrid the number of knots
293c   in one direction (say nu) has reached the value of its upper bound
294c   (nuest), then from that moment on all subsequent knots are added
295c   in the other (v) direction. this may indicate that the value of
296c   nuest is too small. on the other hand, it gives the user the option
297c   of limiting the number of knots the routine locates in any direction
298c   for example, by setting nuest=8 (the lowest allowable value for
299c   nuest), the user can indicate that he wants an approximation which
300c   is a simple cubic polynomial in the variable u.
301c
302c  other subroutines required:
303c    fpspgr,fpchec,fpchep,fpknot,fpopsp,fprati,fpgrsp,fpsysy,fpback,
304c    fpbacp,fpbspl,fpcyt1,fpcyt2,fpdisc,fpgivs,fprota
305c
306c  references:
307c   dierckx p. : fast algorithms for smoothing data over a disc or a
308c                sphere using tensor product splines, in "algorithms
309c                for approximation", ed. j.c.mason and m.g.cox,
310c                clarendon press oxford, 1987, pp. 51-65
311c   dierckx p. : fast algorithms for smoothing data over a disc or a
312c                sphere using tensor product splines, report tw73, dept.
313c                computer science,k.u.leuven, 1985.
314c   dierckx p. : curve and surface fitting with splines, monographs on
315c                numerical analysis, oxford university press, 1993.
316c
317c  author:
318c    p.dierckx
319c    dept. computer science, k.u. leuven
320c    celestijnenlaan 200a, b-3001 heverlee, belgium.
321c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
322c
323c  creation date : july 1985
324c  latest update : march 1989
325c
326c  ..
327c  ..scalar arguments..
328      real*8 r0,r1,s,fp
329      integer mu,mv,nuest,nvest,nu,nv,lwrk,kwrk,ier
330c  ..array arguments..
331      integer iopt(3),ider(4),iwrk(kwrk)
332      real*8 u(mu),v(mv),r(mu*mv),c((nuest-4)*(nvest-4)),tu(nuest),
333     * tv(nvest),wrk(lwrk)
334c  ..local scalars..
335      real*8 per,pi,tol,uu,ve,rmax,rmin,one,half,rn,rb,re
336      integer i,i1,i2,j,jwrk,j1,j2,kndu,kndv,knru,knrv,kwest,l,
337     * ldr,lfpu,lfpv,lwest,lww,m,maxit,mumin,muu,nc
338c  ..function references..
339      real*8 datan2
340      integer max0
341c  ..subroutine references..
342c    fpchec,fpchep,fpspgr
343c  ..
344c  set constants
345      one = 1d0
346      half = 0.5e0
347      pi = datan2(0d0,-one)
348      per = pi+pi
349      ve = v(1)+per
350c  we set up the parameters tol and maxit.
351      maxit = 20
352      tol = 0.1e-02
353c  before starting computations, a data check is made. if the input data
354c  are invalid, control is immediately repassed to the calling program.
355      ier = 10
356      if(iopt(1).lt.(-1) .or. iopt(1).gt.1) go to 200
357      if(iopt(2).lt.0 .or. iopt(2).gt.1) go to 200
358      if(iopt(3).lt.0 .or. iopt(3).gt.1) go to 200
359      if(ider(1).lt.(-1) .or. ider(1).gt.1) go to 200
360      if(ider(2).lt.0 .or. ider(2).gt.1) go to 200
361      if(ider(2).eq.1 .and. iopt(2).eq.0) go to 200
362      if(ider(3).lt.(-1) .or. ider(3).gt.1) go to 200
363      if(ider(4).lt.0 .or. ider(4).gt.1) go to 200
364      if(ider(4).eq.1 .and. iopt(3).eq.0) go to 200
365      mumin = 4
366      if(ider(1).ge.0) mumin = mumin-1
367      if(iopt(2).eq.1 .and. ider(2).eq.1) mumin = mumin-1
368      if(ider(3).ge.0) mumin = mumin-1
369      if(iopt(3).eq.1 .and. ider(4).eq.1) mumin = mumin-1
370      if(mumin.eq.0) mumin = 1
371      if(mu.lt.mumin .or. mv.lt.4) go to 200
372      if(nuest.lt.8 .or. nvest.lt.8) go to 200
373      m = mu*mv
374      nc = (nuest-4)*(nvest-4)
375      lwest = 12+nuest*(mv+nvest+3)+24*nvest+4*mu+8*mv+
376     * max0(nuest,mv+nvest)
377      kwest = 5+mu+mv+nuest+nvest
378      if(lwrk.lt.lwest .or. kwrk.lt.kwest) go to 200
379      if(u(1).le.0. .or. u(mu).ge.pi) go to 200
380      if(mu.eq.1) go to 30
381      do 20 i=2,mu
382        if(u(i-1).ge.u(i)) go to 200
383  20  continue
384  30  if(v(1).lt. (-pi) .or. v(1).ge.pi ) go to 200
385      if(v(mv).ge.v(1)+per) go to 200
386      do 40 i=2,mv
387        if(v(i-1).ge.v(i)) go to 200
388  40  continue
389      if(iopt(1).gt.0) go to 140
390c  if not given, we compute an estimate for r0.
391      rn = mv
392      if(ider(1).lt.0) go to 45
393      rb = r0
394      go to 55
395  45  rb = 0.
396      do 50 i=1,mv
397         rb = rb+r(i)
398  50  continue
399      rb = rb/rn
400c  if not given, we compute an estimate for r1.
401  55  if(ider(3).lt.0) go to 60
402      re = r1
403      go to 70
404  60  re = 0.
405      j = m
406      do 65 i=1,mv
407         re = re+r(j)
408         j = j-1
409  65  continue
410      re = re/rn
411c  we determine the range of r-values.
412  70  rmin = rb
413      rmax = re
414      do 80 i=1,m
415         if(r(i).lt.rmin) rmin = r(i)
416         if(r(i).gt.rmax) rmax = r(i)
417  80  continue
418      wrk(5) = rb
419      wrk(6) = 0.
420      wrk(7) = 0.
421      wrk(8) = re
422      wrk(9) = 0.
423      wrk(10) = 0.
424      wrk(11) = rmax -rmin
425      wrk(12) = wrk(11)
426      iwrk(4) = mu
427      iwrk(5) = mu
428      if(iopt(1).eq.0) go to 140
429      if(nu.lt.8 .or. nu.gt.nuest) go to 200
430      if(nv.lt.11 .or. nv.gt.nvest) go to 200
431      j = nu
432      do 90 i=1,4
433        tu(i) = 0.
434        tu(j) = pi
435        j = j-1
436  90  continue
437      l = 13
438      wrk(l) = 0.
439      if(iopt(2).eq.0) go to 100
440      l = l+1
441      uu = u(1)
442      if(uu.gt.tu(5)) uu = tu(5)
443      wrk(l) = uu*half
444 100  do 110 i=1,mu
445        l = l+1
446        wrk(l) = u(i)
447 110  continue
448      if(iopt(3).eq.0) go to 120
449      l = l+1
450      uu = u(mu)
451      if(uu.lt.tu(nu-4)) uu = tu(nu-4)
452      wrk(l) = uu+(pi-uu)*half
453 120  l = l+1
454      wrk(l) = pi
455      muu = l-12
456      call fpchec(wrk(13),muu,tu,nu,3,ier)
457      if(ier.ne.0) go to 200
458      j1 = 4
459      tv(j1) = v(1)
460      i1 = nv-3
461      tv(i1) = ve
462      j2 = j1
463      i2 = i1
464      do 130 i=1,3
465        i1 = i1+1
466        i2 = i2-1
467        j1 = j1+1
468        j2 = j2-1
469        tv(j2) = tv(i2)-per
470        tv(i1) = tv(j1)+per
471 130  continue
472      l = 13
473      do 135 i=1,mv
474        wrk(l) = v(i)
475        l = l+1
476 135  continue
477      wrk(l) = ve
478      call fpchep(wrk(13),mv+1,tv,nv,3,ier)
479      if (ier.eq.0) go to 150
480      go to 200
481 140  if(s.lt.0.) go to 200
482      if(s.eq.0. .and. (nuest.lt.(mu+6+iopt(2)+iopt(3)) .or.
483     * nvest.lt.(mv+7)) ) go to 200
484c  we partition the working space and determine the spline approximation
485 150  ldr = 5
486      lfpu = 13
487      lfpv = lfpu+nuest
488      lww = lfpv+nvest
489      jwrk = lwrk-12-nuest-nvest
490      knru = 6
491      knrv = knru+mu
492      kndu = knrv+mv
493      kndv = kndu+nuest
494      call fpspgr(iopt,ider,u,mu,v,mv,r,m,rb,re,s,nuest,nvest,tol,maxit,
495     *
496     * nc,nu,tu,nv,tv,c,fp,wrk(1),wrk(2),wrk(3),wrk(4),wrk(lfpu),
497     * wrk(lfpv),wrk(ldr),wrk(11),iwrk(1),iwrk(2),iwrk(3),iwrk(4),
498     * iwrk(5),iwrk(knru),iwrk(knrv),iwrk(kndu),iwrk(kndv),wrk(lww),
499     * jwrk,ier)
500 200  return
501      end
502