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