V3FIT
r8tcspeval.f
1 c tcspeval -- eval tricubic spline function and/or derivatives
2 c
3  subroutine r8tcspeval(xget,yget,zget,iselect,fval,
4  > x,nx,y,ny,z,nz,ilinx,iliny,ilinz,f,inf4,inf5,
5  > ier)
6 c
7  IMPLICIT NONE
8  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
9  integer iselect(10)
10  integer ilinx,iliny,ilinz,nx,ny,nz,inf4,inf5,ier
11 c
12  real*8 xget,yget,zget
13  real*8 fval(*)
14  real*8 x(nx),y(ny),z(nz),f(4,4,4,inf4,inf5,nz)
15 c
16 c modification -- dmc 11 Jan 1999 -- remove SAVE stmts; break routine
17 C into these parts:
18 C
19 C tcspevxyz -- find grid cell of target pt.
20 C tcspevfn -- evaluate function using output of tcspevxyz
21 C
22 C in cases where multiple functions are defined on the same grid,
23 C time can be saved by using tcspevxyz once and then tcspevfn
24 C multiple times.
25 c
26 c input:
27 c (xget,yget,zget) location where interpolated value is desired
28 c x(1).le.xget.le.x(nx) expected
29 c y(1).le.yget.le.y(ny) expected
30 c z(1).le.zget.le.z(nz) expected
31 c
32 c iselect select desired output
33 c
34 c iselect(1)=1 -- want function value (f) itself
35 c iselect(2)=1 -- want df/dx
36 c iselect(3)=1 -- want df/dy
37 c iselect(4)=1 -- want df/dz
38 c iselect(5)=1 -- want d2f/dx2
39 c iselect(6)=1 -- want d2f/dy2
40 c iselect(7)=1 -- want d2f/dz2
41 c iselect(8)=1 -- want d2f/dxdy
42 c iselect(9)=1 -- want d2f/dxdz
43 c iselect(10)=1 -- want d2f/dydz
44 c
45 c
46 c example: iselect(1)=iselect(2)=iselect(3)=iselect(4)=1
47 c f, df/dx, df/dy, and df/dz all evaluated
48 c iselect(5)=iselect(6)=iselect(7)=0
49 c iselect(8)=iselect(9)=iselect(10)=0
50 c 2nd derivatives not evaluated.
51 c
52 c (new dmc Dec 2005 -- higher derivatives available)
53 c iselect(1)=3 --> 3rd derivative, .le.2 diff. in any coordinate
54 c iselect(2:8) select: fxxy fxxz fxyy fxyz fxzz fyyz fyzz
55 c ->note iselect(1)=3, iselect(5)=1 gives fxyz = d3f/dxdydz
56 c iselect(1)=-3 --> 3rd derivative, 3 in one coordinate
57 c iselect(2:4) select: fxxx fyyy fzzz
58 c iselect(1)=4 --> 3rd derivative, .le.2 diff. in any coordinate
59 c iselect(2:7) select: fxxyy fxxyz fxxzz fxyyz fxyzz fyyzz
60 c iselect(1)=-4 --> 3rd derivative, 3 in one coordinate
61 c iselect(2:7) select: fxxxy fxxxz fxyyy fxzzz fyyyz fyzzz
62 c iselect(1)=5 --> 3rd derivative, .le.2 diff. in any coordinate
63 c iselect(2:4) select: fxxyyz fxxyzz fxyyzz
64 c iselect(1)=-5 --> 3rd derivative, 3 in one coordinate
65 c iselect(2:10) select: fxxxyy fxxxyz fxxxzz fxxyyy fxxzzz
66 c fxyyyz fxyzzz fyyyzz fzzzyy
67 c iselect(1)=6 --> 3rd derivative, .le.2 diff. in any coordinate
68 c fxxyyzz
69 c iselect(1)=-6 --> 3rd derivative, 3 in one coordinate
70 c iselect(2:10) select: fxxxyyy fxxxyyz fxxxyzz fxxxyyz
71 c fxxyyyz fxxyzzz fxyyyzz fxyyzzz fyyyzzz
72 c iselect(1)=-7 --> 7th derivative
73 c iselect(2:7) select: fxxxyyyz fxxxyyzz fxxxyzzz
74 c fxxyyyzz fxxyyzzz fxyyyzzz
75 c iselect(1)=-8 --> 8th derivative
76 c iselect(2:4) select: fxxxyyyzz fxxxyyzzz fxxyyyzzz
77 c iselect(1)=-9 --> 9th derivative: fxxxyyyzzz
78 c
79 c-------
80 c
81 c x(1...nx) independent coordinate x, strict ascending
82 c y(1...ny) independent coordinate y, strict ascending
83 c z(1...nz) independent coordinate y, strict ascending
84 c
85 c ilinx -- =1: flag that x is linearly spaced
86 c
87 c see fval (output) description.
88 c
89 c x(1...nx) independent coordinate x, strict ascending
90 c y(1...ny) independent coordinate y, strict ascending
91 c z(1...nz) independent coordinate y, strict ascending
92 c
93 c ilinx -- =1: flag that x is linearly spaced (avoid search for speed)
94 c iliny -- =1: flag that y is linearly spaced (avoid search for speed)
95 c ilinz -- =1: flat that z is linearly spaced (avoid search for speed)
96 c
97 c **CAUTION** actual even spacing of x, y, z is NOT CHECKED HERE!
98 c
99 c
100 c f the function values (at grid points) and spline coefs
101 c
102 c evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i)
103 c and y btw y(j) and y(j+1), dy=y-y(j),
104 c and z btw z(k) and z(k+1), dz=z-z(k)
105 c
106 c do m=1,4
107 c p(m) =
108 c f(1,1,m,i,j,k)+dx*f(2,1,m,i,j,k)+dx**2*f(3,1,m,i,j,k)+dx**3*f(4,1,m,i,j,k)
109 c +dy*(
110 c f(1,2,m,i,j,k)+dx*f(2,2,m,i,j,k)+dx**2*f(3,2,m,i,j,k)+dx**3*f(4,2,m,i,j,k))
111 c +dy**2*(
112 c f(1,3,m,i,j,k)+dx*f(2,3,m,i,j,k)+dx**2*f(3,3,m,i,j,k)+dx**3*f(4,3,m,i,j,k))
113 c +dy**3*(
114 c f(1,4,m,i,j,k)+dx*f(2,4,m,i,j,k)+dx**2*f(3,4,m,i,j,k)+dx**3*f(4,4,m,i,j,k))
115 c enddo
116 c answer = p(1)+dz*p(2)+dz**2*p(3)+dz**3*p(4)
117 c
118 c where d2=dy**2 and d3=dy**3.
119 c
120 c nb dmc Feb 1999 -- p loops unrolled, by hand, to aid vector compilers
121 c
122 c output:
123 c up to 10 elements of fval, ordered as follows:
124 c fval(1)=function value or lowest order derivative requested
125 c fval(2)=next order derivative
126 c etc
127 c the ordering is a subset of the sequence given under the "iselect"
128 c description; the first M elements of fval are used, where M = the
129 c number of non-zero elements of iselect.
130 c
131 c ier = 0 -- successful completion; = 1 -- an error occurred.
132 c
133 c-------------------------------------------------------------------
134 c local
135 c
136  integer :: i(1),j(1),k(1)
137 c
138  real*8 dx(1),dy(1),dz(1)
139 c
140 c--------------------------
141 c
142  i(1)=0
143  j(1)=0
144  k(1)=0
145 c
146  call r8tcspevxyz(xget,yget,zget,x,nx,y,ny,z,nz,ilinx,iliny,ilinz,
147  > i(1),j(1),k(1),dx(1),dy(1),dz(1),ier)
148  if(ier.ne.0) return
149 c
150  call r8tcspevfn(iselect,1,1,fval,i,j,k,dx,dy,dz,f,inf4,inf5,nz)
151 c
152  return
153  end
154 c
155 c-------------------------------------------------------------------------
156 c tcspevxyz -- look up x-y zone
157 c
158 c this is the "first part" of tcspeval, see comments, above.
159 c
160  subroutine r8tcspevxyz(xget,yget,zget,x,nx,y,ny,z,nz,
161  > ilinx,iliny,ilinz,
162  > i,j,k,dx,dy,dz,ier)
163 c
164 !============
165 ! idecl: explicitize implicit INTEGER declarations:
166  IMPLICIT NONE
167  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
168  INTEGER nxm,nym,nzm,ii,jj,kk
169 !============
170 ! idecl: explicitize implicit REAL declarations:
171  real*8 zxget,zyget,zzget,zxtol,zytol,zztol
172 !============
173  integer nx,ny,nz ! array dimensions
174 c
175  real*8 xget,yget,zget ! target point
176  real*8 x(nx),y(ny),z(nz) ! indep. coords.
177 c
178  integer ilinx ! =1: assume x evenly spaced
179  integer iliny ! =1: assume y evenly spaced
180  integer ilinz ! =1: assume z evenly spaced
181 c
182 c output of tcspevxyz
183 c
184  integer i,j,k ! index to cell containing target pt
185  REAL*8 dx,dy,dz ! displacement of target pt w/in cell
186  ! dx=x-x(i) dy=y-y(j) dz=z-z(k)
187 C
188  integer ier ! return ier.ne.0 on error
189 c
190 c------------------------------------
191 c
192  ier=0
193 c
194 c range check
195 c
196  zxget=xget
197  zyget=yget
198  zzget=zget
199 c
200  if((xget.lt.x(1)).or.(xget.gt.x(nx))) then
201  zxtol=4.0e-7_r8*max(abs(x(1)),abs(x(nx)))
202  if((xget.lt.x(1)-zxtol).or.(xget.gt.x(nx)+zxtol)) then
203  ier=1
204  write(6,1001) xget,x(1),x(nx)
205  1001 format(' ?tcspeval: xget=',1pe11.4,' out of range ',
206  > 1pe11.4,' to ',1pe11.4)
207  else
208  if((xget.lt.x(1)-0.5_r8*zxtol).or.
209  > (xget.gt.x(nx)+0.5_r8*zxtol))
210  > write(6,1011) xget,x(1),x(nx)
211  1011 format(' %tcspeval: xget=',1pe15.8,' beyond range ',
212  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
213  if(xget.lt.x(1)) then
214  zxget=x(1)
215  else
216  zxget=x(nx)
217  endif
218  endif
219  endif
220  if((yget.lt.y(1)).or.(yget.gt.y(ny))) then
221  zytol=4.0e-7_r8*max(abs(y(1)),abs(y(ny)))
222  if((yget.lt.y(1)-zytol).or.(yget.gt.y(ny)+zytol)) then
223  ier=1
224  write(6,1002) yget,y(1),y(ny)
225  1002 format(' ?tcspeval: yget=',1pe11.4,' out of range ',
226  > 1pe11.4,' to ',1pe11.4)
227  else
228  if((yget.lt.y(1)-0.5_r8*zytol).or.
229  > (yget.gt.y(ny)+0.5_r8*zytol))
230  > write(6,1012) yget,y(1),y(ny)
231  1012 format(' %tcspeval: yget=',1pe15.8,' beyond range ',
232  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
233  if(yget.lt.y(1)) then
234  zyget=y(1)
235  else
236  zyget=y(ny)
237  endif
238  endif
239  endif
240  if((zget.lt.z(1)).or.(zget.gt.z(nz))) then
241  zztol=4.0e-7_r8*max(abs(z(1)),abs(z(nz)))
242  if((zget.lt.z(1)-zztol).or.(zget.gt.z(nz)+zztol)) then
243  ier=1
244  write(6,1003) zget,z(1),z(nz)
245  1003 format(' ?tcspeval: zget=',1pe11.4,' out of range ',
246  > 1pe11.4,' to ',1pe11.4)
247  else
248  if((zget.lt.z(1)-0.5_r8*zztol).or.
249  > (zget.gt.z(nz)+0.5_r8*zztol))
250  > write(6,1013) zget,z(1),z(nz)
251  1013 format(' %tcspeval: zget=',1pe15.8,' beyond range ',
252  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
253  if(zget.lt.z(1)) then
254  zzget=z(1)
255  else
256  zzget=z(nz)
257  endif
258  endif
259  endif
260  if(ier.ne.0) return
261 c
262 c now find interval in which target point lies..
263 c
264  nxm=nx-1
265  nym=ny-1
266  nzm=nz-1
267 c
268  if(ilinx.eq.1) then
269  ii=1+nxm*(zxget-x(1))/(x(nx)-x(1))
270  i=min(nxm, ii)
271  if(zxget.lt.x(i)) then
272  i=i-1
273  else if(zxget.gt.x(i+1)) then
274  i=i+1
275  endif
276  else
277  if((1.le.i).and.(i.lt.nxm)) then
278  if((x(i).le.zxget).and.(zxget.le.x(i+1))) then
279  continue ! already have the zone
280  else
281  call r8zonfind(x,nx,zxget,i)
282  endif
283  else
284  i=nx/2
285  call r8zonfind(x,nx,zxget,i)
286  endif
287  endif
288 c
289  if(iliny.eq.1) then
290  jj=1+nym*(zyget-y(1))/(y(ny)-y(1))
291  j=min(nym, jj)
292  if(zyget.lt.y(j)) then
293  j=j-1
294  else if(zyget.gt.y(j+1)) then
295  j=j+1
296  endif
297  else
298  if((1.le.j).and.(j.lt.nym)) then
299  if((y(j).le.zyget).and.(zyget.le.y(j+1))) then
300  continue ! already have the zone
301  else
302  call r8zonfind(y,ny,zyget,j)
303  endif
304  else
305  j=ny/2
306  call r8zonfind(y,ny,zyget,j)
307  endif
308  endif
309 c
310  if(ilinz.eq.1) then
311  kk=1+nzm*(zzget-z(1))/(z(nz)-z(1))
312  k=min(nzm, kk)
313  if(zzget.lt.z(k)) then
314  k=k-1
315  else if(zzget.gt.z(k+1)) then
316  k=k+1
317  endif
318  else
319  if((1.le.k).and.(k.lt.nzm)) then
320  if((z(k).le.zzget).and.(zzget.le.z(k+1))) then
321  continue ! already have the zone
322  else
323  call r8zonfind(z,nz,zzget,k)
324  endif
325  else
326  k=nz/2
327  call r8zonfind(z,nz,zzget,k)
328  endif
329  endif
330 c
331  dx=zxget-x(i)
332  dy=zyget-y(j)
333  dz=zzget-z(k)
334 c
335  return
336  end
337 c------------------------------------------------------------------------
338 c tcspevfn -- OK now evaluate the tricubic spline
339 c evaluate at a vector set of target locations as specified by
340 c input vectors (iv,jv,kv), (dxv,dyv,dzv)
341 c
342  subroutine r8tcspevfn(ict,ivec,ivd,fval,iv,jv,kv,dxv,dyv,dzv,
343  > f,inf4,inf5,nz)
344 c
345 c input:
346 !============
347 ! idecl: explicitize implicit INTEGER declarations:
348  IMPLICIT NONE
349  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
350  INTEGER nz,iaval,i,j,k
351 !============
352 ! idecl: explicitize implicit REAL declarations:
353  real*8 dx,dy,dz,p1,p2,p3,p4
354 !============
355  integer ict(10) ! selector:
356 c ict(1)=1 for f (don't evaluate f if ict(1)=0)
357 c ict(2)=1 for df/dx ""
358 c ict(3)=1 for df/dy ""
359 c ict(4)=1 for df/dz ""
360 c ict(5)=1 for d2f/dx2
361 c ict(6)=1 for d2f/dy2
362 c ict(7)=1 for d2f/dz2
363 c ict(8)=1 for d2f/dxdy
364 c ict(9)=1 for d2f/dxdz
365 c ict(10)=1 for d2f/dydz
366 c
367  integer ivec,ivd ! vector dimensioning
368 c
369 c ivec-- number of vector pts (spline values to look up)
370 c ivd -- 1st dimension of fval, .ge. ivec
371 c
372 c output:
373  real*8 fval(ivd,*) ! output array
374 c
375 c for vector elements v, (iv(v),jv(v),kv(v),dxv(v),dyv(v),dzv(v))
376 c
377 c fval(v,1) = first item requested by ict(...),
378 c fval(v,2) = 2nd item requested, ...etc...
379 c
380 c input:
381  integer iv(ivec),jv(ivec),kv(ivec) ! grid cell indices
382  REAL*8 dxv(ivec),dyv(ivec),dzv(ivec) ! displacements w/in cell
383 c
384  integer inf4 ! 4th dimension of f, .le.nx
385  integer inf5 ! 5th dimension of f, .le.ny
386  real*8 f(4,4,4,inf4,inf5,nz) ! tricubic fcn spline coeffs array
387 c
388 c usage example:
389 c
390 c 1. for each element (xx(v),yy(v),zz(v)) in a vector of (x,y,z)
391 c triples, find the x,y,z zone indices and displacements with
392 c to the "lower left corner" of the zone; store these in vectors
393 c iv,jv,kv and dxv,dyv,dzv
394 c
395 c 2. set ict(1)=0, ict(2)=1, ict(3)=1, ict(4)=1 & the rest zero --
396 c to get only the 1st derivatives.
397 c
398 c 3. ivec is the length of the vector; ivd is the 1st dimension
399 c of the array fval to receive the output
400 c
401 c real fval(ivd,10)
402 c real xv(ivd),yv(ivd),zv(ivd)
403 c integer iv(ivd),jv(ivd),kv(ivd)
404 c real dxv(ivd),dyv(ivd),dzv(ivd)
405 c integer ict(10)
406 c
407 c real fspline(4,4,4,nx,ny,nz) ! spline coeffs
408 c data ict/0,1,1,1,0,0,0,0,0,0/ ! this call: want 1st
409 c ! derivatives only ... these will
410 c ! be output to
411 c ! fval(*,1) fval(*,2) fval(*,3)
412 c ...
413 c do iv=1,ivec
414 c ... ! find indices and displacements
415 c enddo
416 c call tcspevfn(ict,ivec,ivd,fval,iv,jv,kv,dxv,dyv,dzv,
417 c > fspline,nx,ny,nz)
418 c
419 c-------------------------------------------------------------------
420 c
421 c local --
422 c
423 ccc real p(4) use p1,p2,p3,p4 now
424 c
425  integer v ! vector index
426 c
427 c OK can now do evaluations
428 c
429  iaval=0 ! fval addressing
430 c
431  if(abs(ict(1)).le.2) then
432  if((ict(1).gt.0).or.(ict(1).eq.-1)) then
433 c evaluate f
434  iaval=iaval+1
435  do v=1,ivec
436  i=iv(v)
437  j=jv(v)
438  k=kv(v)
439  dx=dxv(v)
440  dy=dyv(v)
441  dz=dzv(v)
442  p1=
443  > f(1,1,1,i,j,k)+dy*(f(1,2,1,i,j,k)+dy*(f(1,3,1,i,j,k)+
444  > dy*f(1,4,1,i,j,k)))
445  > +dx*(f(2,1,1,i,j,k)+dy*(f(2,2,1,i,j,k)+dy*(f(2,3,1,i,j,k)+
446  > dy*f(2,4,1,i,j,k)))
447  > +dx*(f(3,1,1,i,j,k)+dy*(f(3,2,1,i,j,k)+dy*(f(3,3,1,i,j,k)+
448  > dy*f(3,4,1,i,j,k)))
449  > +dx*(f(4,1,1,i,j,k)+dy*(f(4,2,1,i,j,k)+dy*(f(4,3,1,i,j,k)+
450  > dy*f(4,4,1,i,j,k)))
451  > )))
452  p2=
453  > f(1,1,2,i,j,k)+dy*(f(1,2,2,i,j,k)+dy*(f(1,3,2,i,j,k)+
454  > dy*f(1,4,2,i,j,k)))
455  > +dx*(f(2,1,2,i,j,k)+dy*(f(2,2,2,i,j,k)+dy*(f(2,3,2,i,j,k)+
456  > dy*f(2,4,2,i,j,k)))
457  > +dx*(f(3,1,2,i,j,k)+dy*(f(3,2,2,i,j,k)+dy*(f(3,3,2,i,j,k)+
458  > dy*f(3,4,2,i,j,k)))
459  > +dx*(f(4,1,2,i,j,k)+dy*(f(4,2,2,i,j,k)+dy*(f(4,3,2,i,j,k)+
460  > dy*f(4,4,2,i,j,k)))
461  > )))
462  p3=
463  > f(1,1,3,i,j,k)+dy*(f(1,2,3,i,j,k)+dy*(f(1,3,3,i,j,k)+
464  > dy*f(1,4,3,i,j,k)))
465  > +dx*(f(2,1,3,i,j,k)+dy*(f(2,2,3,i,j,k)+dy*(f(2,3,3,i,j,k)+
466  > dy*f(2,4,3,i,j,k)))
467  > +dx*(f(3,1,3,i,j,k)+dy*(f(3,2,3,i,j,k)+dy*(f(3,3,3,i,j,k)+
468  > dy*f(3,4,3,i,j,k)))
469  > +dx*(f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
470  > dy*f(4,4,3,i,j,k)))
471  > )))
472  p4=
473  > f(1,1,4,i,j,k)+dy*(f(1,2,4,i,j,k)+dy*(f(1,3,4,i,j,k)+
474  > dy*f(1,4,4,i,j,k)))
475  > +dx*(f(2,1,4,i,j,k)+dy*(f(2,2,4,i,j,k)+dy*(f(2,3,4,i,j,k)+
476  > dy*f(2,4,4,i,j,k)))
477  > +dx*(f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
478  > dy*f(3,4,4,i,j,k)))
479  > +dx*(f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
480  > dy*f(4,4,4,i,j,k)))
481  > )))
482  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
483  enddo
484  endif
485 c
486  if((ict(2).gt.0).and.(ict(1).ne.-1)) then
487 c evaluate df/dx
488  iaval=iaval+1
489  do v=1,ivec
490  i=iv(v)
491  j=jv(v)
492  k=kv(v)
493  dx=dxv(v)
494  dy=dyv(v)
495  dz=dzv(v)
496  p1=
497  > f(2,1,1,i,j,k)+dy*(f(2,2,1,i,j,k)+dy*(f(2,3,1,i,j,k)+
498  > dy*f(2,4,1,i,j,k)))
499  > +2._r8*dx*
500  > (f(3,1,1,i,j,k)+dy*(f(3,2,1,i,j,k)+dy*(f(3,3,1,i,j,k)+
501  > dy*f(3,4,1,i,j,k)))
502  > +1.5_r8*dx*
503  > (f(4,1,1,i,j,k)+dy*(f(4,2,1,i,j,k)+dy*(f(4,3,1,i,j,k)+
504  > dy*f(4,4,1,i,j,k)))
505  > ))
506  p2=
507  > f(2,1,2,i,j,k)+dy*(f(2,2,2,i,j,k)+dy*(f(2,3,2,i,j,k)+
508  > dy*f(2,4,2,i,j,k)))
509  > +2._r8*dx*
510  > (f(3,1,2,i,j,k)+dy*(f(3,2,2,i,j,k)+dy*(f(3,3,2,i,j,k)+
511  > dy*f(3,4,2,i,j,k)))
512  > +1.5_r8*dx*
513  > (f(4,1,2,i,j,k)+dy*(f(4,2,2,i,j,k)+dy*(f(4,3,2,i,j,k)+
514  > dy*f(4,4,2,i,j,k)))
515  > ))
516  p3=
517  > f(2,1,3,i,j,k)+dy*(f(2,2,3,i,j,k)+dy*(f(2,3,3,i,j,k)+
518  > dy*f(2,4,3,i,j,k)))
519  > +2._r8*dx*
520  > (f(3,1,3,i,j,k)+dy*(f(3,2,3,i,j,k)+dy*(f(3,3,3,i,j,k)+
521  > dy*f(3,4,3,i,j,k)))
522  > +1.5_r8*dx*
523  > (f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
524  > dy*f(4,4,3,i,j,k)))
525  > ))
526  p4=
527  > f(2,1,4,i,j,k)+dy*(f(2,2,4,i,j,k)+dy*(f(2,3,4,i,j,k)+
528  > dy*f(2,4,4,i,j,k)))
529  > +2._r8*dx*
530  > (f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
531  > dy*f(3,4,4,i,j,k)))
532  > +1.5_r8*dx*
533  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
534  > dy*f(4,4,4,i,j,k)))
535  > ))
536  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
537  enddo
538  endif
539 c
540  if((ict(3).gt.0).and.(ict(1).ne.-1)) then
541 c evaluate df/dy
542  iaval=iaval+1
543  do v=1,ivec
544  i=iv(v)
545  j=jv(v)
546  k=kv(v)
547  dx=dxv(v)
548  dy=dyv(v)
549  dz=dzv(v)
550  p1=
551  > f(1,2,1,i,j,k)+
552  > dy*(2.0_r8*f(1,3,1,i,j,k)+dy*3.0_r8*f(1,4,1,i,j,k))
553  >+dx*(f(2,2,1,i,j,k)+
554  > dy*(2.0_r8*f(2,3,1,i,j,k)+dy*3.0_r8*f(2,4,1,i,j,k))
555  >+dx*(f(3,2,1,i,j,k)+
556  > dy*(2.0_r8*f(3,3,1,i,j,k)+dy*3.0_r8*f(3,4,1,i,j,k))
557  >+dx*(f(4,2,1,i,j,k)+
558  > dy*(2.0_r8*f(4,3,1,i,j,k)+dy*3.0_r8*f(4,4,1,i,j,k))
559  > )))
560  p2=
561  > f(1,2,2,i,j,k)+
562  > dy*(2.0_r8*f(1,3,2,i,j,k)+dy*3.0_r8*f(1,4,2,i,j,k))
563  >+dx*(f(2,2,2,i,j,k)+
564  > dy*(2.0_r8*f(2,3,2,i,j,k)+dy*3.0_r8*f(2,4,2,i,j,k))
565  >+dx*(f(3,2,2,i,j,k)+
566  > dy*(2.0_r8*f(3,3,2,i,j,k)+dy*3.0_r8*f(3,4,2,i,j,k))
567  >+dx*(f(4,2,2,i,j,k)+
568  > dy*(2.0_r8*f(4,3,2,i,j,k)+dy*3.0_r8*f(4,4,2,i,j,k))
569  > )))
570  p3=
571  > f(1,2,3,i,j,k)+
572  > dy*(2.0_r8*f(1,3,3,i,j,k)+dy*3.0_r8*f(1,4,3,i,j,k))
573  >+dx*(f(2,2,3,i,j,k)+
574  > dy*(2.0_r8*f(2,3,3,i,j,k)+dy*3.0_r8*f(2,4,3,i,j,k))
575  >+dx*(f(3,2,3,i,j,k)+
576  > dy*(2.0_r8*f(3,3,3,i,j,k)+dy*3.0_r8*f(3,4,3,i,j,k))
577  >+dx*(f(4,2,3,i,j,k)+
578  > dy*(2.0_r8*f(4,3,3,i,j,k)+dy*3.0_r8*f(4,4,3,i,j,k))
579  > )))
580  p4=
581  > f(1,2,4,i,j,k)+
582  > dy*(2.0_r8*f(1,3,4,i,j,k)+dy*3.0_r8*f(1,4,4,i,j,k))
583  >+dx*(f(2,2,4,i,j,k)+
584  > dy*(2.0_r8*f(2,3,4,i,j,k)+dy*3.0_r8*f(2,4,4,i,j,k))
585  >+dx*(f(3,2,4,i,j,k)+
586  > dy*(2.0_r8*f(3,3,4,i,j,k)+dy*3.0_r8*f(3,4,4,i,j,k))
587  >+dx*(f(4,2,4,i,j,k)+
588  > dy*(2.0_r8*f(4,3,4,i,j,k)+dy*3.0_r8*f(4,4,4,i,j,k))
589  > )))
590  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
591  enddo
592  endif
593 c
594  if((ict(4).gt.0).and.(ict(1).ne.-1)) then
595 c evaluate df/dz
596  iaval=iaval+1
597  do v=1,ivec
598  i=iv(v)
599  j=jv(v)
600  k=kv(v)
601  dx=dxv(v)
602  dy=dyv(v)
603  dz=dzv(v)
604  p2=
605  > f(1,1,2,i,j,k)+dy*(f(1,2,2,i,j,k)+dy*(f(1,3,2,i,j,k)+
606  > dy*f(1,4,2,i,j,k)))
607  > +dx*(f(2,1,2,i,j,k)+dy*(f(2,2,2,i,j,k)+dy*(f(2,3,2,i,j,k)+
608  > dy*f(2,4,2,i,j,k)))
609  > +dx*(f(3,1,2,i,j,k)+dy*(f(3,2,2,i,j,k)+dy*(f(3,3,2,i,j,k)+
610  > dy*f(3,4,2,i,j,k)))
611  > +dx*(f(4,1,2,i,j,k)+dy*(f(4,2,2,i,j,k)+dy*(f(4,3,2,i,j,k)+
612  > dy*f(4,4,2,i,j,k)))
613  > )))
614  p3=
615  > f(1,1,3,i,j,k)+dy*(f(1,2,3,i,j,k)+dy*(f(1,3,3,i,j,k)+
616  > dy*f(1,4,3,i,j,k)))
617  > +dx*(f(2,1,3,i,j,k)+dy*(f(2,2,3,i,j,k)+dy*(f(2,3,3,i,j,k)+
618  > dy*f(2,4,3,i,j,k)))
619  > +dx*(f(3,1,3,i,j,k)+dy*(f(3,2,3,i,j,k)+dy*(f(3,3,3,i,j,k)+
620  > dy*f(3,4,3,i,j,k)))
621  > +dx*(f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
622  > dy*f(4,4,3,i,j,k)))
623  > )))
624  p4=
625  > f(1,1,4,i,j,k)+dy*(f(1,2,4,i,j,k)+dy*(f(1,3,4,i,j,k)+
626  > dy*f(1,4,4,i,j,k)))
627  > +dx*(f(2,1,4,i,j,k)+dy*(f(2,2,4,i,j,k)+dy*(f(2,3,4,i,j,k)+
628  > dy*f(2,4,4,i,j,k)))
629  > +dx*(f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
630  > dy*f(3,4,4,i,j,k)))
631  > +dx*(f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
632  > dy*f(4,4,4,i,j,k)))
633  > )))
634  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
635  enddo
636  endif
637 c
638  if((ict(5).gt.0).or.(ict(1).eq.-1)) then
639 c evaluate d2f/dx2
640  iaval=iaval+1
641  do v=1,ivec
642  i=iv(v)
643  j=jv(v)
644  k=kv(v)
645  dx=dxv(v)
646  dy=dyv(v)
647  dz=dzv(v)
648  p1= 2.0_r8*
649  > (f(3,1,1,i,j,k)+dy*(f(3,2,1,i,j,k)+dy*(f(3,3,1,i,j,k)+
650  > dy*f(3,4,1,i,j,k))))
651  > +6.0_r8*dx*
652  > (f(4,1,1,i,j,k)+dy*(f(4,2,1,i,j,k)+dy*(f(4,3,1,i,j,k)+
653  > dy*f(4,4,1,i,j,k))))
654  p2= 2.0_r8*
655  > (f(3,1,2,i,j,k)+dy*(f(3,2,2,i,j,k)+dy*(f(3,3,2,i,j,k)+
656  > dy*f(3,4,2,i,j,k))))
657  > +6.0_r8*dx*
658  > (f(4,1,2,i,j,k)+dy*(f(4,2,2,i,j,k)+dy*(f(4,3,2,i,j,k)+
659  > dy*f(4,4,2,i,j,k))))
660  p3= 2.0_r8*
661  > (f(3,1,3,i,j,k)+dy*(f(3,2,3,i,j,k)+dy*(f(3,3,3,i,j,k)+
662  > dy*f(3,4,3,i,j,k))))
663  > +6.0_r8*dx*
664  > (f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
665  > dy*f(4,4,3,i,j,k))))
666  p4= 2.0_r8*
667  > (f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
668  > dy*f(3,4,4,i,j,k))))
669  > +6.0_r8*dx*
670  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
671  > dy*f(4,4,4,i,j,k))))
672  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
673  enddo
674  endif
675 c
676  if((ict(6).gt.0).or.(ict(1).eq.-1)) then
677 c evaluate d2f/dy2
678  iaval=iaval+1
679  do v=1,ivec
680  i=iv(v)
681  j=jv(v)
682  k=kv(v)
683  dx=dxv(v)
684  dy=dyv(v)
685  dz=dzv(v)
686  p1=
687  > 2.0_r8*f(1,3,1,i,j,k)+6.0_r8*dy*f(1,4,1,i,j,k)
688  > +dx*(2.0_r8*f(2,3,1,i,j,k)+6.0_r8*dy*f(2,4,1,i,j,k)
689  > +dx*(2.0_r8*f(3,3,1,i,j,k)+6.0_r8*dy*f(3,4,1,i,j,k)
690  > +dx*(2.0_r8*f(4,3,1,i,j,k)+6.0_r8*dy*f(4,4,1,i,j,k))))
691  p2=
692  > 2.0_r8*f(1,3,2,i,j,k)+6.0_r8*dy*f(1,4,2,i,j,k)
693  > +dx*(2.0_r8*f(2,3,2,i,j,k)+6.0_r8*dy*f(2,4,2,i,j,k)
694  > +dx*(2.0_r8*f(3,3,2,i,j,k)+6.0_r8*dy*f(3,4,2,i,j,k)
695  > +dx*(2.0_r8*f(4,3,2,i,j,k)+6.0_r8*dy*f(4,4,2,i,j,k))))
696  p3=
697  > 2.0_r8*f(1,3,3,i,j,k)+6.0_r8*dy*f(1,4,3,i,j,k)
698  > +dx*(2.0_r8*f(2,3,3,i,j,k)+6.0_r8*dy*f(2,4,3,i,j,k)
699  > +dx*(2.0_r8*f(3,3,3,i,j,k)+6.0_r8*dy*f(3,4,3,i,j,k)
700  > +dx*(2.0_r8*f(4,3,3,i,j,k)+6.0_r8*dy*f(4,4,3,i,j,k))))
701  p4=
702  > 2.0_r8*f(1,3,4,i,j,k)+6.0_r8*dy*f(1,4,4,i,j,k)
703  > +dx*(2.0_r8*f(2,3,4,i,j,k)+6.0_r8*dy*f(2,4,4,i,j,k)
704  > +dx*(2.0_r8*f(3,3,4,i,j,k)+6.0_r8*dy*f(3,4,4,i,j,k)
705  > +dx*(2.0_r8*f(4,3,4,i,j,k)+6.0_r8*dy*f(4,4,4,i,j,k))))
706  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
707  enddo
708  endif
709 c
710  if((ict(7).gt.0).or.(ict(1).eq.-1)) then
711 c evaluate df2/dz2
712  iaval=iaval+1
713  do v=1,ivec
714  i=iv(v)
715  j=jv(v)
716  k=kv(v)
717  dx=dxv(v)
718  dy=dyv(v)
719  dz=dzv(v)
720  p3=
721  > f(1,1,3,i,j,k)+dy*(f(1,2,3,i,j,k)+dy*(f(1,3,3,i,j,k)+
722  > dy*f(1,4,3,i,j,k)))
723  > +dx*(f(2,1,3,i,j,k)+dy*(f(2,2,3,i,j,k)+dy*(f(2,3,3,i,j,k)+
724  > dy*f(2,4,3,i,j,k)))
725  > +dx*(f(3,1,3,i,j,k)+dy*(f(3,2,3,i,j,k)+dy*(f(3,3,3,i,j,k)+
726  > dy*f(3,4,3,i,j,k)))
727  > +dx*(f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
728  > dy*f(4,4,3,i,j,k)))
729  > )))
730  p4=
731  > f(1,1,4,i,j,k)+dy*(f(1,2,4,i,j,k)+dy*(f(1,3,4,i,j,k)+
732  > dy*f(1,4,4,i,j,k)))
733  > +dx*(f(2,1,4,i,j,k)+dy*(f(2,2,4,i,j,k)+dy*(f(2,3,4,i,j,k)+
734  > dy*f(2,4,4,i,j,k)))
735  > +dx*(f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
736  > dy*f(3,4,4,i,j,k)))
737  > +dx*(f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
738  > dy*f(4,4,4,i,j,k)))
739  > )))
740  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
741  enddo
742  endif
743 c
744  if((ict(8).gt.0).and.(ict(1).ne.-1)) then
745 c evaluate d2f/dxdy
746  iaval=iaval+1
747  do v=1,ivec
748  i=iv(v)
749  j=jv(v)
750  k=kv(v)
751  dx=dxv(v)
752  dy=dyv(v)
753  dz=dzv(v)
754  p1=
755  > f(2,2,1,i,j,k)+
756  > dy*(2.0_r8*f(2,3,1,i,j,k)+dy*3.0_r8*f(2,4,1,i,j,k))
757  > +2._r8*dx*(
758  > f(3,2,1,i,j,k)+
759  > dy*(2.0_r8*f(3,3,1,i,j,k)+dy*3.0_r8*f(3,4,1,i,j,k))
760  > +1.5_r8*dx*(
761  > f(4,2,1,i,j,k)+
762  > dy*(2.0_r8*f(4,3,1,i,j,k)+dy*3.0_r8*f(4,4,1,i,j,k))
763  > ))
764  p2=
765  > f(2,2,2,i,j,k)+
766  > dy*(2.0_r8*f(2,3,2,i,j,k)+dy*3.0_r8*f(2,4,2,i,j,k))
767  > +2._r8*dx*(
768  > f(3,2,2,i,j,k)+
769  > dy*(2.0_r8*f(3,3,2,i,j,k)+dy*3.0_r8*f(3,4,2,i,j,k))
770  > +1.5_r8*dx*(
771  > f(4,2,2,i,j,k)+
772  > dy*(2.0_r8*f(4,3,2,i,j,k)+dy*3.0_r8*f(4,4,2,i,j,k))
773  > ))
774  p3=
775  > f(2,2,3,i,j,k)+
776  > dy*(2.0_r8*f(2,3,3,i,j,k)+dy*3.0_r8*f(2,4,3,i,j,k))
777  > +2._r8*dx*(
778  > f(3,2,3,i,j,k)+
779  > dy*(2.0_r8*f(3,3,3,i,j,k)+dy*3.0_r8*f(3,4,3,i,j,k))
780  > +1.5_r8*dx*(
781  > f(4,2,3,i,j,k)+
782  > dy*(2.0_r8*f(4,3,3,i,j,k)+dy*3.0_r8*f(4,4,3,i,j,k))
783  > ))
784  p4=
785  > f(2,2,4,i,j,k)+
786  > dy*(2.0_r8*f(2,3,4,i,j,k)+dy*3.0_r8*f(2,4,4,i,j,k))
787  > +2._r8*dx*(
788  > f(3,2,4,i,j,k)+
789  > dy*(2.0_r8*f(3,3,4,i,j,k)+dy*3.0_r8*f(3,4,4,i,j,k))
790  > +1.5_r8*dx*(
791  > f(4,2,4,i,j,k)+
792  > dy*(2.0_r8*f(4,3,4,i,j,k)+dy*3.0_r8*f(4,4,4,i,j,k))
793  > ))
794  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
795  enddo
796  endif
797 c
798  if((ict(9).gt.0).and.(ict(1).ne.-1)) then
799 c evaluate d2f/dxdz
800  iaval=iaval+1
801  do v=1,ivec
802  i=iv(v)
803  j=jv(v)
804  k=kv(v)
805  dx=dxv(v)
806  dy=dyv(v)
807  dz=dzv(v)
808  p2=
809  > f(2,1,2,i,j,k)+dy*(f(2,2,2,i,j,k)+dy*(f(2,3,2,i,j,k)+
810  > dy*f(2,4,2,i,j,k)))
811  > +2._r8*dx*
812  > (f(3,1,2,i,j,k)+dy*(f(3,2,2,i,j,k)+dy*(f(3,3,2,i,j,k)+
813  > dy*f(3,4,2,i,j,k)))
814  > +1.5_r8*dx*
815  > (f(4,1,2,i,j,k)+dy*(f(4,2,2,i,j,k)+dy*(f(4,3,2,i,j,k)+
816  > dy*f(4,4,2,i,j,k)))
817  > ))
818  p3=
819  > f(2,1,3,i,j,k)+dy*(f(2,2,3,i,j,k)+dy*(f(2,3,3,i,j,k)+
820  > dy*f(2,4,3,i,j,k)))
821  > +2._r8*dx*
822  > (f(3,1,3,i,j,k)+dy*(f(3,2,3,i,j,k)+dy*(f(3,3,3,i,j,k)+
823  > dy*f(3,4,3,i,j,k)))
824  > +1.5_r8*dx*
825  > (f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
826  > dy*f(4,4,3,i,j,k)))
827  > ))
828  p4=
829  > f(2,1,4,i,j,k)+dy*(f(2,2,4,i,j,k)+dy*(f(2,3,4,i,j,k)+
830  > dy*f(2,4,4,i,j,k)))
831  > +2._r8*dx*
832  > (f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
833  > dy*f(3,4,4,i,j,k)))
834  > +1.5_r8*dx*
835  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
836  > dy*f(4,4,4,i,j,k)))
837  > ))
838  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
839  enddo
840  endif
841 c
842  if((ict(10).gt.0).and.(ict(1).ne.-1)) then
843 c evaluate d2f/dydz
844  iaval=iaval+1
845  do v=1,ivec
846  i=iv(v)
847  j=jv(v)
848  k=kv(v)
849  dx=dxv(v)
850  dy=dyv(v)
851  dz=dzv(v)
852  p2=
853  > f(1,2,2,i,j,k)+
854  > dy*(2.0_r8*f(1,3,2,i,j,k)+dy*3.0_r8*f(1,4,2,i,j,k))
855  > +dx*(
856  > f(2,2,2,i,j,k)+
857  > dy*(2.0_r8*f(2,3,2,i,j,k)+dy*3.0_r8*f(2,4,2,i,j,k))
858  > +dx*(
859  > f(3,2,2,i,j,k)+
860  > dy*(2.0_r8*f(3,3,2,i,j,k)+dy*3.0_r8*f(3,4,2,i,j,k))
861  > +dx*(
862  > f(4,2,2,i,j,k)+
863  > dy*(2.0_r8*f(4,3,2,i,j,k)+dy*3.0_r8*f(4,4,2,i,j,k))
864  > )))
865  p3=
866  > f(1,2,3,i,j,k)+
867  > dy*(2.0_r8*f(1,3,3,i,j,k)+dy*3.0_r8*f(1,4,3,i,j,k))
868  > +dx*(
869  > f(2,2,3,i,j,k)+
870  > dy*(2.0_r8*f(2,3,3,i,j,k)+dy*3.0_r8*f(2,4,3,i,j,k))
871  > +dx*(
872  > f(3,2,3,i,j,k)+
873  > dy*(2.0_r8*f(3,3,3,i,j,k)+dy*3.0_r8*f(3,4,3,i,j,k))
874  > +dx*(
875  > f(4,2,3,i,j,k)+
876  > dy*(2.0_r8*f(4,3,3,i,j,k)+dy*3.0_r8*f(4,4,3,i,j,k))
877  > )))
878  p4=
879  > f(1,2,4,i,j,k)+
880  > dy*(2.0_r8*f(1,3,4,i,j,k)+dy*3.0_r8*f(1,4,4,i,j,k))
881  > +dx*(
882  > f(2,2,4,i,j,k)+
883  > dy*(2.0_r8*f(2,3,4,i,j,k)+dy*3.0_r8*f(2,4,4,i,j,k))
884  > +dx*(
885  > f(3,2,4,i,j,k)+
886  > dy*(2.0_r8*f(3,3,4,i,j,k)+dy*3.0_r8*f(3,4,4,i,j,k))
887  > +dx*(
888  > f(4,2,4,i,j,k)+
889  > dy*(2.0_r8*f(4,3,4,i,j,k)+dy*3.0_r8*f(4,4,4,i,j,k))
890  > )))
891  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
892  enddo
893  endif
894 c
895  if(ict(1).eq.-1) then
896 c evaluate d4f/dx2dy2
897  iaval=iaval+1
898  do v=1,ivec
899  i=iv(v)
900  j=jv(v)
901  k=kv(v)
902  dx=dxv(v)
903  dy=dyv(v)
904  dz=dzv(v)
905  p1=
906  > 4.0_r8*f(3,3,1,i,j,k)+12.0_r8*dy*f(3,4,1,i,j,k)
907  > +dx*(12.0_r8*f(4,3,1,i,j,k)+36.0_r8*dy*f(4,4,1,i,j,k))
908  p2=
909  > 4.0_r8*f(3,3,2,i,j,k)+12.0_r8*dy*f(3,4,2,i,j,k)
910  > +dx*(12.0_r8*f(4,3,2,i,j,k)+36.0_r8*dy*f(4,4,2,i,j,k))
911  p3=
912  > 4.0_r8*f(3,3,3,i,j,k)+12.0_r8*dy*f(3,4,3,i,j,k)
913  > +dx*(12.0_r8*f(4,3,3,i,j,k)+36.0_r8*dy*f(4,4,3,i,j,k))
914  p4=
915  > 4.0_r8*f(3,3,4,i,j,k)+12.0_r8*dy*f(3,4,4,i,j,k)
916  > +dx*(12.0_r8*f(4,3,4,i,j,k)+36.0_r8*dy*f(4,4,4,i,j,k))
917  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
918  enddo
919  endif
920 c
921  if(ict(1).eq.-1) then
922 c evaluate d4f/dx2dz2
923  iaval=iaval+1
924  do v=1,ivec
925  i=iv(v)
926  j=jv(v)
927  k=kv(v)
928  dx=dxv(v)
929  dy=dyv(v)
930  dz=dzv(v)
931  p3= 2.0_r8*
932  > (f(3,1,3,i,j,k)+dy*(f(3,2,3,i,j,k)+dy*(f(3,3,3,i,j,k)+
933  > dy*f(3,4,3,i,j,k))))
934  > +6.0_r8*dx*
935  > (f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
936  > dy*f(4,4,3,i,j,k))))
937  p4= 2.0_r8*
938  > (f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
939  > dy*f(3,4,4,i,j,k))))
940  > +6.0_r8*dx*
941  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
942  > dy*f(4,4,4,i,j,k))))
943  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
944  enddo
945  endif
946 c
947  if(ict(1).eq.-1) then
948 c evaluate d4f/dy2dz2
949  iaval=iaval+1
950  do v=1,ivec
951  i=iv(v)
952  j=jv(v)
953  k=kv(v)
954  dx=dxv(v)
955  dy=dyv(v)
956  dz=dzv(v)
957  p3=
958  > 2.0_r8*f(1,3,3,i,j,k)+6.0_r8*dy*f(1,4,3,i,j,k)
959  > +dx*(2.0_r8*f(2,3,3,i,j,k)+6.0_r8*dy*f(2,4,3,i,j,k)
960  > +dx*(2.0_r8*f(3,3,3,i,j,k)+6.0_r8*dy*f(3,4,3,i,j,k)
961  > +dx*(2.0_r8*f(4,3,3,i,j,k)+6.0_r8*dy*f(4,4,3,i,j,k))))
962  p4=
963  > 2.0_r8*f(1,3,4,i,j,k)+6.0_r8*dy*f(1,4,4,i,j,k)
964  > +dx*(2.0_r8*f(2,3,4,i,j,k)+6.0_r8*dy*f(2,4,4,i,j,k)
965  > +dx*(2.0_r8*f(3,3,4,i,j,k)+6.0_r8*dy*f(3,4,4,i,j,k)
966  > +dx*(2.0_r8*f(4,3,4,i,j,k)+6.0_r8*dy*f(4,4,4,i,j,k))))
967  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
968  enddo
969  endif
970 c
971  if(ict(1).eq.-1) then
972 c evaluate d6f/dx2dy2dz2
973  iaval=iaval+1
974  do v=1,ivec
975  i=iv(v)
976  j=jv(v)
977  k=kv(v)
978  dx=dxv(v)
979  dy=dyv(v)
980  dz=dzv(v)
981  p3=
982  > 4.0_r8*f(3,3,3,i,j,k)+12.0_r8*dy*f(3,4,3,i,j,k)
983  > +dx*(12.0_r8*f(4,3,3,i,j,k)+36.0_r8*dy*f(4,4,3,i,j,k))
984  p4=
985  > 4.0_r8*f(3,3,4,i,j,k)+12.0_r8*dy*f(3,4,4,i,j,k)
986  > +dx*(12.0_r8*f(4,3,4,i,j,k)+36.0_r8*dy*f(4,4,4,i,j,k))
987  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
988  enddo
989  endif
990 c
991  endif
992 C
993 C----------------------------------
994 C 3rd derivatives (.le.2 in each coordinate)
995 C
996  if(ict(1).eq.3) then
997  if(ict(2).eq.1) then
998 C ! fxxy
999  iaval=iaval+1
1000  do v=1,ivec
1001  i=iv(v)
1002  j=jv(v)
1003  k=kv(v)
1004  dx=dxv(v)
1005  dy=dyv(v)
1006  dz=dzv(v)
1007  p1= 2.0_r8*
1008  > (f(3,2,1,i,j,k)+2.0_r8*dy*(f(3,3,1,i,j,k)+
1009  > 1.5_r8*dy*f(3,4,1,i,j,k)))
1010  > +6.0_r8*dx*
1011  > (f(4,2,1,i,j,k)+2.0_r8*dy*(f(4,3,1,i,j,k)+
1012  > 1.5_r8*dy*f(4,4,1,i,j,k)))
1013  p2= 2.0_r8*
1014  > (f(3,2,2,i,j,k)+2.0_r8*dy*(f(3,3,2,i,j,k)+
1015  > 1.5_r8*dy*f(3,4,2,i,j,k)))
1016  > +6.0_r8*dx*
1017  > (f(4,2,2,i,j,k)+2.0_r8*dy*(f(4,3,2,i,j,k)+
1018  > 1.5_r8*dy*f(4,4,2,i,j,k)))
1019  p3= 2.0_r8*
1020  > (f(3,2,3,i,j,k)+2.0_r8*dy*(f(3,3,3,i,j,k)+
1021  > 1.5_r8*dy*f(3,4,3,i,j,k)))
1022  > +6.0_r8*dx*
1023  > (f(4,2,3,i,j,k)+2.0_r8*dy*(f(4,3,3,i,j,k)+
1024  > 1.5_r8*dy*f(4,4,3,i,j,k)))
1025  p4= 2.0_r8*
1026  > (f(3,2,4,i,j,k)+2.0_r8*dy*(f(3,3,4,i,j,k)+
1027  > 1.5_r8*dy*f(3,4,4,i,j,k)))
1028  > +6.0_r8*dx*
1029  > (f(4,2,4,i,j,k)+2.0_r8*dy*(f(4,3,4,i,j,k)+
1030  > 1.5_r8*dy*f(4,4,4,i,j,k)))
1031 C
1032  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
1033  enddo
1034  endif
1035 C
1036  if(ict(3).eq.1) then
1037 C ! fxxz
1038  iaval=iaval+1
1039  do v=1,ivec
1040  i=iv(v)
1041  j=jv(v)
1042  k=kv(v)
1043  dx=dxv(v)
1044  dy=dyv(v)
1045  dz=dzv(v)
1046 C
1047  p2= 2.0_r8*
1048  > (f(3,1,2,i,j,k)+dy*(f(3,2,2,i,j,k)+dy*(f(3,3,2,i,j,k)+
1049  > dy*f(3,4,2,i,j,k))))
1050  > +6.0_r8*dx*
1051  > (f(4,1,2,i,j,k)+dy*(f(4,2,2,i,j,k)+dy*(f(4,3,2,i,j,k)+
1052  > dy*f(4,4,2,i,j,k))))
1053  p3= 2.0_r8*
1054  > (f(3,1,3,i,j,k)+dy*(f(3,2,3,i,j,k)+dy*(f(3,3,3,i,j,k)+
1055  > dy*f(3,4,3,i,j,k))))
1056  > +6.0_r8*dx*
1057  > (f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
1058  > dy*f(4,4,3,i,j,k))))
1059  p4= 2.0_r8*
1060  > (f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
1061  > dy*f(3,4,4,i,j,k))))
1062  > +6.0_r8*dx*
1063  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
1064  > dy*f(4,4,4,i,j,k))))
1065  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
1066  enddo
1067  endif
1068 C
1069  if(ict(4).eq.1) then
1070 C ! fxyy
1071  iaval=iaval+1
1072  do v=1,ivec
1073  i=iv(v)
1074  j=jv(v)
1075  k=kv(v)
1076  dx=dxv(v)
1077  dy=dyv(v)
1078  dz=dzv(v)
1079 C
1080  p1=
1081  > 2.0_r8*f(2,3,1,i,j,k)+6.0_r8*dy*f(2,4,1,i,j,k)
1082  > +2.0_r8*dx*(
1083  > 2.0_r8*f(3,3,1,i,j,k)+6.0_r8*dy*f(3,4,1,i,j,k)
1084  > +1.5_r8*dx*(
1085  > 2.0_r8*f(4,3,1,i,j,k)+6.0_r8*dy*f(4,4,1,i,j,k)))
1086  p2=
1087  > 2.0_r8*f(2,3,2,i,j,k)+6.0_r8*dy*f(2,4,2,i,j,k)
1088  > +2.0_r8*dx*(
1089  > 2.0_r8*f(3,3,2,i,j,k)+6.0_r8*dy*f(3,4,2,i,j,k)
1090  > +1.5_r8*dx*(
1091  > 2.0_r8*f(4,3,2,i,j,k)+6.0_r8*dy*f(4,4,2,i,j,k)))
1092  p3=
1093  > 2.0_r8*f(2,3,3,i,j,k)+6.0_r8*dy*f(2,4,3,i,j,k)
1094  > +2.0_r8*dx*(
1095  > 2.0_r8*f(3,3,3,i,j,k)+6.0_r8*dy*f(3,4,3,i,j,k)
1096  > +1.5_r8*dx*(
1097  > 2.0_r8*f(4,3,3,i,j,k)+6.0_r8*dy*f(4,4,3,i,j,k)))
1098  p4=
1099  > 2.0_r8*f(2,3,4,i,j,k)+6.0_r8*dy*f(2,4,4,i,j,k)
1100  > +2.0_r8*dx*(
1101  > 2.0_r8*f(3,3,4,i,j,k)+6.0_r8*dy*f(3,4,4,i,j,k)
1102  > +1.5_r8*dx*(
1103  > 2.0_r8*f(4,3,4,i,j,k)+6.0_r8*dy*f(4,4,4,i,j,k)))
1104  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
1105  enddo
1106  endif
1107 C
1108  if(ict(5).eq.1) then
1109 C ! fxyz
1110  iaval=iaval+1
1111  do v=1,ivec
1112  i=iv(v)
1113  j=jv(v)
1114  k=kv(v)
1115  dx=dxv(v)
1116  dy=dyv(v)
1117  dz=dzv(v)
1118 C
1119  p2=
1120  > f(2,2,2,i,j,k)+
1121  > dy*(2.0_r8*f(2,3,2,i,j,k)+dy*3.0_r8*f(2,4,2,i,j,k))
1122  > +2._r8*dx*(
1123  > f(3,2,2,i,j,k)+
1124  > dy*(2.0_r8*f(3,3,2,i,j,k)+dy*3.0_r8*f(3,4,2,i,j,k))
1125  > +1.5_r8*dx*(
1126  > f(4,2,2,i,j,k)+
1127  > dy*(2.0_r8*f(4,3,2,i,j,k)+dy*3.0_r8*f(4,4,2,i,j,k))
1128  > ))
1129  p3=
1130  > f(2,2,3,i,j,k)+
1131  > dy*(2.0_r8*f(2,3,3,i,j,k)+dy*3.0_r8*f(2,4,3,i,j,k))
1132  > +2._r8*dx*(
1133  > f(3,2,3,i,j,k)+
1134  > dy*(2.0_r8*f(3,3,3,i,j,k)+dy*3.0_r8*f(3,4,3,i,j,k))
1135  > +1.5_r8*dx*(
1136  > f(4,2,3,i,j,k)+
1137  > dy*(2.0_r8*f(4,3,3,i,j,k)+dy*3.0_r8*f(4,4,3,i,j,k))
1138  > ))
1139  p4=
1140  > f(2,2,4,i,j,k)+
1141  > dy*(2.0_r8*f(2,3,4,i,j,k)+dy*3.0_r8*f(2,4,4,i,j,k))
1142  > +2._r8*dx*(
1143  > f(3,2,4,i,j,k)+
1144  > dy*(2.0_r8*f(3,3,4,i,j,k)+dy*3.0_r8*f(3,4,4,i,j,k))
1145  > +1.5_r8*dx*(
1146  > f(4,2,4,i,j,k)+
1147  > dy*(2.0_r8*f(4,3,4,i,j,k)+dy*3.0_r8*f(4,4,4,i,j,k))
1148  > ))
1149  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
1150  enddo
1151  endif
1152 C
1153  if(ict(6).eq.1) then
1154 C ! fxzz
1155  iaval=iaval+1
1156  do v=1,ivec
1157  i=iv(v)
1158  j=jv(v)
1159  k=kv(v)
1160  dx=dxv(v)
1161  dy=dyv(v)
1162  dz=dzv(v)
1163 C
1164  p3=
1165  > f(2,1,3,i,j,k)+dy*(f(2,2,3,i,j,k)+dy*(f(2,3,3,i,j,k)+
1166  > dy*f(2,4,3,i,j,k)))
1167  > +2._r8*dx*
1168  > (f(3,1,3,i,j,k)+dy*(f(3,2,3,i,j,k)+dy*(f(3,3,3,i,j,k)+
1169  > dy*f(3,4,3,i,j,k)))
1170  > +1.5_r8*dx*
1171  > (f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
1172  > dy*f(4,4,3,i,j,k)))
1173  > ))
1174  p4=
1175  > f(2,1,4,i,j,k)+dy*(f(2,2,4,i,j,k)+dy*(f(2,3,4,i,j,k)+
1176  > dy*f(2,4,4,i,j,k)))
1177  > +2._r8*dx*
1178  > (f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
1179  > dy*f(3,4,4,i,j,k)))
1180  > +1.5_r8*dx*
1181  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
1182  > dy*f(4,4,4,i,j,k)))
1183  > ))
1184  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
1185  enddo
1186  endif
1187 C
1188  if(ict(7).eq.1) then
1189 C ! fyyz
1190  iaval=iaval+1
1191  do v=1,ivec
1192  i=iv(v)
1193  j=jv(v)
1194  k=kv(v)
1195  dx=dxv(v)
1196  dy=dyv(v)
1197  dz=dzv(v)
1198 C
1199  p2=
1200  > 2.0_r8*f(1,3,2,i,j,k)+6.0_r8*dy*f(1,4,2,i,j,k)
1201  > +dx*(2.0_r8*f(2,3,2,i,j,k)+6.0_r8*dy*f(2,4,2,i,j,k)
1202  > +dx*(2.0_r8*f(3,3,2,i,j,k)+6.0_r8*dy*f(3,4,2,i,j,k)
1203  > +dx*(2.0_r8*f(4,3,2,i,j,k)+6.0_r8*dy*f(4,4,2,i,j,k))))
1204  p3=
1205  > 2.0_r8*f(1,3,3,i,j,k)+6.0_r8*dy*f(1,4,3,i,j,k)
1206  > +dx*(2.0_r8*f(2,3,3,i,j,k)+6.0_r8*dy*f(2,4,3,i,j,k)
1207  > +dx*(2.0_r8*f(3,3,3,i,j,k)+6.0_r8*dy*f(3,4,3,i,j,k)
1208  > +dx*(2.0_r8*f(4,3,3,i,j,k)+6.0_r8*dy*f(4,4,3,i,j,k))))
1209  p4=
1210  > 2.0_r8*f(1,3,4,i,j,k)+6.0_r8*dy*f(1,4,4,i,j,k)
1211  > +dx*(2.0_r8*f(2,3,4,i,j,k)+6.0_r8*dy*f(2,4,4,i,j,k)
1212  > +dx*(2.0_r8*f(3,3,4,i,j,k)+6.0_r8*dy*f(3,4,4,i,j,k)
1213  > +dx*(2.0_r8*f(4,3,4,i,j,k)+6.0_r8*dy*f(4,4,4,i,j,k))))
1214  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
1215  enddo
1216  endif
1217 C
1218  if(ict(8).eq.1) then
1219 C ! fyzz
1220  iaval=iaval+1
1221  do v=1,ivec
1222  i=iv(v)
1223  j=jv(v)
1224  k=kv(v)
1225  dx=dxv(v)
1226  dy=dyv(v)
1227  dz=dzv(v)
1228 C
1229  p3=
1230  > f(1,2,3,i,j,k)+
1231  > dy*(2.0_r8*f(1,3,3,i,j,k)+dy*3.0_r8*f(1,4,3,i,j,k))
1232  > +dx*(
1233  > f(2,2,3,i,j,k)+
1234  > dy*(2.0_r8*f(2,3,3,i,j,k)+dy*3.0_r8*f(2,4,3,i,j,k))
1235  > +dx*(
1236  > f(3,2,3,i,j,k)+
1237  > dy*(2.0_r8*f(3,3,3,i,j,k)+dy*3.0_r8*f(3,4,3,i,j,k))
1238  > +dx*(
1239  > f(4,2,3,i,j,k)+
1240  > dy*(2.0_r8*f(4,3,3,i,j,k)+dy*3.0_r8*f(4,4,3,i,j,k))
1241  > )))
1242  p4=
1243  > f(1,2,4,i,j,k)+
1244  > dy*(2.0_r8*f(1,3,4,i,j,k)+dy*3.0_r8*f(1,4,4,i,j,k))
1245  > +dx*(
1246  > f(2,2,4,i,j,k)+
1247  > dy*(2.0_r8*f(2,3,4,i,j,k)+dy*3.0_r8*f(2,4,4,i,j,k))
1248  > +dx*(
1249  > f(3,2,4,i,j,k)+
1250  > dy*(2.0_r8*f(3,3,4,i,j,k)+dy*3.0_r8*f(3,4,4,i,j,k))
1251  > +dx*(
1252  > f(4,2,4,i,j,k)+
1253  > dy*(2.0_r8*f(4,3,4,i,j,k)+dy*3.0_r8*f(4,4,4,i,j,k))
1254  > )))
1255  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
1256  enddo
1257  endif
1258 C
1259  endif
1260 C
1261 C----------------------------------
1262 C 3rd derivatives (3 in each coordinate)
1263 C
1264  if(ict(1).eq.-3) then
1265  if(ict(2).eq.1) then
1266 C ! fxxx
1267  iaval=iaval+1
1268  do v=1,ivec
1269  i=iv(v)
1270  j=jv(v)
1271  k=kv(v)
1272  dx=dxv(v)
1273  dy=dyv(v)
1274  dz=dzv(v)
1275 C
1276  p1= 6.0_r8*
1277  > (f(4,1,1,i,j,k)+dy*(f(4,2,1,i,j,k)+dy*(f(4,3,1,i,j,k)+
1278  > dy*f(4,4,1,i,j,k))))
1279  p2= 6.0_r8*
1280  > (f(4,1,2,i,j,k)+dy*(f(4,2,2,i,j,k)+dy*(f(4,3,2,i,j,k)+
1281  > dy*f(4,4,2,i,j,k))))
1282  p3= 6.0_r8*
1283  > (f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
1284  > dy*f(4,4,3,i,j,k))))
1285  p4= 6.0_r8*
1286  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
1287  > dy*f(4,4,4,i,j,k))))
1288  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
1289  enddo
1290  endif
1291 C
1292  if(ict(3).eq.1) then
1293 C ! fyyy
1294  iaval=iaval+1
1295  do v=1,ivec
1296  i=iv(v)
1297  j=jv(v)
1298  k=kv(v)
1299  dx=dxv(v)
1300  dy=dyv(v)
1301  dz=dzv(v)
1302 C
1303  p1=
1304  > 6.0_r8*(f(1,4,1,i,j,k)
1305  > +dx*(f(2,4,1,i,j,k)
1306  > +dx*(f(3,4,1,i,j,k)
1307  > +dx*f(4,4,1,i,j,k))))
1308  p2=
1309  > 6.0_r8*(f(1,4,2,i,j,k)
1310  > +dx*(f(2,4,2,i,j,k)
1311  > +dx*(f(3,4,2,i,j,k)
1312  > +dx*f(4,4,2,i,j,k))))
1313  p3=
1314  > 6.0_r8*(f(1,4,3,i,j,k)
1315  > +dx*(f(2,4,3,i,j,k)
1316  > +dx*(f(3,4,3,i,j,k)
1317  > +dx*f(4,4,3,i,j,k))))
1318  p4=
1319  > 6.0_r8*(f(1,4,4,i,j,k)
1320  > +dx*(f(2,4,4,i,j,k)
1321  > +dx*(f(3,4,4,i,j,k)
1322  > +dx*f(4,4,4,i,j,k))))
1323  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
1324  enddo
1325  endif
1326 C
1327  if(ict(4).eq.1) then
1328 C ! fzzz
1329  iaval=iaval+1
1330  do v=1,ivec
1331  i=iv(v)
1332  j=jv(v)
1333  k=kv(v)
1334  dx=dxv(v)
1335  dy=dyv(v)
1336  dz=dzv(v)
1337 C
1338  p4=
1339  > f(1,1,4,i,j,k)+dy*(f(1,2,4,i,j,k)+dy*(f(1,3,4,i,j,k)+
1340  > dy*f(1,4,4,i,j,k)))
1341  > +dx*(f(2,1,4,i,j,k)+dy*(f(2,2,4,i,j,k)+dy*(f(2,3,4,i,j,k)+
1342  > dy*f(2,4,4,i,j,k)))
1343  > +dx*(f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
1344  > dy*f(3,4,4,i,j,k)))
1345  > +dx*(f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
1346  > dy*f(4,4,4,i,j,k)))
1347  > )))
1348  fval(v,iaval)=6.0_r8*p4
1349  enddo
1350  endif
1351 C
1352  endif
1353 C
1354 C----------------------------------
1355 C 4th derivatives (.le.2 in each coordinate)
1356 C
1357  if(ict(1).eq.4) then
1358  if(ict(2).eq.1) then
1359 C ! fxxyy
1360  iaval=iaval+1
1361  do v=1,ivec
1362  i=iv(v)
1363  j=jv(v)
1364  k=kv(v)
1365  dx=dxv(v)
1366  dy=dyv(v)
1367  dz=dzv(v)
1368 C
1369  p1=
1370  > 4.0_r8*f(3,3,1,i,j,k)+12.0_r8*dy*f(3,4,1,i,j,k)
1371  > +12.0_r8*dx*(
1372  > f(4,3,1,i,j,k)+3.0_r8*dy*f(4,4,1,i,j,k))
1373  p2=
1374  > 4.0_r8*f(3,3,2,i,j,k)+12.0_r8*dy*f(3,4,2,i,j,k)
1375  > +12.0_r8*dx*(
1376  > f(4,3,2,i,j,k)+3.0_r8*dy*f(4,4,2,i,j,k))
1377  p3=
1378  > 4.0_r8*f(3,3,3,i,j,k)+12.0_r8*dy*f(3,4,3,i,j,k)
1379  > +12.0_r8*dx*(
1380  > f(4,3,3,i,j,k)+3.0_r8*dy*f(4,4,3,i,j,k))
1381  p4=
1382  > 4.0_r8*f(3,3,4,i,j,k)+12.0_r8*dy*f(3,4,4,i,j,k)
1383  > +12.0_r8*dx*(
1384  > f(4,3,4,i,j,k)+3.0_r8*dy*f(4,4,4,i,j,k))
1385 C
1386  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
1387  enddo
1388  endif
1389 C
1390  if(ict(3).eq.1) then
1391 C ! fxxyz
1392  iaval=iaval+1
1393  do v=1,ivec
1394  i=iv(v)
1395  j=jv(v)
1396  k=kv(v)
1397  dx=dxv(v)
1398  dy=dyv(v)
1399  dz=dzv(v)
1400 C
1401  p2= 2.0_r8*
1402  > (f(3,2,2,i,j,k)+2.0_r8*dy*(f(3,3,2,i,j,k)+
1403  > 1.5_r8*dy*f(3,4,2,i,j,k)))
1404  > +6.0_r8*dx*
1405  > (f(4,2,2,i,j,k)+2.0_r8*dy*(f(4,3,2,i,j,k)+
1406  > 1.5_r8*dy*f(4,4,2,i,j,k)))
1407  p3= 2.0_r8*
1408  > (f(3,2,3,i,j,k)+2.0_r8*dy*(f(3,3,3,i,j,k)+
1409  > 1.5_r8*dy*f(3,4,3,i,j,k)))
1410  > +6.0_r8*dx*
1411  > (f(4,2,3,i,j,k)+2.0_r8*dy*(f(4,3,3,i,j,k)+
1412  > 1.5_r8*dy*f(4,4,3,i,j,k)))
1413  p4= 2.0_r8*
1414  > (f(3,2,4,i,j,k)+2.0_r8*dy*(f(3,3,4,i,j,k)+
1415  > 1.5_r8*dy*f(3,4,4,i,j,k)))
1416  > +6.0_r8*dx*
1417  > (f(4,2,4,i,j,k)+2.0_r8*dy*(f(4,3,4,i,j,k)+
1418  > 1.5_r8*dy*f(4,4,4,i,j,k)))
1419 C
1420  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
1421  enddo
1422  endif
1423 C
1424  if(ict(4).eq.1) then
1425 C ! fxxzz
1426  iaval=iaval+1
1427  do v=1,ivec
1428  i=iv(v)
1429  j=jv(v)
1430  k=kv(v)
1431  dx=dxv(v)
1432  dy=dyv(v)
1433  dz=dzv(v)
1434 C
1435  p3= 2.0_r8*
1436  > (f(3,1,3,i,j,k)+dy*(f(3,2,3,i,j,k)+dy*(f(3,3,3,i,j,k)+
1437  > dy*f(3,4,3,i,j,k))))
1438  > +6.0_r8*dx*
1439  > (f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
1440  > dy*f(4,4,3,i,j,k))))
1441  p4= 2.0_r8*
1442  > (f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
1443  > dy*f(3,4,4,i,j,k))))
1444  > +6.0_r8*dx*
1445  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
1446  > dy*f(4,4,4,i,j,k))))
1447 C
1448  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
1449  enddo
1450  endif
1451 C
1452  if(ict(5).eq.1) then
1453 C ! fxyyz
1454  iaval=iaval+1
1455  do v=1,ivec
1456  i=iv(v)
1457  j=jv(v)
1458  k=kv(v)
1459  dx=dxv(v)
1460  dy=dyv(v)
1461  dz=dzv(v)
1462 C
1463  p2=
1464  > 2.0_r8*f(2,3,2,i,j,k)+6.0_r8*dy*f(2,4,2,i,j,k)
1465  > +2.0_r8*dx*(
1466  > 2.0_r8*f(3,3,2,i,j,k)+6.0_r8*dy*f(3,4,2,i,j,k)
1467  > +1.5_r8*dx*(
1468  > 2.0_r8*f(4,3,2,i,j,k)+6.0_r8*dy*f(4,4,2,i,j,k)))
1469  p3=
1470  > 2.0_r8*f(2,3,3,i,j,k)+6.0_r8*dy*f(2,4,3,i,j,k)
1471  > +2.0_r8*dx*(
1472  > 2.0_r8*f(3,3,3,i,j,k)+6.0_r8*dy*f(3,4,3,i,j,k)
1473  > +1.5_r8*dx*(
1474  > 2.0_r8*f(4,3,3,i,j,k)+6.0_r8*dy*f(4,4,3,i,j,k)))
1475  p4=
1476  > 2.0_r8*f(2,3,4,i,j,k)+6.0_r8*dy*f(2,4,4,i,j,k)
1477  > +2.0_r8*dx*(
1478  > 2.0_r8*f(3,3,4,i,j,k)+6.0_r8*dy*f(3,4,4,i,j,k)
1479  > +1.5_r8*dx*(
1480  > 2.0_r8*f(4,3,4,i,j,k)+6.0_r8*dy*f(4,4,4,i,j,k)))
1481 C
1482  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
1483  enddo
1484  endif
1485 C
1486  if(ict(6).eq.1) then
1487 C ! fxyzz
1488  iaval=iaval+1
1489  do v=1,ivec
1490  i=iv(v)
1491  j=jv(v)
1492  k=kv(v)
1493  dx=dxv(v)
1494  dy=dyv(v)
1495  dz=dzv(v)
1496 C
1497  p3=
1498  > f(2,2,3,i,j,k)+
1499  > dy*(2.0_r8*f(2,3,3,i,j,k)+dy*3.0_r8*f(2,4,3,i,j,k))
1500  > +2._r8*dx*(
1501  > f(3,2,3,i,j,k)+
1502  > dy*(2.0_r8*f(3,3,3,i,j,k)+dy*3.0_r8*f(3,4,3,i,j,k))
1503  > +1.5_r8*dx*(
1504  > f(4,2,3,i,j,k)+
1505  > dy*(2.0_r8*f(4,3,3,i,j,k)+dy*3.0_r8*f(4,4,3,i,j,k))
1506  > ))
1507  p4=
1508  > f(2,2,4,i,j,k)+
1509  > dy*(2.0_r8*f(2,3,4,i,j,k)+dy*3.0_r8*f(2,4,4,i,j,k))
1510  > +2._r8*dx*(
1511  > f(3,2,4,i,j,k)+
1512  > dy*(2.0_r8*f(3,3,4,i,j,k)+dy*3.0_r8*f(3,4,4,i,j,k))
1513  > +1.5_r8*dx*(
1514  > f(4,2,4,i,j,k)+
1515  > dy*(2.0_r8*f(4,3,4,i,j,k)+dy*3.0_r8*f(4,4,4,i,j,k))
1516  > ))
1517  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
1518  enddo
1519  endif
1520 C
1521  if(ict(7).eq.1) then
1522 C ! fyyzz
1523  iaval=iaval+1
1524  do v=1,ivec
1525  i=iv(v)
1526  j=jv(v)
1527  k=kv(v)
1528  dx=dxv(v)
1529  dy=dyv(v)
1530  dz=dzv(v)
1531 C
1532  p3=
1533  > 2.0_r8*f(1,3,3,i,j,k)+6.0_r8*dy*f(1,4,3,i,j,k)
1534  > +dx*(2.0_r8*f(2,3,3,i,j,k)+6.0_r8*dy*f(2,4,3,i,j,k)
1535  > +dx*(2.0_r8*f(3,3,3,i,j,k)+6.0_r8*dy*f(3,4,3,i,j,k)
1536  > +dx*(2.0_r8*f(4,3,3,i,j,k)+6.0_r8*dy*f(4,4,3,i,j,k))))
1537  p4=
1538  > 2.0_r8*f(1,3,4,i,j,k)+6.0_r8*dy*f(1,4,4,i,j,k)
1539  > +dx*(2.0_r8*f(2,3,4,i,j,k)+6.0_r8*dy*f(2,4,4,i,j,k)
1540  > +dx*(2.0_r8*f(3,3,4,i,j,k)+6.0_r8*dy*f(3,4,4,i,j,k)
1541  > +dx*(2.0_r8*f(4,3,4,i,j,k)+6.0_r8*dy*f(4,4,4,i,j,k))))
1542 C
1543  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
1544  enddo
1545  endif
1546 C
1547  endif
1548 C
1549 C
1550 C----------------------------------
1551 C 4th derivatives (3 in a coordinate)
1552 C
1553  if(ict(1).eq.-4) then
1554  if(ict(2).eq.1) then
1555 C ! fxxxy
1556  iaval=iaval+1
1557  do v=1,ivec
1558  i=iv(v)
1559  j=jv(v)
1560  k=kv(v)
1561  dx=dxv(v)
1562  dy=dyv(v)
1563  dz=dzv(v)
1564 C
1565  p1= 6.0_r8*
1566  > (f(4,2,1,i,j,k)+2.0_r8*dy*(f(4,3,1,i,j,k)+
1567  > 1.5_r8*dy*f(4,4,1,i,j,k)))
1568  p2= 6.0_r8*
1569  > (f(4,2,2,i,j,k)+2.0_r8*dy*(f(4,3,2,i,j,k)+
1570  > 1.5_r8*dy*f(4,4,2,i,j,k)))
1571  p3= 6.0_r8*
1572  > (f(4,2,3,i,j,k)+2.0_r8*dy*(f(4,3,3,i,j,k)+
1573  > 1.5_r8*dy*f(4,4,3,i,j,k)))
1574  p4= 6.0_r8*
1575  > (f(4,2,4,i,j,k)+2.0_r8*dy*(f(4,3,4,i,j,k)+
1576  > 1.5_r8*dy*f(4,4,4,i,j,k)))
1577 C
1578  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
1579  enddo
1580  endif
1581 C
1582  if(ict(3).eq.1) then
1583 C ! fxxxz
1584  iaval=iaval+1
1585  do v=1,ivec
1586  i=iv(v)
1587  j=jv(v)
1588  k=kv(v)
1589  dx=dxv(v)
1590  dy=dyv(v)
1591  dz=dzv(v)
1592 C
1593  p2= 6.0_r8*
1594  > (f(4,1,2,i,j,k)+dy*(f(4,2,2,i,j,k)+dy*(f(4,3,2,i,j,k)+
1595  > dy*f(4,4,2,i,j,k))))
1596  p3= 6.0_r8*
1597  > (f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
1598  > dy*f(4,4,3,i,j,k))))
1599  p4= 6.0_r8*
1600  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
1601  > dy*f(4,4,4,i,j,k))))
1602 C
1603  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
1604  enddo
1605  endif
1606 C
1607  if(ict(4).eq.1) then
1608 C ! fxyyy
1609  iaval=iaval+1
1610  do v=1,ivec
1611  i=iv(v)
1612  j=jv(v)
1613  k=kv(v)
1614  dx=dxv(v)
1615  dy=dyv(v)
1616  dz=dzv(v)
1617 C
1618  p1=
1619  > 6.0_r8*(f(2,4,1,i,j,k)
1620  > +2.0_r8*dx*(f(3,4,1,i,j,k)
1621  > +1.5_r8*dx*f(4,4,1,i,j,k)))
1622  p2=
1623  > 6.0_r8*(f(2,4,2,i,j,k)
1624  > +2.0_r8*dx*(f(3,4,2,i,j,k)
1625  > +1.5_r8*dx*f(4,4,2,i,j,k)))
1626  p3=
1627  > 6.0_r8*(f(2,4,3,i,j,k)
1628  > +2.0_r8*dx*(f(3,4,3,i,j,k)
1629  > +1.5_r8*dx*f(4,4,3,i,j,k)))
1630  p4=
1631  > 6.0_r8*(f(2,4,4,i,j,k)
1632  > +2.0_r8*dx*(f(3,4,4,i,j,k)
1633  > +1.5_r8*dx*f(4,4,4,i,j,k)))
1634 C
1635  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
1636  enddo
1637  endif
1638 C
1639  if(ict(5).eq.1) then
1640 C ! fxzzz
1641  iaval=iaval+1
1642  do v=1,ivec
1643  i=iv(v)
1644  j=jv(v)
1645  k=kv(v)
1646  dx=dxv(v)
1647  dy=dyv(v)
1648  dz=dzv(v)
1649 C
1650  p4=
1651  > f(2,1,4,i,j,k)+dy*(f(2,2,4,i,j,k)+dy*(f(2,3,4,i,j,k)+
1652  > dy*f(2,4,4,i,j,k)))
1653  > +2._r8*dx*
1654  > (f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
1655  > dy*f(3,4,4,i,j,k)))
1656  > +1.5_r8*dx*
1657  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
1658  > dy*f(4,4,4,i,j,k)))
1659  > ))
1660  fval(v,iaval)=6.0_r8*p4
1661  enddo
1662  endif
1663 C
1664  if(ict(6).eq.1) then
1665 C ! fyyyz
1666  iaval=iaval+1
1667  do v=1,ivec
1668  i=iv(v)
1669  j=jv(v)
1670  k=kv(v)
1671  dx=dxv(v)
1672  dy=dyv(v)
1673  dz=dzv(v)
1674 C
1675  p2=
1676  > 6.0_r8*(f(1,4,2,i,j,k)
1677  > +dx*(f(2,4,2,i,j,k)
1678  > +dx*(f(3,4,2,i,j,k)
1679  > +dx*f(4,4,2,i,j,k))))
1680  p3=
1681  > 6.0_r8*(f(1,4,3,i,j,k)
1682  > +dx*(f(2,4,3,i,j,k)
1683  > +dx*(f(3,4,3,i,j,k)
1684  > +dx*f(4,4,3,i,j,k))))
1685  p4=
1686  > 6.0_r8*(f(1,4,4,i,j,k)
1687  > +dx*(f(2,4,4,i,j,k)
1688  > +dx*(f(3,4,4,i,j,k)
1689  > +dx*f(4,4,4,i,j,k))))
1690  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
1691 C
1692  enddo
1693  endif
1694 C
1695  if(ict(7).eq.1) then
1696 C ! fyzzz
1697  iaval=iaval+1
1698  do v=1,ivec
1699  i=iv(v)
1700  j=jv(v)
1701  k=kv(v)
1702  dx=dxv(v)
1703  dy=dyv(v)
1704  dz=dzv(v)
1705 C
1706  p4=
1707  > f(1,2,4,i,j,k)+
1708  > dy*(2.0_r8*f(1,3,4,i,j,k)+dy*3.0_r8*f(1,4,4,i,j,k))
1709  > +dx*(
1710  > f(2,2,4,i,j,k)+
1711  > dy*(2.0_r8*f(2,3,4,i,j,k)+dy*3.0_r8*f(2,4,4,i,j,k))
1712  > +dx*(
1713  > f(3,2,4,i,j,k)+
1714  > dy*(2.0_r8*f(3,3,4,i,j,k)+dy*3.0_r8*f(3,4,4,i,j,k))
1715  > +dx*(
1716  > f(4,2,4,i,j,k)+
1717  > dy*(2.0_r8*f(4,3,4,i,j,k)+dy*3.0_r8*f(4,4,4,i,j,k))
1718  > )))
1719  fval(v,iaval)=6.0_r8*p4
1720  enddo
1721  endif
1722 C
1723  endif
1724 C
1725 C----------------------------------
1726 C 5th derivatives (.le.2 in each coordinate)
1727 C
1728  if(ict(1).eq.5) then
1729  if(ict(2).eq.1) then
1730 C ! fxxyyz
1731  iaval=iaval+1
1732  do v=1,ivec
1733  i=iv(v)
1734  j=jv(v)
1735  k=kv(v)
1736  dx=dxv(v)
1737  dy=dyv(v)
1738  dz=dzv(v)
1739 C
1740  p2=
1741  > 4.0_r8*f(3,3,2,i,j,k)+12.0_r8*dy*f(3,4,2,i,j,k)
1742  > +12.0_r8*dx*(
1743  > f(4,3,2,i,j,k)+3.0_r8*dy*f(4,4,2,i,j,k))
1744  p3=
1745  > 4.0_r8*f(3,3,3,i,j,k)+12.0_r8*dy*f(3,4,3,i,j,k)
1746  > +12.0_r8*dx*(
1747  > f(4,3,3,i,j,k)+3.0_r8*dy*f(4,4,3,i,j,k))
1748  p4=
1749  > 4.0_r8*f(3,3,4,i,j,k)+12.0_r8*dy*f(3,4,4,i,j,k)
1750  > +12.0_r8*dx*(
1751  > f(4,3,4,i,j,k)+3.0_r8*dy*f(4,4,4,i,j,k))
1752 C
1753  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
1754  enddo
1755  endif
1756 C
1757  if(ict(3).eq.1) then
1758 C ! fxxyzz
1759  iaval=iaval+1
1760  do v=1,ivec
1761  i=iv(v)
1762  j=jv(v)
1763  k=kv(v)
1764  dx=dxv(v)
1765  dy=dyv(v)
1766  dz=dzv(v)
1767 C
1768  p3= 2.0_r8*
1769  > (f(3,2,3,i,j,k)+2.0_r8*dy*(f(3,3,3,i,j,k)+
1770  > 1.5_r8*dy*f(3,4,3,i,j,k)))
1771  > +6.0_r8*dx*
1772  > (f(4,2,3,i,j,k)+2.0_r8*dy*(f(4,3,3,i,j,k)+
1773  > 1.5_r8*dy*f(4,4,3,i,j,k)))
1774  p4= 2.0_r8*
1775  > (f(3,2,4,i,j,k)+2.0_r8*dy*(f(3,3,4,i,j,k)+
1776  > 1.5_r8*dy*f(3,4,4,i,j,k)))
1777  > +6.0_r8*dx*
1778  > (f(4,2,4,i,j,k)+2.0_r8*dy*(f(4,3,4,i,j,k)+
1779  > 1.5_r8*dy*f(4,4,4,i,j,k)))
1780 C
1781  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
1782  enddo
1783  endif
1784 C
1785  if(ict(4).eq.1) then
1786 C ! fxyyzz
1787  iaval=iaval+1
1788  do v=1,ivec
1789  i=iv(v)
1790  j=jv(v)
1791  k=kv(v)
1792  dx=dxv(v)
1793  dy=dyv(v)
1794  dz=dzv(v)
1795 C
1796  p3=
1797  > 2.0_r8*f(2,3,3,i,j,k)+6.0_r8*dy*f(2,4,3,i,j,k)
1798  > +2.0_r8*dx*(
1799  > 2.0_r8*f(3,3,3,i,j,k)+6.0_r8*dy*f(3,4,3,i,j,k)
1800  > +1.5_r8*dx*(
1801  > 2.0_r8*f(4,3,3,i,j,k)+6.0_r8*dy*f(4,4,3,i,j,k)))
1802  p4=
1803  > 2.0_r8*f(2,3,4,i,j,k)+6.0_r8*dy*f(2,4,4,i,j,k)
1804  > +2.0_r8*dx*(
1805  > 2.0_r8*f(3,3,4,i,j,k)+6.0_r8*dy*f(3,4,4,i,j,k)
1806  > +1.5_r8*dx*(
1807  > 2.0_r8*f(4,3,4,i,j,k)+6.0_r8*dy*f(4,4,4,i,j,k)))
1808 C
1809  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
1810  enddo
1811  endif
1812 C
1813  endif
1814 C
1815 C----------------------------------
1816 C 5th derivatives (3 in a coordinate)
1817 C
1818  if(ict(1).eq.-5) then
1819  if(ict(2).eq.1) then
1820 C ! fxxxyy
1821  iaval=iaval+1
1822  do v=1,ivec
1823  i=iv(v)
1824  j=jv(v)
1825  k=kv(v)
1826  dx=dxv(v)
1827  dy=dyv(v)
1828  dz=dzv(v)
1829 C
1830  p1=12.0_r8*(
1831  > f(4,3,1,i,j,k)+3.0_r8*dy*f(4,4,1,i,j,k))
1832  p2=12.0_r8*(
1833  > f(4,3,2,i,j,k)+3.0_r8*dy*f(4,4,2,i,j,k))
1834  p3=12.0_r8*(
1835  > f(4,3,3,i,j,k)+3.0_r8*dy*f(4,4,3,i,j,k))
1836  p4=12.0_r8*(
1837  > f(4,3,4,i,j,k)+3.0_r8*dy*f(4,4,4,i,j,k))
1838 C
1839  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
1840  enddo
1841  endif
1842 C
1843  if(ict(3).eq.1) then
1844 C ! fxxxyz
1845  iaval=iaval+1
1846  do v=1,ivec
1847  i=iv(v)
1848  j=jv(v)
1849  k=kv(v)
1850  dx=dxv(v)
1851  dy=dyv(v)
1852  dz=dzv(v)
1853 C
1854  p2= 6.0_r8*
1855  > (f(4,2,2,i,j,k)+2.0_r8*dy*(f(4,3,2,i,j,k)+
1856  > 1.5_r8*dy*f(4,4,2,i,j,k)))
1857  p3= 6.0_r8*
1858  > (f(4,2,3,i,j,k)+2.0_r8*dy*(f(4,3,3,i,j,k)+
1859  > 1.5_r8*dy*f(4,4,3,i,j,k)))
1860  p4= 6.0_r8*
1861  > (f(4,2,4,i,j,k)+2.0_r8*dy*(f(4,3,4,i,j,k)+
1862  > 1.5_r8*dy*f(4,4,4,i,j,k)))
1863 C
1864  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
1865  enddo
1866  endif
1867 C
1868  if(ict(4).eq.1) then
1869 C ! fxxxzz
1870  iaval=iaval+1
1871  do v=1,ivec
1872  i=iv(v)
1873  j=jv(v)
1874  k=kv(v)
1875  dx=dxv(v)
1876  dy=dyv(v)
1877  dz=dzv(v)
1878 C
1879  p3= 6.0_r8*
1880  > (f(4,1,3,i,j,k)+dy*(f(4,2,3,i,j,k)+dy*(f(4,3,3,i,j,k)+
1881  > dy*f(4,4,3,i,j,k))))
1882  p4= 6.0_r8*
1883  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
1884  > dy*f(4,4,4,i,j,k))))
1885 C
1886  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
1887  enddo
1888  endif
1889 C
1890  if(ict(5).eq.1) then
1891 C ! fxxyyy
1892  iaval=iaval+1
1893  do v=1,ivec
1894  i=iv(v)
1895  j=jv(v)
1896  k=kv(v)
1897  dx=dxv(v)
1898  dy=dyv(v)
1899  dz=dzv(v)
1900 C
1901  p1=12.0_r8*(f(3,4,1,i,j,k)+3.0_r8*dx*f(4,4,1,i,j,k))
1902  p2=12.0_r8*(f(3,4,2,i,j,k)+3.0_r8*dx*f(4,4,2,i,j,k))
1903  p3=12.0_r8*(f(3,4,3,i,j,k)+3.0_r8*dx*f(4,4,3,i,j,k))
1904  p4=12.0_r8*(f(3,4,4,i,j,k)+3.0_r8*dx*f(4,4,4,i,j,k))
1905 C
1906  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
1907  enddo
1908  endif
1909 C
1910  if(ict(6).eq.1) then
1911 C ! fxxzzz
1912  iaval=iaval+1
1913  do v=1,ivec
1914  i=iv(v)
1915  j=jv(v)
1916  k=kv(v)
1917  dx=dxv(v)
1918  dy=dyv(v)
1919  dz=dzv(v)
1920 C
1921  p4= 2.0_r8*
1922  > (f(3,1,4,i,j,k)+dy*(f(3,2,4,i,j,k)+dy*(f(3,3,4,i,j,k)+
1923  > dy*f(3,4,4,i,j,k))))
1924  > +6.0_r8*dx*
1925  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
1926  > dy*f(4,4,4,i,j,k))))
1927  fval(v,iaval)=6.0_r8*p4
1928  enddo
1929  endif
1930 C
1931  if(ict(7).eq.1) then
1932 C ! fxyyyz
1933  iaval=iaval+1
1934  do v=1,ivec
1935  i=iv(v)
1936  j=jv(v)
1937  k=kv(v)
1938  dx=dxv(v)
1939  dy=dyv(v)
1940  dz=dzv(v)
1941 C
1942  p2=
1943  > 6.0_r8*(f(2,4,2,i,j,k)
1944  > +2.0_r8*dx*(f(3,4,2,i,j,k)
1945  > +1.5_r8*dx*f(4,4,2,i,j,k)))
1946  p3=
1947  > 6.0_r8*(f(2,4,3,i,j,k)
1948  > +2.0_r8*dx*(f(3,4,3,i,j,k)
1949  > +1.5_r8*dx*f(4,4,3,i,j,k)))
1950  p4=
1951  > 6.0_r8*(f(2,4,4,i,j,k)
1952  > +2.0_r8*dx*(f(3,4,4,i,j,k)
1953  > +1.5_r8*dx*f(4,4,4,i,j,k)))
1954 C
1955  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
1956  enddo
1957  endif
1958 C
1959  if(ict(8).eq.1) then
1960 C ! fxyzzz
1961  iaval=iaval+1
1962  do v=1,ivec
1963  i=iv(v)
1964  j=jv(v)
1965  k=kv(v)
1966  dx=dxv(v)
1967  dy=dyv(v)
1968  dz=dzv(v)
1969 C
1970  p4=
1971  > f(2,2,4,i,j,k)+
1972  > dy*(2.0_r8*f(2,3,4,i,j,k)+dy*3.0_r8*f(2,4,4,i,j,k))
1973  > +2._r8*dx*(
1974  > f(3,2,4,i,j,k)+
1975  > dy*(2.0_r8*f(3,3,4,i,j,k)+dy*3.0_r8*f(3,4,4,i,j,k))
1976  > +1.5_r8*dx*(
1977  > f(4,2,4,i,j,k)+
1978  > dy*(2.0_r8*f(4,3,4,i,j,k)+dy*3.0_r8*f(4,4,4,i,j,k))
1979  > ))
1980  fval(v,iaval)=6.0_r8*p4
1981  enddo
1982  endif
1983 C
1984  if(ict(9).eq.1) then
1985 C ! fyyyzz
1986  iaval=iaval+1
1987  do v=1,ivec
1988  i=iv(v)
1989  j=jv(v)
1990  k=kv(v)
1991  dx=dxv(v)
1992  dy=dyv(v)
1993  dz=dzv(v)
1994 C
1995  p3=
1996  > 6.0_r8*(f(1,4,3,i,j,k)
1997  > +dx*(f(2,4,3,i,j,k)
1998  > +dx*(f(3,4,3,i,j,k)
1999  > +dx*f(4,4,3,i,j,k))))
2000  p4=
2001  > 6.0_r8*(f(1,4,4,i,j,k)
2002  > +dx*(f(2,4,4,i,j,k)
2003  > +dx*(f(3,4,4,i,j,k)
2004  > +dx*f(4,4,4,i,j,k))))
2005  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
2006  enddo
2007  endif
2008 C
2009  if(ict(10).eq.1) then
2010 C ! fyyzzz
2011  iaval=iaval+1
2012  do v=1,ivec
2013  i=iv(v)
2014  j=jv(v)
2015  k=kv(v)
2016  dx=dxv(v)
2017  dy=dyv(v)
2018  dz=dzv(v)
2019 C
2020  p4=
2021  > 2.0_r8*f(1,3,4,i,j,k)+6.0_r8*dy*f(1,4,4,i,j,k)
2022  > +dx*(2.0_r8*f(2,3,4,i,j,k)+6.0_r8*dy*f(2,4,4,i,j,k)
2023  > +dx*(2.0_r8*f(3,3,4,i,j,k)+6.0_r8*dy*f(3,4,4,i,j,k)
2024  > +dx*(2.0_r8*f(4,3,4,i,j,k)+6.0_r8*dy*f(4,4,4,i,j,k))))
2025  fval(v,iaval)=6.0_r8*p4
2026  enddo
2027  endif
2028 C
2029  endif
2030 C
2031 C----------------------------------
2032 C 6th derivatives (2 in each coordinate)
2033 C
2034  if(ict(1).eq.6) then
2035 C ! fxxyyzz
2036  iaval=iaval+1
2037  do v=1,ivec
2038  i=iv(v)
2039  j=jv(v)
2040  k=kv(v)
2041  dx=dxv(v)
2042  dy=dyv(v)
2043  dz=dzv(v)
2044 C
2045  p3=
2046  > 4.0_r8*f(3,3,3,i,j,k)+12.0_r8*dy*f(3,4,3,i,j,k)
2047  > +12.0_r8*dx*(
2048  > f(4,3,3,i,j,k)+3.0_r8*dy*f(4,4,3,i,j,k))
2049  p4=
2050  > 4.0_r8*f(3,3,4,i,j,k)+12.0_r8*dy*f(3,4,4,i,j,k)
2051  > +12.0_r8*dx*(
2052  > f(4,3,4,i,j,k)+3.0_r8*dy*f(4,4,4,i,j,k))
2053 C
2054  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
2055  enddo
2056  endif
2057 C
2058 C----------------------------------
2059 C 6th derivatives (3 in a coordinate)
2060 C
2061  if(ict(1).eq.-6) then
2062  if(ict(2).eq.1) then
2063 C ! fxxxyyy
2064  iaval=iaval+1
2065  do v=1,ivec
2066  i=iv(v)
2067  j=jv(v)
2068  k=kv(v)
2069  dx=dxv(v)
2070  dy=dyv(v)
2071  dz=dzv(v)
2072 C
2073  p1=36.0_r8*f(4,4,1,i,j,k)
2074  p2=36.0_r8*f(4,4,2,i,j,k)
2075  p3=36.0_r8*f(4,4,3,i,j,k)
2076  p4=36.0_r8*f(4,4,4,i,j,k)
2077 C
2078  fval(v,iaval)=p1+dz*(p2+dz*(p3+dz*p4))
2079  enddo
2080  endif
2081 C
2082  if(ict(3).eq.1) then
2083 C ! fxxxyyz
2084  iaval=iaval+1
2085  do v=1,ivec
2086  i=iv(v)
2087  j=jv(v)
2088  k=kv(v)
2089  dx=dxv(v)
2090  dy=dyv(v)
2091  dz=dzv(v)
2092 C
2093  p2=12.0_r8*(
2094  > f(4,3,2,i,j,k)+3.0_r8*dy*f(4,4,2,i,j,k))
2095  p3=12.0_r8*(
2096  > f(4,3,3,i,j,k)+3.0_r8*dy*f(4,4,3,i,j,k))
2097  p4=12.0_r8*(
2098  > f(4,3,4,i,j,k)+3.0_r8*dy*f(4,4,4,i,j,k))
2099 C
2100  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
2101  enddo
2102  endif
2103 C
2104  if(ict(4).eq.1) then
2105 C ! fxxxyzz
2106  iaval=iaval+1
2107  do v=1,ivec
2108  i=iv(v)
2109  j=jv(v)
2110  k=kv(v)
2111  dx=dxv(v)
2112  dy=dyv(v)
2113  dz=dzv(v)
2114 C
2115  p3= 6.0_r8*
2116  > (f(4,2,3,i,j,k)+2.0_r8*dy*(f(4,3,3,i,j,k)+
2117  > 1.5_r8*dy*f(4,4,3,i,j,k)))
2118  p4= 6.0_r8*
2119  > (f(4,2,4,i,j,k)+2.0_r8*dy*(f(4,3,4,i,j,k)+
2120  > 1.5_r8*dy*f(4,4,4,i,j,k)))
2121 C
2122  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
2123  enddo
2124  endif
2125 C
2126  if(ict(5).eq.1) then
2127 C ! fxxxzzz
2128  iaval=iaval+1
2129  do v=1,ivec
2130  i=iv(v)
2131  j=jv(v)
2132  k=kv(v)
2133  dx=dxv(v)
2134  dy=dyv(v)
2135  dz=dzv(v)
2136 C
2137  p4= 6.0_r8*
2138  > (f(4,1,4,i,j,k)+dy*(f(4,2,4,i,j,k)+dy*(f(4,3,4,i,j,k)+
2139  > dy*f(4,4,4,i,j,k))))
2140  fval(v,iaval)=6.0_r8*p4
2141  enddo
2142  endif
2143 C
2144  if(ict(6).eq.1) then
2145 C ! fxxyyyz
2146  iaval=iaval+1
2147  do v=1,ivec
2148  i=iv(v)
2149  j=jv(v)
2150  k=kv(v)
2151  dx=dxv(v)
2152  dy=dyv(v)
2153  dz=dzv(v)
2154 C
2155  p2=12.0_r8*(f(3,4,2,i,j,k)+3.0_r8*dx*f(4,4,2,i,j,k))
2156  p3=12.0_r8*(f(3,4,3,i,j,k)+3.0_r8*dx*f(4,4,3,i,j,k))
2157  p4=12.0_r8*(f(3,4,4,i,j,k)+3.0_r8*dx*f(4,4,4,i,j,k))
2158 C
2159  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
2160  enddo
2161  endif
2162 C
2163  if(ict(7).eq.1) then
2164 C ! fxxyzzz
2165  iaval=iaval+1
2166  do v=1,ivec
2167  i=iv(v)
2168  j=jv(v)
2169  k=kv(v)
2170  dx=dxv(v)
2171  dy=dyv(v)
2172  dz=dzv(v)
2173 C
2174  p4= 2.0_r8*
2175  > (f(3,2,4,i,j,k)+2.0_r8*dy*(f(3,3,4,i,j,k)+
2176  > 1.5_r8*dy*f(3,4,4,i,j,k)))
2177  > +6.0_r8*dx*
2178  > (f(4,2,4,i,j,k)+2.0_r8*dy*(f(4,3,4,i,j,k)+
2179  > 1.5_r8*dy*f(4,4,4,i,j,k)))
2180 C
2181  fval(v,iaval)=6.0_r8*p4
2182  enddo
2183  endif
2184 C
2185  if(ict(8).eq.1) then
2186 C ! fxyyyzz
2187  iaval=iaval+1
2188  do v=1,ivec
2189  i=iv(v)
2190  j=jv(v)
2191  k=kv(v)
2192  dx=dxv(v)
2193  dy=dyv(v)
2194  dz=dzv(v)
2195 C
2196  p3=
2197  > 6.0_r8*(f(2,4,3,i,j,k)
2198  > +2.0_r8*dx*(f(3,4,3,i,j,k)
2199  > +1.5_r8*dx*f(4,4,3,i,j,k)))
2200  p4=
2201  > 6.0_r8*(f(2,4,4,i,j,k)
2202  > +2.0_r8*dx*(f(3,4,4,i,j,k)
2203  > +1.5_r8*dx*f(4,4,4,i,j,k)))
2204 C
2205  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
2206  enddo
2207  endif
2208 C
2209  if(ict(9).eq.1) then
2210 C ! fxyyzzz
2211  iaval=iaval+1
2212  do v=1,ivec
2213  i=iv(v)
2214  j=jv(v)
2215  k=kv(v)
2216  dx=dxv(v)
2217  dy=dyv(v)
2218  dz=dzv(v)
2219 C
2220  p4=
2221  > 2.0_r8*f(2,3,4,i,j,k)+6.0_r8*dy*f(2,4,4,i,j,k)
2222  > +2.0_r8*dx*(
2223  > 2.0_r8*f(3,3,4,i,j,k)+6.0_r8*dy*f(3,4,4,i,j,k)
2224  > +1.5_r8*dx*(
2225  > 2.0_r8*f(4,3,4,i,j,k)+6.0_r8*dy*f(4,4,4,i,j,k)))
2226 C
2227  fval(v,iaval)=6.0_r8*p4
2228  enddo
2229  endif
2230 C
2231  if(ict(10).eq.1) then
2232 C ! fyyyzzz
2233  iaval=iaval+1
2234  do v=1,ivec
2235  i=iv(v)
2236  j=jv(v)
2237  k=kv(v)
2238  dx=dxv(v)
2239  dy=dyv(v)
2240  dz=dzv(v)
2241 C
2242  p4=
2243  > 6.0_r8*(f(1,4,4,i,j,k)
2244  > +dx*(f(2,4,4,i,j,k)
2245  > +dx*(f(3,4,4,i,j,k)
2246  > +dx*f(4,4,4,i,j,k))))
2247  fval(v,iaval)=6.0_r8*p4
2248  enddo
2249  endif
2250 C
2251  endif
2252 C
2253 C----------------------------------
2254 C 7th derivatives
2255 C
2256  if(abs(ict(1)).eq.7) then
2257  if(ict(2).eq.1) then
2258 C ! fxxxyyyz
2259  iaval=iaval+1
2260  do v=1,ivec
2261  i=iv(v)
2262  j=jv(v)
2263  k=kv(v)
2264  dz=dzv(v)
2265 C
2266  p2=36.0_r8*f(4,4,2,i,j,k)
2267  p3=36.0_r8*f(4,4,3,i,j,k)
2268  p4=36.0_r8*f(4,4,4,i,j,k)
2269  fval(v,iaval)=p2+dz*(2.0_r8*p3+dz*3.0_r8*p4)
2270  enddo
2271  endif
2272 C
2273  if(ict(3).eq.1) then
2274 C ! fxxxyyzz
2275  iaval=iaval+1
2276  do v=1,ivec
2277  i=iv(v)
2278  j=jv(v)
2279  k=kv(v)
2280  dx=dxv(v)
2281  dy=dyv(v)
2282  dz=dzv(v)
2283 C
2284  p3=12.0_r8*(
2285  > f(4,3,3,i,j,k)+3.0_r8*dy*f(4,4,3,i,j,k))
2286  p4=12.0_r8*(
2287  > f(4,3,4,i,j,k)+3.0_r8*dy*f(4,4,4,i,j,k))
2288  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
2289  enddo
2290  endif
2291 C
2292  if(ict(4).eq.1) then
2293 C ! fxxxyzzz
2294  iaval=iaval+1
2295  do v=1,ivec
2296  i=iv(v)
2297  j=jv(v)
2298  k=kv(v)
2299  dx=dxv(v)
2300  dy=dyv(v)
2301  dz=dzv(v)
2302 C
2303  p4= 6.0_r8*
2304  > (f(4,2,4,i,j,k)+2.0_r8*dy*(f(4,3,4,i,j,k)+
2305  > 1.5_r8*dy*f(4,4,4,i,j,k)))
2306 C
2307  fval(v,iaval)=6.0_r8*p4
2308  enddo
2309  endif
2310 C
2311  if(ict(5).eq.1) then
2312 C ! fxxyyyzz
2313  iaval=iaval+1
2314  do v=1,ivec
2315  i=iv(v)
2316  j=jv(v)
2317  k=kv(v)
2318  dx=dxv(v)
2319  dy=dyv(v)
2320  dz=dzv(v)
2321 C
2322  p3=12.0_r8*(f(3,4,3,i,j,k)+3.0_r8*dx*f(4,4,3,i,j,k))
2323  p4=12.0_r8*(f(3,4,4,i,j,k)+3.0_r8*dx*f(4,4,4,i,j,k))
2324 C
2325  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
2326  enddo
2327  endif
2328 C
2329  if(ict(6).eq.1) then
2330 C ! fxxyyzzz
2331  iaval=iaval+1
2332  do v=1,ivec
2333  i=iv(v)
2334  j=jv(v)
2335  k=kv(v)
2336  dx=dxv(v)
2337  dy=dyv(v)
2338  dz=dzv(v)
2339 C
2340  p4=
2341  > 4.0_r8*f(3,3,4,i,j,k)+12.0_r8*dy*f(3,4,4,i,j,k)
2342  > +12.0_r8*dx*(
2343  > f(4,3,4,i,j,k)+3.0_r8*dy*f(4,4,4,i,j,k))
2344 C
2345  fval(v,iaval)=6.0_r8*p4
2346  enddo
2347  endif
2348 C
2349  if(ict(7).eq.1) then
2350 C ! fxyyyzzz
2351  iaval=iaval+1
2352  do v=1,ivec
2353  i=iv(v)
2354  j=jv(v)
2355  k=kv(v)
2356  dx=dxv(v)
2357  dy=dyv(v)
2358  dz=dzv(v)
2359 C
2360  p4=
2361  > 6.0_r8*(f(2,4,4,i,j,k)
2362  > +2.0_r8*dx*(f(3,4,4,i,j,k)
2363  > +1.5_r8*dx*f(4,4,4,i,j,k)))
2364 C
2365  fval(v,iaval)=6.0_r8*p4
2366  enddo
2367  endif
2368 C
2369  endif
2370 C
2371 C----------------------------------
2372 C 8th derivatives
2373 C
2374  if(abs(ict(1)).eq.8) then
2375  if(ict(2).eq.1) then
2376 C ! fxxxyyyzz
2377  iaval=iaval+1
2378  do v=1,ivec
2379  i=iv(v)
2380  j=jv(v)
2381  k=kv(v)
2382  dz=dzv(v)
2383 C
2384  p3=36.0_r8*f(4,4,3,i,j,k)
2385  p4=36.0_r8*f(4,4,4,i,j,k)
2386  fval(v,iaval)=2.0_r8*p3+6.0_r8*dz*p4
2387  enddo
2388  endif
2389 C
2390  if(ict(3).eq.1) then
2391 C ! fxxxyyzzz
2392  iaval=iaval+1
2393  do v=1,ivec
2394  i=iv(v)
2395  j=jv(v)
2396  k=kv(v)
2397  dy=dyv(v)
2398 C
2399  p4=12.0_r8*(
2400  > f(4,3,4,i,j,k)+3.0_r8*dy*f(4,4,4,i,j,k))
2401  fval(v,iaval)=6.0_r8*p4
2402  enddo
2403  endif
2404 C
2405  if(ict(4).eq.1) then
2406 C ! fxxyyyzzz
2407  iaval=iaval+1
2408  do v=1,ivec
2409  i=iv(v)
2410  j=jv(v)
2411  k=kv(v)
2412  dx=dxv(v)
2413 C
2414  p4=12.0_r8*(f(3,4,4,i,j,k)+3.0_r8*dx*f(4,4,4,i,j,k))
2415  fval(v,iaval)=6.0_r8*p4
2416  enddo
2417  endif
2418 C
2419  endif
2420 C
2421 C----------------------------------
2422 C 9th derivative
2423 C
2424  if(abs(ict(1)).eq.9) then
2425 C ! fxxxyyyzzz
2426  iaval=iaval+1
2427  do v=1,ivec
2428  i=iv(v)
2429  j=jv(v)
2430  k=kv(v)
2431 C
2432  p4=36.0_r8*f(4,4,4,i,j,k)
2433  fval(v,iaval)=6.0_r8*p4
2434  enddo
2435  endif
2436 C
2437  return
2438  end
2439 
2440 
2441 
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17