1 subroutine evintrp3d(xget,yget,zget,x,nx,y,ny,z,nz,jspline,
2 > f,icoeff,ixdim,iydim,izdim,ict,fval,ier)
30 integer :: ixdim,iydim,izdim
35 real f(icoeff,ixdim,iydim,izdim)
122 real xparam,yparam,zparam
136 call vecin3d_argchk(
'evintrp3d',jspline,
137 > icoeff,nx,ny,nz,ixdim,iydim,izdim,ier)
142 call herm3xyz(xget,yget,zget,x,nx,y,ny,z,nz,0,0,0,
143 > i,j,k,xparam,yparam,zparam,hx,hxi,hy,hyi,hz,hzi,ier)
146 call fvintrp3d(ict,1,1,
147 > fval,i,j,k,xparam,yparam,zparam,
148 > hx,hxi,hy,hyi,hz,hzi,
149 > jspline,f,icoeff,ixdim,iydim,izdim)
158 subroutine fvintrp3d(ict,ivec,ivecd,
159 > fval,ii,jj,kk,xparam,yparam,zparam,
160 > hx,hxi,hy,hyi,hz,hzi,
161 > jspline,fin,icoeff,ixdim,iydim,izdim)
171 integer ii(ivec),jj(ivec),kk(ivec)
172 real xparam(ivec),yparam(ivec),zparam(ivec)
175 real hx(ivec),hy(ivec),hz(ivec)
176 real hxi(ivec),hyi(ivec),hzi(ivec)
179 integer :: jspline(3)
181 integer :: ixdim,iydim,izdim
182 real fin(icoeff,ixdim,iydim,izdim)
200 integer :: i,j,k1,k2,imaxx,imaxy,imaxz,linrank,cubrank,zonrank
201 logical :: splin_flag
204 integer :: iderivs(3,10),ict8(8),maxd_lin(8),maxd_cub(8)
205 integer :: idcub,idlin,idcub2(2),idlin2(2)
207 real :: xp,yp,hhx,hhy,hhxi,hhyi,fac1,fac2(2)
208 real :: h,hi,hhi,h2(2),h2i(2)
209 real :: f22(2,2),f22b(2,2),f2222(2,2,2,2),ans22(2,2),ans1,ans2
210 real :: f422(4,2,2),f422b(4,2,2)
212 real,
parameter :: ONE = 1.0d0
213 real,
parameter :: sixth = 0.166666666666666667
214 real,
parameter :: z36th = sixth*sixth
223 if(jspline(i).eq.-1)
then
224 zonrank = zonrank + 1
225 else if(jspline(i).eq.0)
then
226 linrank = linrank + 1
228 splin_flag = jspline(i).eq.2
229 cubrank = cubrank + 1
233 if(cubrank.eq.3)
then
236 call fvtricub(ict,ivec,ivecd,
237 > fval,ii,jj,kk,xparam,yparam,zparam,
238 > hx,hxi,hy,hyi,hz,hzi,
239 > fin,ixdim,iydim,izdim)
245 if(maxval(iderivs(1:3,i)).le.1)
then
248 call herm3fcn(ict8,ivec,ivecd,
249 > fval(1,i),ii,jj,kk,xparam,yparam,zparam,
250 > hx,hxi,hy,hyi,hz,hzi,
251 > fin,ixdim,iydim,izdim)
264 if(cubrank.eq.0)
then
265 if(jspline(1).eq.-1)
then
270 if(jspline(2).eq.-1)
then
275 if(jspline(3).eq.-1)
then
282 if((iderivs(1,i).le.imaxx).and.(iderivs(2,i).le.imaxy)
283 > .and.(iderivs(3,i).le.imaxz))
then
284 if(linrank.eq.3)
then
288 call pc3fcn(ict8,ivec,ivecd,
289 > fval(1,i), ii,jj,kk, xparam,yparam,zparam,
290 > hx,hxi, hy,hyi, hz,hzi,
291 > fin, ixdim,iydim,izdim)
292 else if(zonrank.eq.3)
then
295 fval(j,i)=fin(1,ii(j),jj(j),kk(j))
297 else if(linrank.eq.1)
then
298 if(jspline(1).eq.0)
then
301 if(iderivs(1,i).eq.0)
then
303 > (one-xparam(j))*fin(1,ii(j),jj(j),kk(j))
304 > +xparam(j)*fin(1,ii(j)+1,jj(j),kk(j))
307 > (fin(1,ii(j)+1,jj(j),kk(j))-
308 > fin(1,ii(j),jj(j),kk(j)))
311 else if(jspline(2).eq.0)
then
314 if(iderivs(2,i).eq.0)
then
316 > (one-yparam(j))*fin(1,ii(j),jj(j),kk(j))
317 > +yparam(j)*fin(1,ii(j),jj(j)+1,kk(j))
320 > (fin(1,ii(j),jj(j)+1,kk(j))-
321 > fin(1,ii(j),jj(j),kk(j)))
327 if(iderivs(3,i).eq.0)
then
329 > (one-zparam(j))*fin(1,ii(j),jj(j),kk(j))
330 > +zparam(j)*fin(1,ii(j),jj(j),kk(j)+1)
333 > (fin(1,ii(j),jj(j),kk(j)+1)-
334 > fin(1,ii(j),jj(j),kk(j)))
339 if(jspline(1).eq.-1)
then
342 if((iderivs(2,i).eq.0).and.
343 > (iderivs(3,i).eq.0))
then
344 fval(j,i) = (one-zparam(j))*
346 > fin(1,ii(j),jj(j),kk(j))
347 > +yparam(j)*fin(1,ii(j),jj(j)+1,kk(j)))
350 > fin(1,ii(j),jj(j),kk(j)+1)
351 > +yparam(j)*fin(1,ii(j),jj(j)+1,kk(j)+1))
353 else if((iderivs(2,i).eq.1).and.
354 > (iderivs(3,i).eq.0))
then
356 > (one-zparam(j))*fin(1,ii(j),jj(j)+1,kk(j))
357 > +zparam(j)*fin(1,ii(j),jj(j)+1,kk(j)+1)
358 > -(one-zparam(j))*fin(1,ii(j),jj(j),kk(j))
359 > -zparam(j)*fin(1,ii(j),jj(j),kk(j)+1) )
361 else if((iderivs(2,i).eq.0).and.
362 > (iderivs(3,i).eq.1))
then
364 > (one-yparam(j))*fin(1,ii(j),jj(j),kk(j)+1)
365 > +yparam(j)*fin(1,ii(j),jj(j)+1,kk(j)+1)
366 > -(one-yparam(j))*fin(1,ii(j),jj(j),kk(j))
367 > -yparam(j)*fin(1,ii(j),jj(j)+1,kk(j)) )
370 fval(j,i) = hyi(j)*hzi(j)*(
371 > fin(1,ii(j),jj(j),kk(j))
372 > -fin(1,ii(j),jj(j)+1,kk(j))
373 > -fin(1,ii(j),jj(j),kk(j)+1)
374 > +fin(1,ii(j),jj(j)+1,kk(j)+1))
378 else if(jspline(2).eq.-1)
then
381 if((iderivs(1,i).eq.0).and.
382 > (iderivs(3,i).eq.0))
then
383 fval(j,i) = (one-zparam(j))*
385 > fin(1,ii(j),jj(j),kk(j))
386 > +xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)))
389 > fin(1,ii(j),jj(j),kk(j)+1)
390 > +xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)+1))
392 else if((iderivs(1,i).eq.1).and.
393 > (iderivs(3,i).eq.0))
then
395 > (one-zparam(j))*fin(1,ii(j)+1,jj(j),kk(j))
396 > +zparam(j)*fin(1,ii(j)+1,jj(j),kk(j)+1)
397 > -(one-zparam(j))*fin(1,ii(j),jj(j),kk(j))
398 > -zparam(j)*fin(1,ii(j),jj(j),kk(j)+1) )
400 else if((iderivs(1,i).eq.0).and.
401 > (iderivs(3,i).eq.1))
then
403 > (one-xparam(j))*fin(1,ii(j),jj(j),kk(j)+1)
404 > +xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)+1)
405 > -(one-xparam(j))*fin(1,ii(j),jj(j),kk(j))
406 > -xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)) )
409 fval(j,i) = hxi(j)*hzi(j)*(
410 > fin(1,ii(j),jj(j),kk(j))
411 > -fin(1,ii(j)+1,jj(j),kk(j))
412 > -fin(1,ii(j),jj(j),kk(j)+1)
413 > +fin(1,ii(j)+1,jj(j),kk(j)+1))
420 if((iderivs(1,i).eq.0).and.
421 > (iderivs(2,i).eq.0))
then
422 fval(j,i) = (one-yparam(j))*
424 > fin(1,ii(j),jj(j),kk(j))
425 > +xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)))
428 > fin(1,ii(j),jj(j)+1,kk(j))
429 > +xparam(j)*fin(1,ii(j)+1,jj(j)+1,kk(j)))
431 else if((iderivs(1,i).eq.1).and.
432 > (iderivs(2,i).eq.0))
then
434 > (one-yparam(j))*fin(1,ii(j)+1,jj(j),kk(j))
435 > +yparam(j)*fin(1,ii(j)+1,jj(j)+1,kk(j))
436 > -(one-yparam(j))*fin(1,ii(j),jj(j),kk(j))
437 > -yparam(j)*fin(1,ii(j),jj(j)+1,kk(j)) )
439 else if((iderivs(1,i).eq.0).and.
440 > (iderivs(2,i).eq.1))
then
442 > (one-xparam(j))*fin(1,ii(j),jj(j)+1,kk(j))
443 > +xparam(j)*fin(1,ii(j)+1,jj(j)+1,kk(j))
444 > -(one-xparam(j))*fin(1,ii(j),jj(j),kk(j))
445 > -xparam(j)*fin(1,ii(j)+1,jj(j),kk(j)) )
448 fval(j,i) = hxi(j)*hyi(j)*(
449 > fin(1,ii(j),jj(j),kk(j))
450 > -fin(1,ii(j)+1,jj(j),kk(j))
451 > -fin(1,ii(j),jj(j)+1,kk(j))
452 > +fin(1,ii(j)+1,jj(j)+1,kk(j)))
469 if(jspline(1).eq.2)
then
471 else if(jspline(1).eq.-1)
then
477 if(jspline(2).eq.2)
then
479 else if(jspline(2).eq.-1)
then
485 if(jspline(3).eq.2)
then
487 else if(jspline(3).eq.-1)
then
493 if(cubrank.eq.1)
then
497 if(jspline(1).ge.1)
then
499 idlin2(1)=iderivs(2,i)
500 idlin2(2)=iderivs(3,i)
501 else if(jspline(2).ge.1)
then
503 idlin2(1)=iderivs(1,i)
504 idlin2(2)=iderivs(3,i)
505 else if(jspline(3).ge.1)
then
507 idlin2(1)=iderivs(1,i)
508 idlin2(2)=iderivs(2,i)
511 if((iderivs(1,i).le.imaxx).and.(iderivs(2,i).le.imaxy)
512 > .and.(iderivs(3,i).le.imaxz))
then
513 if(linrank.eq.2)
then
515 if(jspline(1).ge.1)
then
523 f2222(1:2,1,1:2,1:2) =
524 > fin(1:2,ii(j),jj(j):jj(j)+1,kk(j):kk(j)+1)
525 f2222(1:2,2,1:2,1:2) =
526 > fin(1:2,ii(j)+1,jj(j):jj(j)+1,kk(j):kk(j)+1)
527 else if(jspline(2).ge.1)
then
535 f2222(1:2,1,1:2,1:2) =
536 > fin(1:2,ii(j):ii(j)+1,jj(j),kk(j):kk(j)+1)
537 f2222(1:2,2,1:2,1:2) =
538 > fin(1:2,ii(j):ii(j)+1,jj(j)+1,kk(j):kk(j)+1)
539 else if(jspline(3).ge.1)
then
547 f2222(1:2,1,1:2,1:2) =
548 > fin(1:2,ii(j):ii(j)+1,jj(j):jj(j)+1,kk(j))
549 f2222(1:2,2,1:2,1:2) =
550 > fin(1:2,ii(j):ii(j)+1,jj(j):jj(j)+1,kk(j)+1)
557 call sp1d(ans22(k1,k2),
558 > f2222(1:2,1:2,k1,k2))
560 call hm1d(ans22(k1,k2),
561 > f2222(1:2,1:2,k1,k2))
567 if((idlin2(1).eq.0).and.(idlin2(2).eq.0))
then
569 ans1 = (one-fac2(2))*
570 > ((one-fac2(1))*ans22(1,1)+fac2(1)*ans22(2,1))
572 > ((one-fac2(1))*ans22(1,2)+fac2(1)*ans22(2,2))
573 else if((idlin2(1).eq.1).and.(idlin2(2).eq.0))
then
575 ans1 = (one-fac2(2))*(ans22(2,1)-ans22(1,1))
576 > +fac2(2)*(ans22(2,2)-ans22(1,2))
578 else if((idlin2(1).eq.0).and.(idlin2(2).eq.1))
then
580 ans1 = (one-fac2(1))*(ans22(1,2)-ans22(1,1))
581 > +fac2(1)*(ans22(2,2)-ans22(2,1))
585 ans1 = ans22(2,2)-ans22(1,2)
586 > -ans22(2,1)+ans22(1,1)
587 ans1 = ans1*h2i(1)*h2i(2)
592 else if(linrank.eq.0)
then
594 if(jspline(1).ge.1)
then
598 f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
599 f22(1:2,2)=fin(1:2,ii(j)+1,jj(j),kk(j))
600 else if(jspline(2).ge.1)
then
604 f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
605 f22(1:2,2)=fin(1:2,ii(j),jj(j)+1,kk(j))
606 else if(jspline(3).ge.1)
then
610 f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
611 f22(1:2,2)=fin(1:2,ii(j),jj(j),kk(j)+1)
621 else if(linrank.eq.1)
then
623 if(jspline(1).ge.1)
then
627 if(jspline(2).eq.0)
then
631 f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
632 f22(1:2,2)=fin(1:2,ii(j)+1,jj(j),kk(j))
633 f22b(1:2,1)=fin(1:2,ii(j),jj(j)+1,kk(j))
634 f22b(1:2,2)=fin(1:2,ii(j)+1,jj(j)+1,kk(j))
639 f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
640 f22(1:2,2)=fin(1:2,ii(j)+1,jj(j),kk(j))
641 f22b(1:2,1)=fin(1:2,ii(j),jj(j),kk(j)+1)
642 f22b(1:2,2)=fin(1:2,ii(j)+1,jj(j),kk(j)+1)
644 else if(jspline(2).ge.1)
then
648 if(jspline(1).eq.0)
then
652 f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
653 f22(1:2,2)=fin(1:2,ii(j),jj(j)+1,kk(j))
654 f22b(1:2,1)=fin(1:2,ii(j)+1,jj(j),kk(j))
655 f22b(1:2,2)=fin(1:2,ii(j)+1,jj(j)+1,kk(j))
660 f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
661 f22(1:2,2)=fin(1:2,ii(j),jj(j)+1,kk(j))
662 f22b(1:2,1)=fin(1:2,ii(j),jj(j),kk(j)+1)
663 f22b(1:2,2)=fin(1:2,ii(j),jj(j)+1,kk(j)+1)
665 else if(jspline(3).ge.1)
then
669 if(jspline(1).eq.0)
then
673 f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
674 f22(1:2,2)=fin(1:2,ii(j),jj(j),kk(j)+1)
675 f22b(1:2,1)=fin(1:2,ii(j)+1,jj(j),kk(j))
676 f22b(1:2,2)=fin(1:2,ii(j)+1,jj(j),kk(j)+1)
681 f22(1:2,1)=fin(1:2,ii(j),jj(j),kk(j))
682 f22(1:2,2)=fin(1:2,ii(j),jj(j),kk(j)+1)
683 f22b(1:2,1)=fin(1:2,ii(j),jj(j)+1,kk(j))
684 f22b(1:2,2)=fin(1:2,ii(j),jj(j)+1,kk(j)+1)
695 fval(j,i)=(one-fac1)*ans1 + fac1*ans2
697 fval(j,i)=(ans2-ans1)*hhi
712 if(jspline(1).le.0)
then
714 idcub2(1)=iderivs(2,i)
715 idcub2(2)=iderivs(3,i)
716 else if(jspline(2).le.0)
then
718 idcub2(1)=iderivs(1,i)
719 idcub2(2)=iderivs(3,i)
720 else if(jspline(3).le.0)
then
722 idcub2(1)=iderivs(1,i)
723 idcub2(2)=iderivs(2,i)
726 if((iderivs(1,i).le.imaxx).and.(iderivs(2,i).le.imaxy)
727 > .and.(iderivs(3,i).le.imaxz))
then
729 if(linrank.eq.1)
then
732 if(jspline(1).le.0)
then
741 f422 = fin(1:4,ii(j),jj(j):jj(j)+1,kk(j):kk(j)+1)
742 f422b= fin(1:4,ii(j)+1,jj(j):jj(j)+1,kk(j):kk(j)+1)
743 else if(jspline(2).le.0)
then
752 f422 = fin(1:4,ii(j):ii(j)+1,jj(j),kk(j):kk(j)+1)
753 f422b= fin(1:4,ii(j):ii(j)+1,jj(j)+1,kk(j):kk(j)+1)
754 else if(jspline(3).le.0)
then
763 f422 = fin(1:4,ii(j):ii(j)+1,jj(j):jj(j)+1,kk(j))
764 f422b= fin(1:4,ii(j):ii(j)+1,jj(j):jj(j)+1,kk(j)+1)
768 call sp2d(ans2,f422b)
771 call hm2d(ans2,f422b)
774 fval(j,i) = (ans2-ans1)*hi
776 fval(j,i) = (one-fac1)*ans1 + fac1*ans2
782 if(jspline(1).le.0)
then
789 f422 = fin(1:4,ii(j),jj(j):jj(j)+1,kk(j):kk(j)+1)
790 else if(jspline(2).le.0)
then
797 f422 = fin(1:4,ii(j):ii(j)+1,jj(j),kk(j):kk(j)+1)
798 else if(jspline(3).le.0)
then
805 f422 = fin(1:4,ii(j):ii(j)+1,jj(j):jj(j)+1,kk(j))
822 if(iderivs(1,i).eq.0)
then
823 if(iderivs(2,i).eq.0)
then
824 if(iderivs(3,i).eq.0)
then
830 if(iderivs(3,i).eq.0)
then
837 if(iderivs(2,i).eq.0)
then
838 if(iderivs(3,i).eq.0)
then
844 if(iderivs(3,i).eq.0)
then
851 end subroutine gen_ict8
866 if(abs(ict(1)).le.(2))
then
894 if(ict(10).eq.1)
then
898 else if(ict(1).eq.3)
then
921 else if(ict(1).eq.-3)
then
932 else if(ict(1).eq.4)
then
952 else if(ict(1).eq.-4)
then
972 else if(ict(1).eq.5)
then
983 else if(ict(1).eq.-5)
then
1002 if(ict(8).eq.1)
then
1005 if(ict(9).eq.1)
then
1008 if(ict(10).eq.1)
then
1012 else if(ict(1).eq.6)
then
1015 else if(ict(1).eq.-6)
then
1016 if(ict(2).eq.1)
then
1019 if(ict(3).eq.1)
then
1022 if(ict(4).eq.1)
then
1025 if(ict(5).eq.1)
then
1028 if(ict(6).eq.1)
then
1031 if(ict(7).eq.1)
then
1034 if(ict(8).eq.1)
then
1037 if(ict(9).eq.1)
then
1040 if(ict(10).eq.1)
then
1044 else if(abs(ict(1)).eq.7)
then
1045 if(ict(2).eq.1)
then
1048 if(ict(3).eq.1)
then
1051 if(ict(4).eq.1)
then
1054 if(ict(5).eq.1)
then
1057 if(ict(6).eq.1)
then
1060 if(ict(7).eq.1)
then
1064 else if(abs(ict(1)).eq.8)
then
1065 if(ict(2).eq.1)
then
1068 if(ict(3).eq.1)
then
1071 if(ict(4).eq.1)
then
1075 else if(abs(ict(1)).eq.9)
then
1080 end subroutine dtrans
1082 subroutine add1(idx,idy,idz)
1085 integer,
intent(in) :: idx,idy,idz
1090 if(jspline(1).le.0)
then
1099 if(jspline(2).le.0)
then
1100 maxd_lin(inum)=max(maxd_lin(inum),idy)
1102 maxd_cub(inum)=max(maxd_cub(inum),idy)
1107 if(jspline(3).le.0)
then
1108 maxd_lin(inum)=max(maxd_lin(inum),idz)
1110 maxd_cub(inum)=max(maxd_cub(inum),idz)
1116 subroutine sp1d(ans,f22)
1117 real,
intent(in) :: f22(2,2)
1118 real,
intent(out) :: ans
1124 real :: xpi,xp2,xpi2,cx,cxi,hx2,cxd,cxdi
1139 ans=xpi*f22(1,1) +xp*f22(1,2)
1140 ans=ans+sixth*hx2*(cxi*f22(2,1) +cx*f22(2,2))
1142 else if(idcub.eq.1)
then
1148 ans=hi*(f22(1,2)-f22(1,1))
1149 ans=ans+sixth*h*(cxdi*f22(2,1) +cxd*f22(2,2))
1151 else if(idcub.eq.2)
then
1155 ans=xpi*f22(2,1) + xp*f22(2,2)
1160 ans = hi*(f22(2,2)-f22(2,1))
1167 subroutine hm1d(ans,f22)
1168 real,
intent(in) :: f22(2,2)
1169 real,
intent(out) :: ans
1175 real :: xpi,xp2,xpi2,ax,axbar,bx,bxbar
1176 real :: axp,axbarp,bxp,bxbarp
1190 ans=axbar*f22(1,1) + ax*f22(1,2)
1191 ans=ans + h*(bxbar*f22(2,1) + bx*f22(2,2))
1199 bxbarp=xpi*(3.0*xpi-2.0)
1201 ans=hi*(axbarp*f22(1,1) +axp*f22(1,2))
1202 ans=ans + bxbarp*f22(2,1) + bxp*f22(2,2)
1209 subroutine sp2d(ans,ff)
1210 real,
intent(in) :: ff(0:3,2,2)
1211 real,
intent(out) :: ans
1217 real xpi,xp2,xpi2,cx,cxi,hx2,cxd,cxdi
1218 real ypi,yp2,ypi2,cy,cyi,hy2,cyd,cydi
1232 if((idcub2(1).eq.0).and.(idcub2(2).eq.0))
then
1244 ans=xpi*(ypi*ff(0,1,1) +yp*ff(0,1,2))+
1245 > xp*(ypi*ff(0,2,1)+yp*ff(0,2,2))
1248 > cxi*(ypi*ff(1,1,1) +yp*ff(1,1,2))+
1249 > cx*(ypi*ff(1,2,1)+yp*ff(1,2,2)))
1252 > xpi*(cyi*ff(2,1,1) +cy*ff(2,1,2))+
1253 > xp*(cyi*ff(2,2,1)+cy*ff(2,2,2)))
1255 ans=ans+z36th*hx2*hy2*(
1256 > cxi*(cyi*ff(3,1,1) +cy*ff(3,1,2))+
1257 > cx*(cyi*ff(3,2,1)+cy*ff(3,2,2)))
1260 else if((idcub2(1).eq.1).and.(idcub2(2).eq.0))
then
1272 > -(ypi*ff(0,1,1) +yp*ff(0,1,2))
1273 > +(ypi*ff(0,2,1)+yp*ff(0,2,2)))
1276 > cxdi*(ypi*ff(1,1,1) +yp*ff(1,1,2))+
1277 > cxd*(ypi*ff(1,2,1)+yp*ff(1,2,2)))
1279 ans=ans+sixth*hhxi*hy2*(
1280 > -(cyi*ff(2,1,1) +cy*ff(2,1,2))
1281 > +(cyi*ff(2,2,1)+cy*ff(2,2,2)))
1283 ans=ans+z36th*hhx*hy2*(
1284 > cxdi*(cyi*ff(3,1,1) +cy*ff(3,1,2))+
1285 > cxd*(cyi*ff(3,2,1)+cy*ff(3,2,2)))
1288 else if((idcub2(1).eq.0).and.(idcub2(2).eq.1))
then
1300 > xpi*(-ff(0,1,1) +ff(0,1,2))+
1301 > xp*(-ff(0,2,1)+ff(0,2,2)))
1303 ans=ans+sixth*hx2*hhyi*(
1304 > cxi*(-ff(1,1,1) +ff(1,1,2))+
1305 > cx*(-ff(1,2,1)+ff(1,2,2)))
1308 > xpi*(cydi*ff(2,1,1) +cyd*ff(2,1,2))+
1309 > xp*(cydi*ff(2,2,1)+cyd*ff(2,2,2)))
1311 ans=ans+z36th*hx2*hhy*(
1312 > cxi*(cydi*ff(3,1,1) +cyd*ff(3,1,2))+
1313 > cx*(cydi*ff(3,2,1)+cyd*ff(3,2,2)))
1316 else if((idcub2(1).eq.2).and.(idcub2(2).eq.0))
then
1325 > xpi*(ypi*ff(1,1,1) +yp*ff(1,1,2))+
1326 > xp*(ypi*ff(1,2,1)+yp*ff(1,2,2)))
1329 > xpi*(cyi*ff(3,1,1) +cy*ff(3,1,2))+
1330 > xp*(cyi*ff(3,2,1)+cy*ff(3,2,2)))
1333 else if((idcub2(1).eq.0).and.(idcub2(2).eq.2))
then
1342 > xpi*(ypi*ff(2,1,1) +yp*ff(2,1,2))+
1343 > xp*(ypi*ff(2,2,1)+yp*ff(2,2,2)))
1346 > cxi*(ypi*ff(3,1,1) +yp*ff(3,1,2))+
1347 > cx*(ypi*ff(3,2,1)+yp*ff(3,2,2)))
1350 else if((idcub2(1).eq.1).and.(idcub2(2).eq.1))
then
1361 > ff(0,1,1) -ff(0,1,2)
1362 > -ff(0,2,1)+ff(0,2,2))
1364 ans=ans+sixth*hhx*hhyi*(
1365 > cxdi*(-ff(1,1,1) +ff(1,1,2))+
1366 > cxd*(-ff(1,2,1)+ff(1,2,2)))
1368 ans=ans+sixth*hhxi*hhy*(
1369 > -(cydi*ff(2,1,1) +cyd*ff(2,1,2))
1370 > +(cydi*ff(2,2,1)+cyd*ff(2,2,2)))
1372 ans=ans+z36th*hhx*hhy*(
1373 > cxdi*(cydi*ff(3,1,1) +cyd*ff(3,1,2))+
1374 > cxd*(cydi*ff(3,2,1)+cyd*ff(3,2,2)))
1376 else if((idcub2(1).eq.3).and.(idcub2(2).eq.0))
then
1385 > -(ypi*ff(1,1,1) +yp*ff(1,1,2))
1386 > +(ypi*ff(1,2,1)+yp*ff(1,2,2)))
1388 ans=ans+sixth*hy2*hhxi*(
1389 > -(cyi*ff(3,1,1) +cy*ff(3,1,2))
1390 > +(cyi*ff(3,2,1)+cy*ff(3,2,2)))
1393 else if((idcub2(1).eq.2).and.(idcub2(2).eq.1))
then
1401 > xpi*(-ff(1,1,1) +ff(1,1,2))+
1402 > xp*(-ff(1,2,1) +ff(1,2,2)))
1405 > xpi*(cydi*ff(3,1,1) +cyd*ff(3,1,2))+
1406 > xp*(cydi*ff(3,2,1)+cyd*ff(3,2,2)))
1409 else if((idcub2(1).eq.1).and.(idcub2(2).eq.2))
then
1417 > -(ypi*ff(2,1,1) +yp*ff(2,1,2))
1418 > +(ypi*ff(2,2,1)+yp*ff(2,2,2)))
1421 > cxdi*(ypi*ff(3,1,1) +yp*ff(3,1,2))+
1422 > cxd*(ypi*ff(3,2,1)+yp*ff(3,2,2)))
1424 else if((idcub2(1).eq.0).and.(idcub2(2).eq.3))
then
1433 > xpi*(-ff(2,1,1) +ff(2,1,2))+
1434 > xp*(-ff(2,2,1) +ff(2,2,2)))
1436 ans=ans+sixth*hx2*hhyi*(
1437 > cxi*(-ff(3,1,1) +ff(3,1,2))+
1438 > cx*(-ff(3,2,1) +ff(3,2,2)))
1441 else if((idcub2(1).eq.3).and.(idcub2(2).eq.1))
then
1449 > +( ff(1,1,1) -ff(1,1,2))
1450 > +(-ff(1,2,1)+ff(1,2,2)))
1452 ans=ans+sixth*hhy*hhxi*(
1453 > -(cydi*ff(3,1,1) +cyd*ff(3,1,2))
1454 > +(cydi*ff(3,2,1)+cyd*ff(3,2,2)))
1457 else if((idcub2(1).eq.2).and.(idcub2(2).eq.2))
then
1461 ans=xpi*(ypi*ff(3,1,1) +yp*ff(3,1,2))+
1462 > xp*(ypi*ff(3,2,1) +yp*ff(3,2,2))
1464 else if((idcub2(1).eq.1).and.(idcub2(2).eq.3))
then
1472 > +( ff(2,1,1) -ff(2,1,2))
1473 > +(-ff(2,2,1)+ff(2,2,2)))
1475 ans=ans+sixth*hhx*hhyi*(
1476 > cxdi*(-ff(3,1,1) +ff(3,1,2))+
1477 > cxd*(-ff(3,2,1) +ff(3,2,2)))
1480 else if((idcub2(1).eq.3).and.(idcub2(2).eq.2))
then
1485 > -(ypi*ff(3,1,1) +yp*ff(3,1,2))
1486 > +(ypi*ff(3,2,1)+yp*ff(3,2,2)))
1488 else if((idcub2(1).eq.2).and.(idcub2(2).eq.3))
then
1493 > xpi*(-ff(3,1,1) +ff(3,1,2))+
1494 > xp*(-ff(3,2,1)+ff(3,2,2)))
1497 else if((idcub2(1).eq.3).and.(idcub2(2).eq.3))
then
1502 > +( ff(3,1,1) -ff(3,1,2))
1503 > +(-ff(3,2,1) +ff(3,2,2)))
1510 subroutine hm2d(ans,ff)
1511 real,
intent(in) :: ff(0:3,2,2)
1512 real,
intent(out) :: ans
1518 real :: xpi,xp2,xpi2,ax,axbar,bx,bxbar
1519 real :: axp,axbarp,bxp,bxbarp
1521 real :: ypi,yp2,ypi2,ay,aybar,by,bybar
1522 real :: ayp,aybarp,byp,bybarp
1532 if((idcub2(1).eq.0).and.(idcub2(2).eq.0))
then
1545 ans=axbar*(aybar*ff(0,1,1) +ay*ff(0,1,2))+
1546 > ax*(aybar*ff(0,2,1)+ay*ff(0,2,2))
1549 > bxbar*(aybar*ff(1,1,1) +ay*ff(1,1,2))+
1550 > bx*(aybar*ff(1,2,1)+ay*ff(1,2,2)))
1553 > axbar*(bybar*ff(2,1,1) +by*ff(2,1,2))+
1554 > ax*(bybar*ff(2,2,1)+by*ff(2,2,2)))
1557 > bxbar*(bybar*ff(3,1,1) +by*ff(3,1,2))+
1558 > bx*(bybar*ff(3,2,1)+by*ff(3,2,2)))
1560 else if((idcub2(1).eq.1).and.(idcub2(2).eq.0))
then
1566 bxbarp=xpi*(3.0*xpi-2.0)
1574 > axbarp*(aybar*ff(0,1,1) +ay*ff(0,1,2))+
1575 > axp*(aybar*ff(0,2,1)+ay*ff(0,2,2)))
1578 > bxbarp*(aybar*ff(1,1,1) +ay*ff(1,1,2))+
1579 > bxp*(aybar*ff(1,2,1)+ay*ff(1,2,2))
1582 > axbarp*(bybar*ff(2,1,1) +by*ff(2,1,2))+
1583 > axp*(bybar*ff(2,2,1)+by*ff(2,2,2)))
1586 > bxbarp*(bybar*ff(3,1,1) +by*ff(3,1,2))+
1587 > bxp*(bybar*ff(3,2,1)+by*ff(3,2,2)))
1589 else if((idcub2(1).eq.0).and.(idcub2(2).eq.1))
then
1600 bybarp=ypi*(3.0*ypi-2.0)
1603 > axbar*(aybarp*ff(0,1,1) +ayp*ff(0,1,2))+
1604 > ax*(aybarp*ff(0,2,1)+ayp*ff(0,2,2)))
1607 > bxbar*(aybarp*ff(1,1,1) +ayp*ff(1,1,2))+
1608 > bx*(aybarp*ff(1,2,1)+ayp*ff(1,2,2)))
1611 > axbar*(bybarp*ff(2,1,1) +byp*ff(2,1,2))+
1612 > ax*(bybarp*ff(2,2,1)+byp*ff(2,2,2))
1615 > bxbar*(bybarp*ff(3,1,1) +byp*ff(3,1,2))+
1616 > bx*(bybarp*ff(3,2,1)+byp*ff(3,2,2)))
1618 else if((idcub2(1).eq.1).and.(idcub2(2).eq.1))
then
1624 bxbarp=xpi*(3.0*xpi-2.0)
1629 bybarp=ypi*(3.0*ypi-2.0)
1632 > axbarp*(aybarp*ff(0,1,1) +ayp*ff(0,1,2))+
1633 > axp*(aybarp*ff(0,2,1)+ayp*ff(0,2,2)))
1636 > bxbarp*(aybarp*ff(1,1,1) +ayp*ff(1,1,2))+
1637 > bxp*(aybarp*ff(1,2,1)+ayp*ff(1,2,2)))
1640 > axbarp*(bybarp*ff(2,1,1) +byp*ff(2,1,2))+
1641 > axp*(bybarp*ff(2,2,1)+byp*ff(2,2,2)))
1644 > bxbarp*(bybarp*ff(3,1,1) +byp*ff(3,1,2))+
1645 > bxp*(bybarp*ff(3,2,1)+byp*ff(3,2,2))