V3FIT
tcspline.f
1 c tcspline -- dmc 20 Jan 1999
2 c
3 c set up coefficients for bicubic spline with following BC's:
4 c * LHS and RHS handled as in cubspl.for for 1st coordinate
5 c * derivatives periodic in second coordinate (use pspline.for)
6 c
7 c workspace:
8 c if phi bdy cond. is periodic, not-a-knot, df/dphi = 0 everywhere,
9 c or d2f/dphi2 = 0 everywhere, then the phi boundary condition is
10 c "linear" and a workspace of size at least:
11 c
12 c nwk = 20*inx*inth + 10*max(inx,inth,inph)
13 c
14 c will suffice.
15 c
16 c if the phi bdy cond. involves specification of df/dphi .ne. 0 or
17 c d2f/dphi .ne. 0 at any (x,theta) grid point, then, the phi boundary
18 c condition is "non-linear", a correction step is needed, and a workspace
19 c of size at least:
20 c
21 c nwk = 16*inx*inth*inph
22 c
23 c is required.
24 c
25  subroutine tcspline(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
26  > ibcxmin,bcxmin,ibcxmax,bcxmax,inb1x,
27  > ibcthmin,bcthmin,ibcthmax,bcthmax,inb1th,
28  > ibcphmin,bcphmin,ibcphmax,bcphmax,inb1ph,
29  > wk,nwk,ilinx,ilinth,ilinph,ier)
30 c
31  real x(inx),th(inth),ph(inph)
32  real fspl(4,4,4,inf4,inf5,inph),wk(nwk)
33  real bcxmin(inb1x,*),bcxmax(inb1x,*) ! inth x inph defined (if used)
34  real bcthmin(inb1th,*),bcthmax(inb1th,*) ! inx x inph defined (if used)
35  real bcphmin(inb1ph,*),bcphmax(inb1ph,*) ! inx x inth defined (if used)
36 c
37 c input:
38 c x(1...inx) -- abscissae, first dimension of data
39 c th(1...inth) -- abscissae, second (periodic) dimension of data
40 c ph(1...inph) -- abscissae, third (periodic) dimension of data
41 c fspl(1,1,1,1..inx,1..inth,1..inph) -- function values
42 c inf4 -- fspl dimensioning, inf4.ge.inx required.
43 c inf5 -- fspl dimensioning, inf5.ge.inth required.
44 c
45 c boundary conditions input:
46 c
47 c bc data at xmin, xmax vs. theta,phi
48 c bc data at thmin, thmax vs. x,phi
49 c bc data at phmin, phmax vs. x,theta
50 c
51 c ibcxmin -- indicator for boundary condition at x(1):
52 c bcxmin(...) -- boundary condition data
53 c =-1 -- use periodic boundary condition
54 c =0 -- use "not a knot", bcxmin(...) ignored
55 c =1 -- match slope, specified at x(1),th(ith),ph(iph) by bcxmin(ith,iph)
56 c =2 -- match 2nd derivative, specified at x(1),th(ith),ph(iph)
57 c by bcxmin(ith,iph
58 c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
59 c =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
60 c =5 -- match 1st derivative to 1st divided difference
61 c =6 -- match 2nd derivative to 2nd divided difference
62 c =7 -- match 3rd derivative to 3rd divided difference
63 c (for more detailed definition of BCs 5-7, see the
64 c comments of subroutine mkspline)
65 c NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
66 c
67 c ibcxmax -- indicator for boundary condition at x(nx):
68 c bcxmax(...) -- boundary condition data
69 c (interpretation as with ibcxmin, bcxmin)
70 c NOTE: if ibcxmin=-1 then the periodic BC applies on both sides
71 c and ibcxmax is ignored.
72 c inb1x -- 1st dimension of bcxmin, bcxmax: if ibcxmin or ibcxmax .gt. 0
73 c this must be .ge. inth:
74 c
75 c interpretation of ibcthmin,bcthmin,ibcthmax,bcthmax,inb1th
76 c is same as with ibcxmin,...
77 c
78 c interpretation of ibcphmin,bcphmin,ibcphmax,bcphmax,inb1ph
79 c is same as with ibcxmin,...
80 c
81 c the explicit bdy condition arrays are referenced only if the
82 c corresponding "ibc" flag values are set to 1 or 2.
83 c
84 c output:
85 c fspl(*,*,*,1..inx,1..inth,1..inph) -- bicubic spline coeffs (4x4)
86 c ...fspl(1,1,1,*,*,*) is not replaced.
87 c
88 c ilinx -- =1 on output if x(inx) pts are nearly evenly spaced (tol=1e-3)
89 c ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
90 c ilinph-- =1 on output if ph(inph) evenly spaced (tol=1e-3)
91 c
92 c ier -- completion code, 0 for normal
93 c
94 c workspace:
95 c wk -- must be at least 5*max(inx,inth,inph) large -- or more, see
96 c comments, above.
97 c nwk -- size of workspace
98 c
99 c---------------------------------
100 c ** in what follows, f is an abbreviation for fspl **
101 c
102 c compute tricubic spline of 3d function, given values at the
103 c grid crossing points, f(1,1,1,i,j,k)=f(x(i),th(j),ph(k)).
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 and ph btw ph(k) and ph(k+1), dp=ph-ph(k),
108 c
109 c spline =
110 c f(1,1,1,i,j,k)+dx*f(2,1,1,i,j,k)+dx2*f(3,1,1,i,j,k)+dx3*f(4,1,1,i,j,k)
111 c +dt*(f(1,2,1,i,j,k)+dx*f(2,2,1,i,j,k)+dx2*f(3,2,1,i,j,k)+dx3*f(4,2,1,i,j,k))
112 c +dt2*(f(1,3,1,i,j,k)+dx*f(2,3,1,i,j,k)+dx2*f(3,3,1,i,j,k)+dx3*f(4,3,1,i,j,k))
113 c +dt3*(f(1,4,1,i,j,k)+dx*f(2,4,1,i,j,k)+dx2*f(3,4,1,i,j,k)+dx3*f(4,4,1,i,j,k))
114 c +dp*(
115 c f(1,1,2,i,j,k)+dx*f(2,1,2,i,j,k)+dx2*f(3,1,2,i,j,k)+dx3*f(4,1,2,i,j,k)
116 c +dt*(f(1,2,2,i,j,k)+dx*f(2,2,2,i,j,k)+dx2*f(3,2,2,i,j,k)+dx3*f(4,2,2,i,j,k))
117 c +dt2*(f(1,3,2,i,j,k)+dx*f(2,3,2,i,j,k)+dx2*f(3,3,2,i,j,k)+dx3*f(4,3,2,i,j,k))
118 c +dt3*(f(1,4,2,i,j,k)+dx*f(2,4,2,i,j,k)+dx2*f(3,4,2,i,j,k)+dx3*f(4,4,2,i,j,k)))
119 c +dp2*(
120 c f(1,1,3,i,j,k)+dx*f(2,1,3,i,j,k)+dx2*f(3,1,3,i,j,k)+dx3*f(4,1,3,i,j,k)
121 c +dt*(f(1,2,3,i,j,k)+dx*f(2,2,3,i,j,k)+dx2*f(3,2,3,i,j,k)+dx3*f(4,2,3,i,j,k))
122 c +dt2*(f(1,3,3,i,j,k)+dx*f(2,3,3,i,j,k)+dx2*f(3,3,3,i,j,k)+dx3*f(4,3,3,i,j,k))
123 c +dt3*(f(1,4,3,i,j,k)+dx*f(2,4,3,i,j,k)+dx2*f(3,4,3,i,j,k)+dx3*f(4,4,3,i,j,k)))
124 c +dp3*(
125 c f(1,1,4,i,j,k)+dx*f(2,1,4,i,j,k)+dx2*f(3,1,4,i,j,k)+dx3*f(4,1,4,i,j,k)
126 c +dt*(f(1,2,4,i,j,k)+dx*f(2,2,4,i,j,k)+dx2*f(3,2,4,i,j,k)+dx3*f(4,2,4,i,j,k))
127 c +dt2*(f(1,3,4,i,j,k)+dx*f(2,3,4,i,j,k)+dx2*f(3,3,4,i,j,k)+dx3*f(4,3,4,i,j,k))
128 c +dt3*(f(1,4,4,i,j,k)+dx*f(2,4,4,i,j,k)+dx2*f(3,4,4,i,j,k)+dx3*f(4,4,4,i,j,k)))
129 c
130 c where dx2=dx**2 and dx3=dx**3.
131 c where dt2=dt**2 and dt3=dt**3.
132 c where dp2=dp**2 and dp3=dp**3.
133 c
134 c---------------------------------
135  integer iselect1(10)
136  integer iselect2(10)
137 c
138  real z0,z1,ztol
139  real zcur(1)
140 c
141  data z0/0.0e0/
142  data z1/1.0e0/
143  data ztol/1.0e-3/
144 c
145 c---------------------------------
146 c
147  ier=0
148 c
149  iflg=0
150 c
151 c check phi bdy condition "linearity"
152 c
153  if(ibcphmin.ne.-1) then
154  if((ibcphmin.eq.1).or.(ibcphmin.eq.2)) then
155  do ith=1,inth
156  do ix=1,inx
157  if(bcphmin(ix,ith).ne.z0) iflg=1
158  enddo
159  enddo
160  endif
161  if((ibcphmax.eq.1).or.(ibcphmax.eq.2)) then
162  do ith=1,inth
163  do ix=1,inx
164  if(bcphmax(ix,ith).ne.z0) iflg=1
165  enddo
166  enddo
167  endif
168  endif
169 c
170  itest=10*max(inx,inth,inph) + 20*inx*inth
171  if(iflg.eq.1) then
172  itest=16*inx*inth*inph
173  endif
174 C
175  if(nwk.lt.itest) then
176  write(6,9901) nwk,itest
177  ier=1
178  9901 format(' ?tcspline: workspace too small.'/
179  > ' user supplied nwk=',i7,'; need at least: ',i7/
180  > ' If no explicit df/dph boundary condition is set,'/
181  > ' nwk = 20*inx*inth + 10*max(inx,inth,inph) can be used.'/
182  > ' If an explicit df/dph or d2f/dph2 boundary condition'/
183  > ' is set, nwk=16*inx*inth*inph is required.')
184  endif
185  if(inx.lt.2) then
186  write(6,'('' ?tcspline: at least 2 x points required.'')')
187  ier=1
188  endif
189  if(inth.lt.2) then
190  write(6,'('' ?tcspline: need at least 2 theta points.'')')
191  ier=1
192  endif
193  if(inph.lt.2) then
194  write(6,'('' ?tcspline: need at least 2 phi points.'')')
195  ier=1
196  endif
197 c
198  if((ibcxmin.eq.1).or.(ibcxmax.eq.1).or.(ibcxmin.eq.2).or.
199  > (ibcxmax.eq.2)) then
200  if(inb1x.lt.inth) then
201  ier=1
202  write(6,
203  >'('.lt.' ?tcspline: 1st dim of bcxmin/max arrays inth'')')
204  endif
205  endif
206 c
207  if((ibcthmin.eq.1).or.(ibcthmax.eq.1).or.(ibcthmin.eq.2).or.
208  > (ibcthmax.eq.2)) then
209  if(inb1th.lt.inx) then
210  ier=1
211  write(6,
212  >'('.lt.' ?tcspline: 1st dim of bcthmin/max arrays inx'')')
213  endif
214  endif
215 c
216  if((ibcphmin.eq.1).or.(ibcphmax.eq.1).or.(ibcphmin.eq.2).or.
217  > (ibcphmax.eq.2)) then
218  if(inb1ph.lt.inx) then
219  ier=1
220  write(6,
221  >'('.lt.' ?tcspline: 1st dim of bphmin/max arrays inx'')')
222  endif
223  endif
224 c
225  call ibc_ck(ibcxmin,'tcspline','xmin',-1,7,ier)
226  if(ibcxmin.ge.0) call ibc_ck(ibcxmax,'tcspline','xmax',0,7,ier)
227 c
228  call ibc_ck(ibcthmin,'tcspline','thmin',-1,7,ier)
229  if(ibcthmin.ge.0) call ibc_ck(ibcthmax,'tcspline','thmax',0,7,ier)
230 c
231  call ibc_ck(ibcphmin,'tcspline','phmin',-1,7,ier)
232  if(ibcphmax.ge.0) call ibc_ck(ibcphmax,'tcspline','phmax',0,7,ier)
233 c
234 c check ilinx & x vector
235 c
236  call splinck(x,inx,ilinx,ztol,ierx)
237  if(ierx.ne.0) ier=2
238 c
239  if(ier.eq.2) then
240  write(6,'('' ?tcspline: x axis not strict ascending'')')
241  endif
242 c
243 c check ilinth & th vector
244 c
245  call splinck(th,inth,ilinth,ztol,ierth)
246  if(ierth.ne.0) ier=3
247 c
248  if(ier.eq.3) then
249  write(6,'('' ?tcspline: theta axis not strict ascending'')')
250  endif
251 c
252 c check ilinth & th vector
253 c
254  call splinck(ph,inph,ilinph,ztol,ierph)
255  if(ierph.ne.0) ier=4
256 c
257  if(ier.eq.4) then
258  write(6,'('' ?tcspline: phi axis not strict ascending'')')
259  endif
260 c
261  if(ier.ne.0) return
262 c
263 c------------------------------------
264 c
265 c part 1. compute (x,theta) spline coeffs via an intermediate
266 c routine that call bcspline
267 c
268 c workspace addresses
269 c
270  iaspl2=1
271  iabcx1=iaspl2+16*inx*inth
272  iabcx2=iabcx1+inth
273  iabcth1=iabcx2+inth
274  iabcth2=iabcth1+inx
275  iawk=iabcth2+inx
276  inwk=nwk-iawk+1
277 c
278  do iph=1,inph
279 c
280 c copy bc data
281 c
282  do ix=1,inx
283  wk(iabcth1+ix-1)=0.0
284  wk(iabcth2+ix-1)=0.0
285  if((ibcthmin.eq.1).or.(ibcthmin.eq.2)) then
286  wk(iabcth1+ix-1)=bcthmin(ix,iph)
287  endif
288  if((ibcthmin.ne.-1).and.
289  > ((ibcthmax.eq.1).or.(ibcthmax.eq.2))) then
290  wk(iabcth2+ix-1)=bcthmax(ix,iph)
291  endif
292  enddo
293  do ith=1,inth
294  wk(iabcx1+ith-1)=0.0
295  wk(iabcx2+ith-1)=0.0
296  if((ibcxmin.eq.1).or.(ibcxmin.eq.2)) then
297  wk(iabcx1+ith-1)=bcxmin(ith,iph)
298  endif
299  if((ibcxmin.ne.-1).and.
300  > ((ibcxmax.eq.1).or.(ibcxmax.eq.2))) then
301  wk(iabcx2+ith-1)=bcxmax(ith,iph)
302  endif
303  enddo
304 c
305 c call 2d spline intermediary routine
306 c
307  call tcsp23(x,inx,th,inth,fspl(1,1,1,1,1,iph),inf4,
308  > ibcxmin,wk(iabcx1),ibcxmax,wk(iabcx2),
309  > ibcthmin,wk(iabcth1),ibcthmax,wk(iabcth2),
310  > wk(iaspl2),wk(iawk),inwk,ilinx,ilinth,ier)
311 c
312  if(ier.ne.0) then
313  write(6,*) ' ?tcspline: error in 2d spline, exiting.'
314  return
315  endif
316 c
317  enddo
318 c
319 c ok now fspl(*,*,1,*,*,*) have been evaluated and C2 in (x,theta)
320 c now need to extend to coeffs in phi direction.
321 c
322  xo2=0.5e0
323  xo6=1.0e0/6.0e0
324 c
325 c spline each (x,th) coeff in the phi direction
326 c
327  inpho=4*(inph-1)
328  do ith=1,inth-1
329  do ix=1,inx-1
330 c
331  do ic1=1,4
332  do ic2=1,4
333 c
334 c copy coeff. ordinates in
335 c
336  do iph=1,inph
337  wk(4*(iph-1)+1)=fspl(ic1,ic2,1,ix,ith,iph)
338  enddo
339 c
340 c use linear BC on this first pass; will correct later if
341 c necessary
342 c
343  wk(2)=0.0
344  wk(3)=0.0
345  wk(inpho+2)=0.0
346  wk(inpho+3)=0.0
347 c
348  ibcphmina=ibcphmin
349  ibcphmaxa=ibcphmax
350  if(iflg.eq.1) then
351  if((ibcphmin.eq.1).or.(ibcphmin.eq.2)) ibcphmina=0
352  if((ibcphmax.eq.1).or.(ibcphmax.eq.2)) ibcphmaxa=0
353  endif
354 c
355  call v_spline(ibcphmina,ibcphmaxa,inph,ph,wk,
356  > wk(4*inph+1))
357 c
358 c copy coeffs out
359 c
360  do iph=1,inph-1
361  fspl(ic1,ic2,2,ix,ith,iph)=wk(4*(iph-1)+2)
362  fspl(ic1,ic2,3,ix,ith,iph)=wk(4*(iph-1)+3)*xo2
363  fspl(ic1,ic2,4,ix,ith,iph)=wk(4*(iph-1)+4)*xo6
364  enddo
365 c
366  enddo ! ic2
367  enddo ! ic1
368 c
369  enddo ! ix
370  enddo ! ith
371 c
372 c if there are "non-linear" BCs requiring correction...
373 c
374 c at each (x(ix),th(ith)) get the d/dph BC's right while preserving C2
375 c everywhere...
376 c
377  if(iflg.eq.1) then
378 c
379 c first get BC correction numbers
380 c
381  iabcph1=1
382  iabcph2=iabcph1+inx*inth
383  iaccoef=iabcph2+inx*inth
384  iawk=iaccoef+12*inx*inth*inph
385  inwk=nwk-iawk+1
386 c
387  do i=1,10
388  iselect1(i)=0
389  iselect2(i)=0
390  enddo
391 c
392 c note because iflg=1, we know at least one of ibcphmin/max = 1 or 2
393 c
394  iskip1=0
395  if(ibcphmin.eq.1) then
396  iselect1(4)=1 ! df/dph
397  else if(ibcphmin.eq.2) then
398  iselect1(7)=1 ! d2f/dph2
399  else
400  iskip1=1
401  endif
402 c
403  iskip2=0
404  if(ibcphmax.eq.1) then
405  iselect2(4)=1 ! df/dph
406  else if(ibcphmax.eq.2) then
407  iselect2(7)=1 ! d2f/dph2
408  else
409  iskip2=1
410  endif
411 c
412  ia1=iabcph1-1
413  ia2=iabcph2-1
414  do ith=1,inth
415  do ix=1,inx
416  ia1=ia1+1
417  ia2=ia2+1
418 c
419  if(iskip1.eq.0) then
420  call tcspeval(x(ix),th(ith),ph(1),iselect1, zcur,
421  > x,inx,th,inth,ph,inph,ilinx,ilinth,ilinph,
422  > fspl,inf4,inf5,ier)
423  if(ier.ne.0) then
424  write(6,*) ' ?? tcspline: error in tcspeval call'
425  return
426  endif
427  wk(ia1)=bcphmin(ix,ith)-zcur(1) ! correction needed
428  else
429  wk(ia1)=z0
430  endif
431 c
432  if(iskip2.eq.0) then
433  call tcspeval(x(ix),th(ith),ph(inph),iselect2, zcur,
434  > x,inx,th,inth,ph,inph,ilinx,ilinth,ilinph,
435  > fspl,inf4,inf5,ier)
436  if(ier.ne.0) then
437  write(6,*) ' ?? tcspline: error in tcspeval call'
438  return
439  endif
440  wk(ia2)=bcphmax(ix,ith)-zcur(1) ! correction needed
441  else
442  wk(ia2)=z0
443  endif
444  enddo
445  enddo
446 c
447  call tcspcorr(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
448  > ibcxmin,ibcxmax,ibcthmin,ibcthmax,
449  > ibcphmin,wk(iabcph1),ibcphmax,wk(iabcph2),
450  > wk(iaccoef),wk(iawk),inwk)
451 c
452  endif
453 c
454  return
455  end
456 c-----------------------------------------
457  subroutine tcspcorr(x,inx,th,inth,ph,inph,fspl,inf4,inf5,
458  > ibcxmin,ibcxmax,ibcthmin,ibcthmax,
459  > ibcphmin,bcph1,ibcphmax,bcph2,ccorr,wk,nwk)
460 c
461 c intermediary routine for tcspline:
462 c do correction needed to get C2 3d spline with phi bdy conditions
463 c matched.
464 c
465 c all input unless noted:
466 c
467  real x(inx) ! x axis
468  real th(inth) ! th axis
469  real ph(inph) ! ph axis
470 c
471  real fspl(4,4,4,inf4,inf5,inph) ! spline coeffs -- adjusted
472 c
473  integer ibcxmin,ibcxmax ! x BC flags
474  integer ibcthmin,ibcthmax ! th BC flags
475  integer ibcphmin,ibcphmax ! ph BC flags
476 c
477  real bcph1(inx,inth) ! ph BC correction array @ phi(1)
478  real bcph2(inx,inth) ! ph BC correction array @ phi(inph)
479 c
480 c workspaces:
481 c
482  real ccorr(3,4,inx,inth,inph) ! correction coefficients (partial)
483 c
484  real wk(nwk) ! workspace
485 c
486 c---------------------------
487 c
488  xo2=0.5e0
489  xo6=1.0e0/6.0e0
490 c
491 c 1. splines in phi -- fcns are zero everywhere but have non-zero BCs
492 c
493  if(nwk.lt.10*max(inx,inth,inph)) then
494  write(6,*) ' ?? programming error in tcspcorr (tcspline)'
495  return
496  endif
497 c
498  z0=0.0
499  iawk2=4*inph+1
500 c
501  inpho=4*(inph-1)
502  do ith=1,inth
503  do ix=1,inx
504 c
505  do iph=1,inph
506  wk(4*(iph-1)+1)=z0
507  enddo
508 c
509 c set BC for this 1d spline correction
510 c
511  if(ibcphmin.eq.1) then
512  wk(2)=bcph1(ix,ith)
513  else if(ibcphmin.eq.2) then
514  wk(3)=bcph1(ix,ith)
515  endif
516 c
517  if(ibcphmax.eq.1) then
518  wk(inpho+2)=bcph2(ix,ith)
519  else if(ibcphmax.eq.2) then
520  wk(inpho+3)=bcph2(ix,ith)
521  endif
522 c
523  call v_spline(ibcphmin,ibcphmax,inph,ph,wk,wk(iawk2))
524 c
525 c copy non-zero coeffs out to ccorr
526 c
527  do iph=1,inph-1
528  ccorr(1,1,ix,ith,iph)=wk(4*(iph-1)+2)
529  ccorr(2,1,ix,ith,iph)=wk(4*(iph-1)+3)*xo2
530  ccorr(3,1,ix,ith,iph)=wk(4*(iph-1)+4)*xo6
531  enddo
532 c
533  enddo
534  enddo
535 c
536 c 2. spline the coeffs in x -- use ibcx flags & zero for derivative
537 c bc if necessary
538 c
539  iawk2=4*inx+1
540 c
541  inxo=4*(inx-1)
542  do iph=1,inph-1
543  do ith=1,inth
544 c
545  do icph=1,3
546 c
547  do ix=1,inx
548  wk(4*(ix-1)+1)=ccorr(icph,1,ix,ith,iph)
549  enddo
550 c
551 c zero BC: correction spline
552 c
553  wk(2)=0.0
554  wk(3)=0.0
555  wk(inxo+2)=0.0
556  wk(inxo+3)=0.0
557 c
558  call v_spline(ibcxmin,ibcxmax,inx,x,wk,wk(iawk2))
559 c
560  do ix=1,inx-1
561  ccorr(icph,2,ix,ith,iph)=wk(4*(ix-1)+2)
562  ccorr(icph,3,ix,ith,iph)=wk(4*(ix-1)+3)*xo2
563  ccorr(icph,4,ix,ith,iph)=wk(4*(ix-1)+4)*xo6
564  enddo
565 c
566  enddo
567 c
568  enddo
569  enddo
570 c
571 c 3. spline all the ccorr coefs in th -- use ibcth flags & zero for
572 c derivative correction BC if necessary
573 c
574 c add the results into fspl
575 c
576  iawk2=4*inth+1
577 c
578  intho=4*(inth-1)
579  do iph=1,inph-1
580  do ix=1,inx-1
581 c
582  do icx=1,4
583  do icph=1,3
584 c
585  do ith=1,inth
586  wk(4*(ith-1)+1)=ccorr(icph,icx,ix,ith,iph)
587  enddo
588 c
589 c zero BC: correction spline
590 c
591  wk(2)=0.0
592  wk(3)=0.0
593  wk(intho+2)=0.0
594  wk(intho+3)=0.0
595 c
596  call v_spline(ibcthmin,ibcthmax,inth,th,wk,
597  > wk(iawk2))
598 c
599  do ith=1,inth-1
600  do icth=1,4
601  zfac=1.0
602  if(icth.eq.3) zfac=xo2
603  if(icth.eq.4) zfac=xo6
604  fspl(icx,icth,icph+1,ix,ith,iph)=
605  > fspl(icx,icth,icph+1,ix,ith,iph)+
606  > wk(4*(ith-1)+icth)*zfac
607  enddo
608  enddo
609 c
610  enddo
611  enddo
612 c
613  enddo
614  enddo
615 c
616  return
617  end
618 c-----------------------------------------
619  subroutine tcsp23(x,inx,th,inth,fspl,inf4,
620  > ibcxmin,bcxmin,ibcxmax,bcxmax,
621  > ibcthmin,bcthmin,ibcthmax,bcthmax,
622  > fspl2,wk,nwk,ilinx,ilinth,ier)
623 c
624 c intermediary routines
625 c call bcspline from tcspline loop
626 c to set up 2d splines in each phi plane
627 c
628  real x(inx) ! x axis
629  real th(inth) ! th axis
630 c
631  real fspl(4,4,4,inf4,inth) ! fspl array at one phi pt.
632  real fspl2(4,4,inx,inth) ! temp fspl array for bcspline
633 c
634  real bcxmin(inth),bcxmax(inth) ! d/dx BC's @ x(1),x(inx), th(*)
635  real bcthmin(inx),bcthmax(inx) ! d/dth BC's @ th(1),th(inth), x(*)
636 c
637  real wk(nwk)
638 c
639 c--------------------
640 c
641 c 1. copy spline data in
642 c
643  do ith=1,inth
644  do ix=1,inx
645  fspl2(1,1,ix,ith)=fspl(1,1,1,ix,ith)
646  enddo
647  enddo
648 c
649 c 2. compute the 2d spline
650 c
651  call bcspline(x,inx,th,inth,fspl2,inx,
652  > ibcxmin,bcxmin,ibcxmax,bcxmax,
653  > ibcthmin,bcthmin,ibcthmax,bcthmax,
654  > wk,nwk,ilinx,ilinth,ier)
655  if(ier.ne.0) return
656 c
657 c 3. copy spline coeff results out
658 c
659  do ith=1,inth-1
660  do ix=1,inx-1
661  do j=1,4
662  do i=1,4
663  fspl(i,j,1,ix,ith)=fspl2(i,j,ix,ith)
664  enddo
665  enddo
666  enddo
667  enddo
668 c
669  return
670  end