V3FIT
r8mktricubw.f
1  subroutine r8mktricubw(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  > wk,nwk,ilinx,iliny,ilinz,ier)
6 c
7 c setup a tricubic spline; store coefficients in compatct 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 input:
13 !============
14 ! idecl: explicitize implicit INTEGER declarations:
15  IMPLICIT NONE
16  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
17  INTEGER itest,iadfp,isiz1,iadfw,inwk
18 !============
19 ! idecl: explicitize implicit REAL declarations:
20  real*8 hxlast,hylast,hzlast
21 !============
22  integer nx ! length of x vector
23  integer ny ! length of y vector
24  integer nz ! length of z vector
25  REAL*8 x(nx) ! x vector, strict ascending
26  real*8 y(ny) ! y vector, strict ascending
27  real*8 z(nz) ! z vector, strict ascending
28 c
29  integer nf2 ! 2nd dim. of f array, nf2.ge.nx
30  integer nf3 ! 3rd dim. of f array, nf3.ge.ny
31 c
32 c input/output:
33 c
34  real*8 f(8,nf2,nf3,nz) ! data and spline coefficients
35 c
36 C on input: f(1,i,j,k) = f(x(i),y(j),z(k))
37 C on output: f(1,i,j,k) unchanged
38 C f(2,i,j,k) = d2f/dx2(x(i),y(j),z(k))
39 C f(3,i,j,k) = d2f/dy2(x(i),y(j),z(k))
40 C f(4,i,j,k) = d2f/dz2(x(i),y(j),z(k))
41 C f(5,i,j,k) = d4f/dx2dy2(x(i),y(j),z(k))
42 C f(6,i,j,k) = d4f/dx2dz2(x(i),y(j),z(k))
43 C f(7,i,j,k) = d4f/dy2dz2(x(i),y(j),z(k))
44 C f(8,i,j,k) = d6f/dx2dy2dz2(x(i),y(j),z(k))
45 C
46 C there is a rather Hermite like interpolation formula to go with
47 C this-- see evtricub.for. Also the bicubic formula is given in
48 C mkbicubw.for; the tricubic formula is precisely analogous.
49 C
50 C boundary condition data
51 C inputs:
52  integer inb1x ! 1st dim of xmin & xmax bc arrays
53  integer inb1y ! 1st dim of ymin & ymax bc arrays
54  integer inb1z ! 1st dim of zmin & zmax bc arrays
55 C
56  integer ibcxmin,ibcxmax ! BC type flag @xmin, xmax
57  integer ibcymin,ibcymax ! BC type flag @ymin, ymax
58  integer ibczmin,ibczmax ! BC type flag @zmin, zmax
59 C
60  real*8 bcxmin(inb1x,nz),bcxmax(inb1x,nz) ! xmin & xmax BC data, ny x nz
61  real*8 bcymin(inb1y,nz),bcymax(inb1y,nz) ! ymin & ymax BC data, nx x nz
62  real*8 bczmin(inb1z,ny),bczmax(inb1z,ny) ! zmin & zmax BC data, nx x ny
63 c
64 c where BC data is not required, dummy scalars may be passed.
65 C the ibc* flags determine whether BC data isneeded.
66 c
67 c BC data: bcxmin & bcxmax: BC vs. y,z @xmin,xmax
68 C bcymin & bcymax: BC vs. x,z @ymin,ymax
69 C bczmin & bczmax: BC vs. x,y @zmin,zmax
70 C
71 c ibcxmin -- indicator for boundary condition at xmin=x(1):
72 c bcxmin(...) -- boundary condition data
73 c =-1 -- use periodic boundary condition
74 c =0 -- use "not a knot"
75 c =1 -- match slope, specified at x(1),y(iy),z(iz) by bcxmin(iy,iz)
76 c =2 -- match 2nd derivative, specified at x(1),y(iy),z(iz)
77 c by bcxmin(iy,iz
78 c =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all y(j)
79 c =4 -- boundary condition is d2f/dx2=0 at x(1), all y(j)
80 c =5 -- match 1st derivative to 1st divided difference
81 c =6 -- match 2nd derivative to 2nd divided difference
82 c =7 -- match 3rd derivative to 3rd divided difference
83 c (for more detailed definition of BCs 5-7, see the
84 c comments of subroutine mkspline)
85 c ***NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
86 c
87 c ibcxmax -- indicator for boundary condition at x(nx):
88 c bcxmax(...) -- boundary condition data
89 c (interpretation as with ibcxmin, bcxmin)
90 c NOTE: if ibcxmin=-1 then the periodic BC applies on both sides
91 c and ibcxmax, bcxmax are ignored.
92 c inb1x -- 1st dimension of bcxmin, bcxmax: if ibcxmin or ibcxmax .gt. 0
93 c this must be .ge. ny.
94 c
95 c interpretation of ibcymin,bcymin,ibcymax,bcymax,inb1y
96 c is same as with ibcxmin,...
97 c
98 c interpretation of ibczmin,bczmin,ibczmax,bczmax,inb1z
99 c is same as with ibcxmin,...
100 c
101 c the explicit bdy condition arrays are referenced only if the
102 c corresponding "ibc" flag values are set to 1 or 2.
103 c
104 c workspace:
105  integer nwk ! size of workspace
106  REAL*8 wk(nwk) ! workspace array
107 c
108 c this version requires a very large workspace, nwk.ge.80*nx*ny*nz
109 c so as to be able to use tcspline to calculate the spline coefficients.
110 c
111 c output:
112  integer ilinx ! x vector equal spacing flag
113  integer iliny ! y vector equal spacing flag
114  integer ilinz ! z vector equal spacing flag
115 c
116 c ilinx -- =1 on output if x(nx) pts are nearly evenly spaced (tol=1e-3)
117 c iliny -- =1 on output if y(ny) evenly spaced (tol=1e-3)
118 c ilinz -- =1 on output if z(nz) evenly spaced (tol=1e-3)
119 c
120  integer ier ! exit code
121 c ier -- completion code, 0 for normal
122 c
123 C-----------------------------------------------------
124 C
125  itest=80*nx*ny*nz
126  if(nwk.lt.itest) then
127  write(6,9901) nwk,itest
128  9901 format(' ?mktricubw: workspace too small:'/
129  > ' user supplied: nwk=',i7,'; need at least: ',i7/
130  > ' nwk = at least 21*nx*ny is required.')
131  ier=1
132  return
133  endif
134 C
135  iadfp=1
136  isiz1=64*nx*ny*nz
137  iadfw=iadfp+isiz1
138  inwk = nwk-isiz1
139 C
140  call r8mktricop(f,nf2,nf3,wk(iadfp),nx,ny,nz)
141 C
142 C evaluate 4x4x4 continuous tricubic spline
143 C
144  call r8tcspline(x,nx,y,ny,z,nz,wk(iadfp),nx,ny,
145  > ibcxmin,bcxmin,ibcxmax,bcxmax,inb1x,
146  > ibcymin,bcymin,ibcymax,bcymax,inb1y,
147  > ibczmin,bczmin,ibczmax,bczmax,inb1z,
148  > wk(iadfw),inwk,ilinx,iliny,ilinz,ier)
149 C
150 C convert to 8-coefficient form
151 C
152  hxlast=x(nx)-x(nx-1)
153  hylast=y(ny)-y(ny-1)
154  hzlast=z(nz)-z(nz-1)
155  call r8mktricon(f,nf2,nf3,wk(iadfp),nx,ny,nz,hxlast,hylast,hzlast)
156 C
157  return
158  end
159 C----------------------------------------------------------------
160 C mktricop -- copy spline function input data
161 C
162  subroutine r8mktricop(fin,nf2,nf3,fwk,nx,ny,nz)
163 C
164 !============
165 ! idecl: explicitize implicit INTEGER declarations:
166  IMPLICIT NONE
167  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
168  INTEGER nf3,nx,ny,nz,nf2,iz,iy,ix
169 !============
170  real*8 fin(8,nf2,nf3,nz)
171  real*8 fwk(4,4,4,nx,ny,nz)
172 C
173  do iz=1,nz
174  do iy=1,ny
175  do ix=1,nx
176  fwk(1,1,1,ix,iy,iz)=fin(1,ix,iy,iz)
177  enddo
178  enddo
179  enddo
180 C
181  return
182  end
183 C----------------------------------------------------------------
184 C mktricon -- create compact spline representation from 4x4
185 C (bcspline) representation
186 C
187  subroutine r8mktricon(fin,nf2,nf3,fwk,nx,ny,nz,
188  > hxlast,hylast,hzlast)
189 C
190 !============
191 ! idecl: explicitize implicit INTEGER declarations:
192  IMPLICIT NONE
193  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
194  INTEGER nf3,nx,ny,nz,nf2,iz,iy,ix,iflag,ixuse,iyuse,izuse,j
195 !============
196 ! idecl: explicitize implicit REAL declarations:
197  real*8 hxlast,hylast,hzlast,dxuse,dyuse,dzuse
198 !============
199  real*8 fin(8,nf2,nf3,nz)
200  real*8 fwk(4,4,4,nx,ny,nz)
201 C-----------------------------------------------------
202 C local arrays
203 C
204  integer iselect(10)
205  real*8 zvalues(10)
206 C
207  data iselect/-1,0,0,0,0,0,0,0,0,0/
208 C
209 C-----------------------------------------------------
210 C
211  do iz=1,nz
212  do iy=1,ny
213  do ix=1,nx
214 C
215 C copy derivatives from result. Special treatment needed for end zones
216 C
217  iflag=0
218  dxuse=0.0_r8
219  dyuse=0.0_r8
220  dzuse=0.0_r8
221  ixuse=ix
222  iyuse=iy
223  izuse=iz
224 C
225  if(ix.eq.nx) then
226  iflag=1
227  dxuse=hxlast
228  ixuse=ix-1
229  endif
230  if(iy.eq.ny) then
231  iflag=1
232  dyuse=hylast
233  iyuse=iy-1
234  endif
235  if(iz.eq.nz) then
236  iflag=1
237  dzuse=hzlast
238  izuse=iz-1
239  endif
240 C
241  if(iflag.eq.1) then
242  call r8tcspevfn(iselect,1,1,zvalues,
243  > ixuse,iyuse,izuse,dxuse,dyuse,dzuse,
244  > fwk,nx,ny,nz)
245  do j=2,8
246  fin(j,ix,iy,iz)=zvalues(j)
247  enddo
248  else
249  fin(2,ix,iy,iz)=2.0_r8*fwk(3,1,1,ix,iy,iz)
250  fin(3,ix,iy,iz)=2.0_r8*fwk(1,3,1,ix,iy,iz)
251  fin(4,ix,iy,iz)=2.0_r8*fwk(1,1,3,ix,iy,iz)
252  fin(5,ix,iy,iz)=4.0_r8*fwk(3,3,1,ix,iy,iz)
253  fin(6,ix,iy,iz)=4.0_r8*fwk(3,1,3,ix,iy,iz)
254  fin(7,ix,iy,iz)=4.0_r8*fwk(1,3,3,ix,iy,iz)
255  fin(8,ix,iy,iz)=8.0_r8*fwk(3,3,3,ix,iy,iz)
256  endif
257 C
258  enddo
259  enddo
260  enddo
261 C
262  return
263  end