1 subroutine evbicub(xget,yget,x,nx,y,ny,ilinx,iliny,
101 real xparam(1),yparam(1)
117 call herm2xy(xget,yget,x,nx,y,ny,ilinx,iliny,
118 > i(1),j(1),xparam(1),yparam(1),hx(1),hxi(1),hy(1),hyi(1),ier)
121 call fvbicub(ict,1,1,
122 > fval,i,j,xparam,yparam,hx,hxi,hy,hyi,f,inf2,ny)
130 subroutine fvbicub(ict,ivec,ivecd,
131 > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
138 integer ii(ivec),jj(ivec)
139 real xparam(ivec),yparam(ivec)
142 real hx(ivec),hy(ivec)
143 real hxi(ivec),hyi(ivec)
146 real fin(0:3,inf2,ny)
171 real,
parameter :: sixth = 0.166666666666666667
214 sum=xpi*(ypi*fin(0,i,j) +yp*fin(0,i,j+1))+
215 > xp*(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1))
218 > cxi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
219 > cx*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
222 > xpi*(cyi*fin(2,i,j) +cy*fin(2,i,j+1))+
223 > xp*(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
225 sum=sum+z36th*hx2*hy2*(
226 > cxi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
227 > cx*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
264 > -(ypi*fin(0,i,j) +yp*fin(0,i,j+1))
265 > +(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1)))
267 sum=sum+sixth*hx(v)*(
268 > cxdi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
269 > cxd*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
271 sum=sum+sixth*hxi(v)*hy2*(
272 > -(cyi*fin(2,i,j) +cy*fin(2,i,j+1))
273 > +(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
275 sum=sum+z36th*hx(v)*hy2*(
276 > cxdi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
277 > cxd*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
314 > xpi*(-fin(0,i,j) +fin(0,i,j+1))+
315 > xp*(-fin(0,i+1,j)+fin(0,i+1,j+1)))
317 sum=sum+sixth*hx2*hyi(v)*(
318 > cxi*(-fin(1,i,j) +fin(1,i,j+1))+
319 > cx*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
321 sum=sum+sixth*hy(v)*(
322 > xpi*(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))+
323 > xp*(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
325 sum=sum+z36th*hx2*hy(v)*(
326 > cxi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
327 > cx*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
359 > xpi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
360 > xp*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
363 > xpi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
364 > xp*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
396 > xpi*(ypi*fin(2,i,j) +yp*fin(2,i,j+1))+
397 > xp*(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
400 > cxi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
401 > cx*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
437 > fin(0,i,j) -fin(0,i,j+1)
438 > -fin(0,i+1,j)+fin(0,i+1,j+1))
440 sum=sum+sixth*hx(v)*hyi(v)*(
441 > cxdi*(-fin(1,i,j) +fin(1,i,j+1))+
442 > cxd*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
444 sum=sum+sixth*hxi(v)*hy(v)*(
445 > -(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))
446 > +(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
448 sum=sum+z36th*hx(v)*hy(v)*(
449 > cxdi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
450 > cxd*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
458 else if(ict(1).eq.3)
then
473 > -(ypi*fin(1,i,j) +yp*fin(1,i,j+1))
474 > +(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
476 sum=sum+sixth*hy2*hxi(v)*(
477 > -(cyi*fin(3,i,j) +cy*fin(3,i,j+1))
478 > +(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
500 > xpi*(-fin(1,i,j) +fin(1,i,j+1))+
501 > xp*(-fin(1,i+1,j) +fin(1,i+1,j+1)))
503 sum=sum+sixth*hy(v)*(
504 > xpi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
505 > xp*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
527 > -(ypi*fin(2,i,j) +yp*fin(2,i,j+1))
528 > +(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
530 sum=sum+sixth*hx(v)*(
531 > cxdi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
532 > cxd*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
555 > xpi*(-fin(2,i,j) +fin(2,i,j+1))+
556 > xp*(-fin(2,i+1,j) +fin(2,i+1,j+1)))
558 sum=sum+sixth*hx2*hyi(v)*(
559 > cxi*(-fin(3,i,j) +fin(3,i,j+1))+
560 > cx*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
569 else if(ict(1).eq.4)
then
584 > +( fin(1,i,j) -fin(1,i,j+1))
585 > +(-fin(1,i+1,j)+fin(1,i+1,j+1)))
587 sum=sum+sixth*hy(v)*hxi(v)*(
588 > -(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))
589 > +(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
607 sum=xpi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
608 > xp*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1))
630 > +( fin(2,i,j) -fin(2,i,j+1))
631 > +(-fin(2,i+1,j)+fin(2,i+1,j+1)))
633 sum=sum+sixth*hx(v)*hyi(v)*(
634 > cxdi*(-fin(3,i,j) +fin(3,i,j+1))+
635 > cxd*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
644 else if(ict(1).eq.5)
then
656 > -(ypi*fin(3,i,j) +yp*fin(3,i,j+1))
657 > +(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
674 > xpi*(-fin(3,i,j) +fin(3,i,j+1))+
675 > xp*(-fin(3,i+1,j)+fin(3,i+1,j+1)))
684 else if(ict(1).eq.6)
then
691 > +( fin(3,i,j) -fin(3,i,j+1))
692 > +(-fin(3,i+1,j)+fin(3,i+1,j+1)))
704 subroutine fvbicubx(ict,ivec,ivecd,
705 > fval,ii,jj,xparam,yparam,hx,hxi,hy,hyi,
713 real xparam(ivec),yparam
720 real fin(0:3,inf2,ny)
745 real,
parameter :: sixth = 0.166666666666666667
789 sum=xpi*(ypi*fin(0,i,j) +yp*fin(0,i,j+1))+
790 > xp*(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1))
793 > cxi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
794 > cx*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
797 > xpi*(cyi*fin(2,i,j) +cy*fin(2,i,j+1))+
798 > xp*(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
800 sum=sum+z36th*hx2*hy2*(
801 > cxi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
802 > cx*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
840 > -(ypi*fin(0,i,j) +yp*fin(0,i,j+1))
841 > +(ypi*fin(0,i+1,j)+yp*fin(0,i+1,j+1)))
843 sum=sum+sixth*hx(v)*(
844 > cxdi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
845 > cxd*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
847 sum=sum+sixth*hxi(v)*hy2*(
848 > -(cyi*fin(2,i,j) +cy*fin(2,i,j+1))
849 > +(cyi*fin(2,i+1,j)+cy*fin(2,i+1,j+1)))
851 sum=sum+z36th*hx(v)*hy2*(
852 > cxdi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
853 > cxd*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
891 > xpi*(-fin(0,i,j) +fin(0,i,j+1))+
892 > xp*(-fin(0,i+1,j)+fin(0,i+1,j+1)))
894 sum=sum+sixth*hx2*hyi*(
895 > cxi*(-fin(1,i,j) +fin(1,i,j+1))+
896 > cx*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
899 > xpi*(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))+
900 > xp*(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
902 sum=sum+z36th*hx2*hy*(
903 > cxi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
904 > cx*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
937 > xpi*(ypi*fin(1,i,j) +yp*fin(1,i,j+1))+
938 > xp*(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
941 > xpi*(cyi*fin(3,i,j) +cy*fin(3,i,j+1))+
942 > xp*(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
975 > xpi*(ypi*fin(2,i,j) +yp*fin(2,i,j+1))+
976 > xp*(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
979 > cxi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
980 > cx*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1017 > fin(0,i,j) -fin(0,i,j+1)
1018 > -fin(0,i+1,j)+fin(0,i+1,j+1))
1020 sum=sum+sixth*hx(v)*hyi*(
1021 > cxdi*(-fin(1,i,j) +fin(1,i,j+1))+
1022 > cxd*(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1024 sum=sum+sixth*hxi(v)*hy*(
1025 > -(cydi*fin(2,i,j) +cyd*fin(2,i,j+1))
1026 > +(cydi*fin(2,i+1,j)+cyd*fin(2,i+1,j+1)))
1028 sum=sum+z36th*hx(v)*hy*(
1029 > cxdi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1030 > cxd*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1038 else if(ict(1).eq.3)
then
1039 if(ict(2).eq.1)
then
1054 > -(ypi*fin(1,i,j) +yp*fin(1,i,j+1))
1055 > +(ypi*fin(1,i+1,j)+yp*fin(1,i+1,j+1)))
1057 sum=sum+sixth*hy2*hxi(v)*(
1058 > -(cyi*fin(3,i,j) +cy*fin(3,i,j+1))
1059 > +(cyi*fin(3,i+1,j)+cy*fin(3,i+1,j+1)))
1065 if(ict(3).eq.1)
then
1083 > xpi*(-fin(1,i,j) +fin(1,i,j+1))+
1084 > xp*(-fin(1,i+1,j) +fin(1,i+1,j+1)))
1087 > xpi*(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))+
1088 > xp*(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1094 if(ict(4).eq.1)
then
1111 > -(ypi*fin(2,i,j) +yp*fin(2,i,j+1))
1112 > +(ypi*fin(2,i+1,j)+yp*fin(2,i+1,j+1)))
1114 sum=sum+sixth*hx(v)*(
1115 > cxdi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1116 > cxd*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1122 if(ict(5).eq.1)
then
1139 > xpi*(-fin(2,i,j) +fin(2,i,j+1))+
1140 > xp*(-fin(2,i+1,j) +fin(2,i+1,j+1)))
1142 sum=sum+sixth*hx2*hyi*(
1143 > cxi*(-fin(3,i,j) +fin(3,i,j+1))+
1144 > cx*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
1153 else if(ict(1).eq.4)
then
1154 if(ict(2).eq.1)
then
1170 > +( fin(1,i,j) -fin(1,i,j+1))
1171 > +(-fin(1,i+1,j)+fin(1,i+1,j+1)))
1173 sum=sum+sixth*hy*hxi(v)*(
1174 > -(cydi*fin(3,i,j) +cyd*fin(3,i,j+1))
1175 > +(cydi*fin(3,i+1,j)+cyd*fin(3,i+1,j+1)))
1181 if(ict(3).eq.1)
then
1193 sum=xpi*(ypi*fin(3,i,j) +yp*fin(3,i,j+1))+
1194 > xp*(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1))
1200 if(ict(4).eq.1)
then
1217 > +( fin(2,i,j) -fin(2,i,j+1))
1218 > +(-fin(2,i+1,j)+fin(2,i+1,j+1)))
1220 sum=sum+sixth*hx(v)*hyi*(
1221 > cxdi*(-fin(3,i,j) +fin(3,i,j+1))+
1222 > cxd*(-fin(3,i+1,j) +fin(3,i+1,j+1)))
1231 else if(ict(1).eq.5)
then
1232 if(ict(2).eq.1)
then
1244 > -(ypi*fin(3,i,j) +yp*fin(3,i,j+1))
1245 > +(ypi*fin(3,i+1,j)+yp*fin(3,i+1,j+1)))
1251 if(ict(3).eq.1)
then
1262 > xpi*(-fin(3,i,j) +fin(3,i,j+1))+
1263 > xp*(-fin(3,i+1,j)+fin(3,i+1,j+1)))
1272 else if(ict(1).eq.6)
then
1279 > +( fin(3,i,j) -fin(3,i,j+1))
1280 > +(-fin(3,i+1,j)+fin(3,i+1,j+1)))