V3FIT
r8bcspeval.f
1 c
2 c bcspeval -- eval bicubic spline function and/or derivatives
3 c
4  subroutine r8bcspeval(xget,yget,iselect,fval,
5  > x,nx,y,ny,ilinx,iliny,f,inf3,ier)
6 c
7  IMPLICIT NONE
8  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
9  integer iselect(6)
10  integer ilinx,iliny,nx,ny,inf3,ier
11 c
12  real*8 xget,yget
13  real*8 fval(*)
14  real*8 x(nx),y(ny),f(4,4,inf3,ny)
15 c
16 c modification -- dmc 11 Jan 1999 -- remove SAVE stmts; break routine
17 C into these parts:
18 C
19 C bcspevxy -- find grid cell of target pt.
20 C bcspevfn -- evaluate function using output of bcpsevxy
21 C
22 C in cases where multiple functions are defined on the same grid,
23 C time can be saved by using bcspevxy once and then bcspevfn
24 C multiple times.
25 c
26 c input:
27 c (xget,yget) 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
31 c iselect select desired output
32 c
33 c iselect(1)=1 -- want function value (f) itself
34 c iselect(2)=1 -- want df/dx
35 c iselect(3)=1 -- want df/dy
36 c iselect(4)=1 -- want d2f/dx2
37 c iselect(5)=1 -- want d2f/dy2
38 c iselect(6)=1 -- want d2f/dxdy
39 c
40 c example: iselect(1)=iselect(2)=iselect(3)=1
41 c f, df/dx, and df/dy all evaluated
42 c iselect(4)=iselect(5)=iselect(6)=0
43 c 2nd derivatives not evaluated.
44 c
45 c the number of non zero values iselect(1:6)
46 c determines the number of outputs...
47 c see fval (output) description.
48 c
49 c new dmc December 2005 -- access to higher derivatives (even if not
50 c continuous-- but can only go up to 3rd derivatives on any one coordinate.
51 c if iselect(1)=3 -- want 3rd derivatives
52 c iselect(2)=1 for d3f/dx3
53 c iselect(3)=1 for d3f/dx2dy
54 c iselect(4)=1 for d3f/dxdy2
55 c iselect(5)=1 for d3f/dy3
56 c number of non-zero values iselect(2:5) gives no. of outputs
57 c if iselect(1)=4 -- want 4th derivatives
58 c iselect(2)=1 for d4f/dx3dy
59 c iselect(3)=1 for d4f/dx2dy2
60 c iselect(4)=1 for d4f/dxdy3
61 c number of non-zero values iselect(2:4) gives no. of outputs
62 c if iselect(1)=5 -- want 5th derivatives
63 c iselect(2)=1 for d5f/dx3dy2
64 c iselect(3)=1 for d5f/dx2dy3
65 c number of non-zero values iselect(2:3) gives no. of outputs
66 c if iselect(1)=6 -- want 6th derivatives
67 c d6f/dx3dy3 -- one value is returned.
68 c
69 c x(1...nx) independent coordinate x, strict ascending
70 c y(1...ny) independent coordinate y, strict ascending
71 c
72 c ilinx -- =1: flag that x is linearly spaced (avoid search for speed)
73 c iliny -- =1: flag that y is linearly spaced (avoid search for speed)
74 c
75 c **CAUTION** actual even spacing of x, y is NOT CHECKED HERE!
76 c
77 c
78 c f the function values (at grid points) and spline coefs
79 c
80 c evaluation formula: for point x btw x(i) and x(i+1), dx=x-x(i)
81 c and y btw y(j) and y(j+1), dy=y-y(j),
82 c
83 c spline value =
84 c f(1,1,i,j) + dx*f(2,1,i,j) + dx**2*f(3,1,i,j) + dx**3*f(4,1,i,j)
85 c +dy*(f(1,2,i,j) + dx*f(2,2,i,j) + dx**2*f(3,2,i,j) + dx**3*f(4,2,i,j))
86 c +d2*(f(1,3,i,j) + dx*f(2,3,i,j) + dx**2*f(3,3,i,j) + dx**3*f(4,3,i,j))
87 c +d3*(f(1,4,i,j) + dx*f(2,4,i,j) + dx**2*f(3,4,i,j) + dx**3*f(4,4,i,j))
88 c
89 c where d2=dy**2 and d3=dy**3.
90 c
91 c output:
92 c up to 6 elements of fval, ordered as follows:
93 c fval(1)=function value or lowest order derivative requested
94 c fval(2)=next order derivative
95 c etc
96 c the ordering is a subset of the sequence given under the "iselect"
97 c description.
98 c
99 c ier = 0 -- successful completion; = 1 -- an error occurred.
100 c
101 c-------------------------------------------------------------------
102 c local
103 c
104  integer :: i(1)
105  integer :: j(1)
106 c
107  real*8 dx(1),dy(1)
108 c
109 c--------------------------
110 c
111  i(1)=0
112  j(1)=0
113 c
114  call r8bcspevxy(xget,yget,x,nx,y,ny,ilinx,iliny,
115  > i(1),j(1),dx(1),dy(1),ier)
116  if(ier.ne.0) return
117 c
118  call r8bcspevfn(iselect,1,1,fval,i,j,dx,dy,f,inf3,ny)
119 c
120  return
121  end
122 c
123 c-------------------------------------------------------------------------
124 c bcspevxy -- look up x-y zone
125 c
126 c this is the "first part" of bcspeval, see comments, above.
127 c
128  subroutine r8bcspevxy(xget,yget,x,nx,y,ny,ilinx,iliny,
129  > i,j,dx,dy,ier)
130 c
131 !============
132 ! idecl: explicitize implicit INTEGER declarations:
133  IMPLICIT NONE
134  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
135  INTEGER nxm,nym,ii,jj
136 !============
137 ! idecl: explicitize implicit REAL declarations:
138  real*8 zxget,zyget,zxtol,zytol
139 !============
140  integer nx,ny ! array dimensions
141 c
142  real*8 xget,yget ! target point
143  real*8 x(nx),y(ny) ! indep. coords.
144 c
145  integer ilinx ! =1: assume x evenly spaced
146  integer iliny ! =1: assume y evenly spaced
147 c
148 c output of bcspevxy
149 c
150  integer i,j ! index to cell containing target pt
151  REAL*8 dx,dy ! displacement of target pt w/in cell
152  ! dx=x-x(i) dy=y-y(j)
153 C
154  integer ier ! return ier.ne.0 on error
155 c
156 c------------------------------------
157 c
158  ier=0
159 c
160 c range check
161 c
162  zxget=xget
163  zyget=yget
164 
165  if((xget.lt.x(1)).or.(xget.gt.x(nx))) then
166  zxtol=4.0e-7_r8*max(abs(x(1)),abs(x(nx)))
167  if((xget.lt.x(1)-zxtol).or.(xget.gt.x(nx)+zxtol)) then
168  ier=1
169  write(6,1001) xget,x(1),x(nx)
170  1001 format(' ?bcspeval: xget=',1pe11.4,' out of range ',
171  > 1pe11.4,' to ',1pe11.4)
172  else
173  if((xget.lt.x(1)-0.5_r8*zxtol).or.
174  > (xget.gt.x(nx)+0.5_r8*zxtol))
175  > write(6,1011) xget,x(1),x(nx)
176  1011 format(' %bcspeval: xget=',1pe15.8,' beyond range ',
177  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
178  if(xget.lt.x(1)) then
179  zxget=x(1)
180  else
181  zxget=x(nx)
182  endif
183  endif
184  endif
185  if((yget.lt.y(1)).or.(yget.gt.y(ny))) then
186  zytol=4.0e-7_r8*max(abs(y(1)),abs(y(ny)))
187  if((yget.lt.y(1)-zytol).or.(yget.gt.y(ny)+zytol)) then
188  ier=1
189  write(6,1002) yget,y(1),y(ny)
190  1002 format(' ?bcspeval: yget=',1pe11.4,' out of range ',
191  > 1pe11.4,' to ',1pe11.4)
192  else
193  if((yget.lt.y(1)-0.5_r8*zytol).or.(yget.gt.y(ny)+0.5_r8*zytol))
194  > write(6,1012) yget,y(1),y(ny)
195  1012 format(' %bcspeval: yget=',1pe15.8,' beyond range ',
196  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
197  if(yget.lt.y(1)) then
198  zyget=y(1)
199  else
200  zyget=y(ny)
201  endif
202  endif
203  endif
204  if(ier.ne.0) return
205 c
206 c now find interval in which target point lies..
207 c
208  nxm=nx-1
209  nym=ny-1
210 c
211  if(ilinx.eq.1) then
212  ii=1+nxm*(zxget-x(1))/(x(nx)-x(1))
213  i=min(nxm, ii)
214  if(zxget.lt.x(i)) then
215  i=i-1
216  else if(zxget.gt.x(i+1)) then
217  i=i+1
218  endif
219  else
220  if((1.le.i).and.(i.lt.nxm)) then
221  if((x(i).le.zxget).and.(zxget.le.x(i+1))) then
222  continue ! already have the zone
223  else
224  call r8zonfind(x,nx,zxget,i)
225  endif
226  else
227  i=nx/2
228  call r8zonfind(x,nx,zxget,i)
229  endif
230  endif
231 c
232  if(iliny.eq.1) then
233  jj=1+nym*(zyget-y(1))/(y(ny)-y(1))
234  j=min(nym, jj)
235  if(zyget.lt.y(j)) then
236  j=j-1
237  else if(zyget.gt.y(j+1)) then
238  j=j+1
239  endif
240  else
241  if((1.le.j).and.(j.lt.nym)) then
242  if((y(j).le.zyget).and.(zyget.le.y(j+1))) then
243  continue ! already have the zone
244  else
245  call r8zonfind(y,ny,zyget,j)
246  endif
247  else
248  j=ny/2
249  call r8zonfind(y,ny,zyget,j)
250  endif
251  endif
252 c
253  dx=zxget-x(i)
254  dy=zyget-y(j)
255 c
256  return
257  end
258 c------------------------------------------------------------------------
259 c bcspevfn -- OK now evaluate the bicubic spline
260 c
261  subroutine r8bcspevfn(ict,ivec,ivd,fval,iv,jv,dxv,dyv,f,inf3,ny)
262 c
263 c input:
264 !============
265 ! idecl: explicitize implicit INTEGER declarations:
266  IMPLICIT NONE
267  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
268  INTEGER ny,iaval,i,j
269 !============
270 ! idecl: explicitize implicit REAL declarations:
271  real*8 dx,dy
272 !============
273  integer ict(6) ! selector:
274 c ict(1)=1 for f (don't evaluate f if ict(1)=0)
275 c ict(2)=1 for df/dx ""
276 c ict(3)=1 for df/dy ""
277 c ict(4)=1 for d2f/dx2
278 c ict(5)=1 for d2f/dy2
279 c ict(6)=1 for d2f/dxdy
280 c
281 c note: if ict(1)=-1, evaluate f,d2f/dx2,d2f/dy2,d4f/dx2dy2
282 c
283 c the number of non zero values ict(1:6)
284 c determines the number of outputs...
285 c
286 c new dmc December 2005 -- access to higher derivatives (even if not
287 c continuous-- but can only go up to 3rd derivatives on any one coordinate.
288 c if ict(1)=3 -- want 3rd derivatives
289 c ict(2)=1 for d3f/dx3
290 c ict(3)=1 for d3f/dx2dy
291 c ict(4)=1 for d3f/dxdy2
292 c ict(5)=1 for d3f/dy3
293 c number of non-zero values ict(2:5) gives no. of outputs
294 c if ict(1)=4 -- want 4th derivatives
295 c ict(2)=1 for d4f/dx3dy
296 c ict(3)=1 for d4f/dx2dy2
297 c ict(4)=1 for d4f/dxdy3
298 c number of non-zero values ict(2:4) gives no. of outputs
299 c if ict(1)=5 -- want 5th derivatives
300 c ict(2)=1 for d5f/dx3dy2
301 c ict(3)=1 for d5f/dx2dy3
302 c number of non-zero values ict(2:3) gives no. of outputs
303 c if ict(1)=6 -- want 6th derivatives
304 c d6f/dx3dy3 -- one value is returned.
305 c
306  integer ivec,ivd ! vector dimensioning
307 c
308 c ivec-- number of vector pts (spline values to look up)
309 c ivd -- 1st dimension of fval, .ge.ivec
310 c
311 c output:
312  real*8 fval(ivd,*) ! output array
313 c
314 c v = index to element in vector;
315 c fval(v,1) = first item requested by ict(...),
316 c fval(v,2) = 2nd item requested, ...etc...
317 c
318 c input:
319  integer iv(ivec),jv(ivec) ! grid cell indices -- vectors
320  REAL*8 dxv(ivec),dyv(ivec) ! displacements w/in cell -- vectors
321 c
322  integer inf3 ! 3rd dimension of f -- .ge. nx
323  REAL*8 f(4,4,inf3,ny) ! bicubic fcn spline coeffs array
324 c
325 c usage example:
326 c
327 c 1. for each element (xx(v),yy(v)) in a vector of (x,y) pairs,
328 c find the x and y zone indices and displacements with respect
329 c to the "lower left corner" of the zone; store these in vectors
330 c iv,jv and dxv,dyv.
331 c
332 c 2. set ict(1)=0, ict(2)=1, ict(3)=1, the rest zero -- get only
333 c the 1st derivatives.
334 c
335 c 3. ivec is the length of the vector; ivd is the 1st dimension
336 c of the array fval to receive the output
337 c
338 c real fval(ivd,6)
339 c real xv(ivd),yv(ivd)
340 c integer iv(ivd),jv(ivd)
341 c real dxv(ivd),dyv(ivd)
342 c integer ict(6)
343 c
344 c real fspline(4,4,nx,ny) ! spline coeffs
345 c data ict/0,1,1,0,0,0/ ! this call: want 1st derivatives
346 c ! only ... these will be output to
347 c ! fval(*,1) fval(*,2)
348 c ...
349 c do iv=1,ivec
350 c ... ! find indices and displacements
351 c enddo
352 c call bcspevfn(ict,ivec,ivd,fval,iv,jv,dxv,dyv,fspline,nx,ny)
353 c
354 c-------------------------------------------------------------------
355 c local:
356 c
357  integer v ! vector element index
358 c
359 c OK can now do evaluations
360 c
361  iaval=0 ! fval addressing
362 c
363  if(ict(1).le.2) then
364  if((ict(1).gt.0).or.(ict(1).eq.-1)) then
365 c evaluate f
366  iaval=iaval+1
367  do v=1,ivec
368  i=iv(v)
369  j=jv(v)
370  dx=dxv(v)
371  dy=dyv(v)
372  fval(v,iaval)=
373  > f(1,1,i,j)+dy*(f(1,2,i,j)+dy*(f(1,3,i,j)+dy*f(1,4,i,j)))
374  > +dx*(f(2,1,i,j)+dy*(f(2,2,i,j)+dy*(f(2,3,i,j)+dy*f(2,4,i,j)))
375  > +dx*(f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j)))
376  > +dx*(f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j))))))
377  enddo
378  endif
379 c
380  if((ict(2).gt.0).and.(ict(1).ne.-1)) then
381 c evaluate df/dx
382  iaval=iaval+1
383  do v=1,ivec
384  i=iv(v)
385  j=jv(v)
386  dx=dxv(v)
387  dy=dyv(v)
388  fval(v,iaval)=
389  > f(2,1,i,j)+dy*(f(2,2,i,j)+dy*(f(2,3,i,j)+dy*f(2,4,i,j)))
390  > +2.0_r8*dx*(
391  > f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j)))
392  > +1.5_r8*dx*(
393  > f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j)))
394  > ))
395  enddo
396  endif
397 c
398  if((ict(3).gt.0).and.(ict(1).ne.-1)) then
399 c evaluate df/dy
400  iaval=iaval+1
401  do v=1,ivec
402  i=iv(v)
403  j=jv(v)
404  dx=dxv(v)
405  dy=dyv(v)
406  fval(v,iaval)=
407  > f(1,2,i,j)+dy*(2.0_r8*f(1,3,i,j)+dy*3.0_r8*f(1,4,i,j))
408  > +dx*(f(2,2,i,j)+dy*(2.0_r8*f(2,3,i,j)+dy*3.0_r8*f(2,4,i,j))
409  > +dx*(f(3,2,i,j)+dy*(2.0_r8*f(3,3,i,j)+dy*3.0_r8*f(3,4,i,j))
410  > +dx*(f(4,2,i,j)+dy*(2.0_r8*f(4,3,i,j)+dy*3.0_r8*f(4,4,i,j))
411  > )))
412  enddo
413  endif
414 c
415  if((ict(4).gt.0).or.(ict(1).eq.-1)) then
416 c evaluate d2f/dx2
417  iaval=iaval+1
418  do v=1,ivec
419  i=iv(v)
420  j=jv(v)
421  dx=dxv(v)
422  dy=dyv(v)
423  fval(v,iaval)=
424  > 2.0_r8*(
425  > f(3,1,i,j)+dy*(f(3,2,i,j)+dy*(f(3,3,i,j)+dy*f(3,4,i,j))))
426  > +6.0_r8*dx*(
427  > f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j))))
428  enddo
429  endif
430 c
431  if((ict(5).gt.0).or.(ict(1).eq.-1)) then
432 c evaluate d2f/dy2
433  iaval=iaval+1
434  do v=1,ivec
435  i=iv(v)
436  j=jv(v)
437  dx=dxv(v)
438  dy=dyv(v)
439  fval(v,iaval)=
440  > 2.0_r8*f(1,3,i,j)+6.0_r8*dy*f(1,4,i,j)
441  > +dx*(2.0_r8*f(2,3,i,j)+6.0_r8*dy*f(2,4,i,j)
442  > +dx*(2.0_r8*f(3,3,i,j)+6.0_r8*dy*f(3,4,i,j)
443  > +dx*(2.0_r8*f(4,3,i,j)+6.0_r8*dy*f(4,4,i,j))))
444  enddo
445  endif
446 c
447  if((ict(6).gt.0).and.(ict(1).ne.-1)) then
448 c evaluate d2f/dxdy
449  iaval=iaval+1
450  do v=1,ivec
451  i=iv(v)
452  j=jv(v)
453  dx=dxv(v)
454  dy=dyv(v)
455  fval(v,iaval)=
456  > f(2,2,i,j)+dy*(2.0_r8*f(2,3,i,j)+dy*3.0_r8*f(2,4,i,j))
457  > +2._r8*dx*(f(3,2,i,j)+dy*(2.0_r8*f(3,3,i,j)+dy*3.0_r8*f(3,4,i,j))
458  >+1.5_r8*dx*(f(4,2,i,j)+dy*(2.0_r8*f(4,3,i,j)+dy*3.0_r8*f(4,4,i,j))
459  > ))
460  enddo
461  endif
462 c
463  if(ict(1).eq.-1) then
464  iaval=iaval+1
465  do v=1,ivec
466  i=iv(v)
467  j=jv(v)
468  dx=dxv(v)
469  dy=dyv(v)
470  fval(v,iaval)=
471  > 4.0_r8*f(3,3,i,j)+12.0_r8*dy*f(3,4,i,j)
472  > +dx*(12.0_r8*f(4,3,i,j)+36.0_r8*dy*f(4,4,i,j))
473  enddo
474  endif
475 c
476 c-----------------------------------
477 c access to 3rd derivatives
478 c
479  else if(ict(1).eq.3) then
480  if(ict(2).eq.1) then
481 c evaluate d3f/dx3 (not continuous)
482  iaval=iaval+1
483  do v=1,ivec
484  i=iv(v)
485  j=jv(v)
486  dy=dyv(v)
487  fval(v,iaval)=
488  > +6.0_r8*(
489  > f(4,1,i,j)+dy*(f(4,2,i,j)+dy*(f(4,3,i,j)+dy*f(4,4,i,j))))
490  enddo
491  endif
492 c
493  if(ict(3).eq.1) then
494 c evaluate d3f/dx2dy
495  iaval=iaval+1
496  do v=1,ivec
497  i=iv(v)
498  j=jv(v)
499  dx=dxv(v)
500  dy=dyv(v)
501  fval(v,iaval)=
502  > 2.0_r8*(
503  > f(3,2,i,j)+dy*(2.0_r8*f(3,3,i,j)+dy*3.0_r8*f(3,4,i,j)))
504  > +6.0_r8*dx*(
505  > f(4,2,i,j)+dy*(2.0_r8*f(4,3,i,j)+dy*3.0_r8*f(4,4,i,j)))
506  enddo
507  endif
508 c
509  if(ict(4).eq.1) then
510 c evaluate d3f/dxdy2
511  iaval=iaval+1
512  do v=1,ivec
513  i=iv(v)
514  j=jv(v)
515  dx=dxv(v)
516  dy=dyv(v)
517  fval(v,iaval)=
518  > (2.0_r8*f(2,3,i,j)+6.0_r8*dy*f(2,4,i,j)
519  > +2.0_r8*dx*(2.0_r8*f(3,3,i,j)+6.0_r8*dy*f(3,4,i,j)
520  > +1.5_r8*dx*(2.0_r8*f(4,3,i,j)+6.0_r8*dy*f(4,4,i,j))
521  > ))
522  enddo
523  endif
524 
525  if(ict(5).eq.1) then
526 c evaluate d3f/dy3 (not continuous)
527  iaval=iaval+1
528  do v=1,ivec
529  i=iv(v)
530  j=jv(v)
531  dx=dxv(v)
532  fval(v,iaval)=6.0_r8*(f(1,4,i,j)+
533  > dx*(f(2,4,i,j)+dx*(f(3,4,i,j)+dx*f(4,4,i,j))))
534  enddo
535  endif
536 c
537 c-----------------------------------
538 c access to 4th derivatives
539 c
540  else if(ict(1).eq.4) then
541  if(ict(2).eq.1) then
542 c evaluate d4f/dx3dy (not continuous)
543  iaval=iaval+1
544  do v=1,ivec
545  i=iv(v)
546  j=jv(v)
547  dy=dyv(v)
548  fval(v,iaval)=
549  > +6.0_r8*(
550  > f(4,2,i,j)+dy*2.0_r8*(f(4,3,i,j)+dy*1.5_r8*f(4,4,i,j)))
551  enddo
552  endif
553 c
554  if(ict(3).eq.1) then
555 c evaluate d4f/dx2dy2
556  iaval=iaval+1
557  do v=1,ivec
558  i=iv(v)
559  j=jv(v)
560  dx=dxv(v)
561  dy=dyv(v)
562  fval(v,iaval)=
563  > 4.0_r8*f(3,3,i,j)+12.0_r8*dy*f(3,4,i,j)
564  > +dx*(12.0_r8*f(4,3,i,j)+36.0_r8*dy*f(4,4,i,j))
565  enddo
566  endif
567 c
568  if(ict(4).eq.1) then
569 c evaluate d4f/dxdy3 (not continuous)
570  iaval=iaval+1
571  do v=1,ivec
572  i=iv(v)
573  j=jv(v)
574  dx=dxv(v)
575  fval(v,iaval)=
576  > 6.0_r8*(f(2,4,i,j)
577  > +2.0_r8*dx*(f(3,4,i,j)+1.5_r8*dx*f(4,4,i,j)))
578  enddo
579  endif
580 c
581 c-----------------------------------
582 c access to 5th derivatives
583 c
584  else if(ict(1).eq.5) then
585  if(ict(2).eq.1) then
586 c evaluate d5f/dx3dy2 (not continuous)
587  iaval=iaval+1
588  do v=1,ivec
589  i=iv(v)
590  j=jv(v)
591  dy=dyv(v)
592  fval(v,iaval)=
593  > +12.0_r8*(f(4,3,i,j)+dy*3.0_r8*f(4,4,i,j))
594  enddo
595  endif
596 c
597  if(ict(3).eq.1) then
598 c evaluate d5f/dx3dy2 (not continuous)
599  iaval=iaval+1
600  do v=1,ivec
601  i=iv(v)
602  j=jv(v)
603  dx=dxv(v)
604  fval(v,iaval)=
605  > 12.0_r8*(f(3,4,i,j)+dx*3.0_r8*f(4,4,i,j))
606  enddo
607  endif
608 c
609 c-----------------------------------
610 c access to 6th derivatives
611 c
612  else if(ict(1).eq.6) then
613 c evaluate d6f/dx3dy3 (not continuous)
614  iaval=iaval+1
615  do v=1,ivec
616  i=iv(v)
617  j=jv(v)
618  fval(v,iaval)=
619  > 36.0_r8*f(4,4,i,j)
620  enddo
621  endif
622 c
623  return
624  end
625 c----------------------
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17