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)
19 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
32 real*8 f(8,nf2,nf3,nz)
54 integer ibcxmin,ibcxmax
55 integer ibcymin,ibcymax
56 integer ibczmin,ibczmax
58 real*8 bcxmin(inb1x,nz),bcxmax(inb1x,nz)
59 real*8 bcymin(inb1y,nz),bcymax(inb1y,nz)
60 real*8 bczmin(inb1z,ny),bczmax(inb1z,ny)
116 real*8,
dimension(:,:,:),
allocatable :: fbicub
117 real*8,
dimension(:,:),
allocatable :: fwk
118 real*8,
dimension(:),
allocatable :: bcx1,bcx2,bcy1,bcy2
120 real*8,
dimension(:,:,:,:),
allocatable :: fcorr
121 real*8,
dimension(:,:),
allocatable :: bcc1,bcc2
123 integer iflg,ierx,iery,ierz
139 if(ibczmin.ne.-1)
then
140 if((ibczmin.eq.1).or.(ibczmin.eq.2))
then
143 if(bczmin(ix,iy).ne.0.0_r8) iflg=1
147 if((ibczmax.eq.1).or.(ibczmax.eq.2))
then
150 if(bczmax(ix,iy).ne.0.0_r8) iflg=1
157 write(6,
'('' ?mktricub: at least 2 x points required.'')')
161 write(6,
'('' ?mktricub: need at least 2 y points.'')')
165 write(6,
'('' ?mktricub: need at least 2 z points.'')')
169 if((ibcxmin.eq.1).or.(ibcxmax.eq.1).or.(ibcxmin.eq.2).or.
170 > (ibcxmax.eq.2))
then
174 >
'('.lt.
' ?mktricub: 1st dim of bcxmin/max arrays ny'')')
178 if((ibcymin.eq.1).or.(ibcymax.eq.1).or.(ibcymin.eq.2).or.
179 > (ibcymax.eq.2))
then
183 >
'('.lt.
' ?mktricub: 1st dim of bcymin/max arrays nx'')')
187 if((ibczmin.eq.1).or.(ibczmax.eq.1).or.(ibczmin.eq.2).or.
188 > (ibczmax.eq.2))
then
192 >
'('.lt.
' ?mktricub: 1st dim of bczmin/max arrays nx'')')
196 call ibc_ck(ibcxmin,
'mktricub',
'xmin',-1,7,ier)
197 if(ibcxmin.ge.0)
call ibc_ck(ibcxmax,
'mktricub',
'xmax',0,7,ier)
199 call ibc_ck(ibcymin,
'mktricub',
'ymin',-1,7,ier)
200 if(ibcymin.ge.0)
call ibc_ck(ibcymax,
'mktricub',
'ymax',0,7,ier)
202 call ibc_ck(ibczmin,
'mktricub',
'zmin',-1,7,ier)
203 if(ibczmax.ge.0)
call ibc_ck(ibczmax,
'mktricub',
'zmax',0,7,ier)
207 call r8splinck(x,nx,ilinx,ztol,ierx)
211 write(6,
'('' ?mktricub: x axis not strict ascending'')')
216 call r8splinck(y,ny,iliny,ztol,iery)
220 write(6,
'('' ?mktricub: y axis not strict ascending'')')
225 call r8splinck(z,nz,ilinz,ztol,ierz)
229 write(6,
'('' ?mktricub: z axis not strict ascending'')')
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
242 if(ibcxmin.ne.-1)
then
243 if((ibcxmin.eq.1).or.(ibcxmin.eq.2))
then
244 bcx1(1:ny)=bcxmin(1:ny,iz)
246 if((ibcxmax.eq.1).or.(ibcxmax.eq.2))
then
247 bcx2(1:ny)=bcxmax(1:ny,iz)
250 if(ibcymin.ne.-1)
then
251 if((ibcymin.eq.1).or.(ibcymin.eq.2))
then
252 bcy1(1:nx)=bcymin(1:nx,iz)
254 if((ibcymax.eq.1).or.(ibcymax.eq.2))
then
255 bcy2(1:nx)=bcymax(1:nx,iz)
259 fbicub(1,1:nx,1:ny) = f(1,1:nx,1:ny,iz)
261 call r8mkbicub(x,nx,y,ny,fbicub,nx,
262 > ibcxmin,bcx1,ibcxmax,bcx2,
263 > ibcymin,bcy1,ibcymax,bcy2,
267 f(2:3,1:nx,1:ny,iz) = fbicub(2:3,1:nx,1:ny)
268 f(5,1:nx,1:ny,iz) = fbicub(4,1:nx,1:ny)
272 deallocate(fbicub,bcx1,bcx2,bcy1,bcy2)
282 if((ibczmin.eq.1).or.(ibczmin.eq.2)) ibc1=0
283 if((ibczmax.eq.1).or.(ibczmax.eq.2)) ibc2=0
291 fwk(1,1:nz) = f(1,ix,iy,1:nz)
292 call r8mkspline(z,nz,fwk,
293 > ibc1,zbc1,ibc2,zbc2,ilinz,ier)
295 f(4,ix,iy,1:nz) = fwk(2,1:nz)
297 fwk(1,1:nz) = f(2,ix,iy,1:nz)
298 call r8mkspline(z,nz,fwk,
299 > ibc1,zbc1,ibc2,zbc2,ilinz,ier)
301 f(6,ix,iy,1:nz) = fwk(2,1:nz)
303 fwk(1,1:nz) = f(3,ix,iy,1:nz)
304 call r8mkspline(z,nz,fwk,
305 > ibc1,zbc1,ibc2,zbc2,ilinz,ier)
307 f(7,ix,iy,1:nz) = fwk(2,1:nz)
309 fwk(1,1:nz) = f(5,ix,iy,1:nz)
310 call r8mkspline(z,nz,fwk,
311 > ibc1,zbc1,ibc2,zbc2,ilinz,ier)
313 f(8,ix,iy,1:nz) = fwk(2,1:nz)
324 allocate(fwk(2,max(nx,ny,nz)))
325 allocate(bcc1(nx,ny),bcc2(nx,ny))
326 allocate(fcorr(4,nx,ny,nz))
333 if(ibczmin.eq.1)
then
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)
347 if(ibczmax.eq.1)
then
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)
361 call r8mkspline(z,nz,fwk,
362 > ibczmin,bcc1(ix,iy),ibczmax,bcc2(ix,iy),ilinz,ier)
364 fcorr(1,ix,iy,1:nz)=fwk(2,1:nz)
375 fwk(1,1:nx)=fcorr(1,1:nx,iy,iz)
376 call r8mkspline(x,nx,fwk,
377 > ibcxmin,zbc1,ibcxmax,zbc2,ilinx,ier)
379 fcorr(2,1:nx,iy,iz)=fwk(2,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)
389 fcorr(3,ix,1:ny,iz)=fwk(2,1:ny)
391 fwk(1,1:ny)=fcorr(2,ix,1:ny,iz)
392 call r8mkspline(y,ny,fwk,
393 > ibcymin,zbc1,ibcymax,zbc2,iliny,ier)
395 fcorr(4,ix,1:ny,iz)=fwk(2,1:ny)
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)
412 deallocate(fwk,fcorr,bcc1,bcc2)