V3FIT
herm2ev.f
1  subroutine herm2ev(xget,yget,x,nx,y,ny,ilinx,iliny,
2  > f,inf2,ict,fval,ier)
3 C
4 C evaluate a 2d cubic Hermite interpolant on a rectilinear
5 C grid -- this is C1 in both directions.
6 C
7 C this subroutine calls two subroutines:
8 C herm2xy -- find cell containing (xget,yget)
9 C herm2fcn -- evaluate interpolant function and (optionally) derivatives
10 C
11 C input arguments:
12 C ================
13 C
14  real xget,yget ! target of this interpolation
15  real x(nx) ! ordered x grid
16  real y(ny) ! ordered y grid
17  integer ilinx ! ilinx=1 => assume x evenly spaced
18  integer iliny ! iliny=1 => assume y evenly spaced
19 C
20  real f(0:3,inf2,ny) ! function data
21 C
22 C f 2nd dimension inf2 must be .ge. nx
23 C contents of f:
24 C
25 C f(0,i,j) = f @ x(i),y(j)
26 C f(1,i,j) = df/dx @ x(i),y(j)
27 C f(2,i,j) = df/dy @ x(i),y(j)
28 C f(3,i,j) = d2f/dxdy @ x(i),y(j)
29 C
30  integer ict(4) ! code specifying output desired
31 C
32 C ict(1)=1 -- return f (0, don't)
33 C ict(2)=1 -- return df/dx (0, don't)
34 C ict(3)=1 -- return df/dy (0, don't)
35 C ict(4)=1 -- return d2f/dxdy (0, don't)
36 C
37 C output arguments:
38 C =================
39 C
40  real fval(*) ! output data
41  integer ier ! error code =0 ==> no error
42 C
43 C fval(1) receives the first output (depends on ict(...) spec)
44 C fval(2) receives the second output (depends on ict(...) spec)
45 C fval(3) receives the third output (depends on ict(...) spec)
46 C fval(4) receives the fourth output (depends on ict(...) spec)
47 C
48 C examples:
49 C on input ict = [1,1,1,1]
50 C on output fval= [f,df/dx,df/dy,d2f/dxdy]
51 C
52 C on input ict = [1,0,0,0]
53 C on output fval= [f] ... elements 2 & 3 & 4 never referenced
54 C
55 C on input ict = [0,1,1,0]
56 C on output fval= [df/dx,df/dy] ... element 3 & 4 never referenced
57 C
58 C on input ict = [0,0,1,0]
59 C on output fval= [df/dy] ... elements 2 & 3 & 4 never referenced.
60 C
61 C ier -- completion code: 0 means OK
62 C-------------------
63 C local:
64 C
65  integer i,j ! cell indices
66 C
67 C normalized displacement from (x(i),y(j)) corner of cell.
68 C xparam=0 @x(i) xparam=1 @x(i+1)
69 C yparam=0 @y(j) yparam=1 @y(j+1)
70 C
71  real xparam,yparam
72 C
73 C cell dimensions and
74 C inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j))
75 C
76  real hx,hy
77  real hxi,hyi
78 C
79 C 0 .le. xparam .le. 1
80 C 0 .le. yparam .le. 1
81 C
82 C---------------------------------------------------------------------
83 C
84  call herm2xy(xget,yget,x,nx,y,ny,ilinx,iliny,
85  > i,j,xparam,yparam,hx,hxi,hy,hyi,ier)
86  if(ier.ne.0) return
87 c
88  call herm2fcn(ict,1,1,
89  > fval,i,j,xparam,yparam,hx,hxi,hy,hyi,f,inf2,ny)
90 C
91  return
92  end
93 C---------------------------------------------------------------------
94 c herm2xy -- look up x-y zone
95 c
96 c this is the "first part" of herm2ev, see comments, above.
97 c
98  subroutine herm2xy(xget,yget,x,nx,y,ny,ilinx,iliny,
99  > i,j,xparam,yparam,hx,hxi,hy,hyi,ier)
100 c
101 c input of herm2xy
102 c ================
103 c
104  integer nx,ny ! array dimensions
105 c
106  real xget,yget ! target point
107  real x(nx),y(ny) ! indep. coords., strict ascending
108 c
109  integer ilinx ! =1: x evenly spaced
110  integer iliny ! =1: y evenly spaced
111 c
112 c output of herm2xy
113 c =================
114  integer i,j ! index to cell containing target pt
115 c on exit: 1.le.i.le.nx-1 1.le.j.le.ny-1
116 c
117 c normalized position w/in (i,j) cell (btw 0 and 1):
118 c
119  real xparam ! (xget-x(i))/(x(i+1)-x(i))
120  real yparam ! (yget-y(j))/(y(j+1)-y(j))
121 c
122 c grid spacing
123 c
124  real hx ! hx = x(i+1)-x(i)
125  real hy ! hy = y(j+1)-y(j)
126 c
127 c inverse grid spacing:
128 c
129  real hxi ! 1/hx = 1/(x(i+1)-x(i))
130  real hyi ! 1/hy = 1/(y(j+1)-y(j))
131 c
132  integer ier ! return ier.ne.0 on error
133 c
134 c------------------------------------
135 c
136  ier=0
137 c
138 c range check
139 c
140  zxget=xget
141  zyget=yget
142  if((xget.lt.x(1)).or.(xget.gt.x(nx))) then
143  zxtol=4.0e-7*max(abs(x(1)),abs(x(nx)))
144  if((xget.lt.x(1)-zxtol).or.(xget.gt.x(nx)+zxtol)) then
145  ier=1
146  write(6,1001) xget,x(1),x(nx)
147  1001 format(' ?herm2ev: xget=',1pe11.4,' out of range ',
148  > 1pe11.4,' to ',1pe11.4)
149  else
150  if((xget.lt.x(1)-0.5*zxtol).or.
151  > (xget.gt.x(nx)+0.5*zxtol))
152  > write(6,1011) xget,x(1),x(nx)
153  1011 format(' %herm2ev: xget=',1pe15.8,' beyond range ',
154  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
155  if(xget.lt.x(1)) then
156  zxget=x(1)
157  else
158  zxget=x(nx)
159  endif
160  endif
161  endif
162  if((yget.lt.y(1)).or.(yget.gt.y(ny))) then
163  zytol=4.0e-7*max(abs(y(1)),abs(y(ny)))
164  if((yget.lt.y(1)-zytol).or.(yget.gt.y(ny)+zytol)) then
165  ier=1
166  write(6,1002) yget,y(1),y(ny)
167  1002 format(' ?herm2ev: yget=',1pe11.4,' out of range ',
168  > 1pe11.4,' to ',1pe11.4)
169  else
170  if((yget.lt.y(1)-0.5*zytol).or.
171  > (yget.gt.y(ny)+0.5*zytol))
172  > write(6,1012) yget,y(1),y(ny)
173  1012 format(' %herm2ev: yget=',1pe15.8,' beyond range ',
174  > 1pe15.8,' to ',1pe15.8,' (fixup applied)')
175  if(yget.lt.y(1)) then
176  zyget=y(1)
177  else
178  zyget=y(ny)
179  endif
180  endif
181  endif
182  if(ier.ne.0) return
183 c
184 c now find interval in which target point lies..
185 c
186  nxm=nx-1
187  nym=ny-1
188 c
189  if(ilinx.eq.1) then
190  ii=1+nxm*(zxget-x(1))/(x(nx)-x(1))
191  i=min(nxm, ii)
192  if(zxget.lt.x(i)) then
193  i=i-1
194  else if(zxget.gt.x(i+1)) then
195  i=i+1
196  endif
197  else
198  if((1.le.i).and.(i.lt.nxm)) then
199  if((x(i).le.zxget).and.(zxget.le.x(i+1))) then
200  continue ! already have the zone
201  else
202  call zonfind(x,nx,zxget,i)
203  endif
204  else
205  i=nx/2
206  call zonfind(x,nx,zxget,i)
207  endif
208  endif
209 c
210  if(iliny.eq.1) then
211  jj=1+nym*(zyget-y(1))/(y(ny)-y(1))
212  j=min(nym, jj)
213  if(zyget.lt.y(j)) then
214  j=j-1
215  else if(zyget.gt.y(j+1)) then
216  j=j+1
217  endif
218  else
219  if((1.le.j).and.(j.lt.nym)) then
220  if((y(j).le.zyget).and.(zyget.le.y(j+1))) then
221  continue ! already have the zone
222  else
223  call zonfind(y,ny,zyget,j)
224  endif
225  else
226  j=ny/2
227  call zonfind(y,ny,zyget,j)
228  endif
229  endif
230 c
231  hx=(x(i+1)-x(i))
232  hy=(y(j+1)-y(j))
233 c
234  hxi=1.0/hx
235  hyi=1.0/hy
236 c
237  xparam=(zxget-x(i))*hxi
238  yparam=(zyget-y(j))*hyi
239 c
240  return
241  end
242 C---------------------------------------------------------------------
243 C evaluate C1 cubic Hermite function interpolation -- 2d fcn
244 C --vectorized-- dmc 10 Feb 1999
245 C
246  subroutine herm2fcn(ict,ivec,ivecd,
247  > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
248  > fin,inf2,ny)
249 C
250  integer ict(4) ! requested output control
251  integer ivec ! vector length
252  integer ivecd ! vector dimension (1st dim of fval)
253 C
254  integer ii(ivec),jj(ivec) ! target cells (i,j)
255  real xparam(ivec),yparam(ivec)
256  ! normalized displacements from (i,j) corners
257 C
258  real hx(ivec),hy(ivec) ! grid spacing, and
259  real hxi(ivec),hyi(ivec) ! inverse grid spacing 1/(x(i+1)-x(i))
260  ! & 1/(y(j+1)-y(j))
261 C
262  real fin(0:3,inf2,ny) ! interpolant data (cf "herm2ev")
263 C
264  real fval(ivecd,*) ! output returned
265 C
266 C for detailed description of fin, ict and fval see subroutine
267 C herm2ev comments. Note ict is not vectorized; the same output
268 C is expected to be returned for all input vector data points.
269 C
270 C note that the index inputs ii,jj and parameter inputs
271 C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
272 C output array fval has a vector ** 1st dimension ** whose
273 C size must be given as a separate argument
274 C
275 C to use this routine in scalar mode, pass in ivec=ivecd=1
276 C
277 C---------------
278 C Hermite cubic basis functions
279 C -->for function value matching
280 C a(0)=0 a(1)=1 a'(0)=0 a'(1)=0
281 C abar(0)=1 abar(1)=0 abar'(0)=0 abar'(1)=0
282 C
283 C a(x)=-2*x**3 + 3*x**2 = x*x*(-2.0*x+3.0)
284 C abar(x)=1-a(x)
285 C a'(x)=-abar'(x) = 6.0*x*(1.0-x)
286 C
287 C -->for derivative matching
288 C b(0)=0 b(1)=0 b'(0)=0 b'(1)=1
289 C bbar(0)=0 bbar(1)=0 bbar'(0)=1 bbar'(1)=0
290 C
291 C b(x)=x**3-x**2 b'(x)=3*x**2-2*x
292 C bbar(x)=x**3-2*x**2+x bbar'(x)=3*x**2-4*x+1
293 C
294  real sum
295  integer v
296 C
297  do v=1,ivec
298  i=ii(v)
299  j=jj(v)
300 C
301 C ...in x direction
302 C
303  xp=xparam(v)
304  xpi=1.0-xp
305  xp2=xp*xp
306  xpi2=xpi*xpi
307  ax=xp2*(3.0-2.0*xp)
308  axbar=1.0-ax
309  bx=-xp2*xpi
310  bxbar=xpi2*xp
311 C
312 C ...in y direction
313 C
314  yp=yparam(v)
315  ypi=1.0-yp
316  yp2=yp*yp
317  ypi2=ypi*ypi
318  ay=yp2*(3.0-2.0*yp)
319  aybar=1.0-ay
320  by=-yp2*ypi
321  bybar=ypi2*yp
322 C
323 C ...derivatives...
324 C
325  axp=6.0*xp*xpi
326  axbarp=-axp
327  bxp=xp*(3.0*xp-2.0)
328  bxbarp=xpi*(3.0*xpi-2.0)
329 C
330  ayp=6.0*yp*ypi
331  aybarp=-ayp
332  byp=yp*(3.0*yp-2.0)
333  bybarp=ypi*(3.0*ypi-2.0)
334 C
335  iadr=0
336 C
337 C get desired values:
338 C
339  if(ict(1).eq.1) then
340 C
341 C function value:
342 C
343  iadr=iadr+1
344  sum=axbar*(aybar*fin(0,i,j) +ay*fin(0,i,j+1))+
345  > ax*(aybar*fin(0,i+1,j)+ay*fin(0,i+1,j+1))
346 C
347  sum=sum+hx(v)*(
348  > bxbar*(aybar*fin(1,i,j) +ay*fin(1,i,j+1))+
349  > bx*(aybar*fin(1,i+1,j)+ay*fin(1,i+1,j+1)))
350 C
351  sum=sum+hy(v)*(
352  > axbar*(bybar*fin(2,i,j) +by*fin(2,i,j+1))+
353  > ax*(bybar*fin(2,i+1,j)+by*fin(2,i+1,j+1)))
354 C
355  sum=sum+hx(v)*hy(v)*(
356  > bxbar*(bybar*fin(3,i,j) +by*fin(3,i,j+1))+
357  > bx*(bybar*fin(3,i+1,j)+by*fin(3,i+1,j+1)))
358 C
359  fval(v,iadr)=sum
360  endif
361 C
362  if(ict(2).eq.1) then
363 C
364 C df/dx:
365 C
366  iadr=iadr+1
367 C
368  sum=hxi(v)*(
369  > axbarp*(aybar*fin(0,i,j) +ay*fin(0,i,j+1))+
370  > axp*(aybar*fin(0,i+1,j)+ay*fin(0,i+1,j+1)))
371 C
372  sum=sum+
373  > bxbarp*(aybar*fin(1,i,j) +ay*fin(1,i,j+1))+
374  > bxp*(aybar*fin(1,i+1,j)+ay*fin(1,i+1,j+1))
375 C
376  sum=sum+hxi(v)*hy(v)*(
377  > axbarp*(bybar*fin(2,i,j) +by*fin(2,i,j+1))+
378  > axp*(bybar*fin(2,i+1,j)+by*fin(2,i+1,j+1)))
379 C
380  sum=sum+hy(v)*(
381  > bxbarp*(bybar*fin(3,i,j) +by*fin(3,i,j+1))+
382  > bxp*(bybar*fin(3,i+1,j)+by*fin(3,i+1,j+1)))
383 C
384  fval(v,iadr)=sum
385  endif
386 C
387  if(ict(3).eq.1) then
388 C
389 C df/dy:
390 C
391  iadr=iadr+1
392 C
393  sum=hyi(v)*(
394  > axbar*(aybarp*fin(0,i,j) +ayp*fin(0,i,j+1))+
395  > ax*(aybarp*fin(0,i+1,j)+ayp*fin(0,i+1,j+1)))
396 C
397  sum=sum+hx(v)*hyi(v)*(
398  > bxbar*(aybarp*fin(1,i,j) +ayp*fin(1,i,j+1))+
399  > bx*(aybarp*fin(1,i+1,j)+ayp*fin(1,i+1,j+1)))
400 C
401  sum=sum+
402  > axbar*(bybarp*fin(2,i,j) +byp*fin(2,i,j+1))+
403  > ax*(bybarp*fin(2,i+1,j)+byp*fin(2,i+1,j+1))
404 C
405  sum=sum+hx(v)*(
406  > bxbar*(bybarp*fin(3,i,j) +byp*fin(3,i,j+1))+
407  > bx*(bybarp*fin(3,i+1,j)+byp*fin(3,i+1,j+1)))
408 C
409  fval(v,iadr)=sum
410  endif
411 C
412  if(ict(4).eq.1) then
413 C
414 C d2f/dxdy:
415 C
416  iadr=iadr+1
417 C
418  sum=hxi(v)*hyi(v)*(
419  > axbarp*(aybarp*fin(0,i,j) +ayp*fin(0,i,j+1))+
420  > axp*(aybarp*fin(0,i+1,j)+ayp*fin(0,i+1,j+1)))
421 C
422  sum=sum+hyi(v)*(
423  > bxbarp*(aybarp*fin(1,i,j) +ayp*fin(1,i,j+1))+
424  > bxp*(aybarp*fin(1,i+1,j)+ayp*fin(1,i+1,j+1)))
425 C
426  sum=sum+hxi(v)*(
427  > axbarp*(bybarp*fin(2,i,j) +byp*fin(2,i,j+1))+
428  > axp*(bybarp*fin(2,i+1,j)+byp*fin(2,i+1,j+1)))
429 C
430  sum=sum+
431  > bxbarp*(bybarp*fin(3,i,j) +byp*fin(3,i,j+1))+
432  > bxp*(bybarp*fin(3,i+1,j)+byp*fin(3,i+1,j+1))
433 C
434  fval(v,iadr)=sum
435  endif
436 C
437  enddo ! vector loop
438 C
439  return
440  end
441 C---------------------------------------------------------------------
442 C evaluate C1 cubic Hermite function interpolation -- 2d fcn
443 C --vectorized-- dmc 10 Feb 1999
444 C --optimized for VARIATION along x axis ONLY--
445 C
446  subroutine herm2fcnx(ict,ivec,ivecd,
447  > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
448  > fin,inf2,ny)
449 C
450  integer ict(4) ! requested output control
451  integer ivec ! vector length
452  integer ivecd ! vector dimension (1st dim of fval)
453 C
454  integer ii(ivec),jj ! target cells (i,j)
455  real xparam(ivec),yparam
456  ! normalized displacements from (i,j) corners
457 C
458  real hx(ivec),hy ! grid spacing, and
459  real hxi(ivec),hyi ! inverse grid spacing 1/(x(i+1)-x(i))
460  ! & 1/(y(j+1)-y(j))
461 C
462  real fin(0:3,inf2,ny) ! interpolant data (cf "herm2ev")
463 C
464  real fval(ivecd,*) ! output returned
465 C
466 C for detailed description of fin, ict and fval see subroutine
467 C herm2ev comments. Note ict is not vectorized; the same output
468 C is expected to be returned for all input vector data points.
469 C
470 C note that the index inputs ii,jj and parameter inputs
471 C xparam,yparam,hx,hxi,hy,hyi are vectorized, and the
472 C output array fval has a vector ** 1st dimension ** whose
473 C size must be given as a separate argument
474 C
475 C to use this routine in scalar mode, pass in ivec=ivecd=1
476 C
477 C---------------
478 C Hermite cubic basis functions
479 C -->for function value matching
480 C a(0)=0 a(1)=1 a'(0)=0 a'(1)=0
481 C abar(0)=1 abar(1)=0 abar'(0)=0 abar'(1)=0
482 C
483 C a(x)=-2*x**3 + 3*x**2 = x*x*(-2.0*x+3.0)
484 C abar(x)=1-a(x)
485 C a'(x)=-abar'(x) = 6.0*x*(1.0-x)
486 C
487 C -->for derivative matching
488 C b(0)=0 b(1)=0 b'(0)=0 b'(1)=1
489 C bbar(0)=0 bbar(1)=0 bbar'(0)=1 bbar'(1)=0
490 C
491 C b(x)=x**3-x**2 b'(x)=3*x**2-2*x
492 C bbar(x)=x**3-2*x**2+x bbar'(x)=3*x**2-4*x+1
493 C
494  real sum
495  integer v
496 C
497  j=jj
498 C
499 C ...in y direction
500 C
501  yp=yparam
502  ypi=1.0-yp
503  yp2=yp*yp
504  ypi2=ypi*ypi
505  ay=yp2*(3.0-2.0*yp)
506  aybar=1.0-ay
507  by=-yp2*ypi
508  bybar=ypi2*yp
509 C
510 C ...derivatives...
511 C
512  ayp=6.0*yp*ypi
513  aybarp=-ayp
514  byp=yp*(3.0*yp-2.0)
515  bybarp=ypi*(3.0*ypi-2.0)
516 C
517  do v=1,ivec
518  i=ii(v)
519 C
520 C ...in x direction
521 C
522  xp=xparam(v)
523  xpi=1.0-xp
524  xp2=xp*xp
525  xpi2=xpi*xpi
526  ax=xp2*(3.0-2.0*xp)
527  axbar=1.0-ax
528  bx=-xp2*xpi
529  bxbar=xpi2*xp
530 C
531 C ...derivatives...
532 C
533  axp=6.0*xp*xpi
534  axbarp=-axp
535  bxp=xp*(3.0*xp-2.0)
536  bxbarp=xpi*(3.0*xpi-2.0)
537 C
538  iadr=0
539 C
540 C get desired values:
541 C
542  if(ict(1).eq.1) then
543 C
544 C function value:
545 C
546  iadr=iadr+1
547  sum=axbar*(aybar*fin(0,i,j) +ay*fin(0,i,j+1))+
548  > ax*(aybar*fin(0,i+1,j)+ay*fin(0,i+1,j+1))
549 C
550  sum=sum+hx(v)*(
551  > bxbar*(aybar*fin(1,i,j) +ay*fin(1,i,j+1))+
552  > bx*(aybar*fin(1,i+1,j)+ay*fin(1,i+1,j+1)))
553 C
554  sum=sum+hy*(
555  > axbar*(bybar*fin(2,i,j) +by*fin(2,i,j+1))+
556  > ax*(bybar*fin(2,i+1,j)+by*fin(2,i+1,j+1)))
557 C
558  sum=sum+hx(v)*hy*(
559  > bxbar*(bybar*fin(3,i,j) +by*fin(3,i,j+1))+
560  > bx*(bybar*fin(3,i+1,j)+by*fin(3,i+1,j+1)))
561 C
562  fval(v,iadr)=sum
563  endif
564 C
565  if(ict(2).eq.1) then
566 C
567 C df/dx:
568 C
569  iadr=iadr+1
570 C
571  sum=hxi(v)*(
572  > axbarp*(aybar*fin(0,i,j) +ay*fin(0,i,j+1))+
573  > axp*(aybar*fin(0,i+1,j)+ay*fin(0,i+1,j+1)))
574 C
575  sum=sum+
576  > bxbarp*(aybar*fin(1,i,j) +ay*fin(1,i,j+1))+
577  > bxp*(aybar*fin(1,i+1,j)+ay*fin(1,i+1,j+1))
578 C
579  sum=sum+hxi(v)*hy*(
580  > axbarp*(bybar*fin(2,i,j) +by*fin(2,i,j+1))+
581  > axp*(bybar*fin(2,i+1,j)+by*fin(2,i+1,j+1)))
582 C
583  sum=sum+hy*(
584  > bxbarp*(bybar*fin(3,i,j) +by*fin(3,i,j+1))+
585  > bxp*(bybar*fin(3,i+1,j)+by*fin(3,i+1,j+1)))
586 C
587  fval(v,iadr)=sum
588  endif
589 C
590  if(ict(3).eq.1) then
591 C
592 C df/dy:
593 C
594  iadr=iadr+1
595 C
596  sum=hyi*(
597  > axbar*(aybarp*fin(0,i,j) +ayp*fin(0,i,j+1))+
598  > ax*(aybarp*fin(0,i+1,j)+ayp*fin(0,i+1,j+1)))
599 C
600  sum=sum+hx(v)*hyi*(
601  > bxbar*(aybarp*fin(1,i,j) +ayp*fin(1,i,j+1))+
602  > bx*(aybarp*fin(1,i+1,j)+ayp*fin(1,i+1,j+1)))
603 C
604  sum=sum+
605  > axbar*(bybarp*fin(2,i,j) +byp*fin(2,i,j+1))+
606  > ax*(bybarp*fin(2,i+1,j)+byp*fin(2,i+1,j+1))
607 C
608  sum=sum+hx(v)*(
609  > bxbar*(bybarp*fin(3,i,j) +byp*fin(3,i,j+1))+
610  > bx*(bybarp*fin(3,i+1,j)+byp*fin(3,i+1,j+1)))
611 C
612  fval(v,iadr)=sum
613  endif
614 C
615  if(ict(4).eq.1) then
616 C
617 C d2f/dxdy:
618 C
619  iadr=iadr+1
620 C
621  sum=hxi(v)*hyi*(
622  > axbarp*(aybarp*fin(0,i,j) +ayp*fin(0,i,j+1))+
623  > axp*(aybarp*fin(0,i+1,j)+ayp*fin(0,i+1,j+1)))
624 C
625  sum=sum+hyi*(
626  > bxbarp*(aybarp*fin(1,i,j) +ayp*fin(1,i,j+1))+
627  > bxp*(aybarp*fin(1,i+1,j)+ayp*fin(1,i+1,j+1)))
628 C
629  sum=sum+hxi(v)*(
630  > axbarp*(bybarp*fin(2,i,j) +byp*fin(2,i,j+1))+
631  > axp*(bybarp*fin(2,i+1,j)+byp*fin(2,i+1,j+1)))
632 C
633  sum=sum+
634  > bxbarp*(bybarp*fin(3,i,j) +byp*fin(3,i,j+1))+
635  > bxp*(bybarp*fin(3,i+1,j)+byp*fin(3,i+1,j+1))
636 C
637  fval(v,iadr)=sum
638  endif
639 C
640  enddo ! vector loop
641 C
642  return
643  end
Generated on Thu Mar 5 2020 15:49:24 for V3FIT by   doxygen 1.8.17