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