V3FIT
r8mktricub.f
1  subroutine r8mktricub(x,nx,y,ny,z,nz,f,nf2,nf3,
2  > ibcxmin,bcxmin,ibcxmax,bcxmax,inb1x,
3  > ibcymin,bcymin,ibcymax,bcymax,inb1y,
4  > ibczmin,bczmin,ibczmax,bczmax,inb1z,
5  > ilinx,iliny,ilinz,ier)
6 c
7 c setup a tricubic spline; store coefficients in compact form
8 c (as per suggestion of L. Zakharov, PPPL, Feb. 1999)
9 C 8 coeffs per (x,y,z) grid point:
10 C f,fxx,fyy,fzz,fxxyy,fxxzz,fyyzz,fxxyyzz
11 C
12 C dmc -- modified Feb 2004 -- rewritten to compute coefficients
13 C directly rather than by conversion from the non-compact representation
14 C (to reduce cpu and memory cost)
15 C
16 C
17 C input:
18  IMPLICIT NONE
19  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
20  integer nx ! length of x vector
21  integer ny ! length of y vector
22  integer nz ! length of z vector
23  real*8 x(nx) ! x vector, strict ascending
24  real*8 y(ny) ! y vector, strict ascending
25  real*8 z(nz) ! z vector, strict ascending
26 c
27  integer nf2 ! 2nd dim. of f array, nf2.ge.nx
28  integer nf3 ! 3rd dim. of f array, nf3.ge.ny
29 c
30 c input/output:
31 c
32  real*8 f(8,nf2,nf3,nz) ! data and spline coefficients
33 c
34 C on input: f(1,i,j,k) = f(x(i),y(j),z(k))
35 C on output: f(1,i,j,k) unchanged
36 C f(2,i,j,k) = d2f/dx2(x(i),y(j),z(k))
37 C f(3,i,j,k) = d2f/dy2(x(i),y(j),z(k))
38 C f(4,i,j,k) = d2f/dz2(x(i),y(j),z(k))
39 C f(5,i,j,k) = d4f/dx2dy2(x(i),y(j),z(k))
40 C f(6,i,j,k) = d4f/dx2dz2(x(i),y(j),z(k))
41 C f(7,i,j,k) = d4f/dy2dz2(x(i),y(j),z(k))
42 C f(8,i,j,k) = d6f/dx2dy2dz2(x(i),y(j),z(k))
43 C
44 C there is a rather Hermite like interpolation formula to go with
45 C this-- see evtricub.for. Also the bicubic formula is given in
46 C mkbicubw.for; the tricubic formula is precisely analogous.
47 C
48 C boundary condition data
49 C inputs:
50  integer inb1x ! 1st dim of xmin & xmax bc arrays
51  integer inb1y ! 1st dim of ymin & ymax bc arrays
52  integer inb1z ! 1st dim of zmin & zmax bc arrays
53 C
54  integer ibcxmin,ibcxmax ! BC type flag @xmin, xmax
55  integer ibcymin,ibcymax ! BC type flag @ymin, ymax
56  integer ibczmin,ibczmax ! BC type flag @zmin, zmax
57 C
58  real*8 bcxmin(inb1x,nz),bcxmax(inb1x,nz) ! xmin & xmax BC data, ny x nz
59  real*8 bcymin(inb1y,nz),bcymax(inb1y,nz) ! ymin & ymax BC data, nx x nz
60  real*8 bczmin(inb1z,ny),bczmax(inb1z,ny) ! zmin & zmax BC data, nx x ny
61 c
62 c where BC data is not required, dummy scalars may be passed.
63 C the ibc* flags determine whether BC data isneeded.
64 c
65 c BC data: bcxmin & bcxmax: BC vs. y,z @xmin,xmax
66 C bcymin & bcymax: BC vs. x,z @ymin,ymax
67 C bczmin & bczmax: BC vs. x,y @zmin,zmax
68 C
69 c ibcxmin -- indicator for boundary condition at xmin=x(1):
70 c bcxmin(...) -- boundary condition data
71 c =-1 -- use periodic boundary condition
72 c =0 -- use "not a knot"
73 c =1 -- match slope, specified at x(1),y(iy),z(iz) by bcxmin(iy,iz)
74 c =2 -- match 2nd derivative, specified at x(1),y(iy),z(iz)
75 c by bcxmin(iy,iz
76 c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all y(j)
77 c =4 -- boundary condition is d2f/dx2=0 at x(1), all y(j)
78 c =5 -- df/dx BC from 1st divided difference
79 c =6 -- d2f/dx2 BC from 2nd divided difference (parabolic fit)
80 c =7 -- d3f/dx3 BC from 3rd divided difference (cubic fit)
81 c ***NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
82 c
83 c ibcxmax -- indicator for boundary condition at x(nx):
84 c bcxmax(...) -- boundary condition data
85 c (interpretation as with ibcxmin, bcxmin)
86 c NOTE: if ibcxmin=-1 then the periodic BC applies on both sides
87 c and ibcxmax, bcxmax are ignored.
88 c inb1x -- 1st dimension of bcxmin, bcxmax: if ibcxmin or ibcxmax .gt. 0
89 c this must be .ge. ny.
90 c
91 c interpretation of ibcymin,bcymin,ibcymax,bcymax,inb1y
92 c is same as with ibcxmin,...
93 c
94 c interpretation of ibczmin,bczmin,ibczmax,bczmax,inb1z
95 c is same as with ibcxmin,...
96 c
97 c the explicit bdy condition arrays are referenced only if the
98 c corresponding "ibc" flag values are set to 1 or 2.
99 c
100 c output:
101  integer ilinx ! x vector equal spacing flag
102  integer iliny ! y vector equal spacing flag
103  integer ilinz ! z vector equal spacing flag
104 c
105 c ilinx -- =1 on output if x(nx) pts are nearly evenly spaced (tol=1e-3)
106 c iliny -- =1 on output if y(ny) evenly spaced (tol=1e-3)
107 c ilinz -- =1 on output if z(nz) evenly spaced (tol=1e-3)
108 c
109  integer ier ! exit code
110 c ier -- completion code, 0 for normal
111 c
112 C-----------------------------------------------------
113 c workspace **dynamic allocation**
114 C f90 dynamic array
115 C
116  real*8, dimension(:,:,:), allocatable :: fbicub ! bicubic subsection
117  real*8, dimension(:,:), allocatable :: fwk ! work array
118  real*8, dimension(:), allocatable :: bcx1,bcx2,bcy1,bcy2 ! BCs for mkbicub
119 c
120  real*8, dimension(:,:,:,:), allocatable :: fcorr ! correction spline
121  real*8, dimension(:,:), allocatable :: bcc1,bcc2 ! correction BCs
122 c
123  integer iflg,ierx,iery,ierz
124  integer ix,iy,iz
125 c
126  real*8 ztol
127  real*8 zbc1,zbc2,hz
128  integer ibc1,ibc2
129 c
130  data ztol/1.0e-3_r8/
131 c-----------------------------------------------------
132 c
133  ier=0
134 c
135  iflg=0
136 c
137 c check z bdy condition "linearity"
138 c
139  if(ibczmin.ne.-1) then
140  if((ibczmin.eq.1).or.(ibczmin.eq.2)) then
141  do iy=1,ny
142  do ix=1,nx
143  if(bczmin(ix,iy).ne.0.0_r8) iflg=1
144  enddo
145  enddo
146  endif
147  if((ibczmax.eq.1).or.(ibczmax.eq.2)) then
148  do iy=1,ny
149  do ix=1,nx
150  if(bczmax(ix,iy).ne.0.0_r8) iflg=1
151  enddo
152  enddo
153  endif
154  endif
155 c
156  if(nx.lt.2) then
157  write(6,'('' ?mktricub: at least 2 x points required.'')')
158  ier=1
159  endif
160  if(ny.lt.2) then
161  write(6,'('' ?mktricub: need at least 2 y points.'')')
162  ier=1
163  endif
164  if(nz.lt.2) then
165  write(6,'('' ?mktricub: need at least 2 z points.'')')
166  ier=1
167  endif
168 c
169  if((ibcxmin.eq.1).or.(ibcxmax.eq.1).or.(ibcxmin.eq.2).or.
170  > (ibcxmax.eq.2)) then
171  if(inb1x.lt.ny) then
172  ier=1
173  write(6,
174  >'('.lt.' ?mktricub: 1st dim of bcxmin/max arrays ny'')')
175  endif
176  endif
177 c
178  if((ibcymin.eq.1).or.(ibcymax.eq.1).or.(ibcymin.eq.2).or.
179  > (ibcymax.eq.2)) then
180  if(inb1y.lt.nx) then
181  ier=1
182  write(6,
183  >'('.lt.' ?mktricub: 1st dim of bcymin/max arrays nx'')')
184  endif
185  endif
186 c
187  if((ibczmin.eq.1).or.(ibczmax.eq.1).or.(ibczmin.eq.2).or.
188  > (ibczmax.eq.2)) then
189  if(inb1z.lt.nx) then
190  ier=1
191  write(6,
192  >'('.lt.' ?mktricub: 1st dim of bczmin/max arrays nx'')')
193  endif
194  endif
195 c
196  call ibc_ck(ibcxmin,'mktricub','xmin',-1,7,ier)
197  if(ibcxmin.ge.0) call ibc_ck(ibcxmax,'mktricub','xmax',0,7,ier)
198 c
199  call ibc_ck(ibcymin,'mktricub','ymin',-1,7,ier)
200  if(ibcymin.ge.0) call ibc_ck(ibcymax,'mktricub','ymax',0,7,ier)
201 c
202  call ibc_ck(ibczmin,'mktricub','zmin',-1,7,ier)
203  if(ibczmax.ge.0) call ibc_ck(ibczmax,'mktricub','zmax',0,7,ier)
204 c
205 c check ilinx & x vector
206 c
207  call r8splinck(x,nx,ilinx,ztol,ierx)
208  if(ierx.ne.0) ier=2
209 c
210  if(ier.eq.2) then
211  write(6,'('' ?mktricub: x axis not strict ascending'')')
212  endif
213 c
214 c check iliny & y vector
215 c
216  call r8splinck(y,ny,iliny,ztol,iery)
217  if(iery.ne.0) ier=3
218 c
219  if(ier.eq.3) then
220  write(6,'('' ?mktricub: y axis not strict ascending'')')
221  endif
222 c
223 c check ilinz & z vector
224 c
225  call r8splinck(z,nz,ilinz,ztol,ierz)
226  if(ierz.ne.0) ier=4
227 c
228  if(ier.eq.4) then
229  write(6,'('' ?mktricub: z axis not strict ascending'')')
230  endif
231 c
232  if(ier.ne.0) return
233 c
234 c------------------------------------
235 c 1. compute (x,y) bicubic splines using mkbicub
236 c
237  allocate(fbicub(4,nx,ny))
238  allocate(bcx1(ny),bcx2(ny),bcy1(nx),bcy2(nx))
239  bcx1=0.0; bcx2=0.0; bcy1=0.0; bcy2=0.0_r8
240 c
241  do iz=1,nz
242  if(ibcxmin.ne.-1) then
243  if((ibcxmin.eq.1).or.(ibcxmin.eq.2)) then
244  bcx1(1:ny)=bcxmin(1:ny,iz)
245  endif
246  if((ibcxmax.eq.1).or.(ibcxmax.eq.2)) then
247  bcx2(1:ny)=bcxmax(1:ny,iz)
248  endif
249  endif
250  if(ibcymin.ne.-1) then
251  if((ibcymin.eq.1).or.(ibcymin.eq.2)) then
252  bcy1(1:nx)=bcymin(1:nx,iz)
253  endif
254  if((ibcymax.eq.1).or.(ibcymax.eq.2)) then
255  bcy2(1:nx)=bcymax(1:nx,iz)
256  endif
257  endif
258 c
259  fbicub(1,1:nx,1:ny) = f(1,1:nx,1:ny,iz)
260 c
261  call r8mkbicub(x,nx,y,ny,fbicub,nx,
262  > ibcxmin,bcx1,ibcxmax,bcx2,
263  > ibcymin,bcy1,ibcymax,bcy2,
264  > ilinx,iliny,ier)
265  if(ier.ne.0) return
266 c
267  f(2:3,1:nx,1:ny,iz) = fbicub(2:3,1:nx,1:ny) ! fxx, fyy
268  f(5,1:nx,1:ny,iz) = fbicub(4,1:nx,1:ny) ! fxxyy
269 c
270  enddo
271 c
272  deallocate(fbicub,bcx1,bcx2,bcy1,bcy2)
273 c
274 c 2. homogeneous spline in z direction; inhomogeneous BC imposed later
275 c if necessary
276 c
277  zbc1=0.0_r8
278  zbc2=0.0_r8
279  ibc1=ibczmin
280  ibc2=ibczmax
281  if(iflg.eq.1) then
282  if((ibczmin.eq.1).or.(ibczmin.eq.2)) ibc1=0
283  if((ibczmax.eq.1).or.(ibczmax.eq.2)) ibc2=0
284  endif
285 c
286  allocate(fwk(2,nz))
287 c
288  do iy=1,ny
289  do ix=1,nx
290 
291  fwk(1,1:nz) = f(1,ix,iy,1:nz)
292  call r8mkspline(z,nz,fwk,
293  > ibc1,zbc1,ibc2,zbc2,ilinz,ier)
294  if(ier.ne.0) return
295  f(4,ix,iy,1:nz) = fwk(2,1:nz) ! fzz
296 
297  fwk(1,1:nz) = f(2,ix,iy,1:nz)
298  call r8mkspline(z,nz,fwk,
299  > ibc1,zbc1,ibc2,zbc2,ilinz,ier)
300  if(ier.ne.0) return
301  f(6,ix,iy,1:nz) = fwk(2,1:nz) ! fxxzz
302 
303  fwk(1,1:nz) = f(3,ix,iy,1:nz)
304  call r8mkspline(z,nz,fwk,
305  > ibc1,zbc1,ibc2,zbc2,ilinz,ier)
306  if(ier.ne.0) return
307  f(7,ix,iy,1:nz) = fwk(2,1:nz) ! fyyzz
308 
309  fwk(1,1:nz) = f(5,ix,iy,1:nz)
310  call r8mkspline(z,nz,fwk,
311  > ibc1,zbc1,ibc2,zbc2,ilinz,ier)
312  if(ier.ne.0) return
313  f(8,ix,iy,1:nz) = fwk(2,1:nz) ! fxxyyzz
314 
315  enddo
316  enddo
317 c
318  deallocate(fwk)
319 c
320  if(iflg.eq.1) then
321 c
322 c 3. inhomogeneous BC correction
323 c
324  allocate(fwk(2,max(nx,ny,nz)))
325  allocate(bcc1(nx,ny),bcc2(nx,ny))
326  allocate(fcorr(4,nx,ny,nz))
327 c
328 c correction BCs
329 c
330  do iy=1,ny
331  do ix=1,nx
332  bcc1(ix,iy)=0.0_r8
333  if(ibczmin.eq.1) then
334  hz=z(2)-z(1)
335  bcc1(ix,iy)=(f(1,ix,iy,2)-f(1,ix,iy,1))/hz +
336  > hz*(-2*f(4,ix,iy,1)-f(4,ix,iy,2))/6
337  bcc1(ix,iy)=bczmin(ix,iy)-bcc1(ix,iy)
338  else if(ibczmin.eq.2) then
339  bcc1(ix,iy)=bczmin(ix,iy)-f(4,ix,iy,1)
340  endif
341  enddo
342  enddo
343 c
344  do iy=1,ny
345  do ix=1,nx
346  bcc2(ix,iy)=0.0_r8
347  if(ibczmax.eq.1) then
348  hz=z(2)-z(1)
349  bcc2(ix,iy)=(f(1,ix,iy,2)-f(1,ix,iy,1))/hz +
350  > hz*(-2*f(4,ix,iy,1)-f(4,ix,iy,2))/6
351  bcc2(ix,iy)=bczmax(ix,iy)-bcc2(ix,iy)
352  else if(ibczmax.eq.2) then
353  bcc2(ix,iy)=bczmax(ix,iy)-f(4,ix,iy,1)
354  endif
355  enddo
356  enddo
357 c
358  fwk(1,1:nz)=0.0_r8 ! values are all zero, only BC is set...
359  do iy=1,ny
360  do ix=1,nx
361  call r8mkspline(z,nz,fwk,
362  > ibczmin,bcc1(ix,iy),ibczmax,bcc2(ix,iy),ilinz,ier)
363  if(ier.ne.0) return
364  fcorr(1,ix,iy,1:nz)=fwk(2,1:nz) ! fzz-correction
365  enddo
366  enddo
367 c
368 c higher order corrections
369 c
370  zbc1=0.0_r8
371  zbc2=0.0_r8
372 c
373  do iz=1,nz
374  do iy=1,ny
375  fwk(1,1:nx)=fcorr(1,1:nx,iy,iz)
376  call r8mkspline(x,nx,fwk,
377  > ibcxmin,zbc1,ibcxmax,zbc2,ilinx,ier)
378  if(ier.ne.0) return
379  fcorr(2,1:nx,iy,iz)=fwk(2,1:nx) ! fxxzz-correction
380  enddo
381  enddo
382 c
383  do iz=1,nz
384  do ix=1,nx
385  fwk(1,1:ny)=fcorr(1,ix,1:ny,iz)
386  call r8mkspline(y,ny,fwk,
387  > ibcymin,zbc1,ibcymax,zbc2,iliny,ier)
388  if(ier.ne.0) return
389  fcorr(3,ix,1:ny,iz)=fwk(2,1:ny) ! fyyzz-correction
390 
391  fwk(1,1:ny)=fcorr(2,ix,1:ny,iz)
392  call r8mkspline(y,ny,fwk,
393  > ibcymin,zbc1,ibcymax,zbc2,iliny,ier)
394  if(ier.ne.0) return
395  fcorr(4,ix,1:ny,iz)=fwk(2,1:ny) ! fxxyyzz-correction
396  enddo
397  enddo
398 c
399 c apply correction
400 c
401  do iz=1,nz
402  do iy=1,ny
403  do ix=1,nx
404  f(4,ix,iy,iz)=f(4,ix,iy,iz)+fcorr(1,ix,iy,iz)
405  f(6,ix,iy,iz)=f(6,ix,iy,iz)+fcorr(2,ix,iy,iz)
406  f(7,ix,iy,iz)=f(7,ix,iy,iz)+fcorr(3,ix,iy,iz)
407  f(8,ix,iy,iz)=f(8,ix,iy,iz)+fcorr(4,ix,iy,iz)
408  enddo
409  enddo
410  enddo
411 c
412  deallocate(fwk,fcorr,bcc1,bcc2)
413 c
414  endif
415 c
416 c that's all
417 c
418  return
419  end
420 
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17