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