1 subroutine r8evtricub(xget,yget,zget,x,nx,y,ny,z,nz,
3 > f,inf2,inf3,ict,fval,ier)
20 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
21 INTEGER ny,nz,inf2,inf3,nx
31 real*8 f(0:7,inf2,inf3,nz)
115 integer i(1),j(1),k(1)
122 real*8 xparam(1),yparam(1),zparam(1)
127 real*8 hx(1),hy(1),hz(1)
128 real*8 hxi(1),hyi(1),hzi(1)
140 call r8herm3xyz(xget,yget,zget,x,nx,y,ny,z,nz,ilinx,iliny,ilinz,
141 > i(1),j(1),k(1),xparam(1),yparam(1),zparam(1),
142 > hx(1),hxi(1),hy(1),hyi(1),hz(1),hzi(1),ier)
145 call r8fvtricub(ict,1,1,
146 > fval,i,j,k,xparam,yparam,zparam,
147 > hx,hxi,hy,hyi,hz,hzi,
156 subroutine r8fvtricub(ict,ivec,ivecd,
157 > fval,ii,jj,kk,xparam,yparam,zparam,
158 > hx,hxi,hy,hyi,hz,hzi,
166 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
167 INTEGER inf3,nz,inf2,iadr,i,j,k
170 real*8 z36th,z216th,xp,xpi,xp2,xpi2,cx,cxi,hx2,yp,ypi,yp2
171 real*8 ypi2,cy,cyi,hy2,zp,zpi,zp2,zpi2,cz,czi,hz2,cxd,cxdi
172 real*8 cyd,cydi,czd,czdi
178 integer ii(ivec),jj(ivec),kk(ivec)
179 real*8 xparam(ivec),yparam(ivec),zparam(ivec)
182 real*8 hx(ivec),hy(ivec),hz(ivec)
183 real*8 hxi(ivec),hyi(ivec),hzi(ivec)
186 real*8 fin(0:7,inf2,inf3,nz)
206 REAL*8,
parameter :: sixth = 0.166666666666666667_r8
211 z216th=sixth*sixth*sixth
214 if(abs(ict(1)).le.2)
then
238 cxi=xpi*(xpi2-1.0_r8)
249 cyi=ypi*(ypi2-1.0_r8)
260 czi=zpi*(zpi2-1.0_r8)
265 > xpi*(ypi*fin(0,i,j,k) +yp*fin(0,i,j+1,k))+
266 > xp*(ypi*fin(0,i+1,j,k)+yp*fin(0,i+1,j+1,k)))
268 > xpi*(ypi*fin(0,i,j,k+1) +yp*fin(0,i,j+1,k+1))+
269 > xp*(ypi*fin(0,i+1,j,k+1)+yp*fin(0,i+1,j+1,k+1))))
273 > cxi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
274 > cx*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
276 > cxi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
277 > cx*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
281 > xpi*(cyi*fin(2,i,j,k) +cy*fin(2,i,j+1,k))+
282 > xp*(cyi*fin(2,i+1,j,k)+cy*fin(2,i+1,j+1,k)))
284 > xpi*(cyi*fin(2,i,j,k+1) +cy*fin(2,i,j+1,k+1))+
285 > xp*(cyi*fin(2,i+1,j,k+1)+cy*fin(2,i+1,j+1,k+1))))
289 > xpi*(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))+
290 > xp*(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
292 > xpi*(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))+
293 > xp*(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
295 sum=sum+z36th*hx2*hy2*(
297 > cxi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
298 > cx*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
300 > cxi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
301 > cx*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
303 sum=sum+z36th*hx2*hz2*(
305 > cxi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
306 > cx*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
308 > cxi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
309 > cx*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
311 sum=sum+z36th*hy2*hz2*(
313 > xpi*(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))+
314 > xp*(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
316 > xpi*(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))+
317 > xp*(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
319 sum=sum+z216th*hx2*hy2*hz2*(
321 > cxi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
322 > cx*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
324 > cxi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
325 > cx*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
348 cxd=3.0_r8*xp2-1.0_r8
349 cxdi=-3.0_r8*xpi2+1.0_r8
359 cyi=ypi*(ypi2-1.0_r8)
370 czi=zpi*(zpi2-1.0_r8)
375 > -(ypi*fin(0,i,j,k) +yp*fin(0,i,j+1,k))
376 > +(ypi*fin(0,i+1,j,k)+yp*fin(0,i+1,j+1,k)))
378 > -(ypi*fin(0,i,j,k+1) +yp*fin(0,i,j+1,k+1))
379 > +(ypi*fin(0,i+1,j,k+1)+yp*fin(0,i+1,j+1,k+1))))
381 sum=sum+sixth*hx(v)*(
383 > cxdi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
384 > cxd*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
386 > cxdi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
387 > cxd*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
389 sum=sum+sixth*hxi(v)*hy2*(
391 > -(cyi*fin(2,i,j,k) +cy*fin(2,i,j+1,k))
392 > +(cyi*fin(2,i+1,j,k)+cy*fin(2,i+1,j+1,k)))
394 > -(cyi*fin(2,i,j,k+1) +cy*fin(2,i,j+1,k+1))
395 > +(cyi*fin(2,i+1,j,k+1)+cy*fin(2,i+1,j+1,k+1))))
397 sum=sum+sixth*hxi(v)*hz2*(
399 > -(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))
400 > +(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
402 > -(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))
403 > +(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
405 sum=sum+z36th*hx(v)*hy2*(
407 > cxdi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
408 > cxd*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
410 > cxdi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
411 > cxd*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
413 sum=sum+z36th*hx(v)*hz2*(
415 > cxdi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
416 > cxd*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
418 > cxdi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
419 > cxd*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
421 sum=sum+z36th*hxi(v)*hy2*hz2*(
423 > -(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))
424 > +(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
426 > -(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))
427 > +(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
429 sum=sum+z216th*hx(v)*hy2*hz2*(
431 > cxdi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
432 > cxd*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
434 > cxdi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
435 > cxd*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
459 cxi=xpi*(xpi2-1.0_r8)
469 cyd=3.0_r8*yp2-1.0_r8
470 cydi=-3.0_r8*ypi2+1.0_r8
480 czi=zpi*(zpi2-1.0_r8)
485 > xpi*(-fin(0,i,j,k) +fin(0,i,j+1,k))+
486 > xp*(-fin(0,i+1,j,k)+fin(0,i+1,j+1,k)))
488 > xpi*(-fin(0,i,j,k+1) +fin(0,i,j+1,k+1))+
489 > xp*(-fin(0,i+1,j,k+1)+fin(0,i+1,j+1,k+1))))
491 sum=sum+sixth*hyi(v)*hx2*(
493 > cxi*(-fin(1,i,j,k) +fin(1,i,j+1,k))+
494 > cx*(-fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
496 > cxi*(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))+
497 > cx*(-fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
499 sum=sum+sixth*hy(v)*(
501 > xpi*(cydi*fin(2,i,j,k) +cyd*fin(2,i,j+1,k))+
502 > xp*(cydi*fin(2,i+1,j,k)+cyd*fin(2,i+1,j+1,k)))
504 > xpi*(cydi*fin(2,i,j,k+1) +cyd*fin(2,i,j+1,k+1))+
505 > xp*(cydi*fin(2,i+1,j,k+1)+cyd*fin(2,i+1,j+1,k+1))))
507 sum=sum+sixth*hyi(v)*hz2*(
509 > xpi*(-fin(3,i,j,k) +fin(3,i,j+1,k))+
510 > xp*(-fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
512 > xpi*(-fin(3,i,j,k+1) +fin(3,i,j+1,k+1))+
513 > xp*(-fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
515 sum=sum+z36th*hx2*hy(v)*(
517 > cxi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
518 > cx*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
520 > cxi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
521 > cx*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
523 sum=sum+z36th*hyi(v)*hx2*hz2*(
525 > cxi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
526 > cx*(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
528 > cxi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
529 > cx*(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
531 sum=sum+z36th*hy(v)*hz2*(
533 > xpi*(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))+
534 > xp*(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
536 > xpi*(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))+
537 > xp*(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
539 sum=sum+z216th*hx2*hy(v)*hz2*(
541 > cxi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
542 > cx*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
544 > cxi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
545 > cx*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
569 cxi=xpi*(xpi2-1.0_r8)
580 cyi=ypi*(ypi2-1.0_r8)
590 czd=3.0_r8*zp2-1.0_r8
591 czdi=-3.0_r8*zpi2+1.0_r8
595 > xpi*(ypi*fin(0,i,j,k) +yp*fin(0,i,j+1,k))+
596 > xp*(ypi*fin(0,i+1,j,k)+yp*fin(0,i+1,j+1,k)))
598 > xpi*(ypi*fin(0,i,j,k+1) +yp*fin(0,i,j+1,k+1))+
599 > xp*(ypi*fin(0,i+1,j,k+1)+yp*fin(0,i+1,j+1,k+1))))
601 sum=sum+sixth*hx2*hzi(v)*(
603 > cxi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
604 > cx*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
606 > cxi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
607 > cx*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
609 sum=sum+sixth*hy2*hzi(v)*(
611 > xpi*(cyi*fin(2,i,j,k) +cy*fin(2,i,j+1,k))+
612 > xp*(cyi*fin(2,i+1,j,k)+cy*fin(2,i+1,j+1,k)))
614 > xpi*(cyi*fin(2,i,j,k+1) +cy*fin(2,i,j+1,k+1))+
615 > xp*(cyi*fin(2,i+1,j,k+1)+cy*fin(2,i+1,j+1,k+1))))
617 sum=sum+sixth*hz(v)*(
619 > xpi*(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))+
620 > xp*(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
622 > xpi*(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))+
623 > xp*(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
625 sum=sum+z36th*hx2*hy2*hzi(v)*(
627 > cxi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
628 > cx*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
630 > cxi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
631 > cx*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
633 sum=sum+z36th*hx2*hz(v)*(
635 > cxi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
636 > cx*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
638 > cxi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
639 > cx*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
641 sum=sum+z36th*hy2*hz(v)*(
643 > xpi*(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))+
644 > xp*(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
646 > xpi*(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))+
647 > xp*(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
649 sum=sum+z216th*hx2*hy2*hz(v)*(
651 > cxi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
652 > cx*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
654 > cxi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
655 > cx*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
684 cyi=ypi*(ypi2-1.0_r8)
695 czi=zpi*(zpi2-1.0_r8)
700 > xpi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
701 > xp*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
703 > xpi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
704 > xp*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
708 > xpi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
709 > xp*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
711 > xpi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
712 > xp*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
716 > xpi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
717 > xp*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
719 > xpi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
720 > xp*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
722 sum=sum+z36th*hy2*hz2*(
724 > xpi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
725 > xp*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
727 > xpi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
728 > xp*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
752 cxi=xpi*(xpi2-1.0_r8)
768 czi=zpi*(zpi2-1.0_r8)
773 > xpi*(ypi*fin(2,i,j,k) +yp*fin(2,i,j+1,k))+
774 > xp*(ypi*fin(2,i+1,j,k)+yp*fin(2,i+1,j+1,k)))
776 > xpi*(ypi*fin(2,i,j,k+1) +yp*fin(2,i,j+1,k+1))+
777 > xp*(ypi*fin(2,i+1,j,k+1)+yp*fin(2,i+1,j+1,k+1))))
781 > cxi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
782 > cx*(ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
784 > cxi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
785 > cx*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
789 > xpi*(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))+
790 > xp*(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
792 > xpi*(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))+
793 > xp*(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
795 sum=sum+z36th*hx2*hz2*(
797 > cxi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
798 > cx*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
800 > cxi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
801 > cx*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
825 cxi=xpi*(xpi2-1.0_r8)
836 cyi=ypi*(ypi2-1.0_r8)
846 > xpi*(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))+
847 > xp*(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
849 > xpi*(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))+
850 > xp*(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
854 > cxi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
855 > cx*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
857 > cxi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
858 > cx*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
862 > xpi*(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))+
863 > xp*(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
865 > xpi*(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))+
866 > xp*(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
868 sum=sum+z36th*hx2*hy2*(
870 > cxi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
871 > cx*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
873 > cxi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
874 > cx*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
897 cxd=3.0_r8*xp2-1.0_r8
898 cxdi=-3.0_r8*xpi2+1.0_r8
907 cyd=3.0_r8*yp2-1.0_r8
908 cydi=-3.0_r8*ypi2+1.0_r8
918 czi=zpi*(zpi2-1.0_r8)
923 > (fin(0,i,j,k) -fin(0,i,j+1,k))-
924 > (fin(0,i+1,j,k)-fin(0,i+1,j+1,k)))
926 > (fin(0,i,j,k+1) -fin(0,i,j+1,k+1))-
927 > (fin(0,i+1,j,k+1)-fin(0,i+1,j+1,k+1))))
929 sum=sum+sixth*hyi(v)*hx(v)*(
931 > cxdi*(-fin(1,i,j,k) +fin(1,i,j+1,k))+
932 > cxd*(-fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
934 > cxdi*(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))+
935 > cxd*(-fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
937 sum=sum+sixth*hxi(v)*hy(v)*(
939 > -(cydi*fin(2,i,j,k) +cyd*fin(2,i,j+1,k))
940 > +(cydi*fin(2,i+1,j,k)+cyd*fin(2,i+1,j+1,k)))
942 > -(cydi*fin(2,i,j,k+1) +cyd*fin(2,i,j+1,k+1))
943 > +(cydi*fin(2,i+1,j,k+1)+cyd*fin(2,i+1,j+1,k+1))))
945 sum=sum+sixth*hxi(v)*hyi(v)*hz2*(
947 > (fin(3,i,j,k) -fin(3,i,j+1,k))-
948 > (fin(3,i+1,j,k)-fin(3,i+1,j+1,k)))
950 > (fin(3,i,j,k+1) -fin(3,i,j+1,k+1))-
951 > (fin(3,i+1,j,k+1)-fin(3,i+1,j+1,k+1))))
953 sum=sum+z36th*hx(v)*hy(v)*(
955 > cxdi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
956 > cxd*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
958 > cxdi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
959 > cxd*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
961 sum=sum+z36th*hyi(v)*hx(v)*hz2*(
963 > cxdi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
964 > cxd*(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
966 > cxdi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
967 > cxd*(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
969 sum=sum+z36th*hxi(v)*hy(v)*hz2*(
971 > -(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))
972 > +(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
974 > -(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))
975 > +(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
977 sum=sum+z216th*hx(v)*hy(v)*hz2*(
979 > cxdi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
980 > cxd*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
982 > cxdi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
983 > cxd*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
1006 cxd=3.0_r8*xp2-1.0_r8
1007 cxdi=-3.0_r8*xpi2+1.0_r8
1017 cyi=ypi*(ypi2-1.0_r8)
1027 czd=3.0_r8*zp2-1.0_r8
1028 czdi=-3.0_r8*zpi2+1.0_r8
1032 > (ypi*fin(0,i,j,k) +yp*fin(0,i,j+1,k)) -
1033 > (ypi*fin(0,i+1,j,k)+yp*fin(0,i+1,j+1,k)))
1035 > (ypi*fin(0,i,j,k+1) +yp*fin(0,i,j+1,k+1)) -
1036 > (ypi*fin(0,i+1,j,k+1)+yp*fin(0,i+1,j+1,k+1))))
1038 sum=sum+sixth*hx(v)*hzi(v)*(
1040 > cxdi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
1041 > cxd*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
1043 > cxdi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
1044 > cxd*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
1046 sum=sum+sixth*hxi(v)*hy2*hzi(v)*(
1048 > (cyi*fin(2,i,j,k) +cy*fin(2,i,j+1,k)) -
1049 > (cyi*fin(2,i+1,j,k)+cy*fin(2,i+1,j+1,k)))
1051 > (cyi*fin(2,i,j,k+1) +cy*fin(2,i,j+1,k+1)) -
1052 > (cyi*fin(2,i+1,j,k+1)+cy*fin(2,i+1,j+1,k+1))))
1054 sum=sum+sixth*hxi(v)*hz(v)*(
1056 > -(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))
1057 > +(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
1059 > -(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))
1060 > +(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
1062 sum=sum+z36th*hx(v)*hy2*hzi(v)*(
1064 > cxdi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
1065 > cxd*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
1067 > cxdi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
1068 > cxd*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
1070 sum=sum+z36th*hx(v)*hz(v)*(
1072 > cxdi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
1073 > cxd*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
1075 > cxdi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
1076 > cxd*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
1078 sum=sum+z36th*hxi(v)*hy2*hz(v)*(
1080 > -(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))
1081 > +(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
1083 > -(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))
1084 > +(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
1086 sum=sum+z216th*hx(v)*hy2*hz(v)*(
1088 > cxdi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
1089 > cxd*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
1091 > cxdi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
1092 > cxd*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
1098 if(ict(10).eq.1)
then
1116 cxi=xpi*(xpi2-1.0_r8)
1126 cyd=3.0_r8*yp2-1.0_r8
1127 cydi=-3.0_r8*ypi2+1.0_r8
1136 czd=3.0_r8*zp2-1.0_r8
1137 czdi=-3.0_r8*zpi2+1.0_r8
1141 > xpi*(fin(0,i,j,k) -fin(0,i,j+1,k))+
1142 > xp*(fin(0,i+1,j,k)-fin(0,i+1,j+1,k)))
1144 > xpi*(fin(0,i,j,k+1) -fin(0,i,j+1,k+1))+
1145 > xp*(fin(0,i+1,j,k+1)-fin(0,i+1,j+1,k+1))))
1147 sum=sum+sixth*hyi(v)*hx2*hzi(v)*(
1149 > cxi*(fin(1,i,j,k) -fin(1,i,j+1,k))+
1150 > cx*(fin(1,i+1,j,k)-fin(1,i+1,j+1,k)))
1152 > cxi*(fin(1,i,j,k+1) -fin(1,i,j+1,k+1))+
1153 > cx*(fin(1,i+1,j,k+1)-fin(1,i+1,j+1,k+1))))
1155 sum=sum+sixth*hy(v)*hzi(v)*(
1157 > xpi*(cydi*fin(2,i,j,k) +cyd*fin(2,i,j+1,k))+
1158 > xp*(cydi*fin(2,i+1,j,k)+cyd*fin(2,i+1,j+1,k)))
1160 > xpi*(cydi*fin(2,i,j,k+1) +cyd*fin(2,i,j+1,k+1))+
1161 > xp*(cydi*fin(2,i+1,j,k+1)+cyd*fin(2,i+1,j+1,k+1))))
1163 sum=sum+sixth*hyi(v)*hz(v)*(
1165 > xpi*(-fin(3,i,j,k) +fin(3,i,j+1,k))+
1166 > xp*(-fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
1168 > xpi*(-fin(3,i,j,k+1) +fin(3,i,j+1,k+1))+
1169 > xp*(-fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
1171 sum=sum+z36th*hx2*hy(v)*hzi(v)*(
1173 > cxi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
1174 > cx*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
1176 > cxi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
1177 > cx*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
1179 sum=sum+z36th*hyi(v)*hx2*hz(v)*(
1181 > cxi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
1182 > cx*(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
1184 > cxi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
1185 > cx*(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
1187 sum=sum+z36th*hy(v)*hz(v)*(
1189 > xpi*(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))+
1190 > xp*(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
1192 > xpi*(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))+
1193 > xp*(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
1195 sum=sum+z216th*hx2*hy(v)*hz(v)*(
1197 > cxi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
1198 > cx*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
1200 > cxi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
1201 > cx*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
1210 else if(ict(1).eq.3)
then
1211 if(ict(2).eq.1)
then
1231 cyd=3.0_r8*yp2-1.0_r8
1232 cydi=-3.0_r8*ypi2+1.0_r8
1242 czi=zpi*(zpi2-1.0_r8)
1247 > xpi*(-fin(1,i,j,k) +fin(1,i,j+1,k))+
1248 > xp*( -fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
1250 > xpi*(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))+
1251 > xp*( -fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
1253 sum=sum+sixth*hy(v)*(
1255 > xpi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
1256 > xp*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
1258 > xpi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
1259 > xp*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
1261 sum=sum+sixth*hz2*hyi(v)*(
1263 > xpi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
1264 > xp*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
1266 > xpi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
1267 > xp*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
1269 sum=sum+z36th*hy(v)*hz2*(
1271 > xpi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
1272 > xp*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
1274 > xpi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
1275 > xp*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
1281 if(ict(3).eq.1)
then
1302 cyi=ypi*(ypi2-1.0_r8)
1312 czd=3.0_r8*zp2-1.0_r8
1313 czdi=-3.0_r8*zpi2+1.0_r8
1317 > xpi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
1318 > xp*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
1320 > xpi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
1321 > xp*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
1323 sum=sum+sixth*hy2*hzi(v)*(
1325 > xpi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
1326 > xp*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
1328 > xpi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
1329 > xp*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
1331 sum=sum+sixth*hz(v)*(
1333 > xpi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
1334 > xp*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
1336 > xpi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
1337 > xp*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
1339 sum=sum+z36th*hy2*hz(v)*(
1341 > xpi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
1342 > xp*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
1344 > xpi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
1345 > xp*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
1351 if(ict(4).eq.1)
then
1366 cxd=3.0_r8*xp2-1.0_r8
1367 cxdi=-3.0_r8*xpi2+1.0_r8
1382 czi=zpi*(zpi2-1.0_r8)
1387 > -(ypi*fin(2,i,j,k) +yp*fin(2,i,j+1,k))
1388 > +(ypi*fin(2,i+1,j,k)+yp*fin(2,i+1,j+1,k)))
1390 > -(ypi*fin(2,i,j,k+1) +yp*fin(2,i,j+1,k+1))
1391 > +(ypi*fin(2,i+1,j,k+1)+yp*fin(2,i+1,j+1,k+1))))
1393 sum=sum+sixth*hx(v)*(
1395 > cxdi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
1396 > cxd*(ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
1398 > cxdi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
1399 > cxd*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
1401 sum=sum+sixth*hz2*hxi(v)*(
1403 > -(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))
1404 > +(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
1406 > -(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))
1407 > +(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
1409 sum=sum+z36th*hx(v)*hz2*(
1411 > cxdi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
1412 > cxd*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
1414 > cxdi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
1415 > cxd*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
1421 if(ict(5).eq.1)
then
1436 cxd=3.0_r8*xp2-1.0_r8
1437 cxdi=-3.0_r8*xpi2+1.0_r8
1446 cyd=3.0_r8*yp2-1.0_r8
1447 cydi=-3.0_r8*ypi2+1.0_r8
1456 czd=3.0_r8*zp2-1.0_r8
1457 czdi=-3.0_r8*zpi2+1.0_r8
1459 sum=hxi(v)*hyi(v)*hzi(v)*(
1461 > (fin(0,i,j,k) -fin(0,i,j+1,k))-
1462 > (fin(0,i+1,j,k)-fin(0,i+1,j+1,k)))
1464 > (fin(0,i,j,k+1) -fin(0,i,j+1,k+1))-
1465 > (fin(0,i+1,j,k+1)-fin(0,i+1,j+1,k+1))))
1467 sum=sum+sixth*hyi(v)*hx(v)*hzi(v)*(
1469 > cxdi*(-fin(1,i,j,k) +fin(1,i,j+1,k))+
1470 > cxd*(-fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
1472 > cxdi*(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))+
1473 > cxd*(-fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
1475 sum=sum+sixth*hxi(v)*hy(v)*hzi(v)*(
1477 > -(cydi*fin(2,i,j,k) +cyd*fin(2,i,j+1,k))
1478 > +(cydi*fin(2,i+1,j,k)+cyd*fin(2,i+1,j+1,k)))
1480 > -(cydi*fin(2,i,j,k+1) +cyd*fin(2,i,j+1,k+1))
1481 > +(cydi*fin(2,i+1,j,k+1)+cyd*fin(2,i+1,j+1,k+1))))
1483 sum=sum+sixth*hxi(v)*hyi(v)*hz(v)*(
1485 > (fin(3,i,j,k) -fin(3,i,j+1,k))-
1486 > (fin(3,i+1,j,k)-fin(3,i+1,j+1,k)))
1488 > (fin(3,i,j,k+1) -fin(3,i,j+1,k+1))-
1489 > (fin(3,i+1,j,k+1)-fin(3,i+1,j+1,k+1))))
1491 sum=sum+z36th*hx(v)*hy(v)*hzi(v)*(
1493 > cxdi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
1494 > cxd*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
1496 > cxdi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
1497 > cxd*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
1499 sum=sum+z36th*hyi(v)*hx(v)*hz(v)*(
1501 > cxdi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
1502 > cxd*(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
1504 > cxdi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
1505 > cxd*(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
1507 sum=sum+z36th*hxi(v)*hy(v)*hz(v)*(
1509 > -(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))
1510 > +(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
1512 > -(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))
1513 > +(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
1515 sum=sum+z216th*hx(v)*hy(v)*hz(v)*(
1517 > cxdi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
1518 > cxd*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
1520 > cxdi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
1521 > cxd*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
1527 if(ict(6).eq.1)
then
1542 cxd=3.0_r8*xp2-1.0_r8
1543 cxdi=-3.0_r8*xpi2+1.0_r8
1553 cyi=ypi*(ypi2-1.0_r8)
1564 > -(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))
1565 > +(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
1567 > -(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))
1568 > +(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
1570 sum=sum+sixth*hx(v)*(
1572 > cxdi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
1573 > cxd*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
1575 > cxdi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
1576 > cxd*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
1578 sum=sum+sixth*hy2*hxi(v)*(
1580 > -(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))
1581 > +(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
1583 > -(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))
1584 > +(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
1586 sum=sum+z36th*hx(v)*hy2*(
1588 > cxdi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
1589 > cxd*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
1591 > cxdi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
1592 > cxd*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
1598 if(ict(7).eq.1)
then
1614 cxi=xpi*(xpi2-1.0_r8)
1629 czd=3.0_r8*zp2-1.0_r8
1630 czdi=-3.0_r8*zpi2+1.0_r8
1634 > xpi*(ypi*fin(2,i,j,k) +yp*fin(2,i,j+1,k))+
1635 > xp*(ypi*fin(2,i+1,j,k)+yp*fin(2,i+1,j+1,k)))
1637 > xpi*(ypi*fin(2,i,j,k+1) +yp*fin(2,i,j+1,k+1))+
1638 > xp*(ypi*fin(2,i+1,j,k+1)+yp*fin(2,i+1,j+1,k+1))))
1640 sum=sum+sixth*hx2*hzi(v)*(
1642 > cxi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
1643 > cx*(ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
1645 > cxi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
1646 > cx*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
1648 sum=sum+sixth*hz(v)*(
1650 > xpi*(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))+
1651 > xp*(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
1653 > xpi*(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))+
1654 > xp*(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
1656 sum=sum+z36th*hx2*hz(v)*(
1658 > cxi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
1659 > cx*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
1661 > cxi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
1662 > cx*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
1668 if(ict(8).eq.1)
then
1684 cxi=xpi*(xpi2-1.0_r8)
1694 cyd=3.0_r8*yp2-1.0_r8
1695 cydi=-3.0_r8*ypi2+1.0_r8
1704 > xpi*(-fin(3,i,j,k) +fin(3,i,j+1,k))+
1705 > xp*( -fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
1707 > xpi*(-fin(3,i,j,k+1) +fin(3,i,j+1,k+1))+
1708 > xp*( -fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
1710 sum=sum+sixth*hx2*hyi(v)*(
1712 > cxi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
1713 > cx*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
1715 > cxi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
1716 > cx*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
1718 sum=sum+sixth*hy(v)*(
1720 > xpi*(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))+
1721 > xp*(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
1723 > xpi*(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))+
1724 > xp*(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
1726 sum=sum+z36th*hx2*hy(v)*(
1728 > cxi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
1729 > cx*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
1731 > cxi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
1732 > cx*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
1740 else if(ict(1).eq.-3)
then
1741 if(ict(2).eq.1)
then
1757 cyi=ypi*(ypi2-1.0_r8)
1768 czi=zpi*(zpi2-1.0_r8)
1773 > -(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))
1774 > +(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
1776 > -(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))
1777 > +(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
1779 sum=sum+sixth*hy2*hxi(v)*(
1781 > -(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))
1782 > +(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
1784 > -(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))
1785 > +(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
1787 sum=sum+sixth*hz2*hxi(v)*(
1789 > -(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))
1790 > +(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
1792 > -(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))
1793 > +(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
1795 sum=sum+z36th*hy2*hz2*hxi(v)*(
1797 > -(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))
1798 > +(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
1800 > -(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))
1801 > +(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
1807 if(ict(3).eq.1)
then
1823 cxi=xpi*(xpi2-1.0_r8)
1834 czi=zpi*(zpi2-1.0_r8)
1839 > xpi*(-fin(2,i,j,k) +fin(2,i,j+1,k))+
1840 > xp*( -fin(2,i+1,j,k)+fin(2,i+1,j+1,k)))
1842 > xpi*(-fin(2,i,j,k+1) +fin(2,i,j+1,k+1))+
1843 > xp*( -fin(2,i+1,j,k+1)+fin(2,i+1,j+1,k+1))))
1845 sum=sum+sixth*hx2*hyi(v)*(
1847 > cxi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
1848 > cx*( -fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
1850 > cxi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
1851 > cx*( -fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
1853 sum=sum+sixth*hz2*hyi(v)*(
1855 > xpi*(-fin(6,i,j,k) +fin(6,i,j+1,k))+
1856 > xp*( -fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
1858 > xpi*(-fin(6,i,j,k+1) +fin(6,i,j+1,k+1))+
1859 > xp*( -fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
1861 sum=sum+z36th*hx2*hz2*hyi(v)*(
1863 > cxi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
1864 > cx*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
1866 > cxi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
1867 > cx*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
1873 if(ict(4).eq.1)
then
1889 cxi=xpi*(xpi2-1.0_r8)
1900 cyi=ypi*(ypi2-1.0_r8)
1905 > xpi*(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))+
1906 > xp*(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
1908 > xpi*(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))+
1909 > xp*(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
1911 sum=sum+sixth*hx2*hzi(v)*(
1913 > cxi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
1914 > cx*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
1916 > cxi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
1917 > cx*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
1919 sum=sum+sixth*hy2*hzi(v)*(
1921 > xpi*(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))+
1922 > xp*(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
1924 > xpi*(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))+
1925 > xp*(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
1927 sum=sum+z36th*hx2*hy2*hzi(v)*(
1929 > cxi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
1930 > cx*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
1932 > cxi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
1933 > cx*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
1942 else if(ict(1).eq.4)
then
1943 if(ict(2).eq.1)
then
1969 czi=zpi*(zpi2-1.0_r8)
1974 > xpi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
1975 > xp*( ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
1977 > xpi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
1978 > xp*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
1982 > xpi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
1983 > xp*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
1985 > xpi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
1986 > xp*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
1992 if(ict(3).eq.1)
then
2012 cyd=3.0_r8*yp2-1.0_r8
2013 cydi=-3.0_r8*ypi2+1.0_r8
2022 czd=3.0_r8*zp2-1.0_r8
2023 czdi=-3.0_r8*zpi2+1.0_r8
2027 > xpi*(-fin(1,i,j,k) +fin(1,i,j+1,k))+
2028 > xp*( -fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
2030 > xpi*(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))+
2031 > xp*( -fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
2033 sum=sum+sixth*hy(v)*hzi(v)*(
2035 > xpi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
2036 > xp*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
2038 > xpi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
2039 > xp*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
2041 sum=sum+sixth*hz(v)*hyi(v)*(
2043 > xpi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
2044 > xp*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
2046 > xpi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
2047 > xp*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
2049 sum=sum+z36th*hy(v)*hz(v)*(
2051 > xpi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
2052 > xp*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
2054 > xpi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
2055 > xp*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
2061 if(ict(4).eq.1)
then
2082 cyi=ypi*(ypi2-1.0_r8)
2092 > xpi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
2093 > xp*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
2095 > xpi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
2096 > xp*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
2100 > xpi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
2101 > xp*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
2103 > xpi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
2104 > xp*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
2110 if(ict(5).eq.1)
then
2125 cxd=3.0_r8*xp2-1.0_r8
2126 cxdi=-3.0_r8*xpi2+1.0_r8
2140 czd=3.0_r8*zp2-1.0_r8
2141 czdi=-3.0_r8*zpi2+1.0_r8
2145 > -(ypi*fin(2,i,j,k) +yp*fin(2,i,j+1,k))
2146 > +(ypi*fin(2,i+1,j,k)+yp*fin(2,i+1,j+1,k)))
2148 > -(ypi*fin(2,i,j,k+1) +yp*fin(2,i,j+1,k+1))
2149 > +(ypi*fin(2,i+1,j,k+1)+yp*fin(2,i+1,j+1,k+1))))
2151 sum=sum+sixth*hx(v)*hzi(v)*(
2153 > cxdi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
2154 > cxd*(ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
2156 > cxdi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
2157 > cxd*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
2159 sum=sum+sixth*hz(v)*hxi(v)*(
2161 > -(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))
2162 > +(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
2164 > -(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))
2165 > +(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
2167 sum=sum+z36th*hx(v)*hz(v)*(
2169 > cxdi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
2170 > cxd*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
2172 > cxdi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
2173 > cxd*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
2179 if(ict(6).eq.1)
then
2194 cxd=3.0_r8*xp2-1.0_r8
2195 cxdi=-3.0_r8*xpi2+1.0_r8
2204 cyd=3.0_r8*yp2-1.0_r8
2205 cydi=-3.0_r8*ypi2+1.0_r8
2214 > ( +fin(3,i,j,k) -fin(3,i,j+1,k))
2215 > +( -fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
2217 > ( +fin(3,i,j,k+1) -fin(3,i,j+1,k+1))
2218 > +( -fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
2220 sum=sum+sixth*hx(v)*hyi(v)*(
2222 > cxdi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
2223 > cxd*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
2225 > cxdi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
2226 > cxd*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
2228 sum=sum+sixth*hy(v)*hxi(v)*(
2230 > -(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))
2231 > +(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
2233 > -(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))
2234 > +(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
2236 sum=sum+z36th*hx(v)*hy(v)*(
2238 > cxdi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
2239 > cxd*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
2241 > cxdi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
2242 > cxd*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
2248 if(ict(7).eq.1)
then
2264 cxi=xpi*(xpi2-1.0_r8)
2279 > xpi*(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))+
2280 > xp*( ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
2282 > xpi*(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))+
2283 > xp*(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
2287 > cxi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
2288 > cx*( ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
2290 > cxi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
2291 > cx*( ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
2300 else if(ict(1).eq.-4)
then
2301 if(ict(2).eq.1)
then
2316 cyd=3.0_r8*yp2-1.0_r8
2317 cydi=-3.0_r8*ypi2+1.0_r8
2327 czi=zpi*(zpi2-1.0_r8)
2332 > ( fin(1,i,j,k) -fin(1,i,j+1,k))+
2333 > ( -fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
2335 > ( fin(1,i,j,k+1) -fin(1,i,j+1,k+1))+
2336 > ( -fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
2338 sum=sum+sixth*hy(v)*hxi(v)*(
2340 > -(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
2341 > (cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
2343 > -(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
2344 > (cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
2346 sum=sum+sixth*hz2*hyi(v)*hxi(v)*(
2348 > ( fin(5,i,j,k) -fin(5,i,j+1,k))+
2349 > ( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
2351 > ( fin(5,i,j,k+1) -fin(5,i,j+1,k+1))+
2352 > ( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
2354 sum=sum+z36th*hy(v)*hz2*hxi(v)*(
2356 > -(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
2357 > (cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
2359 > -(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
2360 > (cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
2366 if(ict(3).eq.1)
then
2382 cyi=ypi*(ypi2-1.0_r8)
2392 czd=3.0_r8*zp2-1.0_r8
2393 czdi=-3.0_r8*zpi2+1.0_r8
2397 > +(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))
2398 > -(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
2400 > -(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))
2401 > +(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
2403 sum=sum+sixth*hy2*hzi(v)*hxi(v)*(
2405 > +(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))
2406 > -(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
2408 > -(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))
2409 > +(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
2411 sum=sum+sixth*hz(v)*hxi(v)*(
2413 > -(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))
2414 > +(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
2416 > -(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))
2417 > +(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
2419 sum=sum+z36th*hy2*hz(v)*hxi(v)*(
2421 > -(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))
2422 > +(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
2424 > -(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))
2425 > +(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
2431 if(ict(4).eq.1)
then
2446 cxd=3.0_r8*xp2-1.0_r8
2447 cxdi=-3.0_r8*xpi2+1.0_r8
2457 czi=zpi*(zpi2-1.0_r8)
2462 > ( fin(2,i,j,k) -fin(2,i,j+1,k))
2463 > +(-fin(2,i+1,j,k)+fin(2,i+1,j+1,k)))
2465 > ( fin(2,i,j,k+1) -fin(2,i,j+1,k+1))
2466 > +(-fin(2,i+1,j,k+1)+fin(2,i+1,j+1,k+1))))
2468 sum=sum+sixth*hx(v)*hyi(v)*(
2470 > cxdi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
2471 > cxd*( -fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
2473 > cxdi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
2474 > cxd*( -fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
2476 sum=sum+sixth*hz2*hxi(v)*hyi(v)*(
2478 > ( fin(6,i,j,k) -fin(6,i,j+1,k))
2479 > +(-fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
2481 > ( fin(6,i,j,k+1) -fin(6,i,j+1,k+1))
2482 > +(-fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
2484 sum=sum+z36th*hx(v)*hz2*hyi(v)*(
2486 > cxdi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
2487 > cxd*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
2489 > cxdi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
2490 > cxd*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
2496 if(ict(5).eq.1)
then
2511 cxd=3.0_r8*xp2-1.0_r8
2512 cxdi=-3.0_r8*xpi2+1.0_r8
2522 cyi=ypi*(ypi2-1.0_r8)
2527 > -(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))
2528 > +(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
2530 > -(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))
2531 > +(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
2533 sum=sum+sixth*hx(v)*hzi(v)*(
2535 > cxdi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
2536 > cxd*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
2538 > cxdi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
2539 > cxd*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
2541 sum=sum+sixth*hy2*hxi(v)*hzi(v)*(
2543 > -(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))
2544 > +(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
2546 > -(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))
2547 > +(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
2549 sum=sum+z36th*hx(v)*hy2*hzi(v)*(
2551 > cxdi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
2552 > cxd*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
2554 > cxdi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
2555 > cxd*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
2561 if(ict(6).eq.1)
then
2577 cxi=xpi*(xpi2-1.0_r8)
2587 czd=3.0_r8*zp2-1.0_r8
2588 czdi=-3.0_r8*zpi2+1.0_r8
2592 > xpi*(-fin(2,i,j,k) +fin(2,i,j+1,k))+
2593 > xp*( -fin(2,i+1,j,k)+fin(2,i+1,j+1,k)))
2595 > xpi*(-fin(2,i,j,k+1) +fin(2,i,j+1,k+1))+
2596 > xp*( -fin(2,i+1,j,k+1)+fin(2,i+1,j+1,k+1))))
2598 sum=sum+sixth*hx2*hyi(v)*hzi(v)*(
2600 > cxi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
2601 > cx*( -fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
2603 > cxi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
2604 > cx*( -fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
2606 sum=sum+sixth*hz(v)*hyi(v)*(
2608 > xpi*(-fin(6,i,j,k) +fin(6,i,j+1,k))+
2609 > xp*( -fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
2611 > xpi*(-fin(6,i,j,k+1) +fin(6,i,j+1,k+1))+
2612 > xp*( -fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
2614 sum=sum+z36th*hx2*hz(v)*hyi(v)*(
2616 > cxi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
2617 > cx*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
2619 > cxi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
2620 > cx*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
2626 if(ict(7).eq.1)
then
2642 cxi=xpi*(xpi2-1.0_r8)
2652 cyd=3.0_r8*yp2-1.0_r8
2653 cydi=-3.0_r8*ypi2+1.0_r8
2657 > xpi*(-fin(3,i,j,k) +fin(3,i,j+1,k))+
2658 > xp*( -fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
2660 > xpi*(-fin(3,i,j,k+1) +fin(3,i,j+1,k+1))+
2661 > xp*( -fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
2663 sum=sum+sixth*hx2*hyi(v)*hzi(v)*(
2665 > cxi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
2666 > cx*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
2668 > cxi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
2669 > cx*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
2671 sum=sum+sixth*hy(v)*hzi(v)*(
2673 > xpi*(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))+
2674 > xp*(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
2676 > xpi*(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))+
2677 > xp*(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
2679 sum=sum+z36th*hx2*hy(v)*hzi(v)*(
2681 > cxi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
2682 > cx*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
2684 > cxi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
2685 > cx*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
2694 else if(ict(1).eq.5)
then
2695 if(ict(2).eq.1)
then
2720 czd=3.0_r8*zp2-1.0_r8
2721 czdi=-3.0_r8*zpi2+1.0_r8
2725 > xpi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
2726 > xp*( ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
2728 > xpi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
2729 > xp*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
2731 sum=sum+sixth*hz(v)*(
2733 > xpi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
2734 > xp*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
2736 > xpi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
2737 > xp*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
2743 if(ict(3).eq.1)
then
2763 cyd=3.0_r8*yp2-1.0_r8
2764 cydi=-3.0_r8*ypi2+1.0_r8
2773 > xpi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
2774 > xp*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
2776 > xpi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
2777 > xp*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
2779 sum=sum+sixth*hy(v)*(
2781 > xpi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
2782 > xp*( cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
2784 > xpi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
2785 > xp*( cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
2791 if(ict(4).eq.1)
then
2806 cxd=3.0_r8*xp2-1.0_r8
2807 cxdi=-3.0_r8*xpi2+1.0_r8
2821 > -(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))
2822 > +(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
2824 > -(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))
2825 > +(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
2827 sum=sum+sixth*hx(v)*(
2829 > cxdi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
2830 > cxd*(ypi*fin(7,i+1,j,k) +yp*fin(7,i+1,j+1,k)))
2832 > cxdi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
2833 > cxd*( ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
2842 else if(ict(1).eq.-5)
then
2843 if(ict(2).eq.1)
then
2864 czi=zpi*(zpi2-1.0_r8)
2869 > -(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))
2870 > +( ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
2872 > -(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))
2873 > +(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
2875 sum=sum+sixth*hz2*hxi(v)*(
2877 > -(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))
2878 > +(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
2880 > -(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))
2881 > +(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
2887 if(ict(3).eq.1)
then
2902 cyd=3.0_r8*yp2-1.0_r8
2903 cydi=-3.0_r8*ypi2+1.0_r8
2912 czd=3.0_r8*zp2-1.0_r8
2913 czdi=-3.0_r8*zpi2+1.0_r8
2915 sum=hyi(v)*hzi(v)*hxi(v)*(
2917 > -(-fin(1,i,j,k) +fin(1,i,j+1,k))
2918 > +( -fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
2920 > -(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))
2921 > +( -fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
2923 sum=sum+sixth*hy(v)*hzi(v)*hxi(v)*(
2925 > -(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))
2926 > +(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
2928 > -(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))
2929 > +(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
2931 sum=sum+sixth*hz(v)*hyi(v)*hxi(v)*(
2933 > -(-fin(5,i,j,k) +fin(5,i,j+1,k))
2934 > +( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
2936 > -(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))
2937 > +( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
2939 sum=sum+z36th*hy(v)*hz(v)*hxi(v)*(
2941 > -(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))
2942 > +(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
2944 > -(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))
2945 > +(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
2951 if(ict(4).eq.1)
then
2967 cyi=ypi*(ypi2-1.0_r8)
2977 > -(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))
2978 > +(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
2980 > -(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))
2981 > +(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
2983 sum=sum+sixth*hy2*hxi(v)*(
2985 > -(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))
2986 > +(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
2988 > -(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))
2989 > +(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
2995 if(ict(5).eq.1)
then
3016 czi=zpi*(zpi2-1.0_r8)
3022 > xpi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
3023 > xp*( -fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
3025 > xpi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
3026 > xp*( -fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
3028 sum=sum+sixth*hz2*hyi(v)*(
3030 > xpi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
3031 > xp*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
3033 > xpi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
3034 > xp*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
3040 if(ict(6).eq.1)
then
3061 cyi=ypi*(ypi2-1.0_r8)
3066 > xpi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
3067 > xp*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
3069 > xpi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
3070 > xp*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
3072 sum=sum+sixth*hy2*hzi(v)*(
3074 > xpi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
3075 > xp*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
3077 > xpi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
3078 > xp*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
3084 if(ict(7).eq.1)
then
3099 cxd=3.0_r8*xp2-1.0_r8
3100 cxdi=-3.0_r8*xpi2+1.0_r8
3109 czd=3.0_r8*zp2-1.0_r8
3110 czdi=-3.0_r8*zpi2+1.0_r8
3112 sum=hxi(v)*hzi(v)*hyi(v)*(
3114 > ( fin(2,i,j,k) -fin(2,i,j+1,k))
3115 > +(-fin(2,i+1,j,k)+fin(2,i+1,j+1,k)))
3117 > ( fin(2,i,j,k+1) -fin(2,i,j+1,k+1))
3118 > +(-fin(2,i+1,j,k+1)+fin(2,i+1,j+1,k+1))))
3120 sum=sum+sixth*hx(v)*hzi(v)*hyi(v)*(
3122 > cxdi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
3123 > cxd*(-fin(4,i+1,j,k) +fin(4,i+1,j+1,k)))
3125 > cxdi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
3126 > cxd*(-fin(4,i+1,j,k+1) +fin(4,i+1,j+1,k+1))))
3128 sum=sum+sixth*hz(v)*hxi(v)*hyi(v)*(
3130 > ( fin(6,i,j,k) -fin(6,i,j+1,k))
3131 > +(-fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
3133 > ( fin(6,i,j,k+1) -fin(6,i,j+1,k+1))
3134 > +(-fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
3136 sum=sum+z36th*hx(v)*hz(v)*hyi(v)*(
3138 > cxdi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
3139 > cxd*(-fin(7,i+1,j,k) +fin(7,i+1,j+1,k)))
3141 > cxdi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
3142 > cxd*(-fin(7,i+1,j,k+1) +fin(7,i+1,j+1,k+1))))
3148 if(ict(8).eq.1)
then
3163 cxd=3.0_r8*xp2-1.0_r8
3164 cxdi=-3.0_r8*xpi2+1.0_r8
3173 cyd=3.0_r8*yp2-1.0_r8
3174 cydi=-3.0_r8*ypi2+1.0_r8
3176 sum=hyi(v)*hxi(v)*hzi(v)*(
3178 > ( +fin(3,i,j,k) -fin(3,i,j+1,k))
3179 > +( -fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
3181 > ( +fin(3,i,j,k+1) -fin(3,i,j+1,k+1))
3182 > +( -fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
3184 sum=sum+sixth*hx(v)*hyi(v)*hzi(v)*(
3186 > cxdi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
3187 > cxd*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
3189 > cxdi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
3190 > cxd*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
3192 sum=sum+sixth*hy(v)*hxi(v)*hzi(v)*(
3194 > -(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))
3195 > +(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
3197 > -(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))
3198 > +(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
3200 sum=sum+z36th*hx(v)*hy(v)*hzi(v)*(
3202 > cxdi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
3203 > cxd*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
3205 > cxdi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
3206 > cxd*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
3212 if(ict(9).eq.1)
then
3228 cxi=xpi*(xpi2-1.0_r8)
3238 > xpi*(-fin(6,i,j,k) +fin(6,i,j+1,k))+
3239 > xp*( -fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
3241 > xpi*(-fin(6,i,j,k+1) +fin(6,i,j+1,k+1))+
3242 > xp*( -fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
3244 sum=sum+sixth*hx2*hyi(v)*(
3246 > cxi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
3247 > cx*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
3249 > cxi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
3250 > cx*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
3256 if(ict(10).eq.1)
then
3272 cxi=xpi*(xpi2-1.0_r8)
3282 > xpi*(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))+
3283 > xp*( ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
3285 > xpi*(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))+
3286 > xp*(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
3288 sum=sum+sixth*hx2*hzi(v)*(
3290 > cxi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
3291 > cx*( ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
3293 > cxi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
3294 > cx*( ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
3303 else if(ict(1).eq.6)
then
3328 > xpi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
3329 > xp*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
3331 > xpi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
3332 > xp*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
3341 if(ict(1).eq.-6)
then
3342 if(ict(2).eq.1)
then
3358 czi=zpi*(zpi2-1.0_r8)
3363 > ( fin(4,i,j,k) -fin(4,i,j+1,k))
3364 > +(-fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
3366 > ( fin(4,i,j,k+1) -fin(4,i,j+1,k+1))
3367 > +(-fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
3369 sum=sum+sixth*hz2*hyi(v)*hxi(v)*(
3371 > ( fin(7,i,j,k) -fin(7,i,j+1,k))
3372 > +(-fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
3374 > ( fin(7,i,j,k+1) -fin(7,i,j+1,k+1))
3375 > +(-fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
3381 if(ict(3).eq.1)
then
3401 czd=3.0_r8*zp2-1.0_r8
3402 czdi=-3.0_r8*zpi2+1.0_r8
3406 > -(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))
3407 > +( ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
3409 > -(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))
3410 > +(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
3412 sum=sum+sixth*hz(v)*hxi(v)*(
3414 > -(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))
3415 > +(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
3417 > -(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))
3418 > +(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
3424 if(ict(4).eq.1)
then
3439 cyd=3.0_r8*yp2-1.0_r8
3440 cydi=-3.0_r8*ypi2+1.0_r8
3449 > ( fin(5,i,j,k) -fin(5,i,j+1,k))
3450 > +(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
3452 > ( fin(5,i,j,k+1) -fin(5,i,j+1,k+1))
3453 > +(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
3455 sum=sum+sixth*hy(v)*hxi(v)*(
3457 > -(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))
3458 > +(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
3460 > -(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))
3461 > +(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
3467 if(ict(5).eq.1)
then
3483 cyi=ypi*(ypi2-1.0_r8)
3488 > -(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))
3489 > +(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
3491 > -(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))
3492 > +(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
3494 sum=sum+sixth*hy2*hxi(v)*hzi(v)*(
3496 > -(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))
3497 > +(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
3499 > -(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))
3500 > +(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
3506 if(ict(6).eq.1)
then
3526 czd=3.0_r8*zp2-1.0_r8
3527 czdi=-3.0_r8*zpi2+1.0_r8
3531 > xpi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
3532 > xp*(-fin(4,i+1,j,k) +fin(4,i+1,j+1,k)))
3534 > xpi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
3535 > xp*(-fin(4,i+1,j,k+1) +fin(4,i+1,j+1,k+1))))
3537 sum=sum+sixth*hz(v)*hyi(v)*(
3539 > xpi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
3540 > xp*(-fin(7,i+1,j,k) +fin(7,i+1,j+1,k)))
3542 > xpi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
3543 > xp*(-fin(7,i+1,j,k+1) +fin(7,i+1,j+1,k+1))))
3549 if(ict(7).eq.1)
then
3569 cyd=3.0_r8*yp2-1.0_r8
3570 cydi=-3.0_r8*ypi2+1.0_r8
3574 > xpi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
3575 > xp*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
3577 > xpi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
3578 > xp*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
3580 sum=sum+sixth*hy(v)*hzi(v)*(
3582 > xpi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
3583 > xp*( cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
3585 > xpi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
3586 > xp*( cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
3592 if(ict(8).eq.1)
then
3607 cxd=3.0_r8*xp2-1.0_r8
3608 cxdi=-3.0_r8*xpi2+1.0_r8
3617 > ( fin(6,i,j,k) -fin(6,i,j+1,k))
3618 > +(-fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
3620 > ( fin(6,i,j,k+1) -fin(6,i,j+1,k+1))
3621 > +(-fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
3623 sum=sum+sixth*hx(v)*hyi(v)*(
3625 > cxdi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
3626 > cxd*(-fin(7,i+1,j,k) +fin(7,i+1,j+1,k)))
3628 > cxdi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
3629 > cxd*(-fin(7,i+1,j,k+1) +fin(7,i+1,j+1,k+1))))
3635 if(ict(9).eq.1)
then
3650 cxd=3.0_r8*xp2-1.0_r8
3651 cxdi=-3.0_r8*xpi2+1.0_r8
3660 > -(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))
3661 > +(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
3663 > -(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))
3664 > +(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
3666 sum=sum+sixth*hx(v)*hzi(v)*(
3668 > cxdi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
3669 > cxd*(ypi*fin(7,i+1,j,k) +yp*fin(7,i+1,j+1,k)))
3671 > cxdi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
3672 > cxd*( ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
3678 if(ict(10).eq.1)
then
3694 cxi=xpi*(xpi2-1.0_r8)
3699 > xpi*(-fin(6,i,j,k) +fin(6,i,j+1,k))+
3700 > xp*( -fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
3702 > xpi*(-fin(6,i,j,k+1) +fin(6,i,j+1,k+1))+
3703 > xp*( -fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
3705 sum=sum+sixth*hx2*hyi(v)*hzi(v)*(
3707 > cxi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
3708 > cx*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
3710 > cxi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
3711 > cx*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
3720 else if(abs(ict(1)).eq.7)
then
3721 if(ict(2).eq.1)
then
3736 czd=3.0_r8*zp2-1.0_r8
3737 czdi=-3.0_r8*zpi2+1.0_r8
3739 sum=hyi(v)*hxi(v)*hzi(v)*(
3741 > ( fin(4,i,j,k) -fin(4,i,j+1,k))
3742 > +(-fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
3744 > ( fin(4,i,j,k+1) -fin(4,i,j+1,k+1))
3745 > +(-fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
3747 sum=sum+sixth*hz(v)*hyi(v)*hxi(v)*(
3749 > ( fin(7,i,j,k) -fin(7,i,j+1,k))
3750 > +(-fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
3752 > ( fin(7,i,j,k+1) -fin(7,i,j+1,k+1))
3753 > +(-fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
3759 if(ict(3).eq.1)
then
3779 > -(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))
3780 > +(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
3782 > -(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))
3783 > +(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
3789 if(ict(4).eq.1)
then
3804 cyd=3.0_r8*yp2-1.0_r8
3805 cydi=-3.0_r8*ypi2+1.0_r8
3807 sum=hxi(v)*hyi(v)*hzi(v)*(
3809 > ( fin(5,i,j,k) -fin(5,i,j+1,k))
3810 > +(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
3812 > ( fin(5,i,j,k+1) -fin(5,i,j+1,k+1))
3813 > +(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
3815 sum=sum+sixth*hy(v)*hxi(v)*hzi(v)*(
3817 > -(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))
3818 > +(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
3820 > -(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))
3821 > +(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
3827 if(ict(5).eq.1)
then
3847 > xpi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
3848 > xp*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
3850 > xpi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
3851 > xp*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
3857 if(ict(6).eq.1)
then
3877 > xpi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
3878 > xp*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
3880 > xpi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
3881 > xp*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
3887 if(ict(7).eq.1)
then
3902 cxd=3.0_r8*xp2-1.0_r8
3903 cxdi=-3.0_r8*xpi2+1.0_r8
3905 sum=hxi(v)*hyi(v)*hzi(v)*(
3907 > ( fin(6,i,j,k) -fin(6,i,j+1,k))
3908 > +(-fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
3910 > ( fin(6,i,j,k+1) -fin(6,i,j+1,k+1))
3911 > +(-fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
3913 sum=sum+sixth*hx(v)*hyi(v)*hzi(v)*(
3915 > cxdi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
3916 > cxd*(-fin(7,i+1,j,k) +fin(7,i+1,j+1,k)))
3918 > cxdi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
3919 > cxd*(-fin(7,i+1,j,k+1) +fin(7,i+1,j+1,k+1))))
3928 else if(abs(ict(1)).eq.8)
then
3929 if(ict(2).eq.1)
then
3944 > ( fin(7,i,j,k) -fin(7,i,j+1,k))
3945 > +(-fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
3947 > ( fin(7,i,j,k+1) -fin(7,i,j+1,k+1))
3948 > +(-fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
3954 if(ict(3).eq.1)
then
3969 > -(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))
3970 > +(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
3972 > -(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))
3973 > +(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
3979 if(ict(4).eq.1)
then
3994 > xpi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
3995 > xp*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
3997 > xpi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
3998 > xp*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
4007 else if(abs(ict(1)).eq.9)
then
4015 sum=hyi(v)*hxi(v)*hzi(v)*(
4017 > ( fin(7,i,j,k) -fin(7,i,j+1,k))
4018 > +(-fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
4020 > ( fin(7,i,j,k+1) -fin(7,i,j+1,k+1))
4021 > +(-fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
4034 subroutine r8fvtricubx(ict,ivec,ivecd,
4035 > fval,ii,jj,kk,xparam,yparam,zparam,
4036 > hx,hxi,hy,hyi,hz,hzi,
4044 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
4045 INTEGER inf3,nz,inf2,iadr,j,k,i
4048 REAL*8 z36th,z216th,yp,ypi,yp2,ypi2,cy,cyi,hy2,zp,zpi,zp2
4049 REAL*8 zpi2,cz,czi,hz2,xp,xpi,xp2,xpi2,cx,cxi,hx2,cxd,cxdi
4050 REAL*8 cyd,cydi,czd,czdi
4056 integer ii(ivec),jj,kk
4057 real*8 xparam(ivec),yparam,zparam
4060 real*8 hx(ivec),hy,hz
4061 real*8 hxi(ivec),hyi,hzi
4064 real*8 fin(0:7,inf2,inf3,nz)
4066 real*8 fval(ivecd,*)
4084 REAL*8,
parameter :: sixth = 0.166666666666666667_r8
4089 z216th=sixth*sixth*sixth
4092 if(abs(ict(1)).le.2)
then
4098 if(ict(1).eq.1)
then
4113 cyi=ypi*(ypi2-1.0_r8)
4124 czi=zpi*(zpi2-1.0_r8)
4139 cxi=xpi*(xpi2-1.0_r8)
4144 > xpi*(ypi*fin(0,i,j,k) +yp*fin(0,i,j+1,k))+
4145 > xp*(ypi*fin(0,i+1,j,k)+yp*fin(0,i+1,j+1,k)))
4147 > xpi*(ypi*fin(0,i,j,k+1) +yp*fin(0,i,j+1,k+1))+
4148 > xp*(ypi*fin(0,i+1,j,k+1)+yp*fin(0,i+1,j+1,k+1))))
4152 > cxi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
4153 > cx*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
4155 > cxi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
4156 > cx*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
4160 > xpi*(cyi*fin(2,i,j,k) +cy*fin(2,i,j+1,k))+
4161 > xp*(cyi*fin(2,i+1,j,k)+cy*fin(2,i+1,j+1,k)))
4163 > xpi*(cyi*fin(2,i,j,k+1) +cy*fin(2,i,j+1,k+1))+
4164 > xp*(cyi*fin(2,i+1,j,k+1)+cy*fin(2,i+1,j+1,k+1))))
4168 > xpi*(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))+
4169 > xp*(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
4171 > xpi*(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))+
4172 > xp*(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
4174 sum=sum+z36th*hx2*hy2*(
4176 > cxi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
4177 > cx*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
4179 > cxi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
4180 > cx*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
4182 sum=sum+z36th*hx2*hz2*(
4184 > cxi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
4185 > cx*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
4187 > cxi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
4188 > cx*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
4190 sum=sum+z36th*hy2*hz2*(
4192 > xpi*(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))+
4193 > xp*(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
4195 > xpi*(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))+
4196 > xp*(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
4198 sum=sum+z216th*hx2*hy2*hz2*(
4200 > cxi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
4201 > cx*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
4203 > cxi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
4204 > cx*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
4210 if(ict(2).eq.1)
then
4225 cyi=ypi*(ypi2-1.0_r8)
4236 czi=zpi*(zpi2-1.0_r8)
4250 cxd=3.0_r8*xp2-1.0_r8
4251 cxdi=-3.0_r8*xpi2+1.0_r8
4255 > -(ypi*fin(0,i,j,k) +yp*fin(0,i,j+1,k))
4256 > +(ypi*fin(0,i+1,j,k)+yp*fin(0,i+1,j+1,k)))
4258 > -(ypi*fin(0,i,j,k+1) +yp*fin(0,i,j+1,k+1))
4259 > +(ypi*fin(0,i+1,j,k+1)+yp*fin(0,i+1,j+1,k+1))))
4261 sum=sum+sixth*hx(v)*(
4263 > cxdi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
4264 > cxd*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
4266 > cxdi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
4267 > cxd*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
4269 sum=sum+sixth*hxi(v)*hy2*(
4271 > -(cyi*fin(2,i,j,k) +cy*fin(2,i,j+1,k))
4272 > +(cyi*fin(2,i+1,j,k)+cy*fin(2,i+1,j+1,k)))
4274 > -(cyi*fin(2,i,j,k+1) +cy*fin(2,i,j+1,k+1))
4275 > +(cyi*fin(2,i+1,j,k+1)+cy*fin(2,i+1,j+1,k+1))))
4277 sum=sum+sixth*hxi(v)*hz2*(
4279 > -(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))
4280 > +(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
4282 > -(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))
4283 > +(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
4285 sum=sum+z36th*hx(v)*hy2*(
4287 > cxdi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
4288 > cxd*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
4290 > cxdi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
4291 > cxd*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
4293 sum=sum+z36th*hx(v)*hz2*(
4295 > cxdi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
4296 > cxd*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
4298 > cxdi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
4299 > cxd*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
4301 sum=sum+z36th*hxi(v)*hy2*hz2*(
4303 > -(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))
4304 > +(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
4306 > -(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))
4307 > +(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
4309 sum=sum+z216th*hx(v)*hy2*hz2*(
4311 > cxdi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
4312 > cxd*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
4314 > cxdi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
4315 > cxd*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
4321 if(ict(3).eq.1)
then
4335 cyd=3.0_r8*yp2-1.0_r8
4336 cydi=-3.0_r8*ypi2+1.0_r8
4346 czi=zpi*(zpi2-1.0_r8)
4361 cxi=xpi*(xpi2-1.0_r8)
4366 > xpi*(-fin(0,i,j,k) +fin(0,i,j+1,k))+
4367 > xp*(-fin(0,i+1,j,k)+fin(0,i+1,j+1,k)))
4369 > xpi*(-fin(0,i,j,k+1) +fin(0,i,j+1,k+1))+
4370 > xp*(-fin(0,i+1,j,k+1)+fin(0,i+1,j+1,k+1))))
4372 sum=sum+sixth*hyi*hx2*(
4374 > cxi*(-fin(1,i,j,k) +fin(1,i,j+1,k))+
4375 > cx*(-fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
4377 > cxi*(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))+
4378 > cx*(-fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
4382 > xpi*(cydi*fin(2,i,j,k) +cyd*fin(2,i,j+1,k))+
4383 > xp*(cydi*fin(2,i+1,j,k)+cyd*fin(2,i+1,j+1,k)))
4385 > xpi*(cydi*fin(2,i,j,k+1) +cyd*fin(2,i,j+1,k+1))+
4386 > xp*(cydi*fin(2,i+1,j,k+1)+cyd*fin(2,i+1,j+1,k+1))))
4388 sum=sum+sixth*hyi*hz2*(
4390 > xpi*(-fin(3,i,j,k) +fin(3,i,j+1,k))+
4391 > xp*(-fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
4393 > xpi*(-fin(3,i,j,k+1) +fin(3,i,j+1,k+1))+
4394 > xp*(-fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
4396 sum=sum+z36th*hx2*hy*(
4398 > cxi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
4399 > cx*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
4401 > cxi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
4402 > cx*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
4404 sum=sum+z36th*hyi*hx2*hz2*(
4406 > cxi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
4407 > cx*(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
4409 > cxi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
4410 > cx*(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
4412 sum=sum+z36th*hy*hz2*(
4414 > xpi*(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))+
4415 > xp*(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
4417 > xpi*(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))+
4418 > xp*(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
4420 sum=sum+z216th*hx2*hy*hz2*(
4422 > cxi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
4423 > cx*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
4425 > cxi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
4426 > cx*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
4432 if(ict(4).eq.1)
then
4447 cyi=ypi*(ypi2-1.0_r8)
4457 czd=3.0_r8*zp2-1.0_r8
4458 czdi=-3.0_r8*zpi2+1.0_r8
4472 cxi=xpi*(xpi2-1.0_r8)
4477 > xpi*(ypi*fin(0,i,j,k) +yp*fin(0,i,j+1,k))+
4478 > xp*(ypi*fin(0,i+1,j,k)+yp*fin(0,i+1,j+1,k)))
4480 > xpi*(ypi*fin(0,i,j,k+1) +yp*fin(0,i,j+1,k+1))+
4481 > xp*(ypi*fin(0,i+1,j,k+1)+yp*fin(0,i+1,j+1,k+1))))
4483 sum=sum+sixth*hx2*hzi*(
4485 > cxi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
4486 > cx*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
4488 > cxi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
4489 > cx*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
4491 sum=sum+sixth*hy2*hzi*(
4493 > xpi*(cyi*fin(2,i,j,k) +cy*fin(2,i,j+1,k))+
4494 > xp*(cyi*fin(2,i+1,j,k)+cy*fin(2,i+1,j+1,k)))
4496 > xpi*(cyi*fin(2,i,j,k+1) +cy*fin(2,i,j+1,k+1))+
4497 > xp*(cyi*fin(2,i+1,j,k+1)+cy*fin(2,i+1,j+1,k+1))))
4501 > xpi*(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))+
4502 > xp*(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
4504 > xpi*(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))+
4505 > xp*(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
4507 sum=sum+z36th*hx2*hy2*hzi*(
4509 > cxi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
4510 > cx*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
4512 > cxi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
4513 > cx*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
4515 sum=sum+z36th*hx2*hz*(
4517 > cxi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
4518 > cx*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
4520 > cxi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
4521 > cx*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
4523 sum=sum+z36th*hy2*hz*(
4525 > xpi*(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))+
4526 > xp*(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
4528 > xpi*(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))+
4529 > xp*(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
4531 sum=sum+z216th*hx2*hy2*hz*(
4533 > cxi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
4534 > cx*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
4536 > cxi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
4537 > cx*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
4543 if(ict(5).eq.1)
then
4558 cyi=ypi*(ypi2-1.0_r8)
4569 czi=zpi*(zpi2-1.0_r8)
4583 > xpi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
4584 > xp*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
4586 > xpi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
4587 > xp*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
4591 > xpi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
4592 > xp*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
4594 > xpi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
4595 > xp*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
4599 > xpi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
4600 > xp*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
4602 > xpi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
4603 > xp*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
4605 sum=sum+z36th*hy2*hz2*(
4607 > xpi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
4608 > xp*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
4610 > xpi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
4611 > xp*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
4617 if(ict(6).eq.1)
then
4637 czi=zpi*(zpi2-1.0_r8)
4652 cxi=xpi*(xpi2-1.0_r8)
4657 > xpi*(ypi*fin(2,i,j,k) +yp*fin(2,i,j+1,k))+
4658 > xp*(ypi*fin(2,i+1,j,k)+yp*fin(2,i+1,j+1,k)))
4660 > xpi*(ypi*fin(2,i,j,k+1) +yp*fin(2,i,j+1,k+1))+
4661 > xp*(ypi*fin(2,i+1,j,k+1)+yp*fin(2,i+1,j+1,k+1))))
4665 > cxi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
4666 > cx*(ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
4668 > cxi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
4669 > cx*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
4673 > xpi*(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))+
4674 > xp*(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
4676 > xpi*(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))+
4677 > xp*(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
4679 sum=sum+z36th*hx2*hz2*(
4681 > cxi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
4682 > cx*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
4684 > cxi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
4685 > cx*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
4691 if(ict(7).eq.1)
then
4706 cyi=ypi*(ypi2-1.0_r8)
4726 cxi=xpi*(xpi2-1.0_r8)
4731 > xpi*(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))+
4732 > xp*(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
4734 > xpi*(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))+
4735 > xp*(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
4739 > cxi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
4740 > cx*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
4742 > cxi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
4743 > cx*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
4747 > xpi*(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))+
4748 > xp*(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
4750 > xpi*(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))+
4751 > xp*(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
4753 sum=sum+z36th*hx2*hy2*(
4755 > cxi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
4756 > cx*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
4758 > cxi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
4759 > cx*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
4765 if(ict(8).eq.1)
then
4779 cyd=3.0_r8*yp2-1.0_r8
4780 cydi=-3.0_r8*ypi2+1.0_r8
4790 czi=zpi*(zpi2-1.0_r8)
4804 cxd=3.0_r8*xp2-1.0_r8
4805 cxdi=-3.0_r8*xpi2+1.0_r8
4809 > (fin(0,i,j,k) -fin(0,i,j+1,k))-
4810 > (fin(0,i+1,j,k)-fin(0,i+1,j+1,k)))
4812 > (fin(0,i,j,k+1) -fin(0,i,j+1,k+1))-
4813 > (fin(0,i+1,j,k+1)-fin(0,i+1,j+1,k+1))))
4815 sum=sum+sixth*hyi*hx(v)*(
4817 > cxdi*(-fin(1,i,j,k) +fin(1,i,j+1,k))+
4818 > cxd*(-fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
4820 > cxdi*(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))+
4821 > cxd*(-fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
4823 sum=sum+sixth*hxi(v)*hy*(
4825 > -(cydi*fin(2,i,j,k) +cyd*fin(2,i,j+1,k))
4826 > +(cydi*fin(2,i+1,j,k)+cyd*fin(2,i+1,j+1,k)))
4828 > -(cydi*fin(2,i,j,k+1) +cyd*fin(2,i,j+1,k+1))
4829 > +(cydi*fin(2,i+1,j,k+1)+cyd*fin(2,i+1,j+1,k+1))))
4831 sum=sum+sixth*hxi(v)*hyi*hz2*(
4833 > (fin(3,i,j,k) -fin(3,i,j+1,k))-
4834 > (fin(3,i+1,j,k)-fin(3,i+1,j+1,k)))
4836 > (fin(3,i,j,k+1) -fin(3,i,j+1,k+1))-
4837 > (fin(3,i+1,j,k+1)-fin(3,i+1,j+1,k+1))))
4839 sum=sum+z36th*hx(v)*hy*(
4841 > cxdi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
4842 > cxd*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
4844 > cxdi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
4845 > cxd*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
4847 sum=sum+z36th*hyi*hx(v)*hz2*(
4849 > cxdi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
4850 > cxd*(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
4852 > cxdi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
4853 > cxd*(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
4855 sum=sum+z36th*hxi(v)*hy*hz2*(
4857 > -(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))
4858 > +(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
4860 > -(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))
4861 > +(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
4863 sum=sum+z216th*hx(v)*hy*hz2*(
4865 > cxdi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
4866 > cxd*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
4868 > cxdi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
4869 > cxd*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
4875 if(ict(9).eq.1)
then
4890 cyi=ypi*(ypi2-1.0_r8)
4900 czd=3.0_r8*zp2-1.0_r8
4901 czdi=-3.0_r8*zpi2+1.0_r8
4914 cxd=3.0_r8*xp2-1.0_r8
4915 cxdi=-3.0_r8*xpi2+1.0_r8
4919 > (ypi*fin(0,i,j,k) +yp*fin(0,i,j+1,k)) -
4920 > (ypi*fin(0,i+1,j,k)+yp*fin(0,i+1,j+1,k)))
4922 > (ypi*fin(0,i,j,k+1) +yp*fin(0,i,j+1,k+1)) -
4923 > (ypi*fin(0,i+1,j,k+1)+yp*fin(0,i+1,j+1,k+1))))
4925 sum=sum+sixth*hx(v)*hzi*(
4927 > cxdi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
4928 > cxd*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
4930 > cxdi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
4931 > cxd*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
4933 sum=sum+sixth*hxi(v)*hy2*hzi*(
4935 > (cyi*fin(2,i,j,k) +cy*fin(2,i,j+1,k)) -
4936 > (cyi*fin(2,i+1,j,k)+cy*fin(2,i+1,j+1,k)))
4938 > (cyi*fin(2,i,j,k+1) +cy*fin(2,i,j+1,k+1)) -
4939 > (cyi*fin(2,i+1,j,k+1)+cy*fin(2,i+1,j+1,k+1))))
4941 sum=sum+sixth*hxi(v)*hz*(
4943 > -(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))
4944 > +(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
4946 > -(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))
4947 > +(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
4949 sum=sum+z36th*hx(v)*hy2*hzi*(
4951 > cxdi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
4952 > cxd*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
4954 > cxdi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
4955 > cxd*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
4957 sum=sum+z36th*hx(v)*hz*(
4959 > cxdi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
4960 > cxd*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
4962 > cxdi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
4963 > cxd*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
4965 sum=sum+z36th*hxi(v)*hy2*hz*(
4967 > -(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))
4968 > +(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
4970 > -(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))
4971 > +(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
4973 sum=sum+z216th*hx(v)*hy2*hz*(
4975 > cxdi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
4976 > cxd*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
4978 > cxdi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
4979 > cxd*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
4985 if(ict(10).eq.1)
then
4999 cyd=3.0_r8*yp2-1.0_r8
5000 cydi=-3.0_r8*ypi2+1.0_r8
5009 czd=3.0_r8*zp2-1.0_r8
5010 czdi=-3.0_r8*zpi2+1.0_r8
5024 cxi=xpi*(xpi2-1.0_r8)
5029 > xpi*(fin(0,i,j,k) -fin(0,i,j+1,k))+
5030 > xp*(fin(0,i+1,j,k)-fin(0,i+1,j+1,k)))
5032 > xpi*(fin(0,i,j,k+1) -fin(0,i,j+1,k+1))+
5033 > xp*(fin(0,i+1,j,k+1)-fin(0,i+1,j+1,k+1))))
5035 sum=sum+sixth*hyi*hx2*hzi*(
5037 > cxi*(fin(1,i,j,k) -fin(1,i,j+1,k))+
5038 > cx*(fin(1,i+1,j,k)-fin(1,i+1,j+1,k)))
5040 > cxi*(fin(1,i,j,k+1) -fin(1,i,j+1,k+1))+
5041 > cx*(fin(1,i+1,j,k+1)-fin(1,i+1,j+1,k+1))))
5043 sum=sum+sixth*hy*hzi*(
5045 > xpi*(cydi*fin(2,i,j,k) +cyd*fin(2,i,j+1,k))+
5046 > xp*(cydi*fin(2,i+1,j,k)+cyd*fin(2,i+1,j+1,k)))
5048 > xpi*(cydi*fin(2,i,j,k+1) +cyd*fin(2,i,j+1,k+1))+
5049 > xp*(cydi*fin(2,i+1,j,k+1)+cyd*fin(2,i+1,j+1,k+1))))
5051 sum=sum+sixth*hyi*hz*(
5053 > xpi*(-fin(3,i,j,k) +fin(3,i,j+1,k))+
5054 > xp*(-fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
5056 > xpi*(-fin(3,i,j,k+1) +fin(3,i,j+1,k+1))+
5057 > xp*(-fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
5059 sum=sum+z36th*hx2*hy*hzi*(
5061 > cxi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
5062 > cx*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
5064 > cxi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
5065 > cx*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
5067 sum=sum+z36th*hyi*hx2*hz*(
5069 > cxi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
5070 > cx*(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
5072 > cxi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
5073 > cx*(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
5075 sum=sum+z36th*hy*hz*(
5077 > xpi*(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))+
5078 > xp*(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
5080 > xpi*(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))+
5081 > xp*(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
5083 sum=sum+z216th*hx2*hy*hz*(
5085 > cxi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
5086 > cx*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
5088 > cxi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
5089 > cx*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
5098 else if(ict(1).eq.3)
then
5099 if(ict(2).eq.1)
then
5111 cyd=3.0_r8*yp2-1.0_r8
5112 cydi=-3.0_r8*ypi2+1.0_r8
5122 czi=zpi*(zpi2-1.0_r8)
5136 > xpi*(-fin(1,i,j,k) +fin(1,i,j+1,k))+
5137 > xp*( -fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
5139 > xpi*(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))+
5140 > xp*( -fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
5144 > xpi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
5145 > xp*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
5147 > xpi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
5148 > xp*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
5150 sum=sum+sixth*hz2*hyi*(
5152 > xpi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
5153 > xp*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
5155 > xpi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
5156 > xp*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
5158 sum=sum+z36th*hy*hz2*(
5160 > xpi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
5161 > xp*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
5163 > xpi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
5164 > xp*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
5170 if(ict(3).eq.1)
then
5183 cyi=ypi*(ypi2-1.0_r8)
5193 czd=3.0_r8*zp2-1.0_r8
5194 czdi=-3.0_r8*zpi2+1.0_r8
5207 > xpi*(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))+
5208 > xp*(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
5210 > xpi*(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))+
5211 > xp*(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
5213 sum=sum+sixth*hy2*hzi*(
5215 > xpi*(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))+
5216 > xp*(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
5218 > xpi*(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))+
5219 > xp*(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
5223 > xpi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
5224 > xp*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
5226 > xpi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
5227 > xp*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
5229 sum=sum+z36th*hy2*hz*(
5231 > xpi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
5232 > xp*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
5234 > xpi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
5235 > xp*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
5241 if(ict(4).eq.1)
then
5259 czi=zpi*(zpi2-1.0_r8)
5273 cxd=3.0_r8*xp2-1.0_r8
5274 cxdi=-3.0_r8*xpi2+1.0_r8
5278 > -(ypi*fin(2,i,j,k) +yp*fin(2,i,j+1,k))
5279 > +(ypi*fin(2,i+1,j,k)+yp*fin(2,i+1,j+1,k)))
5281 > -(ypi*fin(2,i,j,k+1) +yp*fin(2,i,j+1,k+1))
5282 > +(ypi*fin(2,i+1,j,k+1)+yp*fin(2,i+1,j+1,k+1))))
5284 sum=sum+sixth*hx(v)*(
5286 > cxdi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
5287 > cxd*(ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
5289 > cxdi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
5290 > cxd*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
5292 sum=sum+sixth*hz2*hxi(v)*(
5294 > -(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))
5295 > +(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
5297 > -(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))
5298 > +(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
5300 sum=sum+z36th*hx(v)*hz2*(
5302 > cxdi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
5303 > cxd*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
5305 > cxdi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
5306 > cxd*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
5312 if(ict(5).eq.1)
then
5324 cyd=3.0_r8*yp2-1.0_r8
5325 cydi=-3.0_r8*ypi2+1.0_r8
5334 czd=3.0_r8*zp2-1.0_r8
5335 czdi=-3.0_r8*zpi2+1.0_r8
5348 cxd=3.0_r8*xp2-1.0_r8
5349 cxdi=-3.0_r8*xpi2+1.0_r8
5351 sum=hxi(v)*hyi*hzi*(
5353 > (fin(0,i,j,k) -fin(0,i,j+1,k))-
5354 > (fin(0,i+1,j,k)-fin(0,i+1,j+1,k)))
5356 > (fin(0,i,j,k+1) -fin(0,i,j+1,k+1))-
5357 > (fin(0,i+1,j,k+1)-fin(0,i+1,j+1,k+1))))
5359 sum=sum+sixth*hyi*hx(v)*hzi*(
5361 > cxdi*(-fin(1,i,j,k) +fin(1,i,j+1,k))+
5362 > cxd*(-fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
5364 > cxdi*(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))+
5365 > cxd*(-fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
5367 sum=sum+sixth*hxi(v)*hy*hzi*(
5369 > -(cydi*fin(2,i,j,k) +cyd*fin(2,i,j+1,k))
5370 > +(cydi*fin(2,i+1,j,k)+cyd*fin(2,i+1,j+1,k)))
5372 > -(cydi*fin(2,i,j,k+1) +cyd*fin(2,i,j+1,k+1))
5373 > +(cydi*fin(2,i+1,j,k+1)+cyd*fin(2,i+1,j+1,k+1))))
5375 sum=sum+sixth*hxi(v)*hyi*hz*(
5377 > (fin(3,i,j,k) -fin(3,i,j+1,k))-
5378 > (fin(3,i+1,j,k)-fin(3,i+1,j+1,k)))
5380 > (fin(3,i,j,k+1) -fin(3,i,j+1,k+1))-
5381 > (fin(3,i+1,j,k+1)-fin(3,i+1,j+1,k+1))))
5383 sum=sum+z36th*hx(v)*hy*hzi*(
5385 > cxdi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
5386 > cxd*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
5388 > cxdi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
5389 > cxd*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
5391 sum=sum+z36th*hyi*hx(v)*hz*(
5393 > cxdi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
5394 > cxd*(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
5396 > cxdi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
5397 > cxd*(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
5399 sum=sum+z36th*hxi(v)*hy*hz*(
5401 > -(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))
5402 > +(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
5404 > -(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))
5405 > +(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
5407 sum=sum+z216th*hx(v)*hy*hz*(
5409 > cxdi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
5410 > cxd*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
5412 > cxdi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
5413 > cxd*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
5419 if(ict(6).eq.1)
then
5432 cyi=ypi*(ypi2-1.0_r8)
5452 cxd=3.0_r8*xp2-1.0_r8
5453 cxdi=-3.0_r8*xpi2+1.0_r8
5457 > -(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))
5458 > +(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
5460 > -(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))
5461 > +(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
5463 sum=sum+sixth*hx(v)*(
5465 > cxdi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
5466 > cxd*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
5468 > cxdi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
5469 > cxd*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
5471 sum=sum+sixth*hy2*hxi(v)*(
5473 > -(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))
5474 > +(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
5476 > -(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))
5477 > +(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
5479 sum=sum+z36th*hx(v)*hy2*(
5481 > cxdi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
5482 > cxd*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
5484 > cxdi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
5485 > cxd*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
5491 if(ict(7).eq.1)
then
5509 czd=3.0_r8*zp2-1.0_r8
5510 czdi=-3.0_r8*zpi2+1.0_r8
5524 cxi=xpi*(xpi2-1.0_r8)
5529 > xpi*(ypi*fin(2,i,j,k) +yp*fin(2,i,j+1,k))+
5530 > xp*(ypi*fin(2,i+1,j,k)+yp*fin(2,i+1,j+1,k)))
5532 > xpi*(ypi*fin(2,i,j,k+1) +yp*fin(2,i,j+1,k+1))+
5533 > xp*(ypi*fin(2,i+1,j,k+1)+yp*fin(2,i+1,j+1,k+1))))
5535 sum=sum+sixth*hx2*hzi*(
5537 > cxi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
5538 > cx*(ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
5540 > cxi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
5541 > cx*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
5545 > xpi*(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))+
5546 > xp*(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
5548 > xpi*(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))+
5549 > xp*(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
5551 sum=sum+z36th*hx2*hz*(
5553 > cxi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
5554 > cx*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
5556 > cxi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
5557 > cx*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
5563 if(ict(8).eq.1)
then
5575 cyd=3.0_r8*yp2-1.0_r8
5576 cydi=-3.0_r8*ypi2+1.0_r8
5595 cxi=xpi*(xpi2-1.0_r8)
5600 > xpi*(-fin(3,i,j,k) +fin(3,i,j+1,k))+
5601 > xp*( -fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
5603 > xpi*(-fin(3,i,j,k+1) +fin(3,i,j+1,k+1))+
5604 > xp*( -fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
5606 sum=sum+sixth*hx2*hyi*(
5608 > cxi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
5609 > cx*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
5611 > cxi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
5612 > cx*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
5616 > xpi*(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))+
5617 > xp*(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
5619 > xpi*(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))+
5620 > xp*(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
5622 sum=sum+z36th*hx2*hy*(
5624 > cxi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
5625 > cx*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
5627 > cxi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
5628 > cx*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
5636 else if(ict(1).eq.-3)
then
5637 if(ict(2).eq.1)
then
5653 cyi=ypi*(ypi2-1.0_r8)
5664 czi=zpi*(zpi2-1.0_r8)
5669 > -(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))
5670 > +(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
5672 > -(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))
5673 > +(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
5675 sum=sum+sixth*hy2*hxi(v)*(
5677 > -(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))
5678 > +(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
5680 > -(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))
5681 > +(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
5683 sum=sum+sixth*hz2*hxi(v)*(
5685 > -(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))
5686 > +(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
5688 > -(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))
5689 > +(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
5691 sum=sum+z36th*hy2*hz2*hxi(v)*(
5693 > -(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))
5694 > +(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
5696 > -(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))
5697 > +(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
5703 if(ict(3).eq.1)
then
5716 czi=zpi*(zpi2-1.0_r8)
5731 cxi=xpi*(xpi2-1.0_r8)
5736 > xpi*(-fin(2,i,j,k) +fin(2,i,j+1,k))+
5737 > xp*( -fin(2,i+1,j,k)+fin(2,i+1,j+1,k)))
5739 > xpi*(-fin(2,i,j,k+1) +fin(2,i,j+1,k+1))+
5740 > xp*( -fin(2,i+1,j,k+1)+fin(2,i+1,j+1,k+1))))
5742 sum=sum+sixth*hx2*hyi*(
5744 > cxi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
5745 > cx*( -fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
5747 > cxi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
5748 > cx*( -fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
5750 sum=sum+sixth*hz2*hyi*(
5752 > xpi*(-fin(6,i,j,k) +fin(6,i,j+1,k))+
5753 > xp*( -fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
5755 > xpi*(-fin(6,i,j,k+1) +fin(6,i,j+1,k+1))+
5756 > xp*( -fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
5758 sum=sum+z36th*hx2*hz2*hyi*(
5760 > cxi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
5761 > cx*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
5763 > cxi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
5764 > cx*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
5770 if(ict(4).eq.1)
then
5783 cyi=ypi*(ypi2-1.0_r8)
5798 cxi=xpi*(xpi2-1.0_r8)
5803 > xpi*(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))+
5804 > xp*(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
5806 > xpi*(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))+
5807 > xp*(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
5809 sum=sum+sixth*hx2*hzi*(
5811 > cxi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
5812 > cx*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
5814 > cxi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
5815 > cx*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
5817 sum=sum+sixth*hy2*hzi*(
5819 > xpi*(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))+
5820 > xp*(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
5822 > xpi*(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))+
5823 > xp*(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
5825 sum=sum+z36th*hx2*hy2*hzi*(
5827 > cxi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
5828 > cx*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
5830 > cxi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
5831 > cx*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
5840 else if(ict(1).eq.4)
then
5841 if(ict(2).eq.1)
then
5859 czi=zpi*(zpi2-1.0_r8)
5873 > xpi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
5874 > xp*( ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
5876 > xpi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
5877 > xp*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
5881 > xpi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
5882 > xp*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
5884 > xpi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
5885 > xp*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
5891 if(ict(3).eq.1)
then
5903 cyd=3.0_r8*yp2-1.0_r8
5904 cydi=-3.0_r8*ypi2+1.0_r8
5913 czd=3.0_r8*zp2-1.0_r8
5914 czdi=-3.0_r8*zpi2+1.0_r8
5927 > xpi*(-fin(1,i,j,k) +fin(1,i,j+1,k))+
5928 > xp*( -fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
5930 > xpi*(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))+
5931 > xp*( -fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
5933 sum=sum+sixth*hy*hzi*(
5935 > xpi*(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
5936 > xp*(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
5938 > xpi*(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
5939 > xp*(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
5941 sum=sum+sixth*hz*hyi*(
5943 > xpi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
5944 > xp*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
5946 > xpi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
5947 > xp*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
5949 sum=sum+z36th*hy*hz*(
5951 > xpi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
5952 > xp*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
5954 > xpi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
5955 > xp*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
5961 if(ict(4).eq.1)
then
5974 cyi=ypi*(ypi2-1.0_r8)
5993 > xpi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
5994 > xp*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
5996 > xpi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
5997 > xp*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
6001 > xpi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
6002 > xp*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
6004 > xpi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
6005 > xp*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
6011 if(ict(5).eq.1)
then
6028 czd=3.0_r8*zp2-1.0_r8
6029 czdi=-3.0_r8*zpi2+1.0_r8
6042 cxd=3.0_r8*xp2-1.0_r8
6043 cxdi=-3.0_r8*xpi2+1.0_r8
6047 > -(ypi*fin(2,i,j,k) +yp*fin(2,i,j+1,k))
6048 > +(ypi*fin(2,i+1,j,k)+yp*fin(2,i+1,j+1,k)))
6050 > -(ypi*fin(2,i,j,k+1) +yp*fin(2,i,j+1,k+1))
6051 > +(ypi*fin(2,i+1,j,k+1)+yp*fin(2,i+1,j+1,k+1))))
6053 sum=sum+sixth*hx(v)*hzi*(
6055 > cxdi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
6056 > cxd*(ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
6058 > cxdi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
6059 > cxd*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
6061 sum=sum+sixth*hz*hxi(v)*(
6063 > -(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))
6064 > +(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
6066 > -(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))
6067 > +(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
6069 sum=sum+z36th*hx(v)*hz*(
6071 > cxdi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
6072 > cxd*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
6074 > cxdi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
6075 > cxd*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
6081 if(ict(6).eq.1)
then
6093 cyd=3.0_r8*yp2-1.0_r8
6094 cydi=-3.0_r8*ypi2+1.0_r8
6112 cxd=3.0_r8*xp2-1.0_r8
6113 cxdi=-3.0_r8*xpi2+1.0_r8
6117 > ( +fin(3,i,j,k) -fin(3,i,j+1,k))
6118 > +( -fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
6120 > ( +fin(3,i,j,k+1) -fin(3,i,j+1,k+1))
6121 > +( -fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
6123 sum=sum+sixth*hx(v)*hyi*(
6125 > cxdi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
6126 > cxd*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
6128 > cxdi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
6129 > cxd*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
6131 sum=sum+sixth*hy*hxi(v)*(
6133 > -(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))
6134 > +(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
6136 > -(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))
6137 > +(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
6139 sum=sum+z36th*hx(v)*hy*(
6141 > cxdi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
6142 > cxd*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
6144 > cxdi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
6145 > cxd*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
6151 if(ict(7).eq.1)
then
6178 cxi=xpi*(xpi2-1.0_r8)
6183 > xpi*(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))+
6184 > xp*( ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
6186 > xpi*(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))+
6187 > xp*(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
6191 > cxi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
6192 > cx*( ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
6194 > cxi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
6195 > cx*( ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
6204 else if(ict(1).eq.-4)
then
6205 if(ict(2).eq.1)
then
6217 cyd=3.0_r8*yp2-1.0_r8
6218 cydi=-3.0_r8*ypi2+1.0_r8
6228 czi=zpi*(zpi2-1.0_r8)
6237 > ( fin(1,i,j,k) -fin(1,i,j+1,k))+
6238 > ( -fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
6240 > ( fin(1,i,j,k+1) -fin(1,i,j+1,k+1))+
6241 > ( -fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
6243 sum=sum+sixth*hy*hxi(v)*(
6245 > -(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))+
6246 > (cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
6248 > -(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))+
6249 > (cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
6251 sum=sum+sixth*hz2*hyi*hxi(v)*(
6253 > ( fin(5,i,j,k) -fin(5,i,j+1,k))+
6254 > ( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
6256 > ( fin(5,i,j,k+1) -fin(5,i,j+1,k+1))+
6257 > ( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
6259 sum=sum+z36th*hy*hz2*hxi(v)*(
6261 > -(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
6262 > (cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
6264 > -(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
6265 > (cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
6271 if(ict(3).eq.1)
then
6284 cyi=ypi*(ypi2-1.0_r8)
6294 czd=3.0_r8*zp2-1.0_r8
6295 czdi=-3.0_r8*zpi2+1.0_r8
6303 > +(ypi*fin(1,i,j,k) +yp*fin(1,i,j+1,k))
6304 > -(ypi*fin(1,i+1,j,k)+yp*fin(1,i+1,j+1,k)))
6306 > -(ypi*fin(1,i,j,k+1) +yp*fin(1,i,j+1,k+1))
6307 > +(ypi*fin(1,i+1,j,k+1)+yp*fin(1,i+1,j+1,k+1))))
6309 sum=sum+sixth*hy2*hzi*hxi(v)*(
6311 > +(cyi*fin(4,i,j,k) +cy*fin(4,i,j+1,k))
6312 > -(cyi*fin(4,i+1,j,k)+cy*fin(4,i+1,j+1,k)))
6314 > -(cyi*fin(4,i,j,k+1) +cy*fin(4,i,j+1,k+1))
6315 > +(cyi*fin(4,i+1,j,k+1)+cy*fin(4,i+1,j+1,k+1))))
6317 sum=sum+sixth*hz*hxi(v)*(
6319 > -(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))
6320 > +(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
6322 > -(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))
6323 > +(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
6325 sum=sum+z36th*hy2*hz*hxi(v)*(
6327 > -(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))
6328 > +(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
6330 > -(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))
6331 > +(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
6337 if(ict(4).eq.1)
then
6350 czi=zpi*(zpi2-1.0_r8)
6364 cxd=3.0_r8*xp2-1.0_r8
6365 cxdi=-3.0_r8*xpi2+1.0_r8
6369 > ( fin(2,i,j,k) -fin(2,i,j+1,k))
6370 > +(-fin(2,i+1,j,k)+fin(2,i+1,j+1,k)))
6372 > ( fin(2,i,j,k+1) -fin(2,i,j+1,k+1))
6373 > +(-fin(2,i+1,j,k+1)+fin(2,i+1,j+1,k+1))))
6375 sum=sum+sixth*hx(v)*hyi*(
6377 > cxdi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
6378 > cxd*( -fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
6380 > cxdi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
6381 > cxd*( -fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
6383 sum=sum+sixth*hz2*hxi(v)*hyi*(
6385 > ( fin(6,i,j,k) -fin(6,i,j+1,k))
6386 > +(-fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
6388 > ( fin(6,i,j,k+1) -fin(6,i,j+1,k+1))
6389 > +(-fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
6391 sum=sum+z36th*hx(v)*hz2*hyi*(
6393 > cxdi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
6394 > cxd*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
6396 > cxdi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
6397 > cxd*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
6403 if(ict(5).eq.1)
then
6416 cyi=ypi*(ypi2-1.0_r8)
6430 cxd=3.0_r8*xp2-1.0_r8
6431 cxdi=-3.0_r8*xpi2+1.0_r8
6435 > -(ypi*fin(3,i,j,k) +yp*fin(3,i,j+1,k))
6436 > +(ypi*fin(3,i+1,j,k)+yp*fin(3,i+1,j+1,k)))
6438 > -(ypi*fin(3,i,j,k+1) +yp*fin(3,i,j+1,k+1))
6439 > +(ypi*fin(3,i+1,j,k+1)+yp*fin(3,i+1,j+1,k+1))))
6441 sum=sum+sixth*hx(v)*hzi*(
6443 > cxdi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
6444 > cxd*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
6446 > cxdi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
6447 > cxd*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
6449 sum=sum+sixth*hy2*hxi(v)*hzi*(
6451 > -(cyi*fin(6,i,j,k) +cy*fin(6,i,j+1,k))
6452 > +(cyi*fin(6,i+1,j,k)+cy*fin(6,i+1,j+1,k)))
6454 > -(cyi*fin(6,i,j,k+1) +cy*fin(6,i,j+1,k+1))
6455 > +(cyi*fin(6,i+1,j,k+1)+cy*fin(6,i+1,j+1,k+1))))
6457 sum=sum+z36th*hx(v)*hy2*hzi*(
6459 > cxdi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
6460 > cxd*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
6462 > cxdi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
6463 > cxd*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
6469 if(ict(6).eq.1)
then
6481 czd=3.0_r8*zp2-1.0_r8
6482 czdi=-3.0_r8*zpi2+1.0_r8
6496 cxi=xpi*(xpi2-1.0_r8)
6501 > xpi*(-fin(2,i,j,k) +fin(2,i,j+1,k))+
6502 > xp*( -fin(2,i+1,j,k)+fin(2,i+1,j+1,k)))
6504 > xpi*(-fin(2,i,j,k+1) +fin(2,i,j+1,k+1))+
6505 > xp*( -fin(2,i+1,j,k+1)+fin(2,i+1,j+1,k+1))))
6507 sum=sum+sixth*hx2*hyi*hzi*(
6509 > cxi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
6510 > cx*( -fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
6512 > cxi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
6513 > cx*( -fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
6515 sum=sum+sixth*hz*hyi*(
6517 > xpi*(-fin(6,i,j,k) +fin(6,i,j+1,k))+
6518 > xp*( -fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
6520 > xpi*(-fin(6,i,j,k+1) +fin(6,i,j+1,k+1))+
6521 > xp*( -fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
6523 sum=sum+z36th*hx2*hz*hyi*(
6525 > cxi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
6526 > cx*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
6528 > cxi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
6529 > cx*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
6535 if(ict(7).eq.1)
then
6547 cyd=3.0_r8*yp2-1.0_r8
6548 cydi=-3.0_r8*ypi2+1.0_r8
6562 cxi=xpi*(xpi2-1.0_r8)
6567 > xpi*(-fin(3,i,j,k) +fin(3,i,j+1,k))+
6568 > xp*( -fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
6570 > xpi*(-fin(3,i,j,k+1) +fin(3,i,j+1,k+1))+
6571 > xp*( -fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
6573 sum=sum+sixth*hx2*hyi*hzi*(
6575 > cxi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
6576 > cx*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
6578 > cxi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
6579 > cx*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
6581 sum=sum+sixth*hy*hzi*(
6583 > xpi*(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))+
6584 > xp*(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
6586 > xpi*(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))+
6587 > xp*(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
6589 sum=sum+z36th*hx2*hy*hzi*(
6591 > cxi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
6592 > cx*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
6594 > cxi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
6595 > cx*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
6604 else if(ict(1).eq.5)
then
6605 if(ict(2).eq.1)
then
6622 czd=3.0_r8*zp2-1.0_r8
6623 czdi=-3.0_r8*zpi2+1.0_r8
6636 > xpi*(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))+
6637 > xp*( ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
6639 > xpi*(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))+
6640 > xp*(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
6644 > xpi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
6645 > xp*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
6647 > xpi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
6648 > xp*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
6654 if(ict(3).eq.1)
then
6666 cyd=3.0_r8*yp2-1.0_r8
6667 cydi=-3.0_r8*ypi2+1.0_r8
6685 > xpi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
6686 > xp*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
6688 > xpi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
6689 > xp*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
6693 > xpi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
6694 > xp*( cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
6696 > xpi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
6697 > xp*( cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
6703 if(ict(4).eq.1)
then
6729 cxd=3.0_r8*xp2-1.0_r8
6730 cxdi=-3.0_r8*xpi2+1.0_r8
6734 > -(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))
6735 > +(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
6737 > -(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))
6738 > +(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
6740 sum=sum+sixth*hx(v)*(
6742 > cxdi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
6743 > cxd*(ypi*fin(7,i+1,j,k) +yp*fin(7,i+1,j+1,k)))
6745 > cxdi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
6746 > cxd*( ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
6755 else if(ict(1).eq.-5)
then
6756 if(ict(2).eq.1)
then
6774 czi=zpi*(zpi2-1.0_r8)
6783 > -(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))
6784 > +( ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
6786 > -(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))
6787 > +(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
6789 sum=sum+sixth*hz2*hxi(v)*(
6791 > -(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))
6792 > +(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
6794 > -(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))
6795 > +(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
6801 if(ict(3).eq.1)
then
6813 cyd=3.0_r8*yp2-1.0_r8
6814 cydi=-3.0_r8*ypi2+1.0_r8
6823 czd=3.0_r8*zp2-1.0_r8
6824 czdi=-3.0_r8*zpi2+1.0_r8
6830 sum=hyi*hzi*hxi(v)*(
6832 > -(-fin(1,i,j,k) +fin(1,i,j+1,k))
6833 > +( -fin(1,i+1,j,k)+fin(1,i+1,j+1,k)))
6835 > -(-fin(1,i,j,k+1) +fin(1,i,j+1,k+1))
6836 > +( -fin(1,i+1,j,k+1)+fin(1,i+1,j+1,k+1))))
6838 sum=sum+sixth*hy*hzi*hxi(v)*(
6840 > -(cydi*fin(4,i,j,k) +cyd*fin(4,i,j+1,k))
6841 > +(cydi*fin(4,i+1,j,k)+cyd*fin(4,i+1,j+1,k)))
6843 > -(cydi*fin(4,i,j,k+1) +cyd*fin(4,i,j+1,k+1))
6844 > +(cydi*fin(4,i+1,j,k+1)+cyd*fin(4,i+1,j+1,k+1))))
6846 sum=sum+sixth*hz*hyi*hxi(v)*(
6848 > -(-fin(5,i,j,k) +fin(5,i,j+1,k))
6849 > +( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
6851 > -(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))
6852 > +( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
6854 sum=sum+z36th*hy*hz*hxi(v)*(
6856 > -(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))
6857 > +(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
6859 > -(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))
6860 > +(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
6866 if(ict(4).eq.1)
then
6879 cyi=ypi*(ypi2-1.0_r8)
6893 > -(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))
6894 > +(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
6896 > -(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))
6897 > +(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
6899 sum=sum+sixth*hy2*hxi(v)*(
6901 > -(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))
6902 > +(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
6904 > -(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))
6905 > +(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
6911 if(ict(5).eq.1)
then
6924 czi=zpi*(zpi2-1.0_r8)
6938 > xpi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
6939 > xp*( -fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
6941 > xpi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
6942 > xp*( -fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
6944 sum=sum+sixth*hz2*hyi*(
6946 > xpi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
6947 > xp*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
6949 > xpi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
6950 > xp*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
6956 if(ict(6).eq.1)
then
6969 cyi=ypi*(ypi2-1.0_r8)
6983 > xpi*(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))+
6984 > xp*(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
6986 > xpi*(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))+
6987 > xp*(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
6989 sum=sum+sixth*hy2*hzi*(
6991 > xpi*(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))+
6992 > xp*(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
6994 > xpi*(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))+
6995 > xp*(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
7001 if(ict(7).eq.1)
then
7013 czd=3.0_r8*zp2-1.0_r8
7014 czdi=-3.0_r8*zpi2+1.0_r8
7027 cxd=3.0_r8*xp2-1.0_r8
7028 cxdi=-3.0_r8*xpi2+1.0_r8
7030 sum=hxi(v)*hzi*hyi*(
7032 > ( fin(2,i,j,k) -fin(2,i,j+1,k))
7033 > +(-fin(2,i+1,j,k)+fin(2,i+1,j+1,k)))
7035 > ( fin(2,i,j,k+1) -fin(2,i,j+1,k+1))
7036 > +(-fin(2,i+1,j,k+1)+fin(2,i+1,j+1,k+1))))
7038 sum=sum+sixth*hx(v)*hzi*hyi*(
7040 > cxdi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
7041 > cxd*(-fin(4,i+1,j,k) +fin(4,i+1,j+1,k)))
7043 > cxdi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
7044 > cxd*(-fin(4,i+1,j,k+1) +fin(4,i+1,j+1,k+1))))
7046 sum=sum+sixth*hz*hxi(v)*hyi*(
7048 > ( fin(6,i,j,k) -fin(6,i,j+1,k))
7049 > +(-fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
7051 > ( fin(6,i,j,k+1) -fin(6,i,j+1,k+1))
7052 > +(-fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
7054 sum=sum+z36th*hx(v)*hz*hyi*(
7056 > cxdi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
7057 > cxd*(-fin(7,i+1,j,k) +fin(7,i+1,j+1,k)))
7059 > cxdi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
7060 > cxd*(-fin(7,i+1,j,k+1) +fin(7,i+1,j+1,k+1))))
7066 if(ict(8).eq.1)
then
7078 cyd=3.0_r8*yp2-1.0_r8
7079 cydi=-3.0_r8*ypi2+1.0_r8
7092 cxd=3.0_r8*xp2-1.0_r8
7093 cxdi=-3.0_r8*xpi2+1.0_r8
7095 sum=hyi*hxi(v)*hzi*(
7097 > ( +fin(3,i,j,k) -fin(3,i,j+1,k))
7098 > +( -fin(3,i+1,j,k)+fin(3,i+1,j+1,k)))
7100 > ( +fin(3,i,j,k+1) -fin(3,i,j+1,k+1))
7101 > +( -fin(3,i+1,j,k+1)+fin(3,i+1,j+1,k+1))))
7103 sum=sum+sixth*hx(v)*hyi*hzi*(
7105 > cxdi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
7106 > cxd*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
7108 > cxdi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
7109 > cxd*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
7111 sum=sum+sixth*hy*hxi(v)*hzi*(
7113 > -(cydi*fin(6,i,j,k) +cyd*fin(6,i,j+1,k))
7114 > +(cydi*fin(6,i+1,j,k)+cyd*fin(6,i+1,j+1,k)))
7116 > -(cydi*fin(6,i,j,k+1) +cyd*fin(6,i,j+1,k+1))
7117 > +(cydi*fin(6,i+1,j,k+1)+cyd*fin(6,i+1,j+1,k+1))))
7119 sum=sum+z36th*hx(v)*hy*hzi*(
7121 > cxdi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
7122 > cxd*(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
7124 > cxdi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
7125 > cxd*(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
7131 if(ict(9).eq.1)
then
7153 cxi=xpi*(xpi2-1.0_r8)
7158 > xpi*(-fin(6,i,j,k) +fin(6,i,j+1,k))+
7159 > xp*( -fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
7161 > xpi*(-fin(6,i,j,k+1) +fin(6,i,j+1,k+1))+
7162 > xp*( -fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
7164 sum=sum+sixth*hx2*hyi*(
7166 > cxi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
7167 > cx*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
7169 > cxi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
7170 > cx*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
7176 if(ict(10).eq.1)
then
7198 cxi=xpi*(xpi2-1.0_r8)
7203 > xpi*(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))+
7204 > xp*( ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
7206 > xpi*(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))+
7207 > xp*(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
7209 sum=sum+sixth*hx2*hzi*(
7211 > cxi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
7212 > cx*( ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
7214 > cxi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
7215 > cx*( ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
7224 else if(ict(1).eq.6)
then
7250 > xpi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
7251 > xp*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
7253 > xpi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
7254 > xp*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
7263 if(ict(1).eq.-6)
then
7264 if(ict(2).eq.1)
then
7277 czi=zpi*(zpi2-1.0_r8)
7286 > ( fin(4,i,j,k) -fin(4,i,j+1,k))
7287 > +(-fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
7289 > ( fin(4,i,j,k+1) -fin(4,i,j+1,k+1))
7290 > +(-fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
7292 sum=sum+sixth*hz2*hyi*hxi(v)*(
7294 > ( fin(7,i,j,k) -fin(7,i,j+1,k))
7295 > +(-fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
7297 > ( fin(7,i,j,k+1) -fin(7,i,j+1,k+1))
7298 > +(-fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
7304 if(ict(3).eq.1)
then
7321 czd=3.0_r8*zp2-1.0_r8
7322 czdi=-3.0_r8*zpi2+1.0_r8
7330 > -(ypi*fin(4,i,j,k) +yp*fin(4,i,j+1,k))
7331 > +( ypi*fin(4,i+1,j,k)+yp*fin(4,i+1,j+1,k)))
7333 > -(ypi*fin(4,i,j,k+1) +yp*fin(4,i,j+1,k+1))
7334 > +(ypi*fin(4,i+1,j,k+1)+yp*fin(4,i+1,j+1,k+1))))
7336 sum=sum+sixth*hz*hxi(v)*(
7338 > -(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))
7339 > +(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
7341 > -(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))
7342 > +(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
7348 if(ict(4).eq.1)
then
7360 cyd=3.0_r8*yp2-1.0_r8
7361 cydi=-3.0_r8*ypi2+1.0_r8
7374 > ( fin(5,i,j,k) -fin(5,i,j+1,k))
7375 > +(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
7377 > ( fin(5,i,j,k+1) -fin(5,i,j+1,k+1))
7378 > +(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
7380 sum=sum+sixth*hy*hxi(v)*(
7382 > -(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))
7383 > +(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
7385 > -(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))
7386 > +(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
7392 if(ict(5).eq.1)
then
7405 cyi=ypi*(ypi2-1.0_r8)
7414 > -(ypi*fin(5,i,j,k) +yp*fin(5,i,j+1,k))
7415 > +(ypi*fin(5,i+1,j,k)+yp*fin(5,i+1,j+1,k)))
7417 > -(ypi*fin(5,i,j,k+1) +yp*fin(5,i,j+1,k+1))
7418 > +(ypi*fin(5,i+1,j,k+1)+yp*fin(5,i+1,j+1,k+1))))
7420 sum=sum+sixth*hy2*hxi(v)*hzi*(
7422 > -(cyi*fin(7,i,j,k) +cy*fin(7,i,j+1,k))
7423 > +(cyi*fin(7,i+1,j,k)+cy*fin(7,i+1,j+1,k)))
7425 > -(cyi*fin(7,i,j,k+1) +cy*fin(7,i,j+1,k+1))
7426 > +(cyi*fin(7,i+1,j,k+1)+cy*fin(7,i+1,j+1,k+1))))
7432 if(ict(6).eq.1)
then
7444 czd=3.0_r8*zp2-1.0_r8
7445 czdi=-3.0_r8*zpi2+1.0_r8
7458 > xpi*(-fin(4,i,j,k) +fin(4,i,j+1,k))+
7459 > xp*(-fin(4,i+1,j,k) +fin(4,i+1,j+1,k)))
7461 > xpi*(-fin(4,i,j,k+1) +fin(4,i,j+1,k+1))+
7462 > xp*(-fin(4,i+1,j,k+1) +fin(4,i+1,j+1,k+1))))
7464 sum=sum+sixth*hz*hyi*(
7466 > xpi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
7467 > xp*(-fin(7,i+1,j,k) +fin(7,i+1,j+1,k)))
7469 > xpi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
7470 > xp*(-fin(7,i+1,j,k+1) +fin(7,i+1,j+1,k+1))))
7476 if(ict(7).eq.1)
then
7488 cyd=3.0_r8*yp2-1.0_r8
7489 cydi=-3.0_r8*ypi2+1.0_r8
7502 > xpi*(-fin(5,i,j,k) +fin(5,i,j+1,k))+
7503 > xp*( -fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
7505 > xpi*(-fin(5,i,j,k+1) +fin(5,i,j+1,k+1))+
7506 > xp*( -fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
7508 sum=sum+sixth*hy*hzi*(
7510 > xpi*(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))+
7511 > xp*( cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
7513 > xpi*(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))+
7514 > xp*( cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
7520 if(ict(8).eq.1)
then
7541 cxd=3.0_r8*xp2-1.0_r8
7542 cxdi=-3.0_r8*xpi2+1.0_r8
7546 > ( fin(6,i,j,k) -fin(6,i,j+1,k))
7547 > +(-fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
7549 > ( fin(6,i,j,k+1) -fin(6,i,j+1,k+1))
7550 > +(-fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
7552 sum=sum+sixth*hx(v)*hyi*(
7554 > cxdi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
7555 > cxd*(-fin(7,i+1,j,k) +fin(7,i+1,j+1,k)))
7557 > cxdi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
7558 > cxd*(-fin(7,i+1,j,k+1) +fin(7,i+1,j+1,k+1))))
7564 if(ict(9).eq.1)
then
7585 cxd=3.0_r8*xp2-1.0_r8
7586 cxdi=-3.0_r8*xpi2+1.0_r8
7590 > -(ypi*fin(6,i,j,k) +yp*fin(6,i,j+1,k))
7591 > +(ypi*fin(6,i+1,j,k)+yp*fin(6,i+1,j+1,k)))
7593 > -(ypi*fin(6,i,j,k+1) +yp*fin(6,i,j+1,k+1))
7594 > +(ypi*fin(6,i+1,j,k+1)+yp*fin(6,i+1,j+1,k+1))))
7596 sum=sum+sixth*hx(v)*hzi*(
7598 > cxdi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
7599 > cxd*(ypi*fin(7,i+1,j,k) +yp*fin(7,i+1,j+1,k)))
7601 > cxdi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
7602 > cxd*( ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
7608 if(ict(10).eq.1)
then
7625 cxi=xpi*(xpi2-1.0_r8)
7630 > xpi*(-fin(6,i,j,k) +fin(6,i,j+1,k))+
7631 > xp*( -fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
7633 > xpi*(-fin(6,i,j,k+1) +fin(6,i,j+1,k+1))+
7634 > xp*( -fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
7636 sum=sum+sixth*hx2*hyi*hzi*(
7638 > cxi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
7639 > cx*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
7641 > cxi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
7642 > cx*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
7651 else if(abs(ict(1)).eq.7)
then
7652 if(ict(2).eq.1)
then
7664 czd=3.0_r8*zp2-1.0_r8
7665 czdi=-3.0_r8*zpi2+1.0_r8
7671 sum=hyi*hxi(v)*hzi*(
7673 > ( fin(4,i,j,k) -fin(4,i,j+1,k))
7674 > +(-fin(4,i+1,j,k)+fin(4,i+1,j+1,k)))
7676 > ( fin(4,i,j,k+1) -fin(4,i,j+1,k+1))
7677 > +(-fin(4,i+1,j,k+1)+fin(4,i+1,j+1,k+1))))
7679 sum=sum+sixth*hz*hyi*hxi(v)*(
7681 > ( fin(7,i,j,k) -fin(7,i,j+1,k))
7682 > +(-fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
7684 > ( fin(7,i,j,k+1) -fin(7,i,j+1,k+1))
7685 > +(-fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
7691 if(ict(3).eq.1)
then
7712 > -(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))
7713 > +(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
7715 > -(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))
7716 > +(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
7722 if(ict(4).eq.1)
then
7734 cyd=3.0_r8*yp2-1.0_r8
7735 cydi=-3.0_r8*ypi2+1.0_r8
7741 sum=hxi(v)*hyi*hzi*(
7743 > ( fin(5,i,j,k) -fin(5,i,j+1,k))
7744 > +(-fin(5,i+1,j,k)+fin(5,i+1,j+1,k)))
7746 > ( fin(5,i,j,k+1) -fin(5,i,j+1,k+1))
7747 > +(-fin(5,i+1,j,k+1)+fin(5,i+1,j+1,k+1))))
7749 sum=sum+sixth*hy*hxi(v)*hzi*(
7751 > -(cydi*fin(7,i,j,k) +cyd*fin(7,i,j+1,k))
7752 > +(cydi*fin(7,i+1,j,k)+cyd*fin(7,i+1,j+1,k)))
7754 > -(cydi*fin(7,i,j,k+1) +cyd*fin(7,i,j+1,k+1))
7755 > +(cydi*fin(7,i+1,j,k+1)+cyd*fin(7,i+1,j+1,k+1))))
7761 if(ict(5).eq.1)
then
7782 > xpi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
7783 > xp*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
7785 > xpi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
7786 > xp*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
7792 if(ict(6).eq.1)
then
7813 > xpi*(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))+
7814 > xp*(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
7816 > xpi*(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))+
7817 > xp*(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
7823 if(ict(7).eq.1)
then
7839 cxd=3.0_r8*xp2-1.0_r8
7840 cxdi=-3.0_r8*xpi2+1.0_r8
7842 sum=hxi(v)*hyi*hzi*(
7844 > ( fin(6,i,j,k) -fin(6,i,j+1,k))
7845 > +(-fin(6,i+1,j,k)+fin(6,i+1,j+1,k)))
7847 > ( fin(6,i,j,k+1) -fin(6,i,j+1,k+1))
7848 > +(-fin(6,i+1,j,k+1)+fin(6,i+1,j+1,k+1))))
7850 sum=sum+sixth*hx(v)*hyi*hzi*(
7852 > cxdi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
7853 > cxd*(-fin(7,i+1,j,k) +fin(7,i+1,j+1,k)))
7855 > cxdi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
7856 > cxd*(-fin(7,i+1,j,k+1) +fin(7,i+1,j+1,k+1))))
7865 else if(abs(ict(1)).eq.8)
then
7866 if(ict(2).eq.1)
then
7882 > ( fin(7,i,j,k) -fin(7,i,j+1,k))
7883 > +(-fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
7885 > ( fin(7,i,j,k+1) -fin(7,i,j+1,k+1))
7886 > +(-fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
7892 if(ict(3).eq.1)
then
7908 > -(ypi*fin(7,i,j,k) +yp*fin(7,i,j+1,k))
7909 > +(ypi*fin(7,i+1,j,k)+yp*fin(7,i+1,j+1,k)))
7911 > -(ypi*fin(7,i,j,k+1) +yp*fin(7,i,j+1,k+1))
7912 > +(ypi*fin(7,i+1,j,k+1)+yp*fin(7,i+1,j+1,k+1))))
7918 if(ict(4).eq.1)
then
7934 > xpi*(-fin(7,i,j,k) +fin(7,i,j+1,k))+
7935 > xp*( -fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
7937 > xpi*(-fin(7,i,j,k+1) +fin(7,i,j+1,k+1))+
7938 > xp*( -fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))
7947 else if(abs(ict(1)).eq.9)
then
7956 sum=hyi*hxi(v)*hzi*(
7958 > ( fin(7,i,j,k) -fin(7,i,j+1,k))
7959 > +(-fin(7,i+1,j,k)+fin(7,i+1,j+1,k)))
7961 > ( fin(7,i,j,k+1) -fin(7,i,j+1,k+1))
7962 > +(-fin(7,i+1,j,k+1)+fin(7,i+1,j+1,k+1))))