V3FIT
herm3ev.f
1  subroutine herm3ev(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 3d cubic Hermite interpolant on a rectilinear
6 C grid -- this is C1 in all directions.
7 C
8 C this subroutine calls two subroutines:
9 C herm3xyz -- find cell containing (xget,yget,zget)
10 C herm3fcn -- evaluate interpolant function and (optionally) derivatives
11 C
12 C input arguments:
13 C ================
14 C
15  real xget,yget,zget ! target of this interpolation
16  real x(nx) ! ordered x grid
17  real y(ny) ! ordered y grid
18  real z(nz) ! ordered z grid
19  integer ilinx ! ilinx=1 => assume x evenly spaced
20  integer iliny ! iliny=1 => assume y evenly spaced
21  integer ilinz ! ilinz=1 => assume z evenly spaced
22 C
23  real f(0:7,inf2,inf3,nz) ! function data
24 C
25 C f 2nd dimension inf2 must be .ge. nx; 3rd dim inf3 .ge. ny
26 C contents of f:
27 C
28 C f(0,i,j,k) = f @ x(i),y(j),z(k)
29 C f(1,i,j,k) = df/dx @ x(i),y(j),z(k)
30 C f(2,i,j,k) = df/dy @ x(i),y(j),z(k)
31 C f(3,i,j,k) = df/dz @ x(i),y(j),z(k)
32 C f(4,i,j,k) = d2f/dxdy @ x(i),y(j),z(k)
33 C f(5,i,j,k) = d2f/dxdz @ x(i),y(j),z(k)
34 C f(6,i,j,k) = d2f/dydz @ x(i),y(j),z(k)
35 C f(7,i,j,k) = d3f/dxdydz @ x(i),y(j),z(k)
36 C
37  integer ict(8) ! code specifying output desired
38 C
39 C ict(1)=1 -- return f (0, don't)
40 C ict(2)=1 -- return df/dx (0, don't)
41 C ict(3)=1 -- return df/dy (0, don't)
42 C ict(4)=1 -- return df/dz (0, don't)
43 C ict(5)=1 -- return d2f/dxdy (0, don't)
44 C ict(6)=1 -- return d2f/dxdz (0, don't)
45 C ict(7)=1 -- return d2f/dydz (0, don't)
46 C ict(8)=1 -- return d3f/dxdydz (0, don't)
47 C
48 C output arguments:
49 C =================
50 C
51  real fval(*) ! output data
52  integer ier ! error code =0 ==> no error
53 C
54 C fval(1) receives the first output (depends on ict(...) spec)
55 C fval(2) receives the second output (depends on ict(...) spec)
56 C fval(3) receives the third output (depends on ict(...) spec)
57 C fval(4) receives the 4th output (depends on ict(...) spec)
58 C fval(5-8) receive 5th thru 8th outputs (if required by ict(...) spec)
59 C
60 C examples:
61 C on input ict = [1,1,1,1,0,0,0,0]
62 C on output fval= [f,df/dx,df/dy,df/dz]
63 C
64 C on input ict = [1,0,0,0,0,0,0,0]
65 C on output fval= [f] ... elements 2-8 never referenced
66 C
67 C on input ict = [0,1,1,0,0,0,0,0]
68 C on output fval= [df/dx,df/dy] ... elements 3-8 never referenced
69 C
70 C on input ict = [0,0,0,0,1,0,0,0]
71 C on output fval= [d2f/dxdy] ... elements 2-8 never referenced.
72 C
73 C ier -- completion code: 0 means OK
74 C-------------------
75 C local:
76 C
77  integer i,j,k ! cell indices
78 C
79 C normalized displacement from (x(i),y(j)) corner of cell.
80 C xparam=0 @x(i) xparam=1 @x(i+1)
81 C yparam=0 @y(j) yparam=1 @y(j+1)
82 C zparam=0 @z(k) zparam=1 @z(k+1)
83 C
84  real xparam,yparam,zparam
85 C
86 C cell dimensions and
87 C inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
88 C
89  real hx,hy,hz
90  real hxi,hyi,hzi
91 C
92 C 0 .le. xparam .le. 1
93 C 0 .le. yparam .le. 1
94 C 0 .le. zparam .le. 1
95 C
96 C---------------------------------------------------------------------
97 C
98  call herm3xyz(xget,yget,zget,x,nx,y,ny,z,nz,ilinx,iliny,ilinz,
99  > i,j,k,xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi,ier)
100  if(ier.ne.0) return
101 c
102  call herm3fcn(ict,1,1,
103  > fval,i,j,k,xparam,yparam,zparam,
104  > hx,hxi,hy,hyi,hz,hzi,
105  > f,inf2,inf3,nz)
106 C
107  return
108  end
109 C---------------------------------------------------------------------
110 c herm3xyz -- look up x-y-z zone
111 c
112 c this is the "first part" of herm3ev, see comments, above.
113 c
114  subroutine herm3xyz(xget,yget,zget,x,nx,y,ny,z,nz,
115  > ilinx,iliny,ilinz,
116  > i,j,k,xparam,yparam,zparam,
117  > hx,hxi,hy,hyi,hz,hzi,ier)
118 c
119 c input of herm3xyz
120 c ================
121 c
122  integer nx,ny,nz ! coord. grid dimensions
123 c
124  real xget,yget,zget ! target point
125  real x(nx),y(ny),z(nz) ! indep. coords. (ascending order)
126 c
127  integer ilinx ! =1: x evenly spaced
128  integer iliny ! =1: y evenly spaced
129  integer ilinz ! =1: z evenly spaced
130 c
131 c output of herm3xyz
132 c =================
133  integer i,j,k ! index to cell containing target pt
134 c on exit: 1.le.i.le.nx-1 1.le.j.le.ny-1 1.le.k.le.nz-1
135 c
136 c normalized position w/in (i,j) cell (btw 0 and 1):
137 c
138  real xparam ! (xget-x(i))/(x(i+1)-x(i))
139  real yparam ! (yget-y(j))/(y(j+1)-y(j))
140  real zparam ! (zget-z(k))/(z(k+1)-z(k))
141 c
142 c grid spacing
143 c
144  real hx ! hx = x(i+1)-x(i)
145  real hy ! hy = y(j+1)-y(j)
146  real hz ! hz = z(k+1)-z(k)
147 c
148 c inverse grid spacing:
149 c
150  real hxi ! 1/hx = 1/(x(i+1)-x(i))
151  real hyi ! 1/hy = 1/(y(j+1)-y(j))
152  real hzi ! 1/hz = 1/(z(k+1)-z(k))
153 c
154  integer ier ! return ier.ne.0 on error
155 c
156 c------------------------------------
157 c
158  ier=0
159 c
160 c range check
161 c
162  zxget=xget
163  zyget=yget
164  zzget=zget
165  if((xget.lt.x(1)).or.(xget.gt.x(nx))) then
166  zxtol=4.0e-7*max(abs(x(1)),abs(x(nx)))
167  if((xget.lt.x(1)-zxtol).or.(xget.gt.x(nx)+zxtol)) then
168  ier=1
169  write(6,1001) xget,x(1),x(nx)
170  1001 format(' ?herm3ev: xget=',1pe11.4,' out of range ',
171  > 1pe11.4,' to ',1pe11.4)
172  else
173  if((xget.lt.x(1)-0.5*zxtol).or.
174  > (xget.gt.x(nx)+0.5*zxtol))
175  > write(6,1011) xget,x(1),x(nx)
176  1011 format(' %herm3ev: xget=',1pe15.8,' beyond range ',
177  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
178  if(xget.lt.x(1)) then
179  zxget=x(1)
180  else
181  zxget=x(nx)
182  endif
183  endif
184  endif
185  if((yget.lt.y(1)).or.(yget.gt.y(ny))) then
186  zytol=4.0e-7*max(abs(y(1)),abs(y(ny)))
187  if((yget.lt.y(1)-zytol).or.(yget.gt.y(ny)+zytol)) then
188  ier=1
189  write(6,1002) yget,y(1),y(ny)
190  1002 format(' ?herm3ev: yget=',1pe11.4,' out of range ',
191  > 1pe11.4,' to ',1pe11.4)
192  else
193  if((yget.lt.y(1)-0.5*zytol).or.
194  > (yget.gt.y(ny)+0.5*zytol))
195  > write(6,1012) yget,y(1),y(ny)
196  1012 format(' %herm3ev: yget=',1pe15.8,' beyond range ',
197  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
198  if(yget.lt.y(1)) then
199  zyget=y(1)
200  else
201  zyget=y(ny)
202  endif
203  endif
204  endif
205  if((zget.lt.z(1)).or.(zget.gt.z(nz))) then
206  zztol=4.0e-7*max(abs(z(1)),abs(z(nz)))
207  if((zget.lt.z(1)-zztol).or.(zget.gt.z(nz)+zztol)) then
208  ier=1
209  write(6,1003) zget,z(1),z(nz)
210  1003 format(' ?herm3ev: zget=',1pe11.4,' out of range ',
211  > 1pe11.4,' to ',1pe11.4)
212  else
213  if((zget.lt.z(1)-0.5*zztol).or.
214  > (zget.gt.z(nz)+0.5*zztol))
215  > write(6,1013) zget,z(1),z(nz)
216  1013 format(' %herm3ev: zget=',1pe15.8,' beyond range ',
217  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
218  if(zget.lt.z(1)) then
219  zzget=z(1)
220  else
221  zzget=z(nz)
222  endif
223  endif
224  endif
225  if(ier.ne.0) return
226 c
227 c now find interval in which target point lies..
228 c
229  nxm=nx-1
230  nym=ny-1
231  nzm=nz-1
232 c
233  if(ilinx.eq.1) then
234  ii=1+nxm*(zxget-x(1))/(x(nx)-x(1))
235  i=min(nxm, ii)
236  if(zxget.lt.x(i)) then
237  i=i-1
238  else if(zxget.gt.x(i+1)) then
239  i=i+1
240  endif
241  else
242  if((1.le.i).and.(i.lt.nxm)) then
243  if((x(i).le.zxget).and.(zxget.le.x(i+1))) then
244  continue ! already have the zone
245  else
246  call zonfind(x,nx,zxget,i)
247  endif
248  else
249  i=nx/2
250  call zonfind(x,nx,zxget,i)
251  endif
252  endif
253 c
254  if(iliny.eq.1) then
255  jj=1+nym*(zyget-y(1))/(y(ny)-y(1))
256  j=min(nym, jj)
257  if(zyget.lt.y(j)) then
258  j=j-1
259  else if(zyget.gt.y(j+1)) then
260  j=j+1
261  endif
262  else
263  if((1.le.j).and.(j.lt.nym)) then
264  if((y(j).le.zyget).and.(zyget.le.y(j+1))) then
265  continue ! already have the zone
266  else
267  call zonfind(y,ny,zyget,j)
268  endif
269  else
270  j=ny/2
271  call zonfind(y,ny,zyget,j)
272  endif
273  endif
274 c
275  if(ilinz.eq.1) then
276  kk=1+nzm*(zzget-z(1))/(z(nz)-z(1))
277  k=min(nzm, kk)
278  if(zzget.lt.z(k)) then
279  k=k-1
280  else if(zzget.gt.z(k+1)) then
281  k=k+1
282  endif
283  else
284  if((1.le.k).and.(k.lt.nzm)) then
285  if((z(k).le.zzget).and.(zzget.le.z(k+1))) then
286  continue ! already have the zone
287  else
288  call zonfind(z,nz,zzget,k)
289  endif
290  else
291  k=nz/2
292  call zonfind(z,nz,zzget,k)
293  endif
294  endif
295 c
296  hx=(x(i+1)-x(i))
297  hy=(y(j+1)-y(j))
298  hz=(z(k+1)-z(k))
299 c
300  hxi=1.0/hx
301  hyi=1.0/hy
302  hzi=1.0/hz
303 c
304  xparam=(zxget-x(i))*hxi
305  yparam=(zyget-y(j))*hyi
306  zparam=(zzget-z(k))*hzi
307 c
308  return
309  end
310 C---------------------------------------------------------------------
311 C evaluate C1 cubic Hermite function interpolation -- 3d fcn
312 C --vectorized-- dmc 10 Feb 1999
313 C
314  subroutine herm3fcn(ict,ivec,ivecd,
315  > fval,ii,jj,kk,xparam,yparam,zparam,
316  > hx,hxi,hy,hyi,hz,hzi,
317  > fin,inf2,inf3,nz)
318 C
319  integer ict(8) ! requested output control
320  integer ivec ! vector length
321  integer ivecd ! vector dimension (1st dim of fval)
322 C
323  integer ii(ivec),jj(ivec),kk(ivec) ! target cells (i,j,k)
324  real xparam(ivec),yparam(ivec),zparam(ivec)
325  ! normalized displacements from (i,j,k) corners
326 C
327  real hx(ivec),hy(ivec),hz(ivec) ! grid spacing, and
328  real hxi(ivec),hyi(ivec),hzi(ivec) ! inverse grid spacing
329  ! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
330 C
331  real fin(0:7,inf2,inf3,nz) ! interpolant data (cf "herm3ev")
332 C
333  real fval(ivecd,*) ! output returned
334 C
335 C for detailed description of fin, ict and fval see subroutine herm3ev
336 C comments. Note ict is not vectorized; the same output
337 C is expected to be returned for all input vector data points.
338 C
339 C note that the index inputs ii,jj,kk and parameter inputs
340 C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and the
341 C output array fval has a vector ** 1st dimension ** whose
342 C size must be given as a separate argument
343 C
344 C to use this routine in scalar mode, pass in ivec=ivecd=1
345 C
346 C---------------
347 C Hermite cubic basis functions
348 C -->for function value matching
349 C a(0)=0 a(1)=1 a'(0)=0 a'(1)=0
350 C abar(0)=1 abar(1)=0 abar'(0)=0 abar'(1)=0
351 C
352 C a(x)=-2*x**3 + 3*x**2 = x*x*(-2.0*x+3.0)
353 C abar(x)=1-a(x)
354 C a'(x)=-abar'(x) = 6.0*x*(1.0-x)
355 C
356 C -->for derivative matching
357 C b(0)=0 b(1)=0 b'(0)=0 b'(1)=1
358 C bbar(0)=0 bbar(1)=0 bbar'(0)=1 bbar'(1)=0
359 C
360 C b(x)=x**3-x**2 b'(x)=3*x**2-2*x
361 C bbar(x)=x**3-2*x**2+x bbar'(x)=3*x**2-4*x+1
362 C
363  real sum
364  integer v
365 C
366  do v=1,ivec
367  i=ii(v)
368  j=jj(v)
369  k=kk(v)
370 C
371 C ...in x direction
372 C
373  xp=xparam(v)
374  xpi=1.0-xp
375  xp2=xp*xp
376  xpi2=xpi*xpi
377  ax=xp2*(3.0-2.0*xp)
378  axbar=1.0-ax
379  bx=-xp2*xpi
380  bxbar=xpi2*xp
381 C
382 C ...in y direction
383 C
384  yp=yparam(v)
385  ypi=1.0-yp
386  yp2=yp*yp
387  ypi2=ypi*ypi
388  ay=yp2*(3.0-2.0*yp)
389  aybar=1.0-ay
390  by=-yp2*ypi
391  bybar=ypi2*yp
392 C
393 C ...in z direction
394 C
395  zp=zparam(v)
396  zpi=1.0-zp
397  zp2=zp*zp
398  zpi2=zpi*zpi
399  az=zp2*(3.0-2.0*zp)
400  azbar=1.0-az
401  bz=-zp2*zpi
402  bzbar=zpi2*zp
403 C
404  iadr=0
405 C
406 C derivatives:
407 C
408  axp=6.0*xp*xpi
409  axbarp=-axp
410  bxp=xp*(3.0*xp-2.0)
411  bxbarp=xpi*(3.0*xpi-2.0)
412 C
413  ayp=6.0*yp*ypi
414  aybarp=-ayp
415  byp=yp*(3.0*yp-2.0)
416  bybarp=ypi*(3.0*ypi-2.0)
417 C
418  azp=6.0*zp*zpi
419  azbarp=-azp
420  bzp=zp*(3.0*zp-2.0)
421  bzbarp=zpi*(3.0*zpi-2.0)
422 C
423 C get desired values:
424 C
425  if(ict(1).eq.1) then
426 C
427 C function value:
428 C
429  iadr=iadr+1
430  sum=azbar*(
431  > axbar*(aybar*fin(0,i,j,k) +ay*fin(0,i,j+1,k))+
432  > ax*(aybar*fin(0,i+1,j,k)+ay*fin(0,i+1,j+1,k)))
433  > + az*(
434  > axbar*(aybar*fin(0,i,j,k+1) +ay*fin(0,i,j+1,k+1))+
435  > ax*(aybar*fin(0,i+1,j,k+1)+ay*fin(0,i+1,j+1,k+1)))
436 C
437  sum=sum+hx(v)*(
438  > azbar*(
439  > bxbar*(aybar*fin(1,i,j,k) +ay*fin(1,i,j+1,k))+
440  > bx*(aybar*fin(1,i+1,j,k)+ay*fin(1,i+1,j+1,k)))
441  > + az*(
442  > bxbar*(aybar*fin(1,i,j,k+1) +ay*fin(1,i,j+1,k+1))+
443  > bx*(aybar*fin(1,i+1,j,k+1)+ay*fin(1,i+1,j+1,k+1)))
444  > )
445 C
446  sum=sum+hy(v)*(
447  > azbar*(
448  > axbar*(bybar*fin(2,i,j,k) +by*fin(2,i,j+1,k))+
449  > ax*(bybar*fin(2,i+1,j,k)+by*fin(2,i+1,j+1,k)))
450  > + az*(
451  > axbar*(bybar*fin(2,i,j,k+1) +by*fin(2,i,j+1,k+1))+
452  > ax*(bybar*fin(2,i+1,j,k+1)+by*fin(2,i+1,j+1,k+1)))
453  > )
454 C
455  sum=sum+hz(v)*(
456  > bzbar*(
457  > axbar*(aybar*fin(3,i,j,k) +ay*fin(3,i,j+1,k))+
458  > ax*(aybar*fin(3,i+1,j,k)+ay*fin(3,i+1,j+1,k)))
459  > + bz*(
460  > axbar*(aybar*fin(3,i,j,k+1) +ay*fin(3,i,j+1,k+1))+
461  > ax*(aybar*fin(3,i+1,j,k+1)+ay*fin(3,i+1,j+1,k+1)))
462  > )
463 C
464  sum=sum+hx(v)*hy(v)*(
465  > azbar*(
466  > bxbar*(bybar*fin(4,i,j,k) +by*fin(4,i,j+1,k))+
467  > bx*(bybar*fin(4,i+1,j,k)+by*fin(4,i+1,j+1,k)))
468  > + az*(
469  > bxbar*(bybar*fin(4,i,j,k+1) +by*fin(4,i,j+1,k+1))+
470  > bx*(bybar*fin(4,i+1,j,k+1)+by*fin(4,i+1,j+1,k+1)))
471  > )
472 C
473  sum=sum+hx(v)*hz(v)*(
474  > bzbar*(
475  > bxbar*(aybar*fin(5,i,j,k) +ay*fin(5,i,j+1,k))+
476  > bx*(aybar*fin(5,i+1,j,k)+ay*fin(5,i+1,j+1,k)))
477  > + bz*(
478  > bxbar*(aybar*fin(5,i,j,k+1) +ay*fin(5,i,j+1,k+1))+
479  > bx*(aybar*fin(5,i+1,j,k+1)+ay*fin(5,i+1,j+1,k+1)))
480  > )
481 C
482  sum=sum+hy(v)*hz(v)*(
483  > bzbar*(
484  > axbar*(bybar*fin(6,i,j,k) +by*fin(6,i,j+1,k))+
485  > ax*(bybar*fin(6,i+1,j,k)+by*fin(6,i+1,j+1,k)))
486  > + bz*(
487  > axbar*(bybar*fin(6,i,j,k+1) +by*fin(6,i,j+1,k+1))+
488  > ax*(bybar*fin(6,i+1,j,k+1)+by*fin(6,i+1,j+1,k+1)))
489  > )
490 C
491  sum=sum+hx(v)*hy(v)*hz(v)*(
492  > bzbar*(
493  > bxbar*(bybar*fin(7,i,j,k) +by*fin(7,i,j+1,k))+
494  > bx*(bybar*fin(7,i+1,j,k)+by*fin(7,i+1,j+1,k)))
495  > + bz*(
496  > bxbar*(bybar*fin(7,i,j,k+1) +by*fin(7,i,j+1,k+1))+
497  > bx*(bybar*fin(7,i+1,j,k+1)+by*fin(7,i+1,j+1,k+1)))
498  > )
499 C
500  fval(v,iadr)=sum
501  endif
502 C
503  if(ict(2).eq.1) then
504 C
505 C df/dx:
506 C
507  iadr=iadr+1
508 C
509  sum=hxi(v)*(
510  > azbar*(
511  > axbarp*(aybar*fin(0,i,j,k) +ay*fin(0,i,j+1,k))+
512  > axp*(aybar*fin(0,i+1,j,k)+ay*fin(0,i+1,j+1,k)))
513  > + az*(
514  > axbarp*(aybar*fin(0,i,j,k+1) +ay*fin(0,i,j+1,k+1))+
515  > axp*(aybar*fin(0,i+1,j,k+1)+ay*fin(0,i+1,j+1,k+1)))
516  > )
517 C
518  sum=sum+(
519  > azbar*(
520  > bxbarp*(aybar*fin(1,i,j,k) +ay*fin(1,i,j+1,k))+
521  > bxp*(aybar*fin(1,i+1,j,k)+ay*fin(1,i+1,j+1,k)))
522  > + az*(
523  > bxbarp*(aybar*fin(1,i,j,k+1) +ay*fin(1,i,j+1,k+1))+
524  > bxp*(aybar*fin(1,i+1,j,k+1)+ay*fin(1,i+1,j+1,k+1)))
525  > )
526 C
527  sum=sum+hxi(v)*hy(v)*(
528  > azbar*(
529  > axbarp*(bybar*fin(2,i,j,k) +by*fin(2,i,j+1,k))+
530  > axp*(bybar*fin(2,i+1,j,k)+by*fin(2,i+1,j+1,k)))
531  > + az*(
532  > axbarp*(bybar*fin(2,i,j,k+1) +by*fin(2,i,j+1,k+1))+
533  > axp*(bybar*fin(2,i+1,j,k+1)+by*fin(2,i+1,j+1,k+1)))
534  > )
535 C
536  sum=sum+hxi(v)*hz(v)*(
537  > bzbar*(
538  > axbarp*(aybar*fin(3,i,j,k) +ay*fin(3,i,j+1,k))+
539  > axp*(aybar*fin(3,i+1,j,k)+ay*fin(3,i+1,j+1,k)))
540  > + bz*(
541  > axbarp*(aybar*fin(3,i,j,k+1) +ay*fin(3,i,j+1,k+1))+
542  > axp*(aybar*fin(3,i+1,j,k+1)+ay*fin(3,i+1,j+1,k+1)))
543  > )
544 C
545  sum=sum+hy(v)*(
546  > azbar*(
547  > bxbarp*(bybar*fin(4,i,j,k) +by*fin(4,i,j+1,k))+
548  > bxp*(bybar*fin(4,i+1,j,k)+by*fin(4,i+1,j+1,k)))
549  > + az*(
550  > bxbarp*(bybar*fin(4,i,j,k+1) +by*fin(4,i,j+1,k+1))+
551  > bxp*(bybar*fin(4,i+1,j,k+1)+by*fin(4,i+1,j+1,k+1)))
552  > )
553 C
554  sum=sum+hz(v)*(
555  > bzbar*(
556  > bxbarp*(aybar*fin(5,i,j,k) +ay*fin(5,i,j+1,k))+
557  > bxp*(aybar*fin(5,i+1,j,k)+ay*fin(5,i+1,j+1,k)))
558  > + bz*(
559  > bxbarp*(aybar*fin(5,i,j,k+1) +ay*fin(5,i,j+1,k+1))+
560  > bxp*(aybar*fin(5,i+1,j,k+1)+ay*fin(5,i+1,j+1,k+1)))
561  > )
562 C
563  sum=sum+hxi(v)*hy(v)*hz(v)*(
564  > bzbar*(
565  > axbarp*(bybar*fin(6,i,j,k) +by*fin(6,i,j+1,k))+
566  > axp*(bybar*fin(6,i+1,j,k)+by*fin(6,i+1,j+1,k)))
567  > + bz*(
568  > axbarp*(bybar*fin(6,i,j,k+1) +by*fin(6,i,j+1,k+1))+
569  > axp*(bybar*fin(6,i+1,j,k+1)+by*fin(6,i+1,j+1,k+1)))
570  > )
571 C
572  sum=sum+hy(v)*hz(v)*(
573  > bzbar*(
574  > bxbarp*(bybar*fin(7,i,j,k) +by*fin(7,i,j+1,k))+
575  > bxp*(bybar*fin(7,i+1,j,k)+by*fin(7,i+1,j+1,k)))
576  > + bz*(
577  > bxbarp*(bybar*fin(7,i,j,k+1) +by*fin(7,i,j+1,k+1))+
578  > bxp*(bybar*fin(7,i+1,j,k+1)+by*fin(7,i+1,j+1,k+1)))
579  > )
580 C
581  fval(v,iadr)=sum
582  endif
583 C
584  if(ict(3).eq.1) then
585 C
586 C df/dy:
587 C
588  iadr=iadr+1
589 C
590  sum=hyi(v)*(
591  > azbar*(
592  > axbar*(aybarp*fin(0,i,j,k) +ayp*fin(0,i,j+1,k))+
593  > ax*(aybarp*fin(0,i+1,j,k)+ayp*fin(0,i+1,j+1,k)))
594  > + az*(
595  > axbar*(aybarp*fin(0,i,j,k+1) +ayp*fin(0,i,j+1,k+1))+
596  > ax*(aybarp*fin(0,i+1,j,k+1)+ayp*fin(0,i+1,j+1,k+1)))
597  > )
598 C
599  sum=sum+hyi(v)*hx(v)*(
600  > azbar*(
601  > bxbar*(aybarp*fin(1,i,j,k) +ayp*fin(1,i,j+1,k))+
602  > bx*(aybarp*fin(1,i+1,j,k)+ayp*fin(1,i+1,j+1,k)))
603  > + az*(
604  > bxbar*(aybarp*fin(1,i,j,k+1) +ayp*fin(1,i,j+1,k+1))+
605  > bx*(aybarp*fin(1,i+1,j,k+1)+ayp*fin(1,i+1,j+1,k+1)))
606  > )
607 C
608  sum=sum+(
609  > azbar*(
610  > axbar*(bybarp*fin(2,i,j,k) +byp*fin(2,i,j+1,k))+
611  > ax*(bybarp*fin(2,i+1,j,k)+byp*fin(2,i+1,j+1,k)))
612  > + az*(
613  > axbar*(bybarp*fin(2,i,j,k+1) +byp*fin(2,i,j+1,k+1))+
614  > ax*(bybarp*fin(2,i+1,j,k+1)+byp*fin(2,i+1,j+1,k+1)))
615  > )
616 C
617  sum=sum+hyi(v)*hz(v)*(
618  > bzbar*(
619  > axbar*(aybarp*fin(3,i,j,k) +ayp*fin(3,i,j+1,k))+
620  > ax*(aybarp*fin(3,i+1,j,k)+ayp*fin(3,i+1,j+1,k)))
621  > + bz*(
622  > axbar*(aybarp*fin(3,i,j,k+1) +ayp*fin(3,i,j+1,k+1))+
623  > ax*(aybarp*fin(3,i+1,j,k+1)+ayp*fin(3,i+1,j+1,k+1)))
624  > )
625 C
626  sum=sum+hx(v)*(
627  > azbar*(
628  > bxbar*(bybarp*fin(4,i,j,k) +byp*fin(4,i,j+1,k))+
629  > bx*(bybarp*fin(4,i+1,j,k)+byp*fin(4,i+1,j+1,k)))
630  > + az*(
631  > bxbar*(bybarp*fin(4,i,j,k+1) +byp*fin(4,i,j+1,k+1))+
632  > bx*(bybarp*fin(4,i+1,j,k+1)+byp*fin(4,i+1,j+1,k+1)))
633  > )
634 C
635  sum=sum+hx(v)*hyi(v)*hz(v)*(
636  > bzbar*(
637  > bxbar*(aybarp*fin(5,i,j,k) +ayp*fin(5,i,j+1,k))+
638  > bx*(aybarp*fin(5,i+1,j,k)+ayp*fin(5,i+1,j+1,k)))
639  > + bz*(
640  > bxbar*(aybarp*fin(5,i,j,k+1) +ayp*fin(5,i,j+1,k+1))+
641  > bx*(aybarp*fin(5,i+1,j,k+1)+ayp*fin(5,i+1,j+1,k+1)))
642  > )
643 C
644  sum=sum+hz(v)*(
645  > bzbar*(
646  > axbar*(bybarp*fin(6,i,j,k) +byp*fin(6,i,j+1,k))+
647  > ax*(bybarp*fin(6,i+1,j,k)+byp*fin(6,i+1,j+1,k)))
648  > + bz*(
649  > axbar*(bybarp*fin(6,i,j,k+1) +byp*fin(6,i,j+1,k+1))+
650  > ax*(bybarp*fin(6,i+1,j,k+1)+byp*fin(6,i+1,j+1,k+1)))
651  > )
652 C
653  sum=sum+hx(v)*hz(v)*(
654  > bzbar*(
655  > bxbar*(bybarp*fin(7,i,j,k) +byp*fin(7,i,j+1,k))+
656  > bx*(bybarp*fin(7,i+1,j,k)+byp*fin(7,i+1,j+1,k)))
657  > + bz*(
658  > bxbar*(bybarp*fin(7,i,j,k+1) +byp*fin(7,i,j+1,k+1))+
659  > bx*(bybarp*fin(7,i+1,j,k+1)+byp*fin(7,i+1,j+1,k+1)))
660  > )
661 C
662  fval(v,iadr)=sum
663  endif
664 C
665  if(ict(4).eq.1) then
666 C
667 C df/dz:
668 C
669  iadr=iadr+1
670 C
671  sum=hzi(v)*(
672  > azbarp*(
673  > axbar*(aybar*fin(0,i,j,k) +ay*fin(0,i,j+1,k))+
674  > ax*(aybar*fin(0,i+1,j,k)+ay*fin(0,i+1,j+1,k)))
675  > + azp*(
676  > axbar*(aybar*fin(0,i,j,k+1) +ay*fin(0,i,j+1,k+1))+
677  > ax*(aybar*fin(0,i+1,j,k+1)+ay*fin(0,i+1,j+1,k+1)))
678  > )
679 C
680  sum=sum+hzi(v)*hx(v)*(
681  > azbarp*(
682  > bxbar*(aybar*fin(1,i,j,k) +ay*fin(1,i,j+1,k))+
683  > bx*(aybar*fin(1,i+1,j,k)+ay*fin(1,i+1,j+1,k)))
684  > + azp*(
685  > bxbar*(aybar*fin(1,i,j,k+1) +ay*fin(1,i,j+1,k+1))+
686  > bx*(aybar*fin(1,i+1,j,k+1)+ay*fin(1,i+1,j+1,k+1)))
687  > )
688 C
689  sum=sum+hzi(v)*hy(v)*(
690  > azbarp*(
691  > axbar*(bybar*fin(2,i,j,k) +by*fin(2,i,j+1,k))+
692  > ax*(bybar*fin(2,i+1,j,k)+by*fin(2,i+1,j+1,k)))
693  > + azp*(
694  > axbar*(bybar*fin(2,i,j,k+1) +by*fin(2,i,j+1,k+1))+
695  > ax*(bybar*fin(2,i+1,j,k+1)+by*fin(2,i+1,j+1,k+1)))
696  > )
697 C
698  sum=sum+(
699  > bzbarp*(
700  > axbar*(aybar*fin(3,i,j,k) +ay*fin(3,i,j+1,k))+
701  > ax*(aybar*fin(3,i+1,j,k)+ay*fin(3,i+1,j+1,k)))
702  > + bzp*(
703  > axbar*(aybar*fin(3,i,j,k+1) +ay*fin(3,i,j+1,k+1))+
704  > ax*(aybar*fin(3,i+1,j,k+1)+ay*fin(3,i+1,j+1,k+1)))
705  > )
706 C
707  sum=sum+hzi(v)*hx(v)*hy(v)*(
708  > azbarp*(
709  > bxbar*(bybar*fin(4,i,j,k) +by*fin(4,i,j+1,k))+
710  > bx*(bybar*fin(4,i+1,j,k)+by*fin(4,i+1,j+1,k)))
711  > + azp*(
712  > bxbar*(bybar*fin(4,i,j,k+1) +by*fin(4,i,j+1,k+1))+
713  > bx*(bybar*fin(4,i+1,j,k+1)+by*fin(4,i+1,j+1,k+1)))
714  > )
715 C
716  sum=sum+hx(v)*(
717  > bzbarp*(
718  > bxbar*(aybar*fin(5,i,j,k) +ay*fin(5,i,j+1,k))+
719  > bx*(aybar*fin(5,i+1,j,k)+ay*fin(5,i+1,j+1,k)))
720  > + bzp*(
721  > bxbar*(aybar*fin(5,i,j,k+1) +ay*fin(5,i,j+1,k+1))+
722  > bx*(aybar*fin(5,i+1,j,k+1)+ay*fin(5,i+1,j+1,k+1)))
723  > )
724 C
725  sum=sum+hy(v)*(
726  > bzbarp*(
727  > axbar*(bybar*fin(6,i,j,k) +by*fin(6,i,j+1,k))+
728  > ax*(bybar*fin(6,i+1,j,k)+by*fin(6,i+1,j+1,k)))
729  > + bzp*(
730  > axbar*(bybar*fin(6,i,j,k+1) +by*fin(6,i,j+1,k+1))+
731  > ax*(bybar*fin(6,i+1,j,k+1)+by*fin(6,i+1,j+1,k+1)))
732  > )
733 C
734  sum=sum+hx(v)*hy(v)*(
735  > bzbarp*(
736  > bxbar*(bybar*fin(7,i,j,k) +by*fin(7,i,j+1,k))+
737  > bx*(bybar*fin(7,i+1,j,k)+by*fin(7,i+1,j+1,k)))
738  > + bzp*(
739  > bxbar*(bybar*fin(7,i,j,k+1) +by*fin(7,i,j+1,k+1))+
740  > bx*(bybar*fin(7,i+1,j,k+1)+by*fin(7,i+1,j+1,k+1)))
741  > )
742 C
743  fval(v,iadr)=sum
744  endif
745 C
746  if(ict(5).eq.1) then
747 C
748 C d2f/dxdy:
749 C
750  iadr=iadr+1
751 C
752  sum=hxi(v)*hyi(v)*(
753  > azbar*(
754  > axbarp*(aybarp*fin(0,i,j,k) +ayp*fin(0,i,j+1,k))+
755  > axp*(aybarp*fin(0,i+1,j,k)+ayp*fin(0,i+1,j+1,k)))
756  > + az*(
757  > axbarp*(aybarp*fin(0,i,j,k+1) +ayp*fin(0,i,j+1,k+1))+
758  > axp*(aybarp*fin(0,i+1,j,k+1)+ayp*fin(0,i+1,j+1,k+1)))
759  > )
760 C
761  sum=sum+hyi(v)*(
762  > azbar*(
763  > bxbarp*(aybarp*fin(1,i,j,k) +ayp*fin(1,i,j+1,k))+
764  > bxp*(aybarp*fin(1,i+1,j,k)+ayp*fin(1,i+1,j+1,k)))
765  > + az*(
766  > bxbarp*(aybarp*fin(1,i,j,k+1) +ayp*fin(1,i,j+1,k+1))+
767  > bxp*(aybarp*fin(1,i+1,j,k+1)+ayp*fin(1,i+1,j+1,k+1)))
768  > )
769 C
770  sum=sum+hxi(v)*(
771  > azbar*(
772  > axbarp*(bybarp*fin(2,i,j,k) +byp*fin(2,i,j+1,k))+
773  > axp*(bybarp*fin(2,i+1,j,k)+byp*fin(2,i+1,j+1,k)))
774  > + az*(
775  > axbarp*(bybarp*fin(2,i,j,k+1) +byp*fin(2,i,j+1,k+1))+
776  > axp*(bybarp*fin(2,i+1,j,k+1)+byp*fin(2,i+1,j+1,k+1)))
777  > )
778 C
779  sum=sum+hxi(v)*hyi(v)*hz(v)*(
780  > bzbar*(
781  > axbarp*(aybarp*fin(3,i,j,k) +ayp*fin(3,i,j+1,k))+
782  > axp*(aybarp*fin(3,i+1,j,k)+ayp*fin(3,i+1,j+1,k)))
783  > + bz*(
784  > axbarp*(aybarp*fin(3,i,j,k+1) +ayp*fin(3,i,j+1,k+1))+
785  > axp*(aybarp*fin(3,i+1,j,k+1)+ayp*fin(3,i+1,j+1,k+1)))
786  > )
787 C
788  sum=sum+(
789  > azbar*(
790  > bxbarp*(bybarp*fin(4,i,j,k) +byp*fin(4,i,j+1,k))+
791  > bxp*(bybarp*fin(4,i+1,j,k)+byp*fin(4,i+1,j+1,k)))
792  > + az*(
793  > bxbarp*(bybarp*fin(4,i,j,k+1) +byp*fin(4,i,j+1,k+1))+
794  > bxp*(bybarp*fin(4,i+1,j,k+1)+byp*fin(4,i+1,j+1,k+1)))
795  > )
796 C
797  sum=sum+hyi(v)*hz(v)*(
798  > bzbar*(
799  > bxbarp*(aybarp*fin(5,i,j,k) +ayp*fin(5,i,j+1,k))+
800  > bxp*(aybarp*fin(5,i+1,j,k)+ayp*fin(5,i+1,j+1,k)))
801  > + bz*(
802  > bxbarp*(aybarp*fin(5,i,j,k+1) +ayp*fin(5,i,j+1,k+1))+
803  > bxp*(aybarp*fin(5,i+1,j,k+1)+ayp*fin(5,i+1,j+1,k+1)))
804  > )
805 C
806  sum=sum+hxi(v)*hz(v)*(
807  > bzbar*(
808  > axbarp*(bybarp*fin(6,i,j,k) +byp*fin(6,i,j+1,k))+
809  > axp*(bybarp*fin(6,i+1,j,k)+byp*fin(6,i+1,j+1,k)))
810  > + bz*(
811  > axbarp*(bybarp*fin(6,i,j,k+1) +byp*fin(6,i,j+1,k+1))+
812  > axp*(bybarp*fin(6,i+1,j,k+1)+byp*fin(6,i+1,j+1,k+1)))
813  > )
814 C
815  sum=sum+hz(v)*(
816  > bzbar*(
817  > bxbarp*(bybarp*fin(7,i,j,k) +byp*fin(7,i,j+1,k))+
818  > bxp*(bybarp*fin(7,i+1,j,k)+byp*fin(7,i+1,j+1,k)))
819  > + bz*(
820  > bxbarp*(bybarp*fin(7,i,j,k+1) +byp*fin(7,i,j+1,k+1))+
821  > bxp*(bybarp*fin(7,i+1,j,k+1)+byp*fin(7,i+1,j+1,k+1)))
822  > )
823 C
824  fval(v,iadr)=sum
825  endif
826 C
827  if(ict(6).eq.1) then
828 C
829 C d2f/dxdz:
830 C
831  iadr=iadr+1
832 C
833  sum=hxi(v)*hzi(v)*(
834  > azbarp*(
835  > axbarp*(aybar*fin(0,i,j,k) +ay*fin(0,i,j+1,k))+
836  > axp*(aybar*fin(0,i+1,j,k)+ay*fin(0,i+1,j+1,k)))
837  > + azp*(
838  > axbarp*(aybar*fin(0,i,j,k+1) +ay*fin(0,i,j+1,k+1))+
839  > axp*(aybar*fin(0,i+1,j,k+1)+ay*fin(0,i+1,j+1,k+1)))
840  > )
841 C
842  sum=sum+hzi(v)*(
843  > azbarp*(
844  > bxbarp*(aybar*fin(1,i,j,k) +ay*fin(1,i,j+1,k))+
845  > bxp*(aybar*fin(1,i+1,j,k)+ay*fin(1,i+1,j+1,k)))
846  > + azp*(
847  > bxbarp*(aybar*fin(1,i,j,k+1) +ay*fin(1,i,j+1,k+1))+
848  > bxp*(aybar*fin(1,i+1,j,k+1)+ay*fin(1,i+1,j+1,k+1)))
849  > )
850 C
851  sum=sum+hxi(v)*hy(v)*hzi(v)*(
852  > azbarp*(
853  > axbarp*(bybar*fin(2,i,j,k) +by*fin(2,i,j+1,k))+
854  > axp*(bybar*fin(2,i+1,j,k)+by*fin(2,i+1,j+1,k)))
855  > + azp*(
856  > axbarp*(bybar*fin(2,i,j,k+1) +by*fin(2,i,j+1,k+1))+
857  > axp*(bybar*fin(2,i+1,j,k+1)+by*fin(2,i+1,j+1,k+1)))
858  > )
859 C
860  sum=sum+hxi(v)*(
861  > bzbarp*(
862  > axbarp*(aybar*fin(3,i,j,k) +ay*fin(3,i,j+1,k))+
863  > axp*(aybar*fin(3,i+1,j,k)+ay*fin(3,i+1,j+1,k)))
864  > + bzp*(
865  > axbarp*(aybar*fin(3,i,j,k+1) +ay*fin(3,i,j+1,k+1))+
866  > axp*(aybar*fin(3,i+1,j,k+1)+ay*fin(3,i+1,j+1,k+1)))
867  > )
868 C
869  sum=sum+hy(v)*hzi(v)*(
870  > azbarp*(
871  > bxbarp*(bybar*fin(4,i,j,k) +by*fin(4,i,j+1,k))+
872  > bxp*(bybar*fin(4,i+1,j,k)+by*fin(4,i+1,j+1,k)))
873  > + azp*(
874  > bxbarp*(bybar*fin(4,i,j,k+1) +by*fin(4,i,j+1,k+1))+
875  > bxp*(bybar*fin(4,i+1,j,k+1)+by*fin(4,i+1,j+1,k+1)))
876  > )
877 C
878  sum=sum+(
879  > bzbarp*(
880  > bxbarp*(aybar*fin(5,i,j,k) +ay*fin(5,i,j+1,k))+
881  > bxp*(aybar*fin(5,i+1,j,k)+ay*fin(5,i+1,j+1,k)))
882  > + bzp*(
883  > bxbarp*(aybar*fin(5,i,j,k+1) +ay*fin(5,i,j+1,k+1))+
884  > bxp*(aybar*fin(5,i+1,j,k+1)+ay*fin(5,i+1,j+1,k+1)))
885  > )
886 C
887  sum=sum+hxi(v)*hy(v)*(
888  > bzbarp*(
889  > axbarp*(bybar*fin(6,i,j,k) +by*fin(6,i,j+1,k))+
890  > axp*(bybar*fin(6,i+1,j,k)+by*fin(6,i+1,j+1,k)))
891  > + bzp*(
892  > axbarp*(bybar*fin(6,i,j,k+1) +by*fin(6,i,j+1,k+1))+
893  > axp*(bybar*fin(6,i+1,j,k+1)+by*fin(6,i+1,j+1,k+1)))
894  > )
895 C
896  sum=sum+hy(v)*(
897  > bzbarp*(
898  > bxbarp*(bybar*fin(7,i,j,k) +by*fin(7,i,j+1,k))+
899  > bxp*(bybar*fin(7,i+1,j,k)+by*fin(7,i+1,j+1,k)))
900  > + bzp*(
901  > bxbarp*(bybar*fin(7,i,j,k+1) +by*fin(7,i,j+1,k+1))+
902  > bxp*(bybar*fin(7,i+1,j,k+1)+by*fin(7,i+1,j+1,k+1)))
903  > )
904 C
905  fval(v,iadr)=sum
906  endif
907 C
908  if(ict(7).eq.1) then
909 C
910 C d2f/dydz:
911 C
912  iadr=iadr+1
913 C
914  sum=hyi(v)*hzi(v)*(
915  > azbarp*(
916  > axbar*(aybarp*fin(0,i,j,k) +ayp*fin(0,i,j+1,k))+
917  > ax*(aybarp*fin(0,i+1,j,k)+ayp*fin(0,i+1,j+1,k)))
918  > + azp*(
919  > axbar*(aybarp*fin(0,i,j,k+1) +ayp*fin(0,i,j+1,k+1))+
920  > ax*(aybarp*fin(0,i+1,j,k+1)+ayp*fin(0,i+1,j+1,k+1)))
921  > )
922 C
923  sum=sum+hyi(v)*hzi(v)*hx(v)*(
924  > azbarp*(
925  > bxbar*(aybarp*fin(1,i,j,k) +ayp*fin(1,i,j+1,k))+
926  > bx*(aybarp*fin(1,i+1,j,k)+ayp*fin(1,i+1,j+1,k)))
927  > + azp*(
928  > bxbar*(aybarp*fin(1,i,j,k+1) +ayp*fin(1,i,j+1,k+1))+
929  > bx*(aybarp*fin(1,i+1,j,k+1)+ayp*fin(1,i+1,j+1,k+1)))
930  > )
931 C
932  sum=sum+hzi(v)*(
933  > azbarp*(
934  > axbar*(bybarp*fin(2,i,j,k) +byp*fin(2,i,j+1,k))+
935  > ax*(bybarp*fin(2,i+1,j,k)+byp*fin(2,i+1,j+1,k)))
936  > + azp*(
937  > axbar*(bybarp*fin(2,i,j,k+1) +byp*fin(2,i,j+1,k+1))+
938  > ax*(bybarp*fin(2,i+1,j,k+1)+byp*fin(2,i+1,j+1,k+1)))
939  > )
940 C
941  sum=sum+hyi(v)*(
942  > bzbarp*(
943  > axbar*(aybarp*fin(3,i,j,k) +ayp*fin(3,i,j+1,k))+
944  > ax*(aybarp*fin(3,i+1,j,k)+ayp*fin(3,i+1,j+1,k)))
945  > + bzp*(
946  > axbar*(aybarp*fin(3,i,j,k+1) +ayp*fin(3,i,j+1,k+1))+
947  > ax*(aybarp*fin(3,i+1,j,k+1)+ayp*fin(3,i+1,j+1,k+1)))
948  > )
949 C
950  sum=sum+hx(v)*hzi(v)*(
951  > azbarp*(
952  > bxbar*(bybarp*fin(4,i,j,k) +byp*fin(4,i,j+1,k))+
953  > bx*(bybarp*fin(4,i+1,j,k)+byp*fin(4,i+1,j+1,k)))
954  > + azp*(
955  > bxbar*(bybarp*fin(4,i,j,k+1) +byp*fin(4,i,j+1,k+1))+
956  > bx*(bybarp*fin(4,i+1,j,k+1)+byp*fin(4,i+1,j+1,k+1)))
957  > )
958 C
959  sum=sum+hx(v)*hyi(v)*(
960  > bzbarp*(
961  > bxbar*(aybarp*fin(5,i,j,k) +ayp*fin(5,i,j+1,k))+
962  > bx*(aybarp*fin(5,i+1,j,k)+ayp*fin(5,i+1,j+1,k)))
963  > + bzp*(
964  > bxbar*(aybarp*fin(5,i,j,k+1) +ayp*fin(5,i,j+1,k+1))+
965  > bx*(aybarp*fin(5,i+1,j,k+1)+ayp*fin(5,i+1,j+1,k+1)))
966  > )
967 C
968  sum=sum+(
969  > bzbarp*(
970  > axbar*(bybarp*fin(6,i,j,k) +byp*fin(6,i,j+1,k))+
971  > ax*(bybarp*fin(6,i+1,j,k)+byp*fin(6,i+1,j+1,k)))
972  > + bzp*(
973  > axbar*(bybarp*fin(6,i,j,k+1) +byp*fin(6,i,j+1,k+1))+
974  > ax*(bybarp*fin(6,i+1,j,k+1)+byp*fin(6,i+1,j+1,k+1)))
975  > )
976 C
977  sum=sum+hx(v)*(
978  > bzbarp*(
979  > bxbar*(bybarp*fin(7,i,j,k) +byp*fin(7,i,j+1,k))+
980  > bx*(bybarp*fin(7,i+1,j,k)+byp*fin(7,i+1,j+1,k)))
981  > + bzp*(
982  > bxbar*(bybarp*fin(7,i,j,k+1) +byp*fin(7,i,j+1,k+1))+
983  > bx*(bybarp*fin(7,i+1,j,k+1)+byp*fin(7,i+1,j+1,k+1)))
984  > )
985 C
986  fval(v,iadr)=sum
987  endif
988 C
989  if(ict(8).eq.1) then
990 C
991 C d3f/dxdydz:
992 C
993  iadr=iadr+1
994 C
995  sum=hxi(v)*hyi(v)*hzi(v)*(
996  > azbarp*(
997  > axbarp*(aybarp*fin(0,i,j,k) +ayp*fin(0,i,j+1,k))+
998  > axp*(aybarp*fin(0,i+1,j,k)+ayp*fin(0,i+1,j+1,k)))
999  > + azp*(
1000  > axbarp*(aybarp*fin(0,i,j,k+1) +ayp*fin(0,i,j+1,k+1))+
1001  > axp*(aybarp*fin(0,i+1,j,k+1)+ayp*fin(0,i+1,j+1,k+1)))
1002  > )
1003 C
1004  sum=sum+hyi(v)*hzi(v)*(
1005  > azbarp*(
1006  > bxbarp*(aybarp*fin(1,i,j,k) +ayp*fin(1,i,j+1,k))+
1007  > bxp*(aybarp*fin(1,i+1,j,k)+ayp*fin(1,i+1,j+1,k)))
1008  > + azp*(
1009  > bxbarp*(aybarp*fin(1,i,j,k+1) +ayp*fin(1,i,j+1,k+1))+
1010  > bxp*(aybarp*fin(1,i+1,j,k+1)+ayp*fin(1,i+1,j+1,k+1)))
1011  > )
1012 C
1013  sum=sum+hxi(v)*hzi(v)*(
1014  > azbarp*(
1015  > axbarp*(bybarp*fin(2,i,j,k) +byp*fin(2,i,j+1,k))+
1016  > axp*(bybarp*fin(2,i+1,j,k)+byp*fin(2,i+1,j+1,k)))
1017  > + azp*(
1018  > axbarp*(bybarp*fin(2,i,j,k+1) +byp*fin(2,i,j+1,k+1))+
1019  > axp*(bybarp*fin(2,i+1,j,k+1)+byp*fin(2,i+1,j+1,k+1)))
1020  > )
1021 C
1022  sum=sum+hxi(v)*hyi(v)*(
1023  > bzbarp*(
1024  > axbarp*(aybarp*fin(3,i,j,k) +ayp*fin(3,i,j+1,k))+
1025  > axp*(aybarp*fin(3,i+1,j,k)+ayp*fin(3,i+1,j+1,k)))
1026  > + bzp*(
1027  > axbarp*(aybarp*fin(3,i,j,k+1) +ayp*fin(3,i,j+1,k+1))+
1028  > axp*(aybarp*fin(3,i+1,j,k+1)+ayp*fin(3,i+1,j+1,k+1)))
1029  > )
1030 C
1031  sum=sum+hzi(v)*(
1032  > azbarp*(
1033  > bxbarp*(bybarp*fin(4,i,j,k) +byp*fin(4,i,j+1,k))+
1034  > bxp*(bybarp*fin(4,i+1,j,k)+byp*fin(4,i+1,j+1,k)))
1035  > + azp*(
1036  > bxbarp*(bybarp*fin(4,i,j,k+1) +byp*fin(4,i,j+1,k+1))+
1037  > bxp*(bybarp*fin(4,i+1,j,k+1)+byp*fin(4,i+1,j+1,k+1)))
1038  > )
1039 C
1040  sum=sum+hyi(v)*(
1041  > bzbarp*(
1042  > bxbarp*(aybarp*fin(5,i,j,k) +ayp*fin(5,i,j+1,k))+
1043  > bxp*(aybarp*fin(5,i+1,j,k)+ayp*fin(5,i+1,j+1,k)))
1044  > + bzp*(
1045  > bxbarp*(aybarp*fin(5,i,j,k+1) +ayp*fin(5,i,j+1,k+1))+
1046  > bxp*(aybarp*fin(5,i+1,j,k+1)+ayp*fin(5,i+1,j+1,k+1)))
1047  > )
1048 C
1049  sum=sum+hxi(v)*(
1050  > bzbarp*(
1051  > axbarp*(bybarp*fin(6,i,j,k) +byp*fin(6,i,j+1,k))+
1052  > axp*(bybarp*fin(6,i+1,j,k)+byp*fin(6,i+1,j+1,k)))
1053  > + bzp*(
1054  > axbarp*(bybarp*fin(6,i,j,k+1) +byp*fin(6,i,j+1,k+1))+
1055  > axp*(bybarp*fin(6,i+1,j,k+1)+byp*fin(6,i+1,j+1,k+1)))
1056  > )
1057 C
1058  sum=sum+(
1059  > bzbarp*(
1060  > bxbarp*(bybarp*fin(7,i,j,k) +byp*fin(7,i,j+1,k))+
1061  > bxp*(bybarp*fin(7,i+1,j,k)+byp*fin(7,i+1,j+1,k)))
1062  > + bzp*(
1063  > bxbarp*(bybarp*fin(7,i,j,k+1) +byp*fin(7,i,j+1,k+1))+
1064  > bxp*(bybarp*fin(7,i+1,j,k+1)+byp*fin(7,i+1,j+1,k+1)))
1065  > )
1066 C
1067  fval(v,iadr)=sum
1068  endif
1069 C
1070  enddo ! vector loop
1071 C
1072  return
1073  end
1074 C---------------------------------------------------------------------
1075 C evaluate C1 cubic Hermite function interpolation -- 3d fcn
1076 C --vectorized-- dmc 10 Feb 1999
1077 C --optimized for VARIATION along x axis ONLY--
1078 C
1079  subroutine herm3fcnx(ict,ivec,ivecd,
1080  > fval,ii,jj,kk,xparam,yparam,zparam,
1081  > hx,hxi,hy,hyi,hz,hzi,
1082  > fin,inf2,inf3,nz)
1083 C
1084  integer ict(8) ! requested output control
1085  integer ivec ! vector length
1086  integer ivecd ! vector dimension (1st dim of fval)
1087 C
1088  integer ii(ivec),jj,kk ! target cells (i,j,k)
1089  real xparam(ivec),yparam,zparam
1090  ! normalized displacements from (i,j,k) corners
1091 C
1092  real hx(ivec),hy,hz ! grid spacing, and
1093  real hxi(ivec),hyi,hzi ! inverse grid spacing
1094  ! 1/(x(i+1)-x(i)) & 1/(y(j+1)-y(j)) & 1/(z(k+1)-z(i))
1095 C
1096  real fin(0:7,inf2,inf3,nz) ! interpolant data (cf "herm3ev")
1097 C
1098  real fval(ivecd,*) ! output returned
1099 C
1100 C for detailed description of fin, ict and fval see subroutine herm3ev
1101 C comments. Note ict is not vectorized; the same output
1102 C is expected to be returned for all input vector data points.
1103 C
1104 C note that the index inputs ii,jj,kk and parameter inputs
1105 C xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi are vectorized, and the
1106 C output array fval has a vector ** 1st dimension ** whose
1107 C size must be given as a separate argument
1108 C
1109 C to use this routine in scalar mode, pass in ivec=ivecd=1
1110 C
1111 C---------------
1112 C Hermite cubic basis functions
1113 C -->for function value matching
1114 C a(0)=0 a(1)=1 a'(0)=0 a'(1)=0
1115 C abar(0)=1 abar(1)=0 abar'(0)=0 abar'(1)=0
1116 C
1117 C a(x)=-2*x**3 + 3*x**2 = x*x*(-2.0*x+3.0)
1118 C abar(x)=1-a(x)
1119 C a'(x)=-abar'(x) = 6.0*x*(1.0-x)
1120 C
1121 C -->for derivative matching
1122 C b(0)=0 b(1)=0 b'(0)=0 b'(1)=1
1123 C bbar(0)=0 bbar(1)=0 bbar'(0)=1 bbar'(1)=0
1124 C
1125 C b(x)=x**3-x**2 b'(x)=3*x**2-2*x
1126 C bbar(x)=x**3-2*x**2+x bbar'(x)=3*x**2-4*x+1
1127 C
1128  real sum
1129  integer v
1130 C
1131  j=jj
1132  k=kk
1133 C
1134 C ...in y direction
1135 C
1136  yp=yparam
1137  ypi=1.0-yp
1138  yp2=yp*yp
1139  ypi2=ypi*ypi
1140  ay=yp2*(3.0-2.0*yp)
1141  aybar=1.0-ay
1142  by=-yp2*ypi
1143  bybar=ypi2*yp
1144 C
1145 C ...in z direction
1146 C
1147  zp=zparam
1148  zpi=1.0-zp
1149  zp2=zp*zp
1150  zpi2=zpi*zpi
1151  az=zp2*(3.0-2.0*zp)
1152  azbar=1.0-az
1153  bz=-zp2*zpi
1154  bzbar=zpi2*zp
1155 C
1156 C derivatives
1157 C
1158  ayp=6.0*yp*ypi
1159  aybarp=-ayp
1160  byp=yp*(3.0*yp-2.0)
1161  bybarp=ypi*(3.0*ypi-2.0)
1162 C
1163  azp=6.0*zp*zpi
1164  azbarp=-azp
1165  bzp=zp*(3.0*zp-2.0)
1166  bzbarp=zpi*(3.0*zpi-2.0)
1167 
1168  do v=1,ivec
1169  i=ii(v)
1170 C
1171 C ...in x direction
1172 C
1173  xp=xparam(v)
1174  xpi=1.0-xp
1175  xp2=xp*xp
1176  xpi2=xpi*xpi
1177  ax=xp2*(3.0-2.0*xp)
1178  axbar=1.0-ax
1179  bx=-xp2*xpi
1180  bxbar=xpi2*xp
1181 C
1182  iadr=0
1183 C
1184 C derivatives:
1185 C
1186  axp=6.0*xp*xpi
1187  axbarp=-axp
1188  bxp=xp*(3.0*xp-2.0)
1189  bxbarp=xpi*(3.0*xpi-2.0)
1190 C
1191 C get desired values:
1192 C
1193  if(ict(1).eq.1) then
1194 C
1195 C function value:
1196 C
1197  iadr=iadr+1
1198  sum=azbar*(
1199  > axbar*(aybar*fin(0,i,j,k) +ay*fin(0,i,j+1,k))+
1200  > ax*(aybar*fin(0,i+1,j,k)+ay*fin(0,i+1,j+1,k)))
1201  > + az*(
1202  > axbar*(aybar*fin(0,i,j,k+1) +ay*fin(0,i,j+1,k+1))+
1203  > ax*(aybar*fin(0,i+1,j,k+1)+ay*fin(0,i+1,j+1,k+1)))
1204 C
1205  sum=sum+hx(v)*(
1206  > azbar*(
1207  > bxbar*(aybar*fin(1,i,j,k) +ay*fin(1,i,j+1,k))+
1208  > bx*(aybar*fin(1,i+1,j,k)+ay*fin(1,i+1,j+1,k)))
1209  > + az*(
1210  > bxbar*(aybar*fin(1,i,j,k+1) +ay*fin(1,i,j+1,k+1))+
1211  > bx*(aybar*fin(1,i+1,j,k+1)+ay*fin(1,i+1,j+1,k+1)))
1212  > )
1213 C
1214  sum=sum+hy*(
1215  > azbar*(
1216  > axbar*(bybar*fin(2,i,j,k) +by*fin(2,i,j+1,k))+
1217  > ax*(bybar*fin(2,i+1,j,k)+by*fin(2,i+1,j+1,k)))
1218  > + az*(
1219  > axbar*(bybar*fin(2,i,j,k+1) +by*fin(2,i,j+1,k+1))+
1220  > ax*(bybar*fin(2,i+1,j,k+1)+by*fin(2,i+1,j+1,k+1)))
1221  > )
1222 C
1223  sum=sum+hz*(
1224  > bzbar*(
1225  > axbar*(aybar*fin(3,i,j,k) +ay*fin(3,i,j+1,k))+
1226  > ax*(aybar*fin(3,i+1,j,k)+ay*fin(3,i+1,j+1,k)))
1227  > + bz*(
1228  > axbar*(aybar*fin(3,i,j,k+1) +ay*fin(3,i,j+1,k+1))+
1229  > ax*(aybar*fin(3,i+1,j,k+1)+ay*fin(3,i+1,j+1,k+1)))
1230  > )
1231 C
1232  sum=sum+hx(v)*hy*(
1233  > azbar*(
1234  > bxbar*(bybar*fin(4,i,j,k) +by*fin(4,i,j+1,k))+
1235  > bx*(bybar*fin(4,i+1,j,k)+by*fin(4,i+1,j+1,k)))
1236  > + az*(
1237  > bxbar*(bybar*fin(4,i,j,k+1) +by*fin(4,i,j+1,k+1))+
1238  > bx*(bybar*fin(4,i+1,j,k+1)+by*fin(4,i+1,j+1,k+1)))
1239  > )
1240 C
1241  sum=sum+hx(v)*hz*(
1242  > bzbar*(
1243  > bxbar*(aybar*fin(5,i,j,k) +ay*fin(5,i,j+1,k))+
1244  > bx*(aybar*fin(5,i+1,j,k)+ay*fin(5,i+1,j+1,k)))
1245  > + bz*(
1246  > bxbar*(aybar*fin(5,i,j,k+1) +ay*fin(5,i,j+1,k+1))+
1247  > bx*(aybar*fin(5,i+1,j,k+1)+ay*fin(5,i+1,j+1,k+1)))
1248  > )
1249 C
1250  sum=sum+hy*hz*(
1251  > bzbar*(
1252  > axbar*(bybar*fin(6,i,j,k) +by*fin(6,i,j+1,k))+
1253  > ax*(bybar*fin(6,i+1,j,k)+by*fin(6,i+1,j+1,k)))
1254  > + bz*(
1255  > axbar*(bybar*fin(6,i,j,k+1) +by*fin(6,i,j+1,k+1))+
1256  > ax*(bybar*fin(6,i+1,j,k+1)+by*fin(6,i+1,j+1,k+1)))
1257  > )
1258 C
1259  sum=sum+hx(v)*hy*hz*(
1260  > bzbar*(
1261  > bxbar*(bybar*fin(7,i,j,k) +by*fin(7,i,j+1,k))+
1262  > bx*(bybar*fin(7,i+1,j,k)+by*fin(7,i+1,j+1,k)))
1263  > + bz*(
1264  > bxbar*(bybar*fin(7,i,j,k+1) +by*fin(7,i,j+1,k+1))+
1265  > bx*(bybar*fin(7,i+1,j,k+1)+by*fin(7,i+1,j+1,k+1)))
1266  > )
1267 C
1268  fval(v,iadr)=sum
1269  endif
1270 C
1271  if(ict(2).eq.1) then
1272 C
1273 C df/dx:
1274 C
1275  iadr=iadr+1
1276 C
1277  sum=hxi(v)*(
1278  > azbar*(
1279  > axbarp*(aybar*fin(0,i,j,k) +ay*fin(0,i,j+1,k))+
1280  > axp*(aybar*fin(0,i+1,j,k)+ay*fin(0,i+1,j+1,k)))
1281  > + az*(
1282  > axbarp*(aybar*fin(0,i,j,k+1) +ay*fin(0,i,j+1,k+1))+
1283  > axp*(aybar*fin(0,i+1,j,k+1)+ay*fin(0,i+1,j+1,k+1)))
1284  > )
1285 C
1286  sum=sum+(
1287  > azbar*(
1288  > bxbarp*(aybar*fin(1,i,j,k) +ay*fin(1,i,j+1,k))+
1289  > bxp*(aybar*fin(1,i+1,j,k)+ay*fin(1,i+1,j+1,k)))
1290  > + az*(
1291  > bxbarp*(aybar*fin(1,i,j,k+1) +ay*fin(1,i,j+1,k+1))+
1292  > bxp*(aybar*fin(1,i+1,j,k+1)+ay*fin(1,i+1,j+1,k+1)))
1293  > )
1294 C
1295  sum=sum+hxi(v)*hy*(
1296  > azbar*(
1297  > axbarp*(bybar*fin(2,i,j,k) +by*fin(2,i,j+1,k))+
1298  > axp*(bybar*fin(2,i+1,j,k)+by*fin(2,i+1,j+1,k)))
1299  > + az*(
1300  > axbarp*(bybar*fin(2,i,j,k+1) +by*fin(2,i,j+1,k+1))+
1301  > axp*(bybar*fin(2,i+1,j,k+1)+by*fin(2,i+1,j+1,k+1)))
1302  > )
1303 C
1304  sum=sum+hxi(v)*hz*(
1305  > bzbar*(
1306  > axbarp*(aybar*fin(3,i,j,k) +ay*fin(3,i,j+1,k))+
1307  > axp*(aybar*fin(3,i+1,j,k)+ay*fin(3,i+1,j+1,k)))
1308  > + bz*(
1309  > axbarp*(aybar*fin(3,i,j,k+1) +ay*fin(3,i,j+1,k+1))+
1310  > axp*(aybar*fin(3,i+1,j,k+1)+ay*fin(3,i+1,j+1,k+1)))
1311  > )
1312 C
1313  sum=sum+hy*(
1314  > azbar*(
1315  > bxbarp*(bybar*fin(4,i,j,k) +by*fin(4,i,j+1,k))+
1316  > bxp*(bybar*fin(4,i+1,j,k)+by*fin(4,i+1,j+1,k)))
1317  > + az*(
1318  > bxbarp*(bybar*fin(4,i,j,k+1) +by*fin(4,i,j+1,k+1))+
1319  > bxp*(bybar*fin(4,i+1,j,k+1)+by*fin(4,i+1,j+1,k+1)))
1320  > )
1321 C
1322  sum=sum+hz*(
1323  > bzbar*(
1324  > bxbarp*(aybar*fin(5,i,j,k) +ay*fin(5,i,j+1,k))+
1325  > bxp*(aybar*fin(5,i+1,j,k)+ay*fin(5,i+1,j+1,k)))
1326  > + bz*(
1327  > bxbarp*(aybar*fin(5,i,j,k+1) +ay*fin(5,i,j+1,k+1))+
1328  > bxp*(aybar*fin(5,i+1,j,k+1)+ay*fin(5,i+1,j+1,k+1)))
1329  > )
1330 C
1331  sum=sum+hxi(v)*hy*hz*(
1332  > bzbar*(
1333  > axbarp*(bybar*fin(6,i,j,k) +by*fin(6,i,j+1,k))+
1334  > axp*(bybar*fin(6,i+1,j,k)+by*fin(6,i+1,j+1,k)))
1335  > + bz*(
1336  > axbarp*(bybar*fin(6,i,j,k+1) +by*fin(6,i,j+1,k+1))+
1337  > axp*(bybar*fin(6,i+1,j,k+1)+by*fin(6,i+1,j+1,k+1)))
1338  > )
1339 C
1340  sum=sum+hy*hz*(
1341  > bzbar*(
1342  > bxbarp*(bybar*fin(7,i,j,k) +by*fin(7,i,j+1,k))+
1343  > bxp*(bybar*fin(7,i+1,j,k)+by*fin(7,i+1,j+1,k)))
1344  > + bz*(
1345  > bxbarp*(bybar*fin(7,i,j,k+1) +by*fin(7,i,j+1,k+1))+
1346  > bxp*(bybar*fin(7,i+1,j,k+1)+by*fin(7,i+1,j+1,k+1)))
1347  > )
1348 C
1349  fval(v,iadr)=sum
1350  endif
1351 C
1352  if(ict(3).eq.1) then
1353 C
1354 C df/dy:
1355 C
1356  iadr=iadr+1
1357 C
1358  sum=hyi*(
1359  > azbar*(
1360  > axbar*(aybarp*fin(0,i,j,k) +ayp*fin(0,i,j+1,k))+
1361  > ax*(aybarp*fin(0,i+1,j,k)+ayp*fin(0,i+1,j+1,k)))
1362  > + az*(
1363  > axbar*(aybarp*fin(0,i,j,k+1) +ayp*fin(0,i,j+1,k+1))+
1364  > ax*(aybarp*fin(0,i+1,j,k+1)+ayp*fin(0,i+1,j+1,k+1)))
1365  > )
1366 C
1367  sum=sum+hyi*hx(v)*(
1368  > azbar*(
1369  > bxbar*(aybarp*fin(1,i,j,k) +ayp*fin(1,i,j+1,k))+
1370  > bx*(aybarp*fin(1,i+1,j,k)+ayp*fin(1,i+1,j+1,k)))
1371  > + az*(
1372  > bxbar*(aybarp*fin(1,i,j,k+1) +ayp*fin(1,i,j+1,k+1))+
1373  > bx*(aybarp*fin(1,i+1,j,k+1)+ayp*fin(1,i+1,j+1,k+1)))
1374  > )
1375 C
1376  sum=sum+(
1377  > azbar*(
1378  > axbar*(bybarp*fin(2,i,j,k) +byp*fin(2,i,j+1,k))+
1379  > ax*(bybarp*fin(2,i+1,j,k)+byp*fin(2,i+1,j+1,k)))
1380  > + az*(
1381  > axbar*(bybarp*fin(2,i,j,k+1) +byp*fin(2,i,j+1,k+1))+
1382  > ax*(bybarp*fin(2,i+1,j,k+1)+byp*fin(2,i+1,j+1,k+1)))
1383  > )
1384 C
1385  sum=sum+hyi*hz*(
1386  > bzbar*(
1387  > axbar*(aybarp*fin(3,i,j,k) +ayp*fin(3,i,j+1,k))+
1388  > ax*(aybarp*fin(3,i+1,j,k)+ayp*fin(3,i+1,j+1,k)))
1389  > + bz*(
1390  > axbar*(aybarp*fin(3,i,j,k+1) +ayp*fin(3,i,j+1,k+1))+
1391  > ax*(aybarp*fin(3,i+1,j,k+1)+ayp*fin(3,i+1,j+1,k+1)))
1392  > )
1393 C
1394  sum=sum+hx(v)*(
1395  > azbar*(
1396  > bxbar*(bybarp*fin(4,i,j,k) +byp*fin(4,i,j+1,k))+
1397  > bx*(bybarp*fin(4,i+1,j,k)+byp*fin(4,i+1,j+1,k)))
1398  > + az*(
1399  > bxbar*(bybarp*fin(4,i,j,k+1) +byp*fin(4,i,j+1,k+1))+
1400  > bx*(bybarp*fin(4,i+1,j,k+1)+byp*fin(4,i+1,j+1,k+1)))
1401  > )
1402 C
1403  sum=sum+hx(v)*hyi*hz*(
1404  > bzbar*(
1405  > bxbar*(aybarp*fin(5,i,j,k) +ayp*fin(5,i,j+1,k))+
1406  > bx*(aybarp*fin(5,i+1,j,k)+ayp*fin(5,i+1,j+1,k)))
1407  > + bz*(
1408  > bxbar*(aybarp*fin(5,i,j,k+1) +ayp*fin(5,i,j+1,k+1))+
1409  > bx*(aybarp*fin(5,i+1,j,k+1)+ayp*fin(5,i+1,j+1,k+1)))
1410  > )
1411 C
1412  sum=sum+hz*(
1413  > bzbar*(
1414  > axbar*(bybarp*fin(6,i,j,k) +byp*fin(6,i,j+1,k))+
1415  > ax*(bybarp*fin(6,i+1,j,k)+byp*fin(6,i+1,j+1,k)))
1416  > + bz*(
1417  > axbar*(bybarp*fin(6,i,j,k+1) +byp*fin(6,i,j+1,k+1))+
1418  > ax*(bybarp*fin(6,i+1,j,k+1)+byp*fin(6,i+1,j+1,k+1)))
1419  > )
1420 C
1421  sum=sum+hx(v)*hz*(
1422  > bzbar*(
1423  > bxbar*(bybarp*fin(7,i,j,k) +byp*fin(7,i,j+1,k))+
1424  > bx*(bybarp*fin(7,i+1,j,k)+byp*fin(7,i+1,j+1,k)))
1425  > + bz*(
1426  > bxbar*(bybarp*fin(7,i,j,k+1) +byp*fin(7,i,j+1,k+1))+
1427  > bx*(bybarp*fin(7,i+1,j,k+1)+byp*fin(7,i+1,j+1,k+1)))
1428  > )
1429 C
1430  fval(v,iadr)=sum
1431  endif
1432 C
1433  if(ict(4).eq.1) then
1434 C
1435 C df/dz:
1436 C
1437  iadr=iadr+1
1438 C
1439  sum=hzi*(
1440  > azbarp*(
1441  > axbar*(aybar*fin(0,i,j,k) +ay*fin(0,i,j+1,k))+
1442  > ax*(aybar*fin(0,i+1,j,k)+ay*fin(0,i+1,j+1,k)))
1443  > + azp*(
1444  > axbar*(aybar*fin(0,i,j,k+1) +ay*fin(0,i,j+1,k+1))+
1445  > ax*(aybar*fin(0,i+1,j,k+1)+ay*fin(0,i+1,j+1,k+1)))
1446  > )
1447 C
1448  sum=sum+hzi*hx(v)*(
1449  > azbarp*(
1450  > bxbar*(aybar*fin(1,i,j,k) +ay*fin(1,i,j+1,k))+
1451  > bx*(aybar*fin(1,i+1,j,k)+ay*fin(1,i+1,j+1,k)))
1452  > + azp*(
1453  > bxbar*(aybar*fin(1,i,j,k+1) +ay*fin(1,i,j+1,k+1))+
1454  > bx*(aybar*fin(1,i+1,j,k+1)+ay*fin(1,i+1,j+1,k+1)))
1455  > )
1456 C
1457  sum=sum+hzi*hy*(
1458  > azbarp*(
1459  > axbar*(bybar*fin(2,i,j,k) +by*fin(2,i,j+1,k))+
1460  > ax*(bybar*fin(2,i+1,j,k)+by*fin(2,i+1,j+1,k)))
1461  > + azp*(
1462  > axbar*(bybar*fin(2,i,j,k+1) +by*fin(2,i,j+1,k+1))+
1463  > ax*(bybar*fin(2,i+1,j,k+1)+by*fin(2,i+1,j+1,k+1)))
1464  > )
1465 C
1466  sum=sum+(
1467  > bzbarp*(
1468  > axbar*(aybar*fin(3,i,j,k) +ay*fin(3,i,j+1,k))+
1469  > ax*(aybar*fin(3,i+1,j,k)+ay*fin(3,i+1,j+1,k)))
1470  > + bzp*(
1471  > axbar*(aybar*fin(3,i,j,k+1) +ay*fin(3,i,j+1,k+1))+
1472  > ax*(aybar*fin(3,i+1,j,k+1)+ay*fin(3,i+1,j+1,k+1)))
1473  > )
1474 C
1475  sum=sum+hzi*hx(v)*hy*(
1476  > azbarp*(
1477  > bxbar*(bybar*fin(4,i,j,k) +by*fin(4,i,j+1,k))+
1478  > bx*(bybar*fin(4,i+1,j,k)+by*fin(4,i+1,j+1,k)))
1479  > + azp*(
1480  > bxbar*(bybar*fin(4,i,j,k+1) +by*fin(4,i,j+1,k+1))+
1481  > bx*(bybar*fin(4,i+1,j,k+1)+by*fin(4,i+1,j+1,k+1)))
1482  > )
1483 C
1484  sum=sum+hx(v)*(
1485  > bzbarp*(
1486  > bxbar*(aybar*fin(5,i,j,k) +ay*fin(5,i,j+1,k))+
1487  > bx*(aybar*fin(5,i+1,j,k)+ay*fin(5,i+1,j+1,k)))
1488  > + bzp*(
1489  > bxbar*(aybar*fin(5,i,j,k+1) +ay*fin(5,i,j+1,k+1))+
1490  > bx*(aybar*fin(5,i+1,j,k+1)+ay*fin(5,i+1,j+1,k+1)))
1491  > )
1492 C
1493  sum=sum+hy*(
1494  > bzbarp*(
1495  > axbar*(bybar*fin(6,i,j,k) +by*fin(6,i,j+1,k))+
1496  > ax*(bybar*fin(6,i+1,j,k)+by*fin(6,i+1,j+1,k)))
1497  > + bzp*(
1498  > axbar*(bybar*fin(6,i,j,k+1) +by*fin(6,i,j+1,k+1))+
1499  > ax*(bybar*fin(6,i+1,j,k+1)+by*fin(6,i+1,j+1,k+1)))
1500  > )
1501 C
1502  sum=sum+hx(v)*hy*(
1503  > bzbarp*(
1504  > bxbar*(bybar*fin(7,i,j,k) +by*fin(7,i,j+1,k))+
1505  > bx*(bybar*fin(7,i+1,j,k)+by*fin(7,i+1,j+1,k)))
1506  > + bzp*(
1507  > bxbar*(bybar*fin(7,i,j,k+1) +by*fin(7,i,j+1,k+1))+
1508  > bx*(bybar*fin(7,i+1,j,k+1)+by*fin(7,i+1,j+1,k+1)))
1509  > )
1510 C
1511  fval(v,iadr)=sum
1512  endif
1513 C
1514  if(ict(5).eq.1) then
1515 C
1516 C d2f/dxdy:
1517 C
1518  iadr=iadr+1
1519 C
1520  sum=hxi(v)*hyi*(
1521  > azbar*(
1522  > axbarp*(aybarp*fin(0,i,j,k) +ayp*fin(0,i,j+1,k))+
1523  > axp*(aybarp*fin(0,i+1,j,k)+ayp*fin(0,i+1,j+1,k)))
1524  > + az*(
1525  > axbarp*(aybarp*fin(0,i,j,k+1) +ayp*fin(0,i,j+1,k+1))+
1526  > axp*(aybarp*fin(0,i+1,j,k+1)+ayp*fin(0,i+1,j+1,k+1)))
1527  > )
1528 C
1529  sum=sum+hyi*(
1530  > azbar*(
1531  > bxbarp*(aybarp*fin(1,i,j,k) +ayp*fin(1,i,j+1,k))+
1532  > bxp*(aybarp*fin(1,i+1,j,k)+ayp*fin(1,i+1,j+1,k)))
1533  > + az*(
1534  > bxbarp*(aybarp*fin(1,i,j,k+1) +ayp*fin(1,i,j+1,k+1))+
1535  > bxp*(aybarp*fin(1,i+1,j,k+1)+ayp*fin(1,i+1,j+1,k+1)))
1536  > )
1537 C
1538  sum=sum+hxi(v)*(
1539  > azbar*(
1540  > axbarp*(bybarp*fin(2,i,j,k) +byp*fin(2,i,j+1,k))+
1541  > axp*(bybarp*fin(2,i+1,j,k)+byp*fin(2,i+1,j+1,k)))
1542  > + az*(
1543  > axbarp*(bybarp*fin(2,i,j,k+1) +byp*fin(2,i,j+1,k+1))+
1544  > axp*(bybarp*fin(2,i+1,j,k+1)+byp*fin(2,i+1,j+1,k+1)))
1545  > )
1546 C
1547  sum=sum+hxi(v)*hyi*hz*(
1548  > bzbar*(
1549  > axbarp*(aybarp*fin(3,i,j,k) +ayp*fin(3,i,j+1,k))+
1550  > axp*(aybarp*fin(3,i+1,j,k)+ayp*fin(3,i+1,j+1,k)))
1551  > + bz*(
1552  > axbarp*(aybarp*fin(3,i,j,k+1) +ayp*fin(3,i,j+1,k+1))+
1553  > axp*(aybarp*fin(3,i+1,j,k+1)+ayp*fin(3,i+1,j+1,k+1)))
1554  > )
1555 C
1556  sum=sum+(
1557  > azbar*(
1558  > bxbarp*(bybarp*fin(4,i,j,k) +byp*fin(4,i,j+1,k))+
1559  > bxp*(bybarp*fin(4,i+1,j,k)+byp*fin(4,i+1,j+1,k)))
1560  > + az*(
1561  > bxbarp*(bybarp*fin(4,i,j,k+1) +byp*fin(4,i,j+1,k+1))+
1562  > bxp*(bybarp*fin(4,i+1,j,k+1)+byp*fin(4,i+1,j+1,k+1)))
1563  > )
1564 C
1565  sum=sum+hyi*hz*(
1566  > bzbar*(
1567  > bxbarp*(aybarp*fin(5,i,j,k) +ayp*fin(5,i,j+1,k))+
1568  > bxp*(aybarp*fin(5,i+1,j,k)+ayp*fin(5,i+1,j+1,k)))
1569  > + bz*(
1570  > bxbarp*(aybarp*fin(5,i,j,k+1) +ayp*fin(5,i,j+1,k+1))+
1571  > bxp*(aybarp*fin(5,i+1,j,k+1)+ayp*fin(5,i+1,j+1,k+1)))
1572  > )
1573 C
1574  sum=sum+hxi(v)*hz*(
1575  > bzbar*(
1576  > axbarp*(bybarp*fin(6,i,j,k) +byp*fin(6,i,j+1,k))+
1577  > axp*(bybarp*fin(6,i+1,j,k)+byp*fin(6,i+1,j+1,k)))
1578  > + bz*(
1579  > axbarp*(bybarp*fin(6,i,j,k+1) +byp*fin(6,i,j+1,k+1))+
1580  > axp*(bybarp*fin(6,i+1,j,k+1)+byp*fin(6,i+1,j+1,k+1)))
1581  > )
1582 C
1583  sum=sum+hz*(
1584  > bzbar*(
1585  > bxbarp*(bybarp*fin(7,i,j,k) +byp*fin(7,i,j+1,k))+
1586  > bxp*(bybarp*fin(7,i+1,j,k)+byp*fin(7,i+1,j+1,k)))
1587  > + bz*(
1588  > bxbarp*(bybarp*fin(7,i,j,k+1) +byp*fin(7,i,j+1,k+1))+
1589  > bxp*(bybarp*fin(7,i+1,j,k+1)+byp*fin(7,i+1,j+1,k+1)))
1590  > )
1591 C
1592  fval(v,iadr)=sum
1593  endif
1594 C
1595  if(ict(6).eq.1) then
1596 C
1597 C d2f/dxdz:
1598 C
1599  iadr=iadr+1
1600 C
1601  sum=hxi(v)*hzi*(
1602  > azbarp*(
1603  > axbarp*(aybar*fin(0,i,j,k) +ay*fin(0,i,j+1,k))+
1604  > axp*(aybar*fin(0,i+1,j,k)+ay*fin(0,i+1,j+1,k)))
1605  > + azp*(
1606  > axbarp*(aybar*fin(0,i,j,k+1) +ay*fin(0,i,j+1,k+1))+
1607  > axp*(aybar*fin(0,i+1,j,k+1)+ay*fin(0,i+1,j+1,k+1)))
1608  > )
1609 C
1610  sum=sum+hzi*(
1611  > azbarp*(
1612  > bxbarp*(aybar*fin(1,i,j,k) +ay*fin(1,i,j+1,k))+
1613  > bxp*(aybar*fin(1,i+1,j,k)+ay*fin(1,i+1,j+1,k)))
1614  > + azp*(
1615  > bxbarp*(aybar*fin(1,i,j,k+1) +ay*fin(1,i,j+1,k+1))+
1616  > bxp*(aybar*fin(1,i+1,j,k+1)+ay*fin(1,i+1,j+1,k+1)))
1617  > )
1618 C
1619  sum=sum+hxi(v)*hy*hzi*(
1620  > azbarp*(
1621  > axbarp*(bybar*fin(2,i,j,k) +by*fin(2,i,j+1,k))+
1622  > axp*(bybar*fin(2,i+1,j,k)+by*fin(2,i+1,j+1,k)))
1623  > + azp*(
1624  > axbarp*(bybar*fin(2,i,j,k+1) +by*fin(2,i,j+1,k+1))+
1625  > axp*(bybar*fin(2,i+1,j,k+1)+by*fin(2,i+1,j+1,k+1)))
1626  > )
1627 C
1628  sum=sum+hxi(v)*(
1629  > bzbarp*(
1630  > axbarp*(aybar*fin(3,i,j,k) +ay*fin(3,i,j+1,k))+
1631  > axp*(aybar*fin(3,i+1,j,k)+ay*fin(3,i+1,j+1,k)))
1632  > + bzp*(
1633  > axbarp*(aybar*fin(3,i,j,k+1) +ay*fin(3,i,j+1,k+1))+
1634  > axp*(aybar*fin(3,i+1,j,k+1)+ay*fin(3,i+1,j+1,k+1)))
1635  > )
1636 C
1637  sum=sum+hy*hzi*(
1638  > azbarp*(
1639  > bxbarp*(bybar*fin(4,i,j,k) +by*fin(4,i,j+1,k))+
1640  > bxp*(bybar*fin(4,i+1,j,k)+by*fin(4,i+1,j+1,k)))
1641  > + azp*(
1642  > bxbarp*(bybar*fin(4,i,j,k+1) +by*fin(4,i,j+1,k+1))+
1643  > bxp*(bybar*fin(4,i+1,j,k+1)+by*fin(4,i+1,j+1,k+1)))
1644  > )
1645 C
1646  sum=sum+(
1647  > bzbarp*(
1648  > bxbarp*(aybar*fin(5,i,j,k) +ay*fin(5,i,j+1,k))+
1649  > bxp*(aybar*fin(5,i+1,j,k)+ay*fin(5,i+1,j+1,k)))
1650  > + bzp*(
1651  > bxbarp*(aybar*fin(5,i,j,k+1) +ay*fin(5,i,j+1,k+1))+
1652  > bxp*(aybar*fin(5,i+1,j,k+1)+ay*fin(5,i+1,j+1,k+1)))
1653  > )
1654 C
1655  sum=sum+hxi(v)*hy*(
1656  > bzbarp*(
1657  > axbarp*(bybar*fin(6,i,j,k) +by*fin(6,i,j+1,k))+
1658  > axp*(bybar*fin(6,i+1,j,k)+by*fin(6,i+1,j+1,k)))
1659  > + bzp*(
1660  > axbarp*(bybar*fin(6,i,j,k+1) +by*fin(6,i,j+1,k+1))+
1661  > axp*(bybar*fin(6,i+1,j,k+1)+by*fin(6,i+1,j+1,k+1)))
1662  > )
1663 C
1664  sum=sum+hy*(
1665  > bzbarp*(
1666  > bxbarp*(bybar*fin(7,i,j,k) +by*fin(7,i,j+1,k))+
1667  > bxp*(bybar*fin(7,i+1,j,k)+by*fin(7,i+1,j+1,k)))
1668  > + bzp*(
1669  > bxbarp*(bybar*fin(7,i,j,k+1) +by*fin(7,i,j+1,k+1))+
1670  > bxp*(bybar*fin(7,i+1,j,k+1)+by*fin(7,i+1,j+1,k+1)))
1671  > )
1672 C
1673  fval(v,iadr)=sum
1674  endif
1675 C
1676  if(ict(7).eq.1) then
1677 C
1678 C d2f/dydz:
1679 C
1680  iadr=iadr+1
1681 C
1682  sum=hyi*hzi*(
1683  > azbarp*(
1684  > axbar*(aybarp*fin(0,i,j,k) +ayp*fin(0,i,j+1,k))+
1685  > ax*(aybarp*fin(0,i+1,j,k)+ayp*fin(0,i+1,j+1,k)))
1686  > + azp*(
1687  > axbar*(aybarp*fin(0,i,j,k+1) +ayp*fin(0,i,j+1,k+1))+
1688  > ax*(aybarp*fin(0,i+1,j,k+1)+ayp*fin(0,i+1,j+1,k+1)))
1689  > )
1690 C
1691  sum=sum+hyi*hzi*hx(v)*(
1692  > azbarp*(
1693  > bxbar*(aybarp*fin(1,i,j,k) +ayp*fin(1,i,j+1,k))+
1694  > bx*(aybarp*fin(1,i+1,j,k)+ayp*fin(1,i+1,j+1,k)))
1695  > + azp*(
1696  > bxbar*(aybarp*fin(1,i,j,k+1) +ayp*fin(1,i,j+1,k+1))+
1697  > bx*(aybarp*fin(1,i+1,j,k+1)+ayp*fin(1,i+1,j+1,k+1)))
1698  > )
1699 C
1700  sum=sum+hzi*(
1701  > azbarp*(
1702  > axbar*(bybarp*fin(2,i,j,k) +byp*fin(2,i,j+1,k))+
1703  > ax*(bybarp*fin(2,i+1,j,k)+byp*fin(2,i+1,j+1,k)))
1704  > + azp*(
1705  > axbar*(bybarp*fin(2,i,j,k+1) +byp*fin(2,i,j+1,k+1))+
1706  > ax*(bybarp*fin(2,i+1,j,k+1)+byp*fin(2,i+1,j+1,k+1)))
1707  > )
1708 C
1709  sum=sum+hyi*(
1710  > bzbarp*(
1711  > axbar*(aybarp*fin(3,i,j,k) +ayp*fin(3,i,j+1,k))+
1712  > ax*(aybarp*fin(3,i+1,j,k)+ayp*fin(3,i+1,j+1,k)))
1713  > + bzp*(
1714  > axbar*(aybarp*fin(3,i,j,k+1) +ayp*fin(3,i,j+1,k+1))+
1715  > ax*(aybarp*fin(3,i+1,j,k+1)+ayp*fin(3,i+1,j+1,k+1)))
1716  > )
1717 C
1718  sum=sum+hx(v)*hzi*(
1719  > azbarp*(
1720  > bxbar*(bybarp*fin(4,i,j,k) +byp*fin(4,i,j+1,k))+
1721  > bx*(bybarp*fin(4,i+1,j,k)+byp*fin(4,i+1,j+1,k)))
1722  > + azp*(
1723  > bxbar*(bybarp*fin(4,i,j,k+1) +byp*fin(4,i,j+1,k+1))+
1724  > bx*(bybarp*fin(4,i+1,j,k+1)+byp*fin(4,i+1,j+1,k+1)))
1725  > )
1726 C
1727  sum=sum+hx(v)*hyi*(
1728  > bzbarp*(
1729  > bxbar*(aybarp*fin(5,i,j,k) +ayp*fin(5,i,j+1,k))+
1730  > bx*(aybarp*fin(5,i+1,j,k)+ayp*fin(5,i+1,j+1,k)))
1731  > + bzp*(
1732  > bxbar*(aybarp*fin(5,i,j,k+1) +ayp*fin(5,i,j+1,k+1))+
1733  > bx*(aybarp*fin(5,i+1,j,k+1)+ayp*fin(5,i+1,j+1,k+1)))
1734  > )
1735 C
1736  sum=sum+(
1737  > bzbarp*(
1738  > axbar*(bybarp*fin(6,i,j,k) +byp*fin(6,i,j+1,k))+
1739  > ax*(bybarp*fin(6,i+1,j,k)+byp*fin(6,i+1,j+1,k)))
1740  > + bzp*(
1741  > axbar*(bybarp*fin(6,i,j,k+1) +byp*fin(6,i,j+1,k+1))+
1742  > ax*(bybarp*fin(6,i+1,j,k+1)+byp*fin(6,i+1,j+1,k+1)))
1743  > )
1744 C
1745  sum=sum+hx(v)*(
1746  > bzbarp*(
1747  > bxbar*(bybarp*fin(7,i,j,k) +byp*fin(7,i,j+1,k))+
1748  > bx*(bybarp*fin(7,i+1,j,k)+byp*fin(7,i+1,j+1,k)))
1749  > + bzp*(
1750  > bxbar*(bybarp*fin(7,i,j,k+1) +byp*fin(7,i,j+1,k+1))+
1751  > bx*(bybarp*fin(7,i+1,j,k+1)+byp*fin(7,i+1,j+1,k+1)))
1752  > )
1753 C
1754  fval(v,iadr)=sum
1755  endif
1756 C
1757  if(ict(8).eq.1) then
1758 C
1759 C d3f/dxdydz:
1760 C
1761  iadr=iadr+1
1762 C
1763  sum=hxi(v)*hyi*hzi*(
1764  > azbarp*(
1765  > axbarp*(aybarp*fin(0,i,j,k) +ayp*fin(0,i,j+1,k))+
1766  > axp*(aybarp*fin(0,i+1,j,k)+ayp*fin(0,i+1,j+1,k)))
1767  > + azp*(
1768  > axbarp*(aybarp*fin(0,i,j,k+1) +ayp*fin(0,i,j+1,k+1))+
1769  > axp*(aybarp*fin(0,i+1,j,k+1)+ayp*fin(0,i+1,j+1,k+1)))
1770  > )
1771 C
1772  sum=sum+hyi*hzi*(
1773  > azbarp*(
1774  > bxbarp*(aybarp*fin(1,i,j,k) +ayp*fin(1,i,j+1,k))+
1775  > bxp*(aybarp*fin(1,i+1,j,k)+ayp*fin(1,i+1,j+1,k)))
1776  > + azp*(
1777  > bxbarp*(aybarp*fin(1,i,j,k+1) +ayp*fin(1,i,j+1,k+1))+
1778  > bxp*(aybarp*fin(1,i+1,j,k+1)+ayp*fin(1,i+1,j+1,k+1)))
1779  > )
1780 C
1781  sum=sum+hxi(v)*hzi*(
1782  > azbarp*(
1783  > axbarp*(bybarp*fin(2,i,j,k) +byp*fin(2,i,j+1,k))+
1784  > axp*(bybarp*fin(2,i+1,j,k)+byp*fin(2,i+1,j+1,k)))
1785  > + azp*(
1786  > axbarp*(bybarp*fin(2,i,j,k+1) +byp*fin(2,i,j+1,k+1))+
1787  > axp*(bybarp*fin(2,i+1,j,k+1)+byp*fin(2,i+1,j+1,k+1)))
1788  > )
1789 C
1790  sum=sum+hxi(v)*hyi*(
1791  > bzbarp*(
1792  > axbarp*(aybarp*fin(3,i,j,k) +ayp*fin(3,i,j+1,k))+
1793  > axp*(aybarp*fin(3,i+1,j,k)+ayp*fin(3,i+1,j+1,k)))
1794  > + bzp*(
1795  > axbarp*(aybarp*fin(3,i,j,k+1) +ayp*fin(3,i,j+1,k+1))+
1796  > axp*(aybarp*fin(3,i+1,j,k+1)+ayp*fin(3,i+1,j+1,k+1)))
1797  > )
1798 C
1799  sum=sum+hzi*(
1800  > azbarp*(
1801  > bxbarp*(bybarp*fin(4,i,j,k) +byp*fin(4,i,j+1,k))+
1802  > bxp*(bybarp*fin(4,i+1,j,k)+byp*fin(4,i+1,j+1,k)))
1803  > + azp*(
1804  > bxbarp*(bybarp*fin(4,i,j,k+1) +byp*fin(4,i,j+1,k+1))+
1805  > bxp*(bybarp*fin(4,i+1,j,k+1)+byp*fin(4,i+1,j+1,k+1)))
1806  > )
1807 C
1808  sum=sum+hyi*(
1809  > bzbarp*(
1810  > bxbarp*(aybarp*fin(5,i,j,k) +ayp*fin(5,i,j+1,k))+
1811  > bxp*(aybarp*fin(5,i+1,j,k)+ayp*fin(5,i+1,j+1,k)))
1812  > + bzp*(
1813  > bxbarp*(aybarp*fin(5,i,j,k+1) +ayp*fin(5,i,j+1,k+1))+
1814  > bxp*(aybarp*fin(5,i+1,j,k+1)+ayp*fin(5,i+1,j+1,k+1)))
1815  > )
1816 C
1817  sum=sum+hxi(v)*(
1818  > bzbarp*(
1819  > axbarp*(bybarp*fin(6,i,j,k) +byp*fin(6,i,j+1,k))+
1820  > axp*(bybarp*fin(6,i+1,j,k)+byp*fin(6,i+1,j+1,k)))
1821  > + bzp*(
1822  > axbarp*(bybarp*fin(6,i,j,k+1) +byp*fin(6,i,j+1,k+1))+
1823  > axp*(bybarp*fin(6,i+1,j,k+1)+byp*fin(6,i+1,j+1,k+1)))
1824  > )
1825 C
1826  sum=sum+(
1827  > bzbarp*(
1828  > bxbarp*(bybarp*fin(7,i,j,k) +byp*fin(7,i,j+1,k))+
1829  > bxp*(bybarp*fin(7,i+1,j,k)+byp*fin(7,i+1,j+1,k)))
1830  > + bzp*(
1831  > bxbarp*(bybarp*fin(7,i,j,k+1) +byp*fin(7,i,j+1,k+1))+
1832  > bxp*(bybarp*fin(7,i+1,j,k+1)+byp*fin(7,i+1,j+1,k+1)))
1833  > )
1834 C
1835  fval(v,iadr)=sum
1836  endif
1837 C
1838  enddo ! vector loop
1839 C
1840  return
1841  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17