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