V3FIT
•All Classes Namespaces Files Functions Variables Enumerations Macros Pages
r8pc3ev.f
1  subroutine r8pc3ev(xget,yget,zget,x,nx,y,ny,z,nz,
2  > ilinx,iliny,ilinz,
3  > f,inf2,inf3,ict,fval,ier)
4 C
5 C evaluate a trilinear interpolant on a rectilinear grid
6 C derivatives are available, but, not continuous across grid planes.
7 C
8 C this subroutine calls two subroutines:
9 C herm3xyz -- find cell containing (xget,yget,zget)
10 C pc3fcn -- evaluate interpolant function and (optionally) derivatives
11 C
12 C input arguments:
13 C ================
14 C
15 !============
16 ! idecl: explicitize implicit INTEGER declarations:
17  IMPLICIT NONE
18  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
19  INTEGER ny,nz,inf2,inf3,nx
20 !============
21  real*8 xget,yget,zget ! target of this interpolation
22  real*8 x(nx) ! ordered x grid
23  real*8 y(ny) ! ordered y grid
24  real*8 z(nz) ! ordered z grid
25  integer ilinx ! ilinx=1 => assume x evenly spaced
26  integer iliny ! iliny=1 => assume y evenly spaced
27  integer ilinz ! ilinz=1 => assume z evenly spaced
28 C
29  real*8 f(inf2,inf3,nz) ! function data
30 C
31 C f 2nd dimension inf2 must be .ge. nx; 3rd dim inf3 .ge. ny
32 C contents of f:
33 C
34 C f(i,j,k) = f @ x(i),y(j),z(k)
35 C
36  integer ict(8) ! code specifying output desired
37 C
38 C ict(1)=1 -- return f (0, don't)
39 C ict(2)=1 -- return df/dx (0, don't)
40 C ict(3)=1 -- return df/dy (0, don't)
41 C ict(4)=1 -- return df/dz (0, don't)
42 C ict(5)=1 -- return d2f/dxdy (0, don't)
43 C ict(6)=1 -- return d2f/dxdz (0, don't)
44 C ict(7)=1 -- return d2f/dydz (0, don't)
45 C ict(8)=1 -- return d3f/dxdydz (0, don't)
46 C
47 C output arguments:
48 C =================
49 C
50  real*8 fval(*) ! output data
51  integer ier ! error code =0 ==> no error
52 C
53 C fval(1) receives the first output (depends on ict(...) spec)
54 C fval(2) receives the second output (depends on ict(...) spec)
55 C fval(3) receives the third output (depends on ict(...) spec)
56 C fval(4) receives the 4th output (depends on ict(...) spec)
57 C fval(5-8) receive 5th thru 8th outputs (if required by ict(...) spec)
58 C
59 C examples:
60 C on input ict = [1,1,1,1,0,0,0,0]
61 C on output fval= [f,df/dx,df/dy,df/dz]
62 C
63 C on input ict = [1,0,0,0,0,0,0,0]
64 C on output fval= [f] ... elements 2-8 never referenced
65 C
66 C on input ict = [0,1,1,0,0,0,0,0]
67 C on output fval= [df/dx,df/dy] ... elements 3-8 never referenced
68 C
69 C on input ict = [0,0,0,0,1,0,0,0]
70 C on output fval= [d2f/dxdy] ... elements 2-8 never referenced.
71 C
72 C ier -- completion code: 0 means OK
73 C-------------------
74 C local:
75 C
76  integer i,j,k ! cell indices
77 C
78 C normalized displacement from (x(i),y(j)) corner of cell.
79 C xparam=0 @x(i) xparam=1 @x(i+1)
80 C yparam=0 @y(j) yparam=1 @y(j+1)
81 C zparam=0 @z(k) zparam=1 @z(k+1)
82 C
83  real*8 xparam,yparam,zparam
84 C
85 C cell dimensions and
86 C inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
87 C
88  real*8 hx,hy,hz
89  real*8 hxi,hyi,hzi
90 C
91 C 0 .le. xparam .le. 1
92 C 0 .le. yparam .le. 1
93 C 0 .le. zparam .le. 1
94 C
95 C---------------------------------------------------------------------
96 C
97  call r8herm3xyz(xget,yget,zget,x,nx,y,ny,z,nz,ilinx,iliny,ilinz,
98  > i,j,k,xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi,ier)
99  if(ier.ne.0) return
100 c
101  call r8pc3fcn(ict,1,1,
102  > fval,i,j,k,xparam,yparam,zparam,
103  > hx,hxi,hy,hyi,hz,hzi,
104  > f,inf2,inf3,nz)
105 C
106  return
107  end
108 C---------------------------------------------------------------------
109 C evaluate trilinear function interpolation -- 3d fcn
110 C --vectorized-- dmc 10 Feb 1999
111 C
112  subroutine r8pc3fcn(ict,ivec,ivecd,
113  > fval,ii,jj,kk,xparam,yparam,zparam,
114  > hx,hxi,hy,hyi,hz,hzi,
115  > fin,inf2,inf3,nz)
116 C
117 !============
118 ! idecl: explicitize implicit INTEGER declarations:
119  IMPLICIT NONE
120  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
121  INTEGER inf3,nz,inf2,i,j,k,iadr
122 !============
123 ! idecl: explicitize implicit REAL declarations:
124  real*8 xp,xpi,yp,ypi,zp,zpi
125 !============
126  integer ict(8) ! requested output control
127  integer ivec ! vector length
128  integer ivecd ! vector dimension (1st dim of fval)
129 C
130  integer ii(ivec),jj(ivec),kk(ivec) ! target cells (i,j,k)
131  real*8 xparam(ivec),yparam(ivec),zparam(ivec)
132  ! normalized displacements from (i,j,k) corners
133 C
134  real*8 hx(ivec),hy(ivec),hz(ivec) ! grid spacing, and
135  real*8 hxi(ivec),hyi(ivec),hzi(ivec) ! inverse grid spacing
136  ! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
137 C
138  real*8 fin(inf2,inf3,nz) ! interpolant data (cf "pc3ev")
139 C
140  real*8 fval(ivecd,*) ! output returned
141 C
142 C for detailed description of fin, ict and fval see subroutine pc3ev
143 C comments. Note ict is not vectorized; the same output
144 C is expected to be returned for all input vector data points.
145 C
146 C note that the index inputs ii,jj,kk and parameter inputs
147 C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and the
148 C output array fval has a vector ** 1st dimension ** whose
149 C size must be given as a separate argument
150 C
151 C to use this routine in scalar mode, pass in ivec=ivecd=1
152 C
153 C---------------
154 C
155  real*8 sum
156  integer v
157 C
158  do v=1,ivec
159  i=ii(v)
160  j=jj(v)
161  k=kk(v)
162 C
163 C ...in x direction
164 C
165  xp=xparam(v)
166  xpi=1.0_r8-xp
167 C
168 C ...in y direction
169 C
170  yp=yparam(v)
171  ypi=1.0_r8-yp
172 C
173 C ...in z direction
174 C
175  zp=zparam(v)
176  zpi=1.0_r8-zp
177 C
178  iadr=0
179 C
180 C get desired values:
181 C
182  if(ict(1).eq.1) then
183 C
184 C function value:
185 C
186  iadr=iadr+1
187  sum=zpi*(
188  > xpi*(ypi*fin(i,j,k) +yp*fin(i,j+1,k))+
189  > xp*(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
190  > +zp*(
191  > xpi*(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))+
192  > xp*(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
193 C
194  fval(v,iadr)=sum
195  endif
196 C
197  if(ict(2).eq.1) then
198 C
199 C df/dx:
200 C
201  iadr=iadr+1
202  sum=zpi*(
203  > -(ypi*fin(i,j,k) +yp*fin(i,j+1,k))
204  > +(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
205  > +zp*(
206  > -(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))
207  > +(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
208 C
209  fval(v,iadr)=sum*hxi(v)
210  endif
211 C
212  if(ict(3).eq.1) then
213 C
214 C df/dy:
215 C
216  iadr=iadr+1
217  sum=zpi*(
218  > xpi*(-fin(i,j,k) +fin(i,j+1,k))+
219  > xp*(-fin(i+1,j,k)+fin(i+1,j+1,k)))
220  > +zp*(
221  > xpi*(-fin(i,j,k+1) +fin(i,j+1,k+1))+
222  > xp*(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
223 C
224  fval(v,iadr)=sum*hyi(v)
225  endif
226 C
227  if(ict(4).eq.1) then
228 C
229 C df/dz:
230 C
231  iadr=iadr+1
232  sum= -(
233  > xpi*(ypi*fin(i,j,k) +yp*fin(i,j+1,k))+
234  > xp*(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
235  > +(
236  > xpi*(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))+
237  > xp*(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
238 C
239  fval(v,iadr)=sum*hzi(v)
240  endif
241 C
242  if(ict(5).eq.1) then
243 C
244 C d2f/dxdy:
245 C
246  iadr=iadr+1
247  sum=zpi*(
248  > -(-fin(i,j,k) +fin(i,j+1,k))
249  > +(-fin(i+1,j,k)+fin(i+1,j+1,k)))
250  > +zp*(
251  > -(-fin(i,j,k+1) +fin(i,j+1,k+1))
252  > +(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
253 C
254  fval(v,iadr)=sum*hxi(v)*hyi(v)
255  endif
256 C
257  if(ict(6).eq.1) then
258 C
259 C d2f/dxdz:
260 C
261  iadr=iadr+1
262  sum= -(
263  > -(ypi*fin(i,j,k) +yp*fin(i,j+1,k))
264  > +(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
265  > +(
266  > -(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))
267  > +(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
268 C
269  fval(v,iadr)=sum*hxi(v)*hzi(v)
270  endif
271 C
272  if(ict(7).eq.1) then
273 C
274 C d2f/dydz:
275 C
276  iadr=iadr+1
277  sum= -(
278  > xpi*(-fin(i,j,k) +fin(i,j+1,k))+
279  > xp*(-fin(i+1,j,k)+fin(i+1,j+1,k)))
280  > +(
281  > xpi*(-fin(i,j,k+1) +fin(i,j+1,k+1))+
282  > xp*(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
283 C
284  fval(v,iadr)=sum*hyi(v)*hzi(v)
285  endif
286 C
287  if(ict(8).eq.1) then
288 C
289 C d3f/dxdydz:
290 C
291  iadr=iadr+1
292  sum= -(
293  > -(-fin(i,j,k) +fin(i,j+1,k))
294  > +(-fin(i+1,j,k)+fin(i+1,j+1,k)))
295  > +(
296  > -(-fin(i,j,k+1) +fin(i,j+1,k+1))
297  > +(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
298 C
299  fval(v,iadr)=sum*hxi(v)*hyi(v)*hzi(v)
300  endif
301 C
302  enddo ! vector loop
303 C
304  return
305  end
306 C---------------------------------------------------------------------
307 C evaluate trilinear function interpolation -- 3d fcn
308 C --vectorized-- dmc 10 Feb 1999
309 C --optimized for VARIATION along x axis ONLY--
310 C
311  subroutine r8pc3fcnx(ict,ivec,ivecd,
312  > fval,ii,jj,kk,xparam,yparam,zparam,
313  > hx,hxi,hy,hyi,hz,hzi,
314  > fin,inf2,inf3,nz)
315 C
316 !============
317 ! idecl: explicitize implicit INTEGER declarations:
318  IMPLICIT NONE
319  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
320  INTEGER inf3,nz,inf2,j,k,i,iadr
321 !============
322 ! idecl: explicitize implicit REAL declarations:
323  REAL*8 yp,ypi,zp,zpi,xp,xpi
324 !============
325  integer ict(8) ! requested output control
326  integer ivec ! vector length
327  integer ivecd ! vector dimension (1st dim of fval)
328 C
329  integer ii(ivec),jj,kk ! target cells (i,j,k)
330  REAL*8 xparam(ivec),yparam,zparam
331  ! normalized displacements from (i,j,k) corners
332 C
333  real*8 hx(ivec),hy,hz ! grid spacing, and
334  real*8 hxi(ivec),hyi,hzi ! inverse grid spacing
335  ! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
336 C
337  real*8 fin(inf2,inf3,nz) ! interpolant data (cf "pc3ev")
338 C
339  real*8 fval(ivecd,*) ! output returned
340 C
341 C for detailed description of fin, ict and fval see subroutine pc3ev
342 C comments. Note ict is not vectorized; the same output
343 C is expected to be returned for all input vector data points.
344 C
345 C note that the index inputs ii,jj,kk and parameter inputs
346 C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and the
347 C output array fval has a vector ** 1st dimension ** whose
348 C size must be given as a separate argument
349 C
350 C to use this routine in scalar mode, pass in ivec=ivecd=1
351 C
352 C---------------
353 C
354  real*8 sum
355  integer v
356 C
357  j=jj
358  k=kk
359 C
360 C ...in y direction
361 C
362  yp=yparam
363  ypi=1.0_r8-yp
364 C
365 C ...in z direction
366 C
367  zp=zparam
368  zpi=1.0_r8-zp
369 
370  do v=1,ivec
371  i=ii(v)
372 C
373 C ...in x direction
374 C
375  xp=xparam(v)
376  xpi=1.0_r8-xp
377 C
378  iadr=0
379 C
380 C get desired values:
381 C
382  if(ict(1).eq.1) then
383 C
384 C function value:
385 C
386  iadr=iadr+1
387  sum=zpi*(
388  > xpi*(ypi*fin(i,j,k) +yp*fin(i,j+1,k))+
389  > xp*(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
390  > +zp*(
391  > xpi*(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))+
392  > xp*(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
393 C
394  fval(v,iadr)=sum
395  endif
396 C
397  if(ict(2).eq.1) then
398 C
399 C df/dx:
400 C
401  iadr=iadr+1
402  sum=zpi*(
403  > -(ypi*fin(i,j,k) +yp*fin(i,j+1,k))
404  > +(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
405  > +zp*(
406  > -(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))
407  > +(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
408 C
409  fval(v,iadr)=sum*hxi(v)
410  endif
411 C
412  if(ict(3).eq.1) then
413 C
414 C df/dy:
415 C
416  iadr=iadr+1
417  sum=zpi*(
418  > xpi*(-fin(i,j,k) +fin(i,j+1,k))+
419  > xp*(-fin(i+1,j,k)+fin(i+1,j+1,k)))
420  > +zp*(
421  > xpi*(-fin(i,j,k+1) +fin(i,j+1,k+1))+
422  > xp*(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
423 C
424  fval(v,iadr)=sum*hyi
425  endif
426 C
427  if(ict(4).eq.1) then
428 C
429 C df/dz:
430 C
431  iadr=iadr+1
432  sum= -(
433  > xpi*(ypi*fin(i,j,k) +yp*fin(i,j+1,k))+
434  > xp*(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
435  > +(
436  > xpi*(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))+
437  > xp*(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
438 C
439  fval(v,iadr)=sum*hzi
440  endif
441 C
442  if(ict(5).eq.1) then
443 C
444 C d2f/dxdy:
445 C
446  iadr=iadr+1
447  sum=zpi*(
448  > -(-fin(i,j,k) +fin(i,j+1,k))
449  > +(-fin(i+1,j,k)+fin(i+1,j+1,k)))
450  > +zp*(
451  > -(-fin(i,j,k+1) +fin(i,j+1,k+1))
452  > +(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
453 C
454  fval(v,iadr)=sum*hxi(v)*hyi
455  endif
456 C
457  if(ict(6).eq.1) then
458 C
459 C d2f/dxdz:
460 C
461  iadr=iadr+1
462  sum= -(
463  > -(ypi*fin(i,j,k) +yp*fin(i,j+1,k))
464  > +(ypi*fin(i+1,j,k)+yp*fin(i+1,j+1,k)))
465  > +(
466  > -(ypi*fin(i,j,k+1) +yp*fin(i,j+1,k+1))
467  > +(ypi*fin(i+1,j,k+1)+yp*fin(i+1,j+1,k+1)))
468 C
469  fval(v,iadr)=sum*hxi(v)*hzi
470  endif
471 C
472  if(ict(7).eq.1) then
473 C
474 C d2f/dydz:
475 C
476  iadr=iadr+1
477  sum= -(
478  > xpi*(-fin(i,j,k) +fin(i,j+1,k))+
479  > xp*(-fin(i+1,j,k)+fin(i+1,j+1,k)))
480  > +(
481  > xpi*(-fin(i,j,k+1) +fin(i,j+1,k+1))+
482  > xp*(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
483 C
484  fval(v,iadr)=sum*hyi*hzi
485  endif
486 C
487  if(ict(8).eq.1) then
488 C
489 C d3f/dxdydz:
490 C
491  iadr=iadr+1
492  sum= -(
493  > -(-fin(i,j,k) +fin(i,j+1,k))
494  > +(-fin(i+1,j,k)+fin(i+1,j+1,k)))
495  > +(
496  > -(-fin(i,j,k+1) +fin(i,j+1,k+1))
497  > +(-fin(i+1,j,k+1)+fin(i+1,j+1,k+1)))
498 C
499  fval(v,iadr)=sum*hxi(v)*hyi*hzi
500  endif
501 C
502  enddo ! vector loop
503 C
504  return
505  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17