V3FIT
r8bcspline.f
1 c bcspline -- dmc 30 May 1996
2 c
3 c set up coefficients for bicubic spline with following BC's:
4 c FULL BC CONTROL at all bdys
5 c
6 c inhomogeneous explicit BCs -- this means setting of 1st or 2nd
7 c derivative at boundary to a non-zero value.
8 c
9 c periodic, not-a-knot, zero derivative, and divided-difference based
10 c BCs are "homogeneous"-- i.e. if splines s & t satisfy the BC then
11 c the spline (c*s + t) formed as a linear combination of these two
12 c splines, also satisfies the BC.
13 c
14 c algorithm note -- handling of inhomogeneous explicit BC's while
15 c maintaining full C2 differentiability is delicate. Basic method: use
16 c a fully C2 method based on the "not-a-knot" BC, and then, correct to
17 c meet each user BC by calculating a C2 spline that is zero at all grid
18 c points but satisfies a BC which is the difference btw the user spec
19 c and the not-a-knot result; add the coeffs of this into the original.
20 c
21 c for this more workspace is needed: nwk .ge. 4*inx*inth +5*max(inx,inth)
22 c
23  subroutine r8bcspline(x,inx,th,inth,fspl,inf3,
24  > ibcxmin,bcxmin,ibcxmax,bcxmax,
25  > ibcthmin,bcthmin,ibcthmax,bcthmax,
26  > wk,nwk,ilinx,ilinth,ier)
27 c
28 !============
29 ! idecl: explicitize implicit INTEGER declarations:
30  IMPLICIT NONE
31  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
32  INTEGER inth,inf3,ibcxmin,ibcxmax,ibcthmin,ibcthmax,nwk,ilinx
33  INTEGER ilinth,ier,inx,iflg2,ix,itest,ierx,ierth,inxo,ith
34  INTEGER intho,ic,ibcthmina,ibcthmaxa,iasc,iinc,iawk,jx,jth,ii
35  INTEGER iadr,ia5w,iaspl
36 !============
37 ! idecl: explicitize implicit REAL declarations:
38  real*8 xo2,xo6,zhxn,zhth,zdiff1,zdiff2
39 !============
40  real*8 x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
41  real*8 bcxmin(*),bcxmax(*) ! (inth) if used
42  real*8 bcthmin(*),bcthmax(*) ! (inx) if used
43 c
44 c input:
45 c x(1...inx) -- abscissae, first dimension of data
46 c th(1...inth) -- abscissae, second dimension of data f(x,th)
47 c fspl(1,1,1..inx,1..inth) -- function values
48 c inf3 -- fspl dimensioning, inf3.ge.inx required.
49 c
50 c boundary conditions input:
51 c ibcxmin -- indicator for boundary condition at x(1):
52 c bcxmin(...) -- boundary condition data
53 c =-1 -- periodic boundary condition
54 c =0 -- use "not a knot", bcxmin(...) ignored
55 c =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
56 c =2 -- match 2nd derivative, specified at x(1),th(ith) by bcxmin(ith)
57 c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
58 c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
59 c =5 -- match 1st derivative df/dx to 1st divided difference
60 c =6 -- match 2nd derivative d2f/dx2 to 2nd divided difference
61 c =7 -- match 3rd derivative d3f/dx3 3rd divided difference
62 c (for more detailed definition of BCs 5-7, see the
63 c comments of subroutine mkspline)
64 c NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
65 c
66 c ibcxmax -- indicator for boundary condition at x(nx):
67 c bcxmax(...) -- boundary condition data
68 c (interpretation as with ibcxmin, bcxmin)
69 c NOTE: if ibcxmin=-1, ibcxmax is ignored! ...and the BC is periodic.
70 c
71 c ibcthmin -- indicator for boundary condition at th(1):
72 c bcthmin(...) -- boundary condition data
73 c (interpretation as with ibcxmin, bcxmin)
74 c ibcthmax -- indicator for boundary condition at th(inth):
75 c bcthmax(...) -- boundary condition data
76 c (interpretation as with ibcxmin, bcxmin)
77 c NOTE: if ibcthmin=-1, ibcthmax is ignored! ...and the BC is periodic.
78 c
79 c NOTE the bcxmin,bcxmax,bcthmin,bcthmax arrays are only used if the
80 c corresponding boundary condition flags are set to 1 or 2.
81 c Carefully note the dimensioning of these arrays!
82 c
83 c output:
84 c fspl(*,*,1..inx,1..inth) -- bicubic spline coeffs (4x4)
85 c ...fspl(1,1,*,*) is not replaced.
86 c
87 c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced (tol=1e-3)
88 c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
89 c
90 c ier -- completion code, 0 for normal
91 c
92 c workspace:
93 c wk -- must be at least 5*max(inx,inth) large
94 c 5*max(inx,inth) + 4*inx*inth large
95 c if explicit non-zero d/dth or d2/dth2 BC's
96 c are supplied.
97 c nwk -- size of workspace of workspace provided
98 c
99 c---------------------------------
100 c in what follows, "f" is an abbreviation for "fspl"...
101 c
102 c compute bicubic spline of 2d function, given values at the grid
103 c grid crossing points, f(1,1,i,j)=f(x(i),th(j)).
104 c
105 c on evaluation: for point x btw x(i) and x(i+1), dx=x-x(i)
106 c and th btw th(j) and th(j+1), dt=th-th(j),
107 c
108 c spline =
109 c f(1,1,i,j) + dx*f(2,1,i,j) + dx**2*f(3,1,i,j) + dx**3*f(4,1,i,j)
110 c +dt*(f(1,2,i,j) + dx*f(2,2,i,j) + dx**2*f(3,2,i,j) + dx**3*f(4,2,i,j))
111 c +d2*(f(1,3,i,j) + dx*f(2,3,i,j) + dx**2*f(3,3,i,j) + dx**3*f(4,3,i,j))
112 c +d3*(f(1,4,i,j) + dx*f(2,4,i,j) + dx**2*f(3,4,i,j) + dx**3*f(4,4,i,j))
113 c
114 c where d2=dt**2 and d3=dt**3.
115 c
116  integer iselect1(10)
117  integer iselect2(10)
118 c
119  real*8 zcur(1)
120 c---------------------------------
121 c
122 c see if 2nd pass is needed due to "non-linear" d/dth bdy cond.
123 c
124  iflg2=0
125  if(ibcthmin.ne.-1) then
126  if((ibcthmin.eq.1).or.(ibcthmin.eq.2)) then
127  do ix=1,inx
128  if (bcthmin(ix).ne.0.0_r8) iflg2=1
129  enddo
130  endif
131  if((ibcthmax.eq.1).or.(ibcthmax.eq.2)) then
132  do ix=1,inx
133  if (bcthmax(ix).ne.0.0_r8) iflg2=1
134  enddo
135  endif
136  endif
137 c
138  ier=0
139  itest=5*max(inx,inth)
140  if(iflg2.eq.1) then
141  itest=itest +4*inx*inth
142  endif
143 c
144  if(nwk.lt.itest) then
145  write(6,9901) nwk,itest
146  9901 format(' ?bcspline: workspace too small:'/
147  > ' user supplied: nwk=',i6,'; need at least: ',i6/
148  > ' nwk=4*nx*ny +5*max(nx,ny) will work for any user'/
149  > ' choice of bdy conditions.')
150  ier=1
151  endif
152  if(inx.lt.2) then
153  write(6,'('' ?bcspline: at least 2 x points required.'')')
154  ier=1
155  endif
156  if(inth.lt.2) then
157  write(6,'('' ?bcspline: need at least 2 theta points.'')')
158  ier=1
159  endif
160 c
161  call ibc_ck(ibcxmin,'bcspline','xmin',-1,7,ier)
162  if(ibcxmin.ge.0) call ibc_ck(ibcxmax,'bcspline','xmax',0,7,ier)
163  call ibc_ck(ibcthmin,'bcspline','thmin',-1,7,ier)
164  if(ibcthmin.ge.0) call ibc_ck(ibcthmax,'bcspline','thmax',0,7,ier)
165 c
166 c check ilinx & x vector
167 c
168  call r8splinck(x,inx,ilinx,1.0e-3_r8,ierx)
169  if(ierx.ne.0) ier=2
170 c
171  if(ier.eq.2) then
172  write(6,'('' ?bcspline: x axis not strict ascending'')')
173  endif
174 c
175 c check ilinth & th vector
176 c
177  call r8splinck(th,inth,ilinth,1.0e-3_r8,ierth)
178  if(ierth.ne.0) ier=3
179 c
180  if(ier.eq.3) then
181  write(6,'('' ?bcspline: th axis not strict ascending'')')
182  endif
183 c
184  if(ier.ne.0) return
185 c
186 c------------------------------------
187 c
188  xo2=0.5_r8
189  xo6=1.0_r8/6.0_r8
190 c
191 c spline in x direction
192 c
193  inxo=4*(inx-1)
194  do ith=1,inth
195 c
196 c copy the function in
197 c
198  do ix=1,inx
199  wk(4*(ix-1)+1)=fspl(1,1,ix,ith)
200  enddo
201 c
202  if(ibcxmin.eq.1) then
203  wk(2)=bcxmin(ith)
204  else if(ibcxmin.eq.2) then
205  wk(3)=bcxmin(ith)
206  endif
207 c
208  if(ibcxmax.eq.1) then
209  wk(inxo+2)=bcxmax(ith)
210  else if(ibcxmax.eq.2) then
211  wk(inxo+3)=bcxmax(ith)
212  endif
213 c
214 c use Wayne's routine
215 c
216  call r8v_spline(ibcxmin,ibcxmax,inx,x,wk,wk(4*inx+1))
217 c
218 c copy the coefficients out
219 c
220  do ix=1,inx
221  fspl(2,1,ix,ith)=wk(4*(ix-1)+2)
222  fspl(3,1,ix,ith)=wk(4*(ix-1)+3)*xo2
223  fspl(4,1,ix,ith)=wk(4*(ix-1)+4)*xo6
224  enddo
225 c
226  enddo
227 c
228 c-----------------------------------
229 c
230 c spline in theta direction
231 c
232  intho=4*(inth-1)
233  do ix=1,inx
234 c
235 c spline each x coeff
236 c
237  do ic=1,4
238 c
239 c copy ordinates in
240 c
241  do ith=1,inth
242  wk(4*(ith-1)+1)=fspl(ic,1,ix,ith)
243  enddo
244 c
245 c first pass: use a linear BC -- if flag indicates BC correction
246 c will be needed, it will be done later
247 c
248  wk(2)=0.0_r8
249  wk(3)=0.0_r8
250  wk(intho+2)=0.0_r8
251  wk(intho+3)=0.0_r8
252 c
253  ibcthmina=ibcthmin
254  ibcthmaxa=ibcthmax
255  if(iflg2.eq.1) then
256  if((ibcthmin.eq.1).or.(ibcthmin.eq.2)) ibcthmina=0
257  if((ibcthmax.eq.1).or.(ibcthmax.eq.2)) ibcthmaxa=0
258  endif
259 c
260  call r8v_spline(ibcthmina,ibcthmaxa,inth,th,wk,wk(4*inth+1))
261 c
262 c copy coeffs out
263 c
264  do ith=1,inth
265  fspl(ic,2,ix,ith)=wk(4*(ith-1)+2)
266  fspl(ic,3,ix,ith)=wk(4*(ith-1)+3)*xo2
267  fspl(ic,4,ix,ith)=wk(4*(ith-1)+4)*xo6
268  enddo
269 c
270  enddo
271 c
272  enddo
273 c
274 c now make correction for user BC's if needed
275 c
276  if(iflg2.eq.1) then
277 c
278  iasc=1 ! wk addr for correction splines
279  iinc=4*inth ! spacing btw correction splines
280  iawk=iasc+4*inth*inx
281 c
282 c last grid zone widths
283 c
284  zhxn=x(inx)-x(inx-1)
285  jx=inx-1
286  zhth=th(inth)-th(inth-1)
287  jth=inth-1
288 c
289  do ii=1,10
290  iselect1(ii)=0
291  iselect2(ii)=0
292  enddo
293  if(ibcthmin.eq.1) iselect1(3)=1
294  if(ibcthmin.eq.2) iselect1(5)=1
295  if(ibcthmax.eq.1) iselect2(3)=1
296  if(ibcthmax.eq.2) iselect2(5)=1
297 c
298 c loop over BC's
299 c
300  do ix=1,inx
301 c
302 c (a) d/dth @ th(1) difference btw current BC and user request
303 c
304  if(ibcthmin.eq.1) then
305  if(ix.lt.inx) then
306  zcur(1)=fspl(1,2,ix,1) ! 1st deriv.
307  else
308  zcur(1)=fspl(1,2,jx,1)+zhxn*(fspl(2,2,jx,1)+zhxn*
309  > (fspl(3,2,jx,1)+zhxn*fspl(4,2,jx,1)))
310  endif
311  zdiff1=bcthmin(ix)-zcur(1)
312  else if(ibcthmin.eq.2) then
313  if(ix.lt.inx) then
314  zcur(1)=2.0_r8*fspl(1,3,ix,1) ! 2nd deriv.
315  else
316  zcur(1)=2.0_r8*
317  > (fspl(1,3,jx,1)+zhxn*(fspl(2,3,jx,1)+zhxn*
318  > (fspl(3,3,jx,1)+zhxn*fspl(4,3,jx,1))))
319  endif
320  zdiff1=bcthmin(ix)-zcur(1)
321  else
322  zdiff1=0.0_r8
323  endif
324 c
325 c (b) d/dth @ th(inth) difference btw current BC and user request
326 c
327  if(ibcthmax.eq.1) then
328  if(ix.lt.inx) then
329 c 1st deriv.
330  zcur(1)=
331  > fspl(1,2,ix,jth)+zhth*(2.0_r8*fspl(1,3,ix,jth)+
332  > zhth*3.0_r8*fspl(1,4,ix,jth))
333  else
334  call r8bcspeval(x(inx),th(inth),iselect2, zcur,
335  > x,inx,th,inth,ilinx,ilinth,fspl,inf3,ier)
336  if(ier.ne.0) return
337  endif
338  zdiff2=bcthmax(ix)-zcur(1)
339  else if(ibcthmax.eq.2) then
340  if(ix.lt.inx) then
341 c 2nd deriv.
342  zcur(1)=2.0_r8*fspl(1,3,ix,jth)+
343  > 6.0_r8*zhth*fspl(1,4,ix,jth)
344  else
345  call r8bcspeval(x(inx),th(inth),iselect2, zcur(1),
346  > x,inx,th,inth,ilinx,ilinth,fspl,inf3,ier)
347  if(ier.ne.0) return
348  endif
349  zdiff2=bcthmax(ix)-zcur(1)
350  else
351  zdiff2=0.0_r8
352  endif
353 c
354 c ok compute the theta spline with BC's to span the difference(s)
355 c these theta "correction splines" are zero at all the grid points
356 c but have at least one non-zero 1st or 2nd derivative BC
357 c
358  iadr=iasc+(ix-1)*iinc
359  do ith=1,inth
360  wk(iadr+4*(ith-1))=0.0_r8
361  enddo
362 c
363  wk(iadr+1)=0.0_r8
364  wk(iadr+2)=0.0_r8
365  wk(iadr+intho+1)=0.0_r8
366  wk(iadr+intho+2)=0.0_r8
367 c
368  if(ibcthmin.eq.1) then
369  wk(iadr+1)=zdiff1
370  else if(ibcthmin.eq.2) then
371  wk(iadr+2)=zdiff1
372  endif
373 c
374  if(ibcthmax.eq.1) then
375  wk(iadr+intho+1)=zdiff2
376  else if(ibcthmax.eq.2) then
377  wk(iadr+intho+2)=zdiff2
378  endif
379 c
380  call r8v_spline(ibcthmin,ibcthmax,inth,th,wk(iadr),wk(iawk))
381  enddo
382 c
383 c add in results to main array -- th spline coef corrections
384 c
385  do ix=1,inx
386  iadr=iasc+(ix-1)*iinc
387  do ith=1,inth-1
388  wk(iadr+4*(ith-1)+2)=wk(iadr+4*(ith-1)+2)*xo2
389  wk(iadr+4*(ith-1)+3)=wk(iadr+4*(ith-1)+3)*xo6
390  if(ix.lt.inx) then
391  fspl(1,2,ix,ith)=fspl(1,2,ix,ith)+wk(iadr+4*(ith-1)+1)
392  fspl(1,3,ix,ith)=fspl(1,3,ix,ith)+wk(iadr+4*(ith-1)+2)
393  fspl(1,4,ix,ith)=fspl(1,4,ix,ith)+wk(iadr+4*(ith-1)+3)
394  endif
395  enddo
396  enddo
397 c
398 c compute the x splines of the th spline correction coeffs
399 c
400  ia5w=iawk+4*inx
401 c
402  do ith=1,inth-1
403  do ic=2,4
404  do ix=1,inx
405  iaspl=iasc+iinc*(ix-1)
406  wk(iawk+4*(ix-1))=wk(iaspl+4*(ith-1)+(ic-1))
407  enddo
408 c
409 c use zero BCs for this correction spline
410 c
411  wk(iawk+1)=0.0_r8
412  wk(iawk+2)=0.0_r8
413  wk(iawk+inxo+1)=0.0_r8
414  wk(iawk+inxo+2)=0.0_r8
415 c
416 c periodic spline of correction spline higher coeffs (1st coeffs are
417 c all zero by defn of the correction spline
418 c
419  call r8v_spline(ibcxmin,ibcxmax,inx,x,wk(iawk),wk(ia5w))
420 c
421  do ix=1,inx-1
422  fspl(2,ic,ix,ith)=fspl(2,ic,ix,ith)+
423  > wk(iawk+4*(ix-1)+1)
424  fspl(3,ic,ix,ith)=fspl(3,ic,ix,ith)+
425  > wk(iawk+4*(ix-1)+2)*xo2
426  fspl(4,ic,ix,ith)=fspl(4,ic,ix,ith)+
427  > wk(iawk+4*(ix-1)+3)*xo6
428  enddo
429 c
430  enddo
431  enddo ! ith
432 c
433  endif ! BC correction needs test
434 c
435  return
436  end