V3FIT
mkintrp3d.f
1  subroutine mkintrp3d(x,nx,y,ny,z,nz,jspline,
2  > f,icoeff,ixdim,iydim,izdim,
3  > ibcxmin,bcxmin,ibcxmax,bcxmax,
4  > ibcymin,bcymin,ibcymax,bcymax,
5  > ibczmin,bczmin,ibczmax,bczmax,
6  > ier)
7 c
8 c setup a tricubic spline, or tricubic Hermite, or hybrid linear/zonal 2d or
9 c 1d with 1d or 2d cubic or Hermite spline interpolation
10 c
11  implicit NONE
12 C
13 C input:
14  integer nx ! length of x vector
15  integer ny ! length of y vector
16  integer nz ! length of z vector
17  real x(nx) ! x vector, strict ascending
18  real y(ny) ! y vector, strict ascending
19  real z(nz) ! z vector, strict ascending
20 c
21  integer :: jspline(3) ! interpolation method control
22 C (1) -- 1st dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline
23 C (2) -- 2nd dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline
24 C (3) -- 3rd dimension: -1: zone, 0: pclin, 1: Hermite, 2: spline
25 C
26 C Standard interpolation-- all jspline values match
27 C e.g. jspline(1)=jspline(2)=jspline(3)=2 for tricubic spline
28 C
29 C Hybrid interpolation-- not all jspline values are the same.
30 C
31 C RESTRICTION: if any jspline(...) element has value 1, none can have
32 C value 2. I.e. Spline and Hermite interpolation cannot currently be mixed.
33 C This restriction exists because of technical issues in the
34 C implementation (it could be removed in principle but the work to do
35 C this has not been scheduled).
36 c
37 c coefficient buffer dimensions
38  integer :: icoeff ! #coefficients per data point
39  integer :: ixdim ! nx; nx-1 if jspline(1)==-1
40  integer :: iydim ! ny; ny-1 if jspline(2)==-1
41  integer :: izdim ! nz; nz-1 if jspline(3)==-1
42 c input/output:
43  real f(icoeff,ixdim,iydim,izdim) ! data and spline coefficients
44 c
45 C boundary condition data
46 C Contiguous storage is assumed-- 1st dimension size must match
47 C actual use
48 C
49 C
50  integer ibcxmin,ibcxmax ! BC type flag @xmin, xmax
51  integer ibcymin,ibcymax ! BC type flag @ymin, ymax
52  integer ibczmin,ibczmax ! BC type flag @zmin, zmax
53 C
54  real bcxmin(iydim,izdim),bcxmax(iydim,izdim) ! xmin & xmax BC data
55  real bcymin(ixdim,izdim),bcymax(ixdim,izdim) ! ymin & ymax BC data
56  real bczmin(ixdim,iydim),bczmax(ixdim,iydim) ! zmin & zmax BC data
57 c
58 c where BC data is not required, dummy scalars may be passed.
59 C the ibc* flags determine whether BC data isneeded.
60 c
61 c BC data not required for zonal or piecewise linear interpolation
62 c for Hermite interpolation ibc* values from set {-1,0,1} are accepted.
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
88 c interpretation of ibcymin,bcymin,ibcymax,bcymax
89 c is same as with ibcxmin,...
90 c
91 c interpretation of ibczmin,bczmin,ibczmax,bczmax
92 c is same as with ibcxmin,...
93 c
94 c the explicit bdy condition arrays are referenced only if the
95 c corresponding "ibc" flag values are set > 0.
96 c
97  integer ier ! exit code
98 c ier -- completion code, 0 for normal
99 c
100 C-----------------------
101  integer :: kspline
102  integer :: ii,jj,imul,imin,imax,ickx,icky,ickz,inum
103  integer :: idum1,idum2,idum3,idimcu
104  integer :: ipx,ipy,ipz
105  real :: ztol = 0.0001
106  logical :: ifound(-1:2)
107 
108  real, dimension(:,:), allocatable :: wk2
109  real, dimension(:,:,:), allocatable :: wk3
110 C-----------------------
111 c
112  ier=0
113 c
114  imin=3
115  imax=-2
116  imul=1
117  inum=0
118  idimcu=0
119  ifound = .false.
120 
121  do ii=1,3
122  imin=min(imin,jspline(ii))
123  imax=max(imax,jspline(ii))
124  if(jspline(ii).gt.0) then
125  idimcu=idimcu+1
126  imul=imul*2
127  endif
128  if(.not.ifound(jspline(ii))) then
129  ifound(jspline(ii)) = .true.
130  inum = inum+1
131  endif
132  enddo
133 
134  if((imin.lt.-1).or.(imax.gt.2)) then
135  ier = 1
136  write(6,*)
137  > ' ?mkintrp3d: spline type control out of range -1 to 2: ',
138  > jspline
139  endif
140  if(ier.ne.0) return
141 
142  if(inum.eq.1) then
143  kspline=imin ! same interp type on all dimensions
144  else
145  kspline=-99 ! hybrid
146  if(ifound(1).and.ifound(2)) then
147  ier = 1
148  write(6,*)
149  > ' ?mkintrp3d: spline/Hermite hybrid not supported (',
150  > jspline,')'
151  endif
152  endif
153  if(ier.ne.0) return
154 c
155  if(imul.ne.icoeff) then
156  write(6,*)
157  > ' ?coeff dimension inconsistency for spline type codes ',
158  > jspline
159  write(6,*) ' in mkintrp3d: expected: ',imul,' got: ',icoeff
160  ier=1
161  return
162  endif
163 c
164 c
165 c check dimensioning consistency
166 c
167  if(jspline(1).eq.-1) then
168  ickx=nx-1
169  else
170  ickx=nx
171  endif
172 c
173  if(jspline(2).eq.-1) then
174  icky=ny-1
175  else
176  icky=ny
177  endif
178 c
179  if(jspline(3).eq.-1) then
180  ickz=nz-1
181  else
182  ickz=nz
183  endif
184 c
185  if((ickx.ne.ixdim).or.(icky.ne.iydim).or.(ickz.ne.izdim)) then
186  write(6,*)
187  > ' ?mkintrp2d: dimensioning inconsistent with '//
188  > 'interpolation controls: ',jspline
189  write(6,*) ' expected: ',ickx,icky,ickz,
190  > '; got: ',ixdim,iydim,izdim
191  ier=1
192  return
193  endif
194 c
195  if(jspline(1).le.0) then
196  call splinck(x,nx,idum1,ztol,ier)
197  if(ier.ne.0) then
198  write(6,*) ' ?mkintrp2d: x axis not strict ascending.'
199  return
200  endif
201  endif
202 c
203  if(jspline(2).le.0) then
204  call splinck(y,ny,idum1,ztol,ier)
205  if(ier.ne.0) then
206  write(6,*) ' ?mkintrp2d: y axis not strict ascending.'
207  return
208  endif
209  endif
210 c
211  if(jspline(3).le.0) then
212  call splinck(z,nz,idum1,ztol,ier)
213  if(ier.ne.0) then
214  write(6,*) ' ?mkintrp2d: z axis not strict ascending.'
215  return
216  endif
217  endif
218 c
219 c if no work to be done: exit now
220  if(imul.eq.1) return
221 c
222 c check Hermite BCs if necessary
223 c
224  if(jspline(1).eq.1) then
225  if((min(ibcxmin,ibcxmax).lt.-1).or.
226  > (max(ibcxmin,ibcxmax).gt.1)) then
227  write(6,*) ' ?mkintrp2d: Bdy Cond code out of range for'
228  write(6,*) ' Hermite interpolation; (-1:1) allowed, '//
229  > 'found: ',ibcxmin,ibcxmax
230  ier=1
231  return
232  endif
233  ipx=0
234  if(ibcxmin.eq.-1) then
235  ipx=1
236  else if((ibcxmin.eq.1).or.(ibcxmax.eq.1)) then
237  ipx=2
238  endif
239  endif
240 c
241  if(jspline(2).eq.1) then
242  if((min(ibcymin,ibcymax).lt.-1).or.
243  > (max(ibcymin,ibcymax).gt.1)) then
244  write(6,*) ' ?mkintrp2d: Bdy Cond code out of range for'
245  write(6,*) ' Hermite interpolation; (-1:1) allowed, '//
246  > 'found: ',ibcymin,ibcymax
247  ier=1
248  return
249  endif
250  ipy=0
251  if(ibcymin.eq.-1) then
252  ipy=1
253  else if((ibcymin.eq.1).or.(ibcymax.eq.1)) then
254  ipy=2
255  endif
256  endif
257 c
258  if(jspline(3).eq.1) then
259  if((min(ibczmin,ibczmax).lt.-1).or.
260  > (max(ibczmin,ibczmax).gt.1)) then
261  write(6,*) ' ?mkintrp2d: Bdy Cond code out of range for'
262  write(6,*) ' Hermite interpolation; (-1:1) allowed, '//
263  > 'found: ',ibczmin,ibczmax
264  ier=1
265  return
266  endif
267  ipz=0
268  if(ibczmin.eq.-1) then
269  ipz=1
270  else if((ibczmin.eq.1).or.(ibczmax.eq.1)) then
271  ipz=2
272  endif
273  endif
274 c
275  if(kspline.eq.1) then
276  ! tricubic Hermite
277 
278  ! put the BCs inside the function data at the right locations...
279  call util_bcherm3(f, ixdim, iydim, izdim,
280  > ibcxmin,ibcxmax, ibcymin,ibcymax, ibczmin,ibczmax,
281  > bcxmin, bcxmax, bcymin, bcymax, bczmin, bczmax,
282  > x, y, z)
283 
284  call akherm3p(x,ixdim, y,iydim, z,izdim, f,ixdim,iydim,
285  > idum1,idum2,idum3, ipx,ipy,ipz, ier)
286 
287  else if(kspline.eq.2) then
288  ! tricubic Spline
289 
290  call mktricub(x,nx,y,ny,z,nz,
291  > f,ixdim,iydim,
292  > ibcxmin,bcxmin,ibcxmax,bcxmax,iydim,
293  > ibcymin,bcymin,ibcymax,bcymax,ixdim,
294  > ibczmin,bczmin,ibczmax,bczmax,ixdim,
295  > idum1,idum2,idum3, ier)
296 
297  else
298  ! Hybrid
299  if(idimcu.eq.1) then
300  ! cubic along 1 dimension; other dims are step or pclin
301 
302  if(jspline(1).gt.0) then
303  ! cubic in x direction
304  do jj=1,izdim
305  do ii=1,iydim
306  if(jspline(1).eq.1) then
307  call util_bcherm1(f(1,1,ii,jj), ixdim,
308  > ibcxmin, ibcxmax,
309  > bcxmin(ii,jj), bcxmax(ii,jj), x)
310  call akherm1p(x,ixdim,f(1,1,ii,jj),idum1,ipx,
311  > ier)
312 
313  else if(jspline(1).eq.2) then
314  call mkspline(x,ixdim,f(1,1,ii,jj),
315  > ibcxmin,bcxmin(ii,jj),
316  > ibcxmax,bcxmax(ii,jj),
317  > idum1,ier)
318  endif
319  if(ier.ne.0) exit
320  enddo
321  if(ier.ne.0) exit
322  enddo
323 
324  else if(jspline(2).gt.0) then
325  ! cubic in y direction
326  allocate(wk2(2,ny))
327 
328  do jj=1,izdim
329  do ii=1,ixdim
330  wk2(1,1:iydim) = f(1,ii,1:iydim,jj)
331  wk2(2,1:iydim) = 0.0
332  if(jspline(2).eq.1) then
333  call util_bcherm1(wk2, iydim,
334  > ibcymin, ibcymax,
335  > bcymin(ii,jj), bcymax(ii,jj), y)
336  call akherm1p(y,iydim,wk2,idum2,ipy,ier)
337 
338  else if(jspline(2).eq.2) then
339  call mkspline(y,iydim,wk2,
340  > ibcymin,bcymin(ii,jj),
341  > ibcymax,bcymax(ii,jj),
342  > idum2,ier)
343 
344  endif
345  if(ier.ne.0) exit
346  f(1:2,ii,1:iydim,jj) = wk2(1:2,1:iydim)
347  enddo
348  if(ier.ne.0) exit
349  enddo
350 
351  deallocate(wk2)
352 
353  else
354  ! cubic in z direction
355  allocate(wk2(2,nz))
356 
357  do jj=1,iydim
358  do ii=1,ixdim
359  wk2(1,1:izdim) = f(1,ii,jj,1:izdim)
360  wk2(2,1:izdim) = 0.0
361  if(jspline(3).eq.1) then
362  call util_bcherm1(wk2, izdim,
363  > ibczmin, ibczmax,
364  > bczmin(ii,jj), bczmax(ii,jj), z)
365  call akherm1p(z,izdim,wk2,idum2,ipz,ier)
366 
367  else if(jspline(3).eq.2) then
368  call mkspline(z,izdim,wk2,
369  > ibczmin,bczmin(ii,jj),
370  > ibczmax,bczmax(ii,jj),
371  > idum2,ier)
372 
373  endif
374  if(ier.ne.0) exit
375  f(1:2,ii,jj,1:izdim) = wk2(1:2,1:izdim)
376  enddo
377  if(ier.ne.0) exit
378  enddo
379 
380  deallocate(wk2)
381 
382  endif
383  else
384  ! cubic along 2 dimensions
385  if(jspline(3).le.0) then
386  do ii=1,izdim
387  if(jspline(1).eq.1) then
388  call util_bcherm2(f(1,1,1,ii),ixdim,iydim,
389  > ibcxmin, ibcxmax, ibcymin, ibcymax,
390  > bcxmin(1,ii), bcxmax(1,ii),
391  > bcymin(1,ii), bcymax(1,ii),
392  > x,y)
393  call akherm2p(x,ixdim,y,iydim,f(1,1,1,ii),ixdim,
394  > idum1,idum2,ipx,ipy,ier)
395  else
396  call mkbicub(x,ixdim,y,iydim,f(1,1,1,ii),ixdim,
397  > ibcxmin, bcxmin(1,ii), ibcxmax, bcxmax(1,ii),
398  > ibcymin, bcymin(1,ii), ibcymax, bcymax(1,ii),
399  > idum1,idum2,ier)
400  endif
401  if(ier.ne.0) exit
402  enddo
403 
404  else if(jspline(2).le.0) then
405  allocate(wk3(4,nx,nz))
406 
407  do ii=1,iydim
408  wk3 = f(1:4,1:nx,ii,1:nz)
409  if(jspline(1).eq.1) then
410  call util_bcherm2(wk3,ixdim,izdim,
411  > ibcxmin, ibcxmax, ibczmin, ibczmax,
412  > bcxmin(ii,1:nz), bcxmax(ii,1:nz),
413  > bczmin(1,ii), bczmax(1,ii), x, z)
414  call akherm2p(x,ixdim,z,izdim,wk3,ixdim,
415  > idum1,idum2,ipx,ipz,ier)
416  else
417  call mkbicub(x,ixdim,z,izdim,wk3,ixdim,
418  > ibcxmin,bcxmin(ii,1:nz),
419  > ibcxmax,bcxmax(ii,1:nz),
420  > ibczmin,bczmin(1,ii), ibczmax,bczmax(1,ii),
421  > idum1,idum2,ier)
422  endif
423  if(ier.ne.0) exit
424  f(1:4,1:nx,ii,1:nz) = wk3
425  enddo
426 
427  deallocate(wk3)
428 
429  else
430  allocate(wk3(4,ny,nz))
431 
432  do ii=1,ixdim
433  wk3 = f(1:4,ii,1:ny,1:nz)
434  if(jspline(2).eq.1) then
435  call util_bcherm2(wk3,iydim,izdim,
436  > ibcymin, ibcymax, ibczmin, ibczmax,
437  > bcymin(ii,1:nz), bcymax(ii,1:nz),
438  > bczmin(ii,1:ny), bczmax(ii,1:ny), y, z)
439  call akherm2p(y,iydim,z,izdim,wk3,iydim,
440  > idum1,idum2,ipy,ipz,ier)
441  else
442  call mkbicub(y,iydim,z,izdim,wk3,iydim,
443  > ibcymin,bcymin(ii,1:nz),
444  > ibcymax,bcymax(ii,1:nz),
445  > ibczmin,bczmin(ii,1:ny),
446  > ibczmax,bczmax(ii,1:ny),
447  > idum1,idum2,ier)
448  endif
449  if(ier.ne.0) exit
450  f(1:4,ii,1:ny,1:nz) = wk3
451  enddo
452 
453  deallocate(wk3)
454 
455  endif
456  endif
457 
458  endif
459 
460  return
461  end