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