1 subroutine r8akherm3(x,nx,y,ny,z,nz,fherm,nf2,nf3,
2 > ilinx,iliny,ilinz,ier)
16 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
17 INTEGER iliny,ilinz,ier,ilinx
19 integer nx,ny,nz,nf2,nf3
23 real*8 fherm(0:7,nf2,nf3,nz)
49 call r8akherm3p(x,nx,y,ny,z,nz,fherm,nf2,nf3,
50 > ilinx,iliny,ilinz,0,0,0,ier)
55 subroutine r8akherm3p(x,nx,y,ny,z,nz,fherm,nf2,nf3,
56 > ilinx,iliny,ilinz,ipx,ipy,ipz,ier)
74 INTEGER,
PARAMETER :: R8=selected_real_kind(12,100)
75 INTEGER iliny,ilinz,ier,ilinx,ierx,iery,ierz,ierbc,iz,iy,ix
76 INTEGER icx,incx,ix1,icy,incy,iy1,icz,incz,iz1,izm2,izm1
77 INTEGER izmm2,izmm1,izp2,izp1,izpp2,izpp1,iym2,iym1,iymm2
78 INTEGER iymm1,iyp2,iyp1,iypp2,iypp1,ixm2,ixm1,ixmm2,ixmm1
79 INTEGER iflagx,ixp2,ixp1,ixpp2,ixpp1,iflagy,iflagz,k
82 real*8 cxp,cxm,cxpp,cxmm,cxtrap0,cxtrapn,cyp,cym,cypp,cymm
83 real*8 cytrap0,cytrapn,czp,czm,czpp,czmm,cztrap0,cztrapn
84 real*8 zansr,cxm2,cxp2,cym2,cyp2,fzp,fzm
86 integer nx,ny,nz,nf2,nf3
90 real*8 fherm(0:7,nf2,nf3,nz)
121 real*8 wx(2),wy(2),wz(2),exy(2,2),exz(2,2),eyz(2,2),fxyz(2,2,2)
123 real*8,
dimension(:,:,:),
allocatable :: ftmp
135 call r8splinck(x,nx,ilinx,1.0e-3_r8,ierx)
136 if(ierx.ne.0) ier=ier+1
139 write(6,
'('' ?akherm3: x axis not strict ascending'')')
142 call r8splinck(y,ny,iliny,1.0e-3_r8,iery)
143 if(iery.ne.0) ier=ier+1
146 write(6,
'('' ?akherm3: y axis not strict ascending'')')
149 call r8splinck(z,nz,ilinz,1.0e-3_r8,ierz)
150 if(ierz.ne.0) ier=ier+1
153 write(6,
'('' ?akherm3: z axis not strict ascending'')')
158 write(6,*)
'?akherm3: fherm (x) array dimension too small.'
163 write(6,*)
'?akherm3: fherm (y) array dimension too small.'
167 call ibc_ck(ipx,
'akherm3',
'X Bdy Cond',0,2,ierbc)
171 call ibc_ck(ipy,
'akherm3',
'Y Bdy Cond',0,2,ierbc)
175 call ibc_ck(ipz,
'akherm3',
'Z Bdy Cond',0,2,ierbc)
188 allocate(ftmp(0:nx+1,0:ny+1,0:nz+1))
193 ftmp(ix,iy,iz)=fherm(0,ix,iy,iz)
205 xx(nx+1)=2*x(nx)-x(nx-1)
207 yy(ny+1)=2*y(ny)-y(ny-1)
209 zz(nz+1)=2*z(nz)-z(nz-1)
226 cxp=(ftmp(2,iy,iz)-ftmp(1,iy,iz))/(xx(2)-xx(1))
227 cxm=(ftmp(nx,iy,iz)-ftmp(nx-1,iy,iz))/(xx(nx)-xx(nx-1))
234 cxpp=(ftmp(3,iy,iz)-ftmp(2,iy,iz))/(xx(3)-xx(2))
235 cxmm=(ftmp(nx-1,iy,iz)-ftmp(nx-2,iy,iz))/
236 > (xx(nx-1)-xx(nx-2))
238 call r8akherm0(cxmm,cxm,cxp,cxpp,wx,fherm(1,1,iy,iz))
239 fherm(1,nx,iy,iz)=fherm(1,1,iy,iz)
243 fherm(1,nx,iy,iz)=fherm(1,1,iy,iz)
249 else if(ipx.eq.0)
then
254 cxpp=(ftmp(3,iy,iz)-ftmp(2,iy,iz))/(xx(3)-xx(2))
255 fherm(1,1,iy,iz)=1.5_r8*cxp-0.5_r8*cxpp
257 cxmm=(ftmp(nx-1,iy,iz)-ftmp(nx-2,iy,iz))/
258 > (xx(nx-1)-xx(nx-2))
259 fherm(1,nx,iy,iz)=1.5_r8*cxm-0.5_r8*cxmm
263 fherm(1,nx,iy,iz)=fherm(1,1,iy,iz)
268 cxtrap0=2.0_r8*fherm(1,1,iy,iz)-cxp
269 cxtrapn=2.0_r8*fherm(1,nx,iy,iz)-cxm
275 cxtrap0=2.0_r8*fherm(1,1,iy,iz)-cxp
276 cxtrapn=2.0_r8*fherm(1,nx,iy,iz)-cxm
280 ftmp(0,iy,iz)=ftmp(1,iy,iz)-cxtrap0*(xx(1)-xx(0))
281 ftmp(nx+1,iy,iz)=ftmp(nx,iy,iz)+cxtrapn*(xx(nx+1)-xx(nx))
291 cyp=(ftmp(ix,2,iz)-ftmp(ix,1,iz))/(yy(2)-yy(1))
292 cym=(ftmp(ix,ny,iz)-ftmp(ix,ny-1,iz))/(yy(ny)-yy(ny-1))
299 cypp=(ftmp(ix,3,iz)-ftmp(ix,2,iz))/(yy(3)-yy(2))
300 cymm=(ftmp(ix,ny-1,iz)-ftmp(ix,ny-2,iz))/
301 > (yy(ny-1)-yy(ny-2))
303 call r8akherm0(cymm,cym,cyp,cypp,wy,fherm(2,ix,1,iz))
304 fherm(2,ix,ny,iz)=fherm(2,ix,1,iz)
307 fherm(2,ix,1,iz) = cyp
308 fherm(2,ix,ny,iz)=fherm(2,ix,1,iz)
313 else if(ipy.eq.0)
then
318 cypp=(ftmp(ix,3,iz)-ftmp(ix,2,iz))/(yy(3)-yy(2))
319 fherm(2,ix,1,iz)=1.5_r8*cyp-0.5_r8*cypp
321 cymm=(ftmp(ix,ny-1,iz)-ftmp(ix,ny-2,iz))/
322 > (yy(ny-1)-yy(ny-2))
323 fherm(2,ix,ny,iz)=1.5_r8*cym-0.5_r8*cymm
326 fherm(2,ix,1,iz) = cyp
327 fherm(2,ix,ny,iz)=fherm(2,ix,1,iz)
332 cytrap0=2.0_r8*fherm(2,ix,1,iz)-cyp
333 cytrapn=2.0_r8*fherm(2,ix,ny,iz)-cym
339 cytrap0=2.0_r8*fherm(2,ix,1,iz)-cyp
340 cytrapn=2.0_r8*fherm(2,ix,ny,iz)-cym
344 ftmp(ix,0,iz)=ftmp(ix,1,iz)-cytrap0*(yy(1)-yy(0))
345 ftmp(ix,ny+1,iz)=ftmp(ix,ny,iz)+cytrapn*(yy(ny+1)-yy(ny))
355 czp=(ftmp(ix,iy,2)-ftmp(ix,iy,1))/(zz(2)-zz(1))
356 czm=(ftmp(ix,iy,nz)-ftmp(ix,iy,nz-1))/(zz(nz)-zz(nz-1))
363 czpp=(ftmp(ix,iy,3)-ftmp(ix,iy,2))/(zz(3)-zz(2))
364 czmm=(ftmp(ix,iy,nz-1)-ftmp(ix,iy,nz-2))/
365 > (zz(nz-1)-zz(nz-2))
367 call r8akherm0(czmm,czm,czp,czpp,wz,fherm(3,ix,iy,1))
368 fherm(3,ix,iy,nz)=fherm(3,ix,iy,1)
371 fherm(3,ix,iy,1) = czp
372 fherm(3,ix,iy,nz)=fherm(3,ix,iy,1)
378 else if(ipy.eq.0)
then
383 czpp=(ftmp(ix,iy,3)-ftmp(ix,iy,2))/(zz(3)-zz(2))
384 fherm(3,ix,iy,1)=1.5_r8*czp-0.5_r8*czpp
386 czmm=(ftmp(ix,iy,nz-1)-ftmp(ix,iy,nz-2))/
387 > (zz(nz-1)-zz(nz-2))
388 fherm(3,ix,iy,nz)=1.5_r8*czm-0.5_r8*czmm
391 fherm(3,ix,iy,1) = czp
392 fherm(3,ix,iy,nz)=fherm(3,ix,iy,1)
397 cztrap0=2.0_r8*fherm(3,ix,iy,1)-czp
398 cztrapn=2.0_r8*fherm(3,ix,iy,nz)-czm
404 cztrap0=2.0_r8*fherm(3,ix,iy,1)-czp
405 cztrapn=2.0_r8*fherm(3,ix,iy,nz)-czm
409 ftmp(ix,iy,0)=ftmp(ix,iy,1)-cztrap0*(zz(1)-zz(0))
410 ftmp(ix,iy,nz+1)=ftmp(ix,iy,nz)+cztrapn*(zz(nz+1)-zz(nz))
427 ftmp(icx,icy,iz)=ftmp(icx,iy1,iz)+ftmp(ix1,icy,iz)-
444 ftmp(icx,iy,icz)=ftmp(icx,iy,iz1)+ftmp(ix1,iy,icz)-
461 ftmp(ix,icy,icz)=ftmp(ix,icy,iz1)+ftmp(ix,iy1,icz)-
484 > 2*(ftmp(icx,icy,iz1)+ftmp(icx,iy1,icz)+
485 > ftmp(ix1,icy,icz)) -
486 > (ftmp(icx,iy1,iz1)+ftmp(ix1,icy,iz1)+
487 > ftmp(ix1,iy1,icz)))/3
534 cxm=(ftmp(ixm2,iy,iz)-ftmp(ixm1,iy,iz))/
535 > (xx(ixm2)-xx(ixm1))
537 cxmm=(ftmp(ixmm2,iy,iz)-ftmp(ixmm1,iy,iz))/
538 > (xx(ixmm2)-xx(ixmm1))
541 cxmm=(ftmp(nx-1,iy,iz)-ftmp(nx-2,iy,iz))/
542 > (xx(nx-1)-xx(nx-2))
554 cxp=(ftmp(ixp2,iy,iz)-ftmp(ixp1,iy,iz))/
555 > (xx(ixp2)-xx(ixp1))
557 cxpp=(ftmp(ixpp2,iy,iz)-ftmp(ixpp1,iy,iz))/
558 > (xx(ixpp2)-xx(ixpp1))
561 cxpp=(ftmp(3,iy,iz)-ftmp(2,iy,iz))/
574 call r8akherm0(cxmm,cxm,cxp,cxpp,wx,zansr)
575 if((ix.gt.1).and.(ix.lt.nx)) fherm(1,ix,iy,iz)=zansr
580 cym=(ftmp(ix,iym2,iz)-ftmp(ix,iym1,iz))/
581 > (yy(iym2)-yy(iym1))
583 cymm=(ftmp(ix,iymm2,iz)-ftmp(ix,iymm1,iz))/
584 > (yy(iymm2)-yy(iymm1))
587 cymm=(ftmp(ix,ny-1,iz)-ftmp(ix,ny-2,iz))/
588 > (yy(ny-1)-yy(ny-2))
594 cyp=(ftmp(ix,iyp2,iz)-ftmp(ix,iyp1,iz))/
595 > (yy(iyp2)-yy(iyp1))
597 cypp=(ftmp(ix,iypp2,iz)-ftmp(ix,iypp1,iz))/
598 > (yy(iypp2)-yy(iypp1))
601 cypp=(ftmp(ix,3,iz)-ftmp(ix,2,iz))/(yy(3)-yy(2))
613 call r8akherm0(cymm,cym,cyp,cypp,wy,zansr)
614 if((iy.gt.1).and.(iy.lt.ny)) fherm(2,ix,iy,iz)=zansr
619 czm=(ftmp(ix,iy,izm2)-ftmp(ix,iy,izm1))/
620 > (zz(izm2)-zz(izm1))
622 czmm=(ftmp(ix,iy,izmm2)-ftmp(ix,iy,izmm1))/
623 > (zz(izmm2)-zz(izmm1))
626 czmm=(ftmp(ix,iy,nz-1)-ftmp(ix,iy,nz-2))/
627 > (zz(nz-1)-zz(nz-2))
633 czp=(ftmp(ix,iy,izp2)-ftmp(ix,iy,izp1))/
634 > (zz(izp2)-zz(izp1))
636 czpp=(ftmp(ix,iy,izpp2)-ftmp(ix,iy,izpp1))/
637 > (zz(izpp2)-zz(izpp1))
640 czpp=(ftmp(ix,iy,3)-ftmp(ix,iy,2))/(zz(3)-zz(2))
652 call r8akherm0(czmm,czm,czp,czpp,wz,zansr)
653 if((iz.gt.1).and.(iz.lt.nz)) fherm(3,ix,iy,iz)=zansr
657 cxm2=(ftmp(ixm2,iym1,iz)-ftmp(ixm1,iym1,iz))/
658 > (xx(ixm2)-xx(ixm1))
659 exy(1,1)=(cxm-cxm2)/(yy(iym2)-yy(iym1))
661 cxm2=(ftmp(ixm2,iyp2,iz)-ftmp(ixm1,iyp2,iz))/
662 > (xx(ixm2)-xx(ixm1))
663 exy(1,2)=(cxm2-cxm)/(yy(iyp2)-yy(iyp1))
665 cxp2=(ftmp(ixp2,iym1,iz)-ftmp(ixp1,iym1,iz))/
666 > (xx(ixp2)-xx(ixp1))
667 exy(2,1)=(cxp-cxp2)/(yy(iym2)-yy(iym1))
669 cxp2=(ftmp(ixp2,iyp2,iz)-ftmp(ixp1,iyp2,iz))/
670 > (xx(ixp2)-xx(ixp1))
671 exy(2,2)=(cxp2-cxp)/(yy(iyp2)-yy(iyp1))
675 cxm2=(ftmp(ixm2,iy,izm1)-ftmp(ixm1,iy,izm1))/
676 > (xx(ixm2)-xx(ixm1))
677 exz(1,1)=(cxm-cxm2)/(zz(izm2)-zz(izm1))
679 cxm2=(ftmp(ixm2,iy,izp2)-ftmp(ixm1,iy,izp2))/
680 > (xx(ixm2)-xx(ixm1))
681 exz(1,2)=(cxm2-cxm)/(zz(izp2)-zz(izp1))
683 cxp2=(ftmp(ixp2,iy,izm1)-ftmp(ixp1,iy,izm1))/
684 > (xx(ixp2)-xx(ixp1))
685 exz(2,1)=(cxp-cxp2)/(zz(izm2)-zz(izm1))
687 cxp2=(ftmp(ixp2,iy,izp2)-ftmp(ixp1,iy,izp2))/
688 > (xx(ixp2)-xx(ixp1))
689 exz(2,2)=(cxp2-cxp)/(zz(izp2)-zz(izp1))
693 cym2=(ftmp(ix,iym2,izm1)-ftmp(ix,iym1,izm1))/
694 > (yy(iym2)-yy(iym1))
695 eyz(1,1)=(cym-cym2)/(zz(izm2)-zz(izm1))
697 cym2=(ftmp(ix,iym2,izp2)-ftmp(ix,iym1,izp2))/
698 > (yy(iym2)-yy(iym1))
699 eyz(1,2)=(cym2-cym)/(zz(izp2)-zz(izp1))
701 cyp2=(ftmp(ix,iyp2,izm1)-ftmp(ix,iyp1,izm1))/
702 > (yy(iyp2)-yy(iyp1))
703 eyz(2,1)=(cyp-cyp2)/(zz(izm2)-zz(izm1))
705 cyp2=(ftmp(ix,iyp2,izp2)-ftmp(ix,iyp1,izp2))/
706 > (yy(iyp2)-yy(iyp1))
707 eyz(2,2)=(cyp2-cyp)/(zz(izp2)-zz(izp1))
713 fzp=(ftmp(ixm2,iym2,izm2)-ftmp(ixm2,iym1,izm2)-
714 > ftmp(ixm1,iym2,izm2)+ftmp(ixm1,iym1,izm2))/
715 > ((xx(ixm2)-xx(ixm1))*(yy(iym2)-yy(iym1)))
716 fzm=(ftmp(ixm2,iym2,izm1)-ftmp(ixm2,iym1,izm1)-
717 > ftmp(ixm1,iym2,izm1)+ftmp(ixm1,iym1,izm1))/
718 > ((xx(ixm2)-xx(ixm1))*(yy(iym2)-yy(iym1)))
719 fxyz(1,1,k)=(fzp-fzm)/(zz(izm2)-zz(izm1))
721 fzp=(ftmp(ixm2,iyp2,izm2)-ftmp(ixm2,iyp1,izm2)-
722 > ftmp(ixm1,iyp2,izm2)+ftmp(ixm1,iyp1,izm2))/
723 > ((xx(ixm2)-xx(ixm1))*(yy(iyp2)-yy(iyp1)))
724 fzm=(ftmp(ixm2,iyp2,izm1)-ftmp(ixm2,iyp1,izm1)-
725 > ftmp(ixm1,iyp2,izm1)+ftmp(ixm1,iyp1,izm1))/
726 > ((xx(ixm2)-xx(ixm1))*(yy(iyp2)-yy(iyp1)))
727 fxyz(1,2,k)=(fzp-fzm)/(zz(izm2)-zz(izm1))
729 fzp=(ftmp(ixp2,iym2,izm2)-ftmp(ixp2,iym1,izm2)-
730 > ftmp(ixp1,iym2,izm2)+ftmp(ixp1,iym1,izm2))/
731 > ((xx(ixp2)-xx(ixp1))*(yy(iym2)-yy(iym1)))
732 fzm=(ftmp(ixp2,iym2,izm1)-ftmp(ixp2,iym1,izm1)-
733 > ftmp(ixp1,iym2,izm1)+ftmp(ixp1,iym1,izm1))/
734 > ((xx(ixp2)-xx(ixp1))*(yy(iym2)-yy(iym1)))
735 fxyz(2,1,k)=(fzp-fzm)/(zz(izm2)-zz(izm1))
737 fzp=(ftmp(ixp2,iyp2,izm2)-ftmp(ixp2,iyp1,izm2)-
738 > ftmp(ixp1,iyp2,izm2)+ftmp(ixp1,iyp1,izm2))/
739 > ((xx(ixp2)-xx(ixp1))*(yy(iyp2)-yy(iyp1)))
740 fzm=(ftmp(ixp2,iyp2,izm1)-ftmp(ixp2,iyp1,izm1)-
741 > ftmp(ixp1,iyp2,izm1)+ftmp(ixp1,iyp1,izm1))/
742 > ((xx(ixp2)-xx(ixp1))*(yy(iyp2)-yy(iyp1)))
743 fxyz(2,2,k)=(fzp-fzm)/(zz(izm2)-zz(izm1))
747 fzp=(ftmp(ixm2,iym2,izp2)-ftmp(ixm2,iym1,izp2)-
748 > ftmp(ixm1,iym2,izp2)+ftmp(ixm1,iym1,izp2))/
749 > ((xx(ixm2)-xx(ixm1))*(yy(iym2)-yy(iym1)))
750 fzm=(ftmp(ixm2,iym2,izp1)-ftmp(ixm2,iym1,izp1)-
751 > ftmp(ixm1,iym2,izp1)+ftmp(ixm1,iym1,izp1))/
752 > ((xx(ixm2)-xx(ixm1))*(yy(iym2)-yy(iym1)))
753 fxyz(1,1,k)=(fzp-fzm)/(zz(izp2)-zz(izp1))
755 fzp=(ftmp(ixm2,iyp2,izp2)-ftmp(ixm2,iyp1,izp2)-
756 > ftmp(ixm1,iyp2,izp2)+ftmp(ixm1,iyp1,izp2))/
757 > ((xx(ixm2)-xx(ixm1))*(yy(iyp2)-yy(iyp1)))
758 fzm=(ftmp(ixm2,iyp2,izp1)-ftmp(ixm2,iyp1,izp1)-
759 > ftmp(ixm1,iyp2,izp1)+ftmp(ixm1,iyp1,izp1))/
760 > ((xx(ixm2)-xx(ixm1))*(yy(iyp2)-yy(iyp1)))
761 fxyz(1,2,k)=(fzp-fzm)/(zz(izp2)-zz(izp1))
763 fzp=(ftmp(ixp2,iym2,izp2)-ftmp(ixp2,iym1,izp2)-
764 > ftmp(ixp1,iym2,izp2)+ftmp(ixp1,iym1,izp2))/
765 > ((xx(ixp2)-xx(ixp1))*(yy(iym2)-yy(iym1)))
766 fzm=(ftmp(ixp2,iym2,izp1)-ftmp(ixp2,iym1,izp1)-
767 > ftmp(ixp1,iym2,izp1)+ftmp(ixp1,iym1,izp1))/
768 > ((xx(ixp2)-xx(ixp1))*(yy(iym2)-yy(iym1)))
769 fxyz(2,1,k)=(fzp-fzm)/(zz(izp2)-zz(izp1))
771 fzp=(ftmp(ixp2,iyp2,izp2)-ftmp(ixp2,iyp1,izp2)-
772 > ftmp(ixp1,iyp2,izp2)+ftmp(ixp1,iyp1,izp2))/
773 > ((xx(ixp2)-xx(ixp1))*(yy(iyp2)-yy(iyp1)))
774 fzm=(ftmp(ixp2,iyp2,izp1)-ftmp(ixp2,iyp1,izp1)-
775 > ftmp(ixp1,iyp2,izp1)+ftmp(ixp1,iyp1,izp1))/
776 > ((xx(ixp2)-xx(ixp1))*(yy(iyp2)-yy(iyp1)))
777 fxyz(2,2,k)=(fzp-fzm)/(zz(izp2)-zz(izp1))
782 > (wx(1)*(wy(1)*exy(1,1)+wy(2)*exy(1,2))+
783 > wx(2)*(wy(1)*exy(2,1)+wy(2)*exy(2,2)))/
784 > ((wx(1)+wx(2))*(wy(1)+wy(2)))
787 > (wx(1)*(wz(1)*exz(1,1)+wz(2)*exz(1,2))+
788 > wx(2)*(wz(1)*exz(2,1)+wz(2)*exz(2,2)))/
789 > ((wx(1)+wx(2))*(wz(1)+wz(2)))
792 > (wy(1)*(wz(1)*eyz(1,1)+wz(2)*eyz(1,2))+
793 > wy(2)*(wz(1)*eyz(2,1)+wz(2)*eyz(2,2)))/
794 > ((wy(1)+wy(2))*(wz(1)+wz(2)))
796 fherm(7,ix,iy,iz)=(wz(1)*
797 > (wx(1)*(wy(1)*fxyz(1,1,1)+wy(2)*fxyz(1,2,1))+
798 > wx(2)*(wy(1)*fxyz(2,1,1)+wy(2)*fxyz(2,2,1))) + wz(2)*
799 > (wx(1)*(wy(1)*fxyz(1,1,2)+wy(2)*fxyz(1,2,2))+
800 > wx(2)*(wy(1)*fxyz(2,1,2)+wy(2)*fxyz(2,2,2))) )/
801 > ((wx(1)+wx(2))*(wy(1)+wy(2))*(wz(1)+wz(2)))