V3FIT
r8evbicub.f
1  subroutine r8evbicub(xget,yget,x,nx,y,ny,ilinx,iliny,
2  > f,inf2,ict,fval,ier)
3 C
4 C evaluate a 2d cubic Spline interpolant on a rectilinear
5 C grid -- this is C2 in both directions.
6 C
7 C this subroutine calls two subroutines:
8 C herm2xy -- find cell containing (xget,yget)
9 C fvbicub -- evaluate interpolant function and (optionally) derivatives
10 C
11 C input arguments:
12 C ================
13 C
14 !============
15 ! idecl: explicitize implicit INTEGER declarations:
16  IMPLICIT NONE
17  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
18  INTEGER inf2
19 !============
20  integer nx,ny ! grid sizes
21  real*8 xget,yget ! target of this interpolation
22  real*8 x(nx) ! ordered x grid
23  real*8 y(ny) ! ordered y grid
24  integer ilinx ! ilinx=1 => assume x evenly spaced
25  integer iliny ! iliny=1 => assume y evenly spaced
26 C
27  real*8 f(0:3,inf2,ny) ! function data
28 C
29 C f 2nd dimension inf2 must be .ge. nx
30 C contents of f:
31 C
32 C f(0,i,j) = f @ x(i),y(j)
33 C f(1,i,j) = d2f/dx2 @ x(i),y(j)
34 C f(2,i,j) = d2f/dy2 @ x(i),y(j)
35 C f(3,i,j) = d4f/dx2dy2 @ x(i),y(j)
36 C
37 C (these are spline coefficients selected for continuous 2-
38 C diffentiability, see mkbicub[w].for)
39 C
40  integer ict(6) ! code specifying output desired
41 C
42 C ict(1)=1 -- return f (0, don't)
43 C ict(2)=1 -- return df/dx (0, don't)
44 C ict(3)=1 -- return df/dy (0, don't)
45 C ict(4)=1 -- return d2f/dx2 (0, don't)
46 C ict(5)=1 -- return d2f/dy2 (0, don't)
47 C ict(6)=1 -- return d2f/dxdy (0, don't)
48 c the number of non zero values ict(1:6)
49 c determines the number of outputs...
50 c
51 c new dmc December 2005 -- access to higher derivatives (even if not
52 c continuous-- but can only go up to 3rd derivatives on any one coordinate.
53 c if ict(1)=3 -- want 3rd derivatives
54 c ict(2)=1 for d3f/dx3
55 c ict(3)=1 for d3f/dx2dy
56 c ict(4)=1 for d3f/dxdy2
57 c ict(5)=1 for d3f/dy3
58 c number of non-zero values ict(2:5) gives no. of outputs
59 c if ict(1)=4 -- want 4th derivatives
60 c ict(2)=1 for d4f/dx3dy
61 c ict(3)=1 for d4f/dx2dy2
62 c ict(4)=1 for d4f/dxdy3
63 c number of non-zero values ict(2:4) gives no. of outputs
64 c if ict(1)=5 -- want 5th derivatives
65 c ict(2)=1 for d5f/dx3dy2
66 c ict(3)=1 for d5f/dx2dy3
67 c number of non-zero values ict(2:3) gives no. of outputs
68 c if ict(1)=6 -- want 6th derivatives
69 c d6f/dx3dy3 -- one value is returned.
70 C
71 C output arguments:
72 C =================
73 C
74  real*8 fval(*) ! output data
75  integer ier ! error code =0 ==> no error
76 C
77 C fval(1) receives the first output (depends on ict(...) spec)
78 C fval(2) receives the second output (depends on ict(...) spec)
79 C fval(3) receives the third output (depends on ict(...) spec)
80 C fval(4) receives the fourth output (depends on ict(...) spec)
81 C fval(5) receives the fourth output (depends on ict(...) spec)
82 C fval(6) receives the fourth output (depends on ict(...) spec)
83 C
84 C examples:
85 C on input ict = [1,1,1,0,0,1]
86 C on output fval= [f,df/dx,df/dy,d2f/dxdy], elements 5 & 6 not referenced.
87 C
88 C on input ict = [1,0,0,0,0,0]
89 C on output fval= [f] ... elements 2 -- 6 never referenced.
90 C
91 C on input ict = [0,0,0,1,1,0]
92 C on output fval= [d2f/dx2,d2f/dy2] ... elements 3 -- 6 never referenced.
93 C
94 C on input ict = [0,0,1,0,0,0]
95 C on output fval= [df/dy] ... elements 2 -- 6 never referenced.
96 C
97 C ier -- completion code: 0 means OK
98 C-------------------
99 C local:
100 C
101  integer i(1),j(1) ! cell indices
102 C
103 C normalized displacement from (x(i),y(j)) corner of cell.
104 C xparam=0 @x(i) xparam=1 @x(i+1)
105 C yparam=0 @y(j) yparam=1 @y(j+1)
106 C
107  real*8 xparam(1),yparam(1)
108 C
109 C cell dimensions and
110 C inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
111 C
112  real*8 hx(1),hy(1)
113  real*8 hxi(1),hyi(1)
114 C
115 C 0 .le. xparam .le. 1
116 C 0 .le. yparam .le. 1
117 C
118 C ** the interface is very similar to herm2ev.for; can use herm2xy **
119 C---------------------------------------------------------------------
120 C
121  i(1)=0
122  j(1)=0
123  call r8herm2xy(xget,yget,x,nx,y,ny,ilinx,iliny,
124  > i(1),j(1),xparam(1),yparam(1),hx(1),hxi(1),hy(1),hyi(1),ier)
125  if(ier.ne.0) return
126 c
127  call r8fvbicub(ict,1,1,
128  > fval,i,j,xparam,yparam,hx,hxi,hy,hyi,f,inf2,ny)
129 C
130  return
131  end
132 C---------------------------------------------------------------------
133 C evaluate C1 cubic Hermite function interpolation -- 2d fcn
134 C --vectorized-- dmc 10 Feb 1999
135 C
136  subroutine r8fvbicub(ict,ivec,ivecd,
137  > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
138  > fin,inf2,ny)
139 C
140 !============
141 ! idecl: explicitize implicit INTEGER declarations:
142  IMPLICIT NONE
143  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
144  INTEGER ny,inf2,iadr,i,j
145 !============
146 ! idecl: explicitize implicit REAL declarations:
147  real*8 z36th,xp,xpi,xp2,xpi2,cx,cxi,hx2,yp,ypi,yp2,ypi2,cy
148  real*8 cyi,hy2,cxd,cxdi,cyd,cydi
149 !============
150  integer ict(6) ! requested output control
151  integer ivec ! vector length
152  integer ivecd ! vector dimension (1st dim of fval)
153 C
154  integer ii(ivec),jj(ivec) ! target cells (i,j)
155  real*8 xparam(ivec),yparam(ivec)
156  ! normalized displacements from (i,j) corners
157 C
158  real*8 hx(ivec),hy(ivec) ! grid spacing, and
159  real*8 hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
160  ! & 1/(y(j+1)-y(j))
161 C
162  real*8 fin(0:3,inf2,ny) ! interpolant data (cf "evbicub")
163 C
164  real*8 fval(ivecd,*) ! output returned
165 C
166 C for detailed description of fin, ict and fval see subroutine
167 C evbicub comments. Note ict is not vectorized; the same output
168 C is expected to be returned for all input vector data points.
169 C
170 C note that the index inputs ii,jj and parameter inputs
171 C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
172 C output array fval has a vector ** 1st dimension ** whose
173 C size must be given as a separate argument
174 C
175 C to use this routine in scalar mode, pass in ivec=ivecd=1
176 C
177 C---------------
178 C Spline evaluation consists of a "mixing" of the interpolant
179 C data using the linear functionals xparam, xpi = 1-xparam,
180 C yparam, ypi = 1-yparam, and the cubic functionals
181 C xparam**3-xparam, xpi**3-xpi, yparam**3-yparam, ypi**3-ypi ...
182 C and their derivatives as needed.
183 C
184  integer v
185  REAL*8 sum
186 C
187  real*8, parameter :: sixth = 0.166666666666666667_r8
188 C
189 C---------------
190 C ...in x direction
191 C
192  z36th=sixth*sixth
193  iadr=0
194 C
195  if(ict(1).le.2) then
196 C
197 C get desired values:
198 C
199  if(ict(1).eq.1) then
200 C
201 C function value:
202 C
203  iadr=iadr+1
204  do v=1,ivec
205  i=ii(v)
206  j=jj(v)
207 C
208 C in x direction...
209 C
210  xp=xparam(v)
211  xpi=1.0_r8-xp
212  xp2=xp*xp
213  xpi2=xpi*xpi
214 C
215  cx=xp*(xp2-1.0_r8)
216  cxi=xpi*(xpi2-1.0_r8)
217  hx2=hx(v)*hx(v)
218 C
219 C ...and in y direction
220 C
221  yp=yparam(v)
222  ypi=1.0_r8-yp
223  yp2=yp*yp
224  ypi2=ypi*ypi
225 C
226  cy=yp*(yp2-1.0_r8)
227  cyi=ypi*(ypi2-1.0_r8)
228  hy2=hy(v)*hy(v)
229 C
230  sum=xpi*(ypi*fin(0,i,j) +yp*fin(0,i,j+1))+
231  > xp*(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1))
232 C
233  sum=sum+sixth*hx2*(
234  > cxi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
235  > cx*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
236 C
237  sum=sum+sixth*hy2*(
238  > xpi*(cyi*fin(2,i,j) +cy*fin(2,i,j+1))+
239  > xp*(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
240 C
241  sum=sum+z36th*hx2*hy2*(
242  > cxi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
243  > cx*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
244 C
245  fval(v,iadr)=sum
246  enddo
247  endif
248 C
249  if(ict(2).eq.1) then
250 C
251 C df/dx:
252 C
253  iadr=iadr+1
254  do v=1,ivec
255  i=ii(v)
256  j=jj(v)
257 C
258 C in x direction...
259 C
260  xp=xparam(v)
261  xpi=1.0_r8-xp
262  xp2=xp*xp
263  xpi2=xpi*xpi
264 
265  cxd=3.0_r8*xp2-1.0_r8
266  cxdi=-3.0_r8*xpi2+1.0_r8
267 C
268 C ...and in y direction
269 C
270  yp=yparam(v)
271  ypi=1.0_r8-yp
272  yp2=yp*yp
273  ypi2=ypi*ypi
274 C
275  cy=yp*(yp2-1.0_r8)
276  cyi=ypi*(ypi2-1.0_r8)
277  hy2=hy(v)*hy(v)
278 C
279  sum=hxi(v)*(
280  > -(ypi*fin(0,i,j) +yp*fin(0,i,j+1))
281  > +(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1)))
282 C
283  sum=sum+sixth*hx(v)*(
284  > cxdi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
285  > cxd*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
286 C
287  sum=sum+sixth*hxi(v)*hy2*(
288  > -(cyi*fin(2,i,j) +cy*fin(2,i,j+1))
289  > +(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
290 C
291  sum=sum+z36th*hx(v)*hy2*(
292  > cxdi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
293  > cxd*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
294 C
295  fval(v,iadr)=sum
296  enddo
297  endif
298 C
299  if(ict(3).eq.1) then
300 C
301 C df/dy:
302 C
303  iadr=iadr+1
304  do v=1,ivec
305  i=ii(v)
306  j=jj(v)
307 C
308 C in x direction...
309 C
310  xp=xparam(v)
311  xpi=1.0_r8-xp
312  xp2=xp*xp
313  xpi2=xpi*xpi
314 C
315  cx=xp*(xp2-1.0_r8)
316  cxi=xpi*(xpi2-1.0_r8)
317  hx2=hx(v)*hx(v)
318 C
319 C ...and in y direction
320 C
321  yp=yparam(v)
322  ypi=1.0_r8-yp
323  yp2=yp*yp
324  ypi2=ypi*ypi
325 
326  cyd=3.0_r8*yp2-1.0_r8
327  cydi=-3.0_r8*ypi2+1.0_r8
328 C
329  sum=hyi(v)*(
330  > xpi*(-fin(0,i,j) +fin(0,i,j+1))+
331  > xp*(-fin(0,i+1,j)+fin(0,i+1,j+1)))
332 C
333  sum=sum+sixth*hx2*hyi(v)*(
334  > cxi*(-fin(1,i,j) +fin(1,i,j+1))+
335  > cx*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
336 C
337  sum=sum+sixth*hy(v)*(
338  > xpi*(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))+
339  > xp*(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
340 C
341  sum=sum+z36th*hx2*hy(v)*(
342  > cxi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
343  > cx*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
344 C
345  fval(v,iadr)=sum
346  enddo
347  endif
348 C
349  if(ict(4).eq.1) then
350 C
351 C d2f/dx2:
352 C
353  iadr=iadr+1
354  do v=1,ivec
355  i=ii(v)
356  j=jj(v)
357 C
358 C in x direction...
359 C
360  xp=xparam(v)
361  xpi=1.0_r8-xp
362 C
363 C ...and in y direction
364 C
365  yp=yparam(v)
366  ypi=1.0_r8-yp
367  yp2=yp*yp
368  ypi2=ypi*ypi
369 C
370  cy=yp*(yp2-1.0_r8)
371  cyi=ypi*(ypi2-1.0_r8)
372  hy2=hy(v)*hy(v)
373 C
374  sum=(
375  > xpi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
376  > xp*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
377 C
378  sum=sum+sixth*hy2*(
379  > xpi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
380  > xp*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
381 C
382  fval(v,iadr)=sum
383  enddo
384  endif
385 C
386  if(ict(5).eq.1) then
387 C
388 C d2f/dy2:
389 C
390  iadr=iadr+1
391  do v=1,ivec
392  i=ii(v)
393  j=jj(v)
394 C
395 C in x direction...
396 C
397  xp=xparam(v)
398  xpi=1.0_r8-xp
399  xp2=xp*xp
400  xpi2=xpi*xpi
401 C
402  cx=xp*(xp2-1.0_r8)
403  cxi=xpi*(xpi2-1.0_r8)
404  hx2=hx(v)*hx(v)
405 C
406 C ...and in y direction
407 C
408  yp=yparam(v)
409  ypi=1.0_r8-yp
410 C
411  sum=(
412  > xpi*(ypi*fin(2,i,j) +yp*fin(2,i,j+1))+
413  > xp*(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
414 C
415  sum=sum+sixth*hx2*(
416  > cxi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
417  > cx*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
418 C
419  fval(v,iadr)=sum
420  enddo
421  endif
422 C
423  if(ict(6).eq.1) then
424 C
425 C d2f/dxdy:
426 C
427  iadr=iadr+1
428  do v=1,ivec
429  i=ii(v)
430  j=jj(v)
431 C
432 C in x direction...
433 C
434  xp=xparam(v)
435  xpi=1.0_r8-xp
436  xp2=xp*xp
437  xpi2=xpi*xpi
438 
439  cxd=3.0_r8*xp2-1.0_r8
440  cxdi=-3.0_r8*xpi2+1.0_r8
441 C
442 C ...and in y direction
443 C
444  yp=yparam(v)
445  ypi=1.0_r8-yp
446  yp2=yp*yp
447  ypi2=ypi*ypi
448 
449  cyd=3.0_r8*yp2-1.0_r8
450  cydi=-3.0_r8*ypi2+1.0_r8
451 C
452  sum=hxi(v)*hyi(v)*(
453  > fin(0,i,j) -fin(0,i,j+1)
454  > -fin(0,i+1,j)+fin(0,i+1,j+1))
455 C
456  sum=sum+sixth*hx(v)*hyi(v)*(
457  > cxdi*(-fin(1,i,j) +fin(1,i,j+1))+
458  > cxd*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
459 C
460  sum=sum+sixth*hxi(v)*hy(v)*(
461  > -(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))
462  > +(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
463 C
464  sum=sum+z36th*hx(v)*hy(v)*(
465  > cxdi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
466  > cxd*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
467 C
468  fval(v,iadr)=sum
469  enddo
470  endif
471 C
472 C-------------------------------------------------
473 C
474  else if(ict(1).eq.3) then
475  if(ict(2).eq.1) then
476 c evaluate d3f/dx3 (not continuous)
477  iadr=iadr+1
478  do v=1,ivec
479  i=ii(v)
480  j=jj(v)
481  yp=yparam(v)
482  ypi=1.0_r8-yp
483  yp2=yp*yp
484  ypi2=ypi*ypi
485  cy=yp*(yp2-1.0_r8)
486  cyi=ypi*(ypi2-1.0_r8)
487  hy2=hy(v)*hy(v)
488  sum=hxi(v)*(
489  > -(ypi*fin(1,i,j) +yp*fin(1,i,j+1))
490  > +(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
491 C
492  sum=sum+sixth*hy2*hxi(v)*(
493  > -(cyi*fin(3,i,j) +cy*fin(3,i,j+1))
494  > +(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
495 C
496  fval(v,iadr)=sum
497  enddo
498  endif
499 c
500  if(ict(3).eq.1) then
501 c evaluate d3f/dx2dy
502  iadr=iadr+1
503  do v=1,ivec
504  i=ii(v)
505  j=jj(v)
506  xp=xparam(v)
507  xpi=1.0_r8-xp
508  yp=yparam(v)
509  ypi=1.0_r8-yp
510  yp2=yp*yp
511  ypi2=ypi*ypi
512  cyd=3.0_r8*yp2-1.0_r8
513  cydi=-3.0_r8*ypi2+1.0_r8
514 C
515  sum=hyi(v)*(
516  > xpi*(-fin(1,i,j) +fin(1,i,j+1))+
517  > xp*(-fin(1,i+1,j) +fin(1,i+1,j+1)))
518 C
519  sum=sum+sixth*hy(v)*(
520  > xpi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
521  > xp*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
522 C
523  fval(v,iadr)=sum
524  enddo
525  endif
526 c
527  if(ict(4).eq.1) then
528 c evaluate d3f/dxdy2
529  iadr=iadr+1
530  do v=1,ivec
531  i=ii(v)
532  j=jj(v)
533  xp=xparam(v)
534  xpi=1.0_r8-xp
535  xp2=xp*xp
536  xpi2=xpi*xpi
537  cxd=3.0_r8*xp2-1.0_r8
538  cxdi=-3.0_r8*xpi2+1.0_r8
539  yp=yparam(v)
540  ypi=1.0_r8-yp
541 C
542  sum=hxi(v)*(
543  > -(ypi*fin(2,i,j) +yp*fin(2,i,j+1))
544  > +(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
545 C
546  sum=sum+sixth*hx(v)*(
547  > cxdi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
548  > cxd*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
549 C
550  fval(v,iadr)=sum
551  enddo
552  endif
553 
554  if(ict(5).eq.1) then
555 c evaluate d3f/dy3 (not continuous)
556  iadr=iadr+1
557  do v=1,ivec
558  i=ii(v)
559  j=jj(v)
560 C
561  xp=xparam(v)
562  xpi=1.0_r8-xp
563  xp2=xp*xp
564  xpi2=xpi*xpi
565 C
566  cx=xp*(xp2-1.0_r8)
567  cxi=xpi*(xpi2-1.0_r8)
568  hx2=hx(v)*hx(v)
569 C
570  sum=hyi(v)*(
571  > xpi*(-fin(2,i,j) +fin(2,i,j+1))+
572  > xp*(-fin(2,i+1,j) +fin(2,i+1,j+1)))
573 C
574  sum=sum+sixth*hx2*hyi(v)*(
575  > cxi*(-fin(3,i,j) +fin(3,i,j+1))+
576  > cx*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
577 C
578  fval(v,iadr)=sum
579  enddo
580  endif
581 c
582 c-----------------------------------
583 c access to 4th derivatives
584 c
585  else if(ict(1).eq.4) then
586  if(ict(2).eq.1) then
587 c evaluate d4f/dx3dy (not continuous)
588  iadr=iadr+1
589  do v=1,ivec
590  i=ii(v)
591  j=jj(v)
592  yp=yparam(v)
593  ypi=1.0_r8-yp
594  yp2=yp*yp
595  ypi2=ypi*ypi
596  cyd=3.0_r8*yp2-1.0_r8
597  cydi=-3.0_r8*ypi2+1.0_r8
598 C
599  sum=hxi(v)*hyi(v)*(
600  > +( fin(1,i,j) -fin(1,i,j+1))
601  > +(-fin(1,i+1,j)+fin(1,i+1,j+1)))
602 C
603  sum=sum+sixth*hy(v)*hxi(v)*(
604  > -(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))
605  > +(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
606 C
607  fval(v,iadr)=sum
608  enddo
609  endif
610 c
611  if(ict(3).eq.1) then
612 c evaluate d4f/dx2dy2
613  iadr=iadr+1
614  do v=1,ivec
615  i=ii(v)
616  j=jj(v)
617 C
618  xp=xparam(v)
619  xpi=1.0_r8-xp
620  yp=yparam(v)
621  ypi=1.0_r8-yp
622 C
623  sum=xpi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
624  > xp*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1))
625 C
626  fval(v,iadr)=sum
627  enddo
628  endif
629 c
630  if(ict(4).eq.1) then
631 c evaluate d4f/dxdy3 (not continuous)
632  iadr=iadr+1
633  do v=1,ivec
634  i=ii(v)
635  j=jj(v)
636 C
637  xp=xparam(v)
638  xpi=1.0_r8-xp
639  xp2=xp*xp
640  xpi2=xpi*xpi
641 C
642  cxd=3.0_r8*xp2-1.0_r8
643  cxdi=-3.0_r8*xpi2+1.0_r8
644 C
645  sum=hyi(v)*hxi(v)*(
646  > +( fin(2,i,j) -fin(2,i,j+1))
647  > +(-fin(2,i+1,j)+fin(2,i+1,j+1)))
648 C
649  sum=sum+sixth*hx(v)*hyi(v)*(
650  > cxdi*(-fin(3,i,j) +fin(3,i,j+1))+
651  > cxd*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
652 C
653  fval(v,iadr)=sum
654  enddo
655  endif
656 c
657 c-----------------------------------
658 c access to 5th derivatives
659 c
660  else if(ict(1).eq.5) then
661  if(ict(2).eq.1) then
662 c evaluate d5f/dx3dy2 (not continuous)
663  iadr=iadr+1
664  do v=1,ivec
665  i=ii(v)
666  j=jj(v)
667 C
668  yp=yparam(v)
669  ypi=1.0_r8-yp
670 C
671  sum=hxi(v)*(
672  > -(ypi*fin(3,i,j) +yp*fin(3,i,j+1))
673  > +(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
674 C
675  fval(v,iadr)=sum
676  enddo
677  endif
678 c
679  if(ict(3).eq.1) then
680 c evaluate d5f/dx2dy3 (not continuous)
681  iadr=iadr+1
682  do v=1,ivec
683  i=ii(v)
684  j=jj(v)
685 C
686  xp=xparam(v)
687  xpi=1.0_r8-xp
688 C
689  sum=hyi(v)*(
690  > xpi*(-fin(3,i,j) +fin(3,i,j+1))+
691  > xp*(-fin(3,i+1,j)+fin(3,i+1,j+1)))
692 C
693  fval(v,iadr)=sum
694  enddo
695  endif
696 c
697 c-----------------------------------
698 c access to 6th derivatives
699 c
700  else if(ict(1).eq.6) then
701 c evaluate d6f/dx3dy3 (not continuous)
702  iadr=iadr+1
703  do v=1,ivec
704  i=ii(v)
705  j=jj(v)
706  sum=hxi(v)*hyi(v)*(
707  > +( fin(3,i,j) -fin(3,i,j+1))
708  > +(-fin(3,i+1,j)+fin(3,i+1,j+1)))
709  fval(v,iadr)=sum
710  enddo
711  endif
712 c
713  return
714  end
715 C---------------------------------------------------------------------
716 C evaluate C1 cubic Hermite function interpolation -- 2d fcn
717 C --vectorized-- dmc 10 Feb 1999
718 C --optimized for VARIATION along x axis ONLY--
719 C
720  subroutine r8fvbicubx(ict,ivec,ivecd,
721  > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
722  > fin,inf2,ny)
723 C
724 !============
725 ! idecl: explicitize implicit INTEGER declarations:
726  IMPLICIT NONE
727  INTEGER, PARAMETER :: R8=selected_real_kind(12,100)
728  INTEGER ny,inf2,iadr,j,i
729 !============
730 ! idecl: explicitize implicit REAL declarations:
731  real*8 z36th,yp,ypi,yp2,ypi2,cy,cyi,hy2,xp,xpi,xp2,xpi2,cx
732  real*8 cxi,hx2,cxd,cxdi,cyd,cydi
733 !============
734  integer ict(6) ! requested output control
735  integer ivec ! vector length
736  integer ivecd ! vector dimension (1st dim of fval)
737 C
738  integer ii(ivec),jj ! target cells (i,j)
739  real*8 xparam(ivec),yparam
740  ! normalized displacements from (i,j) corners
741 C
742  real*8 hx(ivec),hy ! grid spacing, and
743  real*8 hxi(ivec),hyi ! inverse grid spacing 1/(x(i+1)-x(i))
744  ! & 1/(y(j+1)-y(j))
745 C
746  real*8 fin(0:3,inf2,ny) ! interpolant data (cf "evbicub")
747 C
748  real*8 fval(ivecd,*) ! output returned
749 C
750 C for detailed description of fin, ict and fval see subroutine
751 C evbicub comments. Note ict is not vectorized; the same output
752 C is expected to be returned for all input vector data points.
753 C
754 C note that the index inputs ii,jj and parameter inputs
755 C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
756 C output array fval has a vector ** 1st dimension ** whose
757 C size must be given as a separate argument
758 C
759 C to use this routine in scalar mode, pass in ivec=ivecd=1
760 C
761 C---------------
762 C Spline evaluation consists of a "mixing" of the interpolant
763 C data using the linear functionals xparam, xpi = 1-xparam,
764 C yparam, ypi = 1-yparam, and the cubic functionals
765 C xparam**3-xparam, xpi**3-xpi, yparam**3-yparam, ypi**3-ypi ...
766 C and their derivatives as needed.
767 C
768  integer v
769  REAL*8 sum
770 C
771  REAL*8, parameter :: sixth = 0.166666666666666667_r8
772 C
773 C---------------
774 C ...in x direction
775 C
776  z36th=sixth*sixth
777  iadr=0
778 C
779  if(ict(1).le.2) then
780 C
781 C get desired values:
782 C
783  if(ict(1).eq.1) then
784 C
785 C function value:
786 C
787  j=jj
788 C
789 C ...and in y direction
790 C
791  yp=yparam
792  ypi=1.0_r8-yp
793  yp2=yp*yp
794  ypi2=ypi*ypi
795 C
796  cy=yp*(yp2-1.0_r8)
797  cyi=ypi*(ypi2-1.0_r8)
798  hy2=hy*hy
799 C
800  iadr=iadr+1
801  do v=1,ivec
802  i=ii(v)
803 C
804 C in x direction...
805 C
806  xp=xparam(v)
807  xpi=1.0_r8-xp
808  xp2=xp*xp
809  xpi2=xpi*xpi
810 C
811  cx=xp*(xp2-1.0_r8)
812  cxi=xpi*(xpi2-1.0_r8)
813  hx2=hx(v)*hx(v)
814 C
815  sum=xpi*(ypi*fin(0,i,j) +yp*fin(0,i,j+1))+
816  > xp*(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1))
817 C
818  sum=sum+sixth*hx2*(
819  > cxi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
820  > cx*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
821 C
822  sum=sum+sixth*hy2*(
823  > xpi*(cyi*fin(2,i,j) +cy*fin(2,i,j+1))+
824  > xp*(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
825 C
826  sum=sum+z36th*hx2*hy2*(
827  > cxi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
828  > cx*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
829 C
830  fval(v,iadr)=sum
831  enddo
832  endif
833 C
834  if(ict(2).eq.1) then
835 C
836 C df/dx:
837 C
838  j=jj
839 C
840 C ...and in y direction
841 C
842  yp=yparam
843  ypi=1.0_r8-yp
844  yp2=yp*yp
845  ypi2=ypi*ypi
846 C
847  cy=yp*(yp2-1.0_r8)
848  cyi=ypi*(ypi2-1.0_r8)
849  hy2=hy*hy
850 C
851  iadr=iadr+1
852  do v=1,ivec
853  i=ii(v)
854 C
855 C in x direction...
856 C
857  xp=xparam(v)
858  xpi=1.0_r8-xp
859  xp2=xp*xp
860  xpi2=xpi*xpi
861 
862  cxd=3.0_r8*xp2-1.0_r8
863  cxdi=-3.0_r8*xpi2+1.0_r8
864 C
865  sum=hxi(v)*(
866  > -(ypi*fin(0,i,j) +yp*fin(0,i,j+1))
867  > +(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1)))
868 C
869  sum=sum+sixth*hx(v)*(
870  > cxdi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
871  > cxd*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
872 C
873  sum=sum+sixth*hxi(v)*hy2*(
874  > -(cyi*fin(2,i,j) +cy*fin(2,i,j+1))
875  > +(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
876 C
877  sum=sum+z36th*hx(v)*hy2*(
878  > cxdi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
879  > cxd*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
880 C
881  fval(v,iadr)=sum
882  enddo
883  endif
884 C
885  if(ict(3).eq.1) then
886 C
887 C df/dy:
888 C
889  j=jj
890 C
891 C ...and in y direction
892 C
893  yp=yparam
894  ypi=1.0_r8-yp
895  yp2=yp*yp
896  ypi2=ypi*ypi
897 
898  cyd=3.0_r8*yp2-1.0_r8
899  cydi=-3.0_r8*ypi2+1.0_r8
900 C
901  iadr=iadr+1
902  do v=1,ivec
903  i=ii(v)
904 C
905 C in x direction...
906 C
907  xp=xparam(v)
908  xpi=1.0_r8-xp
909  xp2=xp*xp
910  xpi2=xpi*xpi
911 C
912  cx=xp*(xp2-1.0_r8)
913  cxi=xpi*(xpi2-1.0_r8)
914  hx2=hx(v)*hx(v)
915 C
916  sum=hyi*(
917  > xpi*(-fin(0,i,j) +fin(0,i,j+1))+
918  > xp*(-fin(0,i+1,j)+fin(0,i+1,j+1)))
919 C
920  sum=sum+sixth*hx2*hyi*(
921  > cxi*(-fin(1,i,j) +fin(1,i,j+1))+
922  > cx*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
923 C
924  sum=sum+sixth*hy*(
925  > xpi*(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))+
926  > xp*(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
927 C
928  sum=sum+z36th*hx2*hy*(
929  > cxi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
930  > cx*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
931 C
932  fval(v,iadr)=sum
933  enddo
934  endif
935 C
936  if(ict(4).eq.1) then
937 C
938 C d2f/dx2:
939 C
940  j=jj
941 C
942 C ...and in y direction
943 C
944  yp=yparam
945  ypi=1.0_r8-yp
946  yp2=yp*yp
947  ypi2=ypi*ypi
948 C
949  cy=yp*(yp2-1.0_r8)
950  cyi=ypi*(ypi2-1.0_r8)
951  hy2=hy*hy
952 C
953  iadr=iadr+1
954  do v=1,ivec
955  i=ii(v)
956 C
957 C in x direction...
958 C
959  xp=xparam(v)
960  xpi=1.0_r8-xp
961 C
962  sum=(
963  > xpi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
964  > xp*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
965 C
966  sum=sum+sixth*hy2*(
967  > xpi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
968  > xp*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
969 C
970  fval(v,iadr)=sum
971  enddo
972  endif
973 C
974  if(ict(5).eq.1) then
975 C
976 C d2f/dy2:
977 C
978  j=jj
979 C
980 C ...and in y direction
981 C
982  yp=yparam
983  ypi=1.0_r8-yp
984 C
985  iadr=iadr+1
986  do v=1,ivec
987  i=ii(v)
988 C
989 C in x direction...
990 C
991  xp=xparam(v)
992  xpi=1.0_r8-xp
993  xp2=xp*xp
994  xpi2=xpi*xpi
995 C
996  cx=xp*(xp2-1.0_r8)
997  cxi=xpi*(xpi2-1.0_r8)
998  hx2=hx(v)*hx(v)
999 C
1000  sum=(
1001  > xpi*(ypi*fin(2,i,j) +yp*fin(2,i,j+1))+
1002  > xp*(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
1003 C
1004  sum=sum+sixth*hx2*(
1005  > cxi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1006  > cx*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1007 C
1008  fval(v,iadr)=sum
1009  enddo
1010  endif
1011 C
1012  if(ict(6).eq.1) then
1013 C
1014 C d2f/dxdy:
1015 C
1016  j=jj
1017 C
1018 C ...and in y direction
1019 C
1020  yp=yparam
1021  ypi=1.0_r8-yp
1022  yp2=yp*yp
1023  ypi2=ypi*ypi
1024 
1025  cyd=3.0_r8*yp2-1.0_r8
1026  cydi=-3.0_r8*ypi2+1.0_r8
1027 C
1028  iadr=iadr+1
1029  do v=1,ivec
1030  i=ii(v)
1031 C
1032 C in x direction...
1033 C
1034  xp=xparam(v)
1035  xpi=1.0_r8-xp
1036  xp2=xp*xp
1037  xpi2=xpi*xpi
1038 
1039  cxd=3.0_r8*xp2-1.0_r8
1040  cxdi=-3.0_r8*xpi2+1.0_r8
1041 C
1042  sum=hxi(v)*hyi*(
1043  > fin(0,i,j) -fin(0,i,j+1)
1044  > -fin(0,i+1,j)+fin(0,i+1,j+1))
1045 C
1046  sum=sum+sixth*hx(v)*hyi*(
1047  > cxdi*(-fin(1,i,j) +fin(1,i,j+1))+
1048  > cxd*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1049 C
1050  sum=sum+sixth*hxi(v)*hy*(
1051  > -(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))
1052  > +(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
1053 C
1054  sum=sum+z36th*hx(v)*hy*(
1055  > cxdi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1056  > cxd*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1057 C
1058  fval(v,iadr)=sum
1059  enddo
1060  endif
1061 C
1062 C-------------------------------------------------
1063 C
1064  else if(ict(1).eq.3) then
1065  if(ict(2).eq.1) then
1066 c evaluate d3f/dx3 (not continuous)
1067  j=jj
1068  yp=yparam
1069  ypi=1.0_r8-yp
1070  yp2=yp*yp
1071  ypi2=ypi*ypi
1072  cy=yp*(yp2-1.0_r8)
1073  cyi=ypi*(ypi2-1.0_r8)
1074  hy2=hy*hy
1075 C
1076  iadr=iadr+1
1077  do v=1,ivec
1078  i=ii(v)
1079  sum=hxi(v)*(
1080  > -(ypi*fin(1,i,j) +yp*fin(1,i,j+1))
1081  > +(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
1082 C
1083  sum=sum+sixth*hy2*hxi(v)*(
1084  > -(cyi*fin(3,i,j) +cy*fin(3,i,j+1))
1085  > +(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
1086 C
1087  fval(v,iadr)=sum
1088  enddo
1089  endif
1090 c
1091  if(ict(3).eq.1) then
1092 c evaluate d3f/dx2dy
1093  j=jj
1094 C
1095  yp=yparam
1096  ypi=1.0_r8-yp
1097  yp2=yp*yp
1098  ypi2=ypi*ypi
1099  cyd=3.0_r8*yp2-1.0_r8
1100  cydi=-3.0_r8*ypi2+1.0_r8
1101 C
1102  iadr=iadr+1
1103  do v=1,ivec
1104  i=ii(v)
1105  xp=xparam(v)
1106  xpi=1.0_r8-xp
1107 C
1108  sum=hyi*(
1109  > xpi*(-fin(1,i,j) +fin(1,i,j+1))+
1110  > xp*(-fin(1,i+1,j) +fin(1,i+1,j+1)))
1111 C
1112  sum=sum+sixth*hy*(
1113  > xpi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1114  > xp*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1115 C
1116  fval(v,iadr)=sum
1117  enddo
1118  endif
1119 c
1120  if(ict(4).eq.1) then
1121 c evaluate d3f/dxdy2
1122  j=jj
1123  yp=yparam
1124  ypi=1.0_r8-yp
1125 C
1126  iadr=iadr+1
1127  do v=1,ivec
1128  i=ii(v)
1129  xp=xparam(v)
1130  xpi=1.0_r8-xp
1131  xp2=xp*xp
1132  xpi2=xpi*xpi
1133  cxd=3.0_r8*xp2-1.0_r8
1134  cxdi=-3.0_r8*xpi2+1.0_r8
1135 C
1136  sum=hxi(v)*(
1137  > -(ypi*fin(2,i,j) +yp*fin(2,i,j+1))
1138  > +(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
1139 C
1140  sum=sum+sixth*hx(v)*(
1141  > cxdi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1142  > cxd*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1143 C
1144  fval(v,iadr)=sum
1145  enddo
1146  endif
1147 
1148  if(ict(5).eq.1) then
1149 c evaluate d3f/dy3 (not continuous)
1150  iadr=iadr+1
1151  j=jj
1152  do v=1,ivec
1153  i=ii(v)
1154 C
1155  xp=xparam(v)
1156  xpi=1.0_r8-xp
1157  xp2=xp*xp
1158  xpi2=xpi*xpi
1159 C
1160  cx=xp*(xp2-1.0_r8)
1161  cxi=xpi*(xpi2-1.0_r8)
1162  hx2=hx(v)*hx(v)
1163 C
1164  sum=hyi*(
1165  > xpi*(-fin(2,i,j) +fin(2,i,j+1))+
1166  > xp*(-fin(2,i+1,j) +fin(2,i+1,j+1)))
1167 C
1168  sum=sum+sixth*hx2*hyi*(
1169  > cxi*(-fin(3,i,j) +fin(3,i,j+1))+
1170  > cx*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
1171 C
1172  fval(v,iadr)=sum
1173  enddo
1174  endif
1175 c
1176 c-----------------------------------
1177 c access to 4th derivatives
1178 c
1179  else if(ict(1).eq.4) then
1180  if(ict(2).eq.1) then
1181 c evaluate d4f/dx3dy (not continuous)
1182  iadr=iadr+1
1183  j=jj
1184 C
1185  yp=yparam
1186  ypi=1.0_r8-yp
1187  yp2=yp*yp
1188  ypi2=ypi*ypi
1189  cyd=3.0_r8*yp2-1.0_r8
1190  cydi=-3.0_r8*ypi2+1.0_r8
1191 C
1192  do v=1,ivec
1193  i=ii(v)
1194 C
1195  sum=hxi(v)*hyi*(
1196  > +( fin(1,i,j) -fin(1,i,j+1))
1197  > +(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1198 C
1199  sum=sum+sixth*hy*hxi(v)*(
1200  > -(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))
1201  > +(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1202 C
1203  fval(v,iadr)=sum
1204  enddo
1205  endif
1206 c
1207  if(ict(3).eq.1) then
1208 c evaluate d4f/dx2dy2
1209  iadr=iadr+1
1210  do v=1,ivec
1211  i=ii(v)
1212  j=jj
1213 C
1214  xp=xparam(v)
1215  xpi=1.0_r8-xp
1216  yp=yparam
1217  ypi=1.0_r8-yp
1218 C
1219  sum=xpi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1220  > xp*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1))
1221 C
1222  fval(v,iadr)=sum
1223  enddo
1224  endif
1225 c
1226  if(ict(4).eq.1) then
1227 c evaluate d4f/dxdy3 (not continuous)
1228  j=jj
1229 C
1230  iadr=iadr+1
1231  do v=1,ivec
1232  i=ii(v)
1233 C
1234  xp=xparam(v)
1235  xpi=1.0_r8-xp
1236  xp2=xp*xp
1237  xpi2=xpi*xpi
1238 C
1239  cxd=3.0_r8*xp2-1.0_r8
1240  cxdi=-3.0_r8*xpi2+1.0_r8
1241 C
1242  sum=hyi*hxi(v)*(
1243  > +( fin(2,i,j) -fin(2,i,j+1))
1244  > +(-fin(2,i+1,j)+fin(2,i+1,j+1)))
1245 C
1246  sum=sum+sixth*hx(v)*hyi*(
1247  > cxdi*(-fin(3,i,j) +fin(3,i,j+1))+
1248  > cxd*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
1249 C
1250  fval(v,iadr)=sum
1251  enddo
1252  endif
1253 c
1254 c-----------------------------------
1255 c access to 5th derivatives
1256 c
1257  else if(ict(1).eq.5) then
1258  if(ict(2).eq.1) then
1259 c evaluate d5f/dx3dy2 (not continuous)
1260  j=jj
1261 C
1262  yp=yparam
1263  ypi=1.0_r8-yp
1264 C
1265  iadr=iadr+1
1266  do v=1,ivec
1267  i=ii(v)
1268 C
1269  sum=hxi(v)*(
1270  > -(ypi*fin(3,i,j) +yp*fin(3,i,j+1))
1271  > +(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1272 C
1273  fval(v,iadr)=sum
1274  enddo
1275  endif
1276 c
1277  if(ict(3).eq.1) then
1278 c evaluate d5f/dx2dy3 (not continuous)
1279  j=jj
1280  iadr=iadr+1
1281  do v=1,ivec
1282  i=ii(v)
1283 C
1284  xp=xparam(v)
1285  xpi=1.0_r8-xp
1286 C
1287  sum=hyi*(
1288  > xpi*(-fin(3,i,j) +fin(3,i,j+1))+
1289  > xp*(-fin(3,i+1,j)+fin(3,i+1,j+1)))
1290 C
1291  fval(v,iadr)=sum
1292  enddo
1293  endif
1294 c
1295 c-----------------------------------
1296 c access to 6th derivatives
1297 c
1298  else if(ict(1).eq.6) then
1299 c evaluate d6f/dx3dy3 (not continuous)
1300  iadr=iadr+1
1301  j=jj
1302  do v=1,ivec
1303  i=ii(v)
1304  sum=hxi(v)*hyi*(
1305  > +( fin(3,i,j) -fin(3,i,j+1))
1306  > +(-fin(3,i+1,j)+fin(3,i+1,j+1)))
1307  fval(v,iadr)=sum
1308  enddo
1309  endif
1310 c
1311  return
1312  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17