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,
43 real f(icoeff,ixdim,iydim,izdim)
50 integer ibcxmin,ibcxmax
51 integer ibcymin,ibcymax
52 integer ibczmin,ibczmax
54 real bcxmin(iydim,izdim),bcxmax(iydim,izdim)
55 real bcymin(ixdim,izdim),bcymax(ixdim,izdim)
56 real bczmin(ixdim,iydim),bczmax(ixdim,iydim)
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)
108 real,
dimension(:,:),
allocatable :: wk2
109 real,
dimension(:,:,:),
allocatable :: wk3
122 imin=min(imin,jspline(ii))
123 imax=max(imax,jspline(ii))
124 if(jspline(ii).gt.0)
then
128 if(.not.ifound(jspline(ii)))
then
129 ifound(jspline(ii)) = .true.
134 if((imin.lt.-1).or.(imax.gt.2))
then
137 >
' ?mkintrp3d: spline type control out of range -1 to 2: ',
146 if(ifound(1).and.ifound(2))
then
149 >
' ?mkintrp3d: spline/Hermite hybrid not supported (',
155 if(imul.ne.icoeff)
then
157 >
' ?coeff dimension inconsistency for spline type codes ',
159 write(6,*)
' in mkintrp3d: expected: ',imul,
' got: ',icoeff
167 if(jspline(1).eq.-1)
then
173 if(jspline(2).eq.-1)
then
179 if(jspline(3).eq.-1)
then
185 if((ickx.ne.ixdim).or.(icky.ne.iydim).or.(ickz.ne.izdim))
then
187 >
' ?mkintrp2d: dimensioning inconsistent with '//
188 >
'interpolation controls: ',jspline
189 write(6,*)
' expected: ',ickx,icky,ickz,
190 >
'; got: ',ixdim,iydim,izdim
195 if(jspline(1).le.0)
then
196 call splinck(x,nx,idum1,ztol,ier)
198 write(6,*)
' ?mkintrp2d: x axis not strict ascending.'
203 if(jspline(2).le.0)
then
204 call splinck(y,ny,idum1,ztol,ier)
206 write(6,*)
' ?mkintrp2d: y axis not strict ascending.'
211 if(jspline(3).le.0)
then
212 call splinck(z,nz,idum1,ztol,ier)
214 write(6,*)
' ?mkintrp2d: z axis not strict ascending.'
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
234 if(ibcxmin.eq.-1)
then
236 else if((ibcxmin.eq.1).or.(ibcxmax.eq.1))
then
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
251 if(ibcymin.eq.-1)
then
253 else if((ibcymin.eq.1).or.(ibcymax.eq.1))
then
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
268 if(ibczmin.eq.-1)
then
270 else if((ibczmin.eq.1).or.(ibczmax.eq.1))
then
275 if(kspline.eq.1)
then
279 call util_bcherm3(f, ixdim, iydim, izdim,
280 > ibcxmin,ibcxmax, ibcymin,ibcymax, ibczmin,ibczmax,
281 > bcxmin, bcxmax, bcymin, bcymax, bczmin, bczmax,
284 call akherm3p(x,ixdim, y,iydim, z,izdim, f,ixdim,iydim,
285 > idum1,idum2,idum3, ipx,ipy,ipz, ier)
287 else if(kspline.eq.2)
then
290 call mktricub(x,nx,y,ny,z,nz,
292 > ibcxmin,bcxmin,ibcxmax,bcxmax,iydim,
293 > ibcymin,bcymin,ibcymax,bcymax,ixdim,
294 > ibczmin,bczmin,ibczmax,bczmax,ixdim,
295 > idum1,idum2,idum3, ier)
302 if(jspline(1).gt.0)
then
306 if(jspline(1).eq.1)
then
307 call util_bcherm1(f(1,1,ii,jj), ixdim,
309 > bcxmin(ii,jj), bcxmax(ii,jj), x)
310 call akherm1p(x,ixdim,f(1,1,ii,jj),idum1,ipx,
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),
324 else if(jspline(2).gt.0)
then
330 wk2(1,1:iydim) = f(1,ii,1:iydim,jj)
332 if(jspline(2).eq.1)
then
333 call util_bcherm1(wk2, iydim,
335 > bcymin(ii,jj), bcymax(ii,jj), y)
336 call akherm1p(y,iydim,wk2,idum2,ipy,ier)
338 else if(jspline(2).eq.2)
then
339 call mkspline(y,iydim,wk2,
340 > ibcymin,bcymin(ii,jj),
341 > ibcymax,bcymax(ii,jj),
346 f(1:2,ii,1:iydim,jj) = wk2(1:2,1:iydim)
359 wk2(1,1:izdim) = f(1,ii,jj,1:izdim)
361 if(jspline(3).eq.1)
then
362 call util_bcherm1(wk2, izdim,
364 > bczmin(ii,jj), bczmax(ii,jj), z)
365 call akherm1p(z,izdim,wk2,idum2,ipz,ier)
367 else if(jspline(3).eq.2)
then
368 call mkspline(z,izdim,wk2,
369 > ibczmin,bczmin(ii,jj),
370 > ibczmax,bczmax(ii,jj),
375 f(1:2,ii,jj,1:izdim) = wk2(1:2,1:izdim)
385 if(jspline(3).le.0)
then
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),
393 call akherm2p(x,ixdim,y,iydim,f(1,1,1,ii),ixdim,
394 > idum1,idum2,ipx,ipy,ier)
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),
404 else if(jspline(2).le.0)
then
405 allocate(wk3(4,nx,nz))
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)
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),
424 f(1:4,1:nx,ii,1:nz) = wk3
430 allocate(wk3(4,ny,nz))
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)
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),
450 f(1:4,ii,1:ny,1:nz) = wk3