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