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