c pfdtd.f c FDTD and Total Tunning by Jianqing Wang, Nagoya Institute of Technology c 2002.2.21 c MPI by Takeshi Shinoda, Nagoya Institute of Technology c 2001.10.30 c c Assignment of data & computation c (numbers means processor ID) c c +---+---+ c / 6 / 7 /| c +---+---+7+ c / 4 / 5 /|/| c +---+---+5+3+ c | 4 | 5 |/|/ c +---+---+1+ c | 0 | 1 |/ c +---+---+ c c k j c | / c |/ c O--i c c Send & recv connection between pc c c 6----7 c /| /| c 4----5 | c | | | | c | 2--|-3 c |/ |/ c 0----1 c c program pfdtd include 'mpif.h' parameter(II=20,JJ=20,KK=20,NSTOP=200) parameter(DX=0.002,DY=0.002,DZ=0.002) parameter(IA=15,JA=15,KGAP=15,KA1=KGAP-1,KA2=KGAP+1) parameter(TAU0=1.0/2.0E9) parameter(C=3.0E8,PI=3.1415926) parameter(EPS0=8.854E-12,XMU0=1.2566306E-6) parameter(LM=4,MM=2,RC0=-60.0, & NPML=(LM+1)*((JJ+1)*(KK+1)+(II+1)*(KK+1)+(II+1)*(JJ+1))) integer n,i,j,k integer istart(0:7),iexstop(0:7),ieystop(0:7),iezstop(0:7), & ihxstop(0:7),ihystop(0:7),ihzstop(0:7), & jstart(0:7),jexstop(0:7),jeystop(0:7),jezstop(0:7), & jhxstop(0:7),jhystop(0:7),jhzstop(0:7), & kstart(0:7),kexstop(0:7),keystop(0:7),kezstop(0:7), & khxstop(0:7),khystop(0:7),khzstop(0:7) double precision t,dt double precision & E(II+1,JJ+1,KK+1,3),H(0:II,0:JJ,0:KK,3), & Cex(II+1,JJ+1,KK+1),Cey(II+1,JJ+1,KK+1),Cez(II+1,JJ+1,KK+1), & Cexl(II+1,JJ+1,KK+1),Ceyl(II+1,JJ+1,KK+1), & Cezl(II+1,JJ+1,KK+1),Chx,Chy,Chz,Vgap(NSTOP),Cgap(NSTOP) integer istepmlx(3,0:7),ispepmlx(3,0:7),jstepmlx(3,0:7), & jspepmlx(3,0:7),kstepmlx(3,0:7),kspepmlx(3,0:7), & istepmly(3,0:7),ispepmly(3,0:7),jstepmly(3,0:7), & jspepmly(3,0:7),kstepmly(3,0:7),kspepmly(3,0:7), & istepmlz(3,0:7),ispepmlz(3,0:7),jstepmlz(3,0:7), & jspepmlz(3,0:7),kstepmlz(3,0:7),kspepmlz(3,0:7), & isthpmlx(3,0:7),isphpmlx(3,0:7),jsthpmlx(3,0:7), & jsphpmlx(3,0:7),ksthpmlx(3,0:7),ksphpmlx(3,0:7), & isthpmly(3,0:7),isphpmly(3,0:7),jsthpmly(3,0:7), & jsphpmly(3,0:7),ksthpmly(3,0:7),ksphpmly(3,0:7), & isthpmlz(3,0:7),isphpmlz(3,0:7),jsthpmlz(3,0:7), & jsphpmlz(3,0:7),ksthpmlz(3,0:7),ksphpmlz(3,0:7), & lpmlst(3) double precision & Exy(NPML),Exz(NPML),Eyx(NPML),Eyz(NPML),Ezx(NPML),Ezy(NPML), & Hxy(NPML),Hxz(NPML),Hyx(NPML),Hyz(NPML),Hzx(NPML),Hzy(NPML), & Cedx(II+1),Cedy(JJ+1),Cedz(KK+1),Cedxl(II+1),Cedyl(JJ+1), & Cedzl(KK+1),Chdx(II+1),Chdy(JJ+1),Chdz(KK+1),Chdxl(II+1), & Chdyl(JJ+1),Chdzl(KK+1) c 出力ファイル名 character fname1*20 character fname2*20 c MPI 送受信のための変数 c プロセッサID integer myrank c 送受信データタグ，エラーコード，通信状態 等 integer tag(3),rquest(3),ierr,status(MPI_STATUS_SIZE) c プロセッサ名，プロセッサ名の長さ c character pcname*MPI_MAX_PROCESSOR_NAME c integer nameln c 送受信データのパッキングのためのデータ型 integer jkplan integer kiplan integer ijplan c 時刻 double precision start c print*,'program start' c MPI 送受信のための準備 c MPI 初期化 call MPI_INIT(ierr) c 自プロセッサIDの取得 call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) c print*,myrank,char(myrank) c 自プロセッサ名の取得 c call MPI_GET_PROCESSOR_NAME(pcname,nameln,ierr) c 送受信データのパッキングのための型の設定 call MPI_TYPE_VECTOR( & 3*(KK+1)*(JJ+1),1,(II+1), & MPI_DOUBLE_PRECISION,jkplan,ierr) call MPI_TYPE_VECTOR( & 3*(KK+1),(II+1)*1,(II+1)*(JJ+1), & MPI_DOUBLE_PRECISION,kiplan,ierr) call MPI_TYPE_VECTOR( & 3,(JJ+1)*(II+1)*1,(II+1)*(JJ+1)*(KK+1), & MPI_DOUBLE_PRECISION,ijplan,ierr) call MPI_TYPE_COMMIT(jkplan,ierr) call MPI_TYPE_COMMIT(kiplan,ierr) call MPI_TYPE_COMMIT(ijplan,ierr) c if (myrank.eq.0) print*,'mpi vector ok' c メッセージタグの設定 tag(1)=1 tag(2)=2 tag(3)=3 c E,H 初期化 call init(E,H,Exy,Exz,Eyx,Eyz,Ezx,Ezy, & Hxy,Hxz,Hyx,Hyz,Hzx,Hzy,II,JJ,KK,NPML,myrank) c 設定 call setup(Cex,Cey,Cez,Cexl,Ceyl,Cezl,Chx,Chy,Chz, & dt,II,JJ,KK,LM,DX,DY,DZ,C,EPS0,XMU0, & istart,iexstop,ieystop,iezstop,ihxstop,ihystop,ihzstop, & jstart,jexstop,jeystop,jezstop,jhxstop,jhystop,jhzstop, & kstart,kexstop,keystop,kezstop,khxstop,khystop,khzstop) call setup_pml(Cedx,Cedy,Cedz,Cedxl,Cedyl,Cedzl, & Chdx,Chdy,Chdz,Chdxl,Chdyl,Chdzl,II,JJ,KK,LM,MM,RC0, & EPS0,XMU0,C,DX,DY,DZ,dt, & istepmlx,ispepmlx,jstepmlx,jspepmlx,kstepmlx,kspepmlx, & istepmly,ispepmly,jstepmly,jspepmly,kstepmly,kspepmly, & istepmlz,ispepmlz,jstepmlz,jspepmlz,kstepmlz,kspepmlz, & isthpmlx,isphpmlx,jsthpmlx,jsphpmlx,ksthpmlx,ksphpmlx, & isthpmly,isphpmly,jsthpmly,jsphpmly,ksthpmly,ksphpmly, & isthpmlz,isphpmlz,jsthpmlz,jsphpmlz,ksthpmlz,ksphpmlz, & lpmlst,myrank) c Open output file c for own processor fname1='01234567' fname1='pvi'//fname1(myrank+1:myrank+1)//'.dat' if (myrank.eq.0) then open (02,FILE=fname1) end if fname2='01234567' fname2='peh'//fname2(myrank+1:myrank+1)//'.dat' if (myrank.eq.7) then open (03,FILE=fname2) end if t=0.0 call MPI_Barrier(MPI_COMM_WORLD,ierr) if (myrank.eq.0) then start=MPI_Wtime() print*,'iteration start' c & ,start,MPI_Wtime(),MPI_Wtick() end if do n=1,NSTOP c E-field update call efld(E,H,Cex,Cey,Cez,Cexl,Ceyl,Cezl,II,JJ,KK,DX,DY,DZ, & istart,iexstop,ieystop,iezstop,jstart,jexstop,jeystop, & jezstop,kstart,kexstop,keystop,kezstop,myrank) c Excitation source if (myrank.eq.0) then do k=KA1,KA2 E(IA,JA,k,3)=0.0 end do if (t.le.(2.0*TAU0)) then Vgap(n)=exp(-(4.0/TAU0)**2*(t-TAU0)**2) else Vgap(n)=0.0 end if E(IA,JA,KGAP,3)=-Vgap(n)/DZ end if c Apply e-field pml absorbing boundary condition call epml(E,H,Exy,Exz,Eyx,Eyz,Ezx,Ezy, & Cedx,Cedy,Cedz,Cedxl,Cedyl,Cedzl,II,JJ,KK,NPML, & istepmlx,ispepmlx,jstepmlx,jspepmlx,kstepmlx,kspepmlx, & istepmly,ispepmly,jstepmly,jspepmly,kstepmly,kspepmly, & istepmlz,ispepmlz,jstepmlz,jspepmlz,kstepmlz,kspepmlz, & lpmlst,myrank) c MPI送受信 E 開始 if (mod(myrank,2).eq.1) then c 1,3,5,7 は 0,2,4,6 に 送信 call MPI_SEND(E( 1,1,1,1),1,jkplan,myrank-1, & tag(1),MPI_COMM_WORLD,ierr) else c 0,2,4,6 は 1,3,5,7 から 受信 call MPI_RECV(E(II+1,1,1,1),1,jkplan,myrank+1, & tag(1),MPI_COMM_WORLD,status,ierr) endif if (mod((myrank/2),2).eq.1) then c 2,3,6,7 は 0,1,4,5 に 送信 call MPI_SEND(E(1, 1,1,1),1,kiplan,myrank-2, & tag(2),MPI_COMM_WORLD,ierr) else c 0,1,4,5 は 2,3,6,7 から 受信 call MPI_RECV(E(1,JJ+1,1,1),1,kiplan,myrank+2, & tag(2),MPI_COMM_WORLD,status,ierr) endif if ((myrank/4).eq.1) then c 4,5,6,7 は 0,1,2,3 へ 送信 call MPI_SEND(E(1,1, 1,1),1,ijplan,myrank-4,tag(3), & MPI_COMM_WORLD,ierr) else c 0,1,2,3 は 4,5,6,7 から 受信 call MPI_RECV(E(1,1,KK+1,1),1,ijplan,myrank+4,tag(3), & MPI_COMM_WORLD,status,ierr) endif c MPI送受信 E 終了 t=t+dt/2.0 c H-field update call hfld(E,H,Chx,Chy,Chz,II,JJ,KK, & istart,ihxstop,ihystop,ihzstop,jstart,jhxstop,jhystop, & jhzstop,kstart,khxstop,khystop,khzstop,myrank) c Apply h-field pml absorbing boundary condition call hpml(E,H,Hxy,Hxz,Hyx,Hyz,Hzx,Hzy, & Chdx,Chdy,Chdz,Chdxl,Chdyl,Chdzl,II,JJ,KK,NPML, & isthpmlx,isphpmlx,jsthpmlx,jsphpmlx,ksthpmlx,ksphpmlx, & isthpmly,isphpmly,jsthpmly,jsphpmly,ksthpmly,ksphpmly, & isthpmlz,isphpmlz,jsthpmlz,jsphpmlz,ksthpmlz,ksphpmlz, & lpmlst,myrank) c MPI 送受信 H 開始 if (mod(myrank,2).eq.0) then c 0,2,4,6 は 1,3,5,7 へ送信 call MPI_SEND(H(II,0,0,1),1,jkplan,myrank+1, & tag(1),MPI_COMM_WORLD,ierr) else c 1,3,5,7 は 0,2,4,6 から 受信 call MPI_RECV(H( 0,0,0,1),1,jkplan,myrank-1, & tag(1),MPI_COMM_WORLD,status,ierr) endif if (mod((myrank/2),2).eq.0) then c 0,1,4,5 は 2,3,6,7 へ 送信 call MPI_SEND(H(0,JJ,0,1),1,kiplan,myrank+2, & tag(2),MPI_COMM_WORLD,ierr) else c 2,3,6,7 は 0,1,4,5 から 受信 call MPI_RECV(H(0, 0,0,1),1,kiplan,myrank-2, & tag(2),MPI_COMM_WORLD,status,ierr) endif if ((myrank/4).eq.0) then c 0,1,2,3 は 4,5,6,7 へ 送信 call MPI_SEND(H(0,0,KK,1),1,ijplan,myrank+4,tag(3), & MPI_COMM_WORLD,ierr) else c 0,1,2,3 は 4,5,6,7 から 受信 call MPI_RECV(H(0,0, 0,1),1,ijplan,myrank-4,tag(3), & MPI_COMM_WORLD,status,ierr) endif c MPI 送受信 H 終了 t=t+dt/2.0 c Save temporary voltage and current at the gap if (myrank.eq.0) then call savvi(H,Vgap,Cgap,n,NSTOP,II,JJ,KK,IA,JA,KGAP,DX,DY) end if c Save temporary E and H fields if (myrank.eq.7) then write(03,*) E(IA-10,JA-10,KGAP-10,3), & H(IA-10,JA-10,KGAP-10,1) end if end do c Close output file close (02) close (03) call MPI_Barrier(MPI_COMM_WORLD,ierr) if (myrank.eq.0) then print*,'iteration finished',MPI_Wtime()-start c & ,(MPI_Wtime()-start),MPI_Wtick() end if c MPI 終了 call MPI_FINALIZE(ierr) stop end c E,H 初期化 subroutine init(E,H,Exy,Exz,Eyx,Eyz,Ezx,Ezy, & Hxy,Hxz,Hyx,Hyz,Hzx,Hzy,II,JJ,KK,NPML,myrank) double precision & E(II+1,JJ+1,KK+1,3),H(0:II,0:JJ,0:KK,3), & Exy(NPML),Exz(NPML),Eyx(NPML),Eyz(NPML),Ezx(NPML),Ezy(NPML), & Hxy(NPML),Hxz(NPML),Hyx(NPML),Hyz(NPML),Hzx(NPML),Hzy(NPML) integer i,j,k,l character fname*20 fname='01234567' fname='const'//fname(myrank+1:myrank+1)//'.dat' C Read material properties c OPEN(01,FILE=fname) c DO I=0,MKIND+2 c READ(01,*) EP(I),SG(I),RHO(I) c EP(I)=EP(I)*EPS0 c END DO c CLOSE(01) c E 初期化 do k=1,KK+1 do j=1,JJ+1 do i=1,II+1 E(i,j,k,1)=0.0 E(i,j,k,2)=0.0 E(i,j,k,3)=0.0 end do end do end do c H 初期化 do k=0,KK do j=0,JJ do i=0,II H(i,j,k,1)=0.0 H(i,j,k,2)=0.0 H(i,j,k,3)=0.0 end do end do end do c PML 初期化 do l=1,NPML Exy(l)=0.0 Exz(l)=0.0 Eyx(l)=0.0 Eyz(l)=0.0 Ezx(l)=0.0 Ezy(l)=0.0 Hxy(l)=0.0 Hxz(l)=0.0 Hyx(l)=0.0 Hyz(l)=0.0 Hzx(l)=0.0 Hzy(l)=0.0 end do return end c パラメータ設定 subroutine setup(Cex,Cey,Cez,Cexl,Ceyl,Cezl,Chx,Chy,Chz, & dt,II,JJ,KK,LM,DX,DY,DZ,C,EPS0,XMU0, & istart,iexstop,ieystop,iezstop,ihxstop,ihystop,ihzstop, & jstart,jexstop,jeystop,jezstop,jhxstop,jhystop,jhzstop, & kstart,kexstop,keystop,kezstop,khxstop,khystop,khzstop) integer i,j,k integer istart(0:7),iexstop(0:7),ieystop(0:7),iezstop(0:7), & ihxstop(0:7),ihystop(0:7),ihzstop(0:7), & jstart(0:7),jexstop(0:7),jeystop(0:7),jezstop(0:7), & jhxstop(0:7),jhystop(0:7),jhzstop(0:7), & kstart(0:7),kexstop(0:7),keystop(0:7),kezstop(0:7), & khxstop(0:7),khystop(0:7),khzstop(0:7) double precision & Cex(II+1,JJ+1,KK+1),Cey(II+1,JJ+1,KK+1),Cez(II+1,JJ+1,KK+1), & Cexl(II+1,JJ+1,KK+1),Ceyl(II+1,JJ+1,KK+1), & Cezl(II+1,JJ+1,KK+1),Chx,Chy,Chz double precision dt dt=0.99/C/sqrt((1.0/DX)**2+(1.0/DY)**2+(1.0/DZ)**2) do k=1,KK+1 do j=1,JJ+1 do i=1,II+1 Cex(i,j,k)=1.0 Cey(i,j,k)=1.0 Cez(i,j,k)=1.0 Cexl(i,j,k)=dt/EPS0 Ceyl(i,j,k)=dt/EPS0 Cezl(i,j,k)=dt/EPS0 end do end do end do Chx=dt/(XMU0*DX) Chy=dt/(XMU0*DY) Chz=dt/(XMU0*DZ) istart(0)=LM+1 iexstop(0)=II ieystop(0)=II iezstop(0)=II ihxstop(0)=II ihystop(0)=II ihzstop(0)=II jstart(0)=LM+1 jexstop(0)=JJ jeystop(0)=JJ jezstop(0)=JJ jhxstop(0)=JJ jhystop(0)=JJ jhzstop(0)=JJ kstart(0)=LM+1 kexstop(0)=KK keystop(0)=KK kezstop(0)=KK khxstop(0)=KK khystop(0)=KK khzstop(0)=KK istart(1)=1 iexstop(1)=II-LM ieystop(1)=II-LM+1 iezstop(1)=II-LM+1 ihxstop(1)=II-LM+1 ihystop(1)=II-LM ihzstop(1)=II-LM jstart(1)=LM+1 jexstop(1)=JJ jeystop(1)=JJ jezstop(1)=JJ jhxstop(1)=JJ jhystop(1)=JJ jhzstop(1)=JJ kstart(1)=LM+1 kexstop(1)=KK keystop(1)=KK kezstop(1)=KK khxstop(1)=KK khystop(1)=KK khzstop(1)=KK istart(2)=LM+1 iexstop(2)=II ieystop(2)=II iezstop(2)=II ihxstop(2)=II ihystop(2)=II ihzstop(2)=II jstart(2)=1 jexstop(2)=JJ-LM+1 jeystop(2)=JJ-LM jezstop(2)=JJ-LM+1 jhxstop(2)=JJ-LM jhystop(2)=JJ-LM+1 jhzstop(2)=JJ-LM kstart(2)=LM+1 kexstop(2)=KK keystop(2)=KK kezstop(2)=KK khxstop(2)=KK khystop(2)=KK khzstop(2)=KK istart(3)=1 iexstop(3)=II-LM ieystop(3)=II-LM+1 iezstop(3)=II-LM+1 ihxstop(3)=II-LM+1 ihystop(3)=II-LM ihzstop(3)=II-LM jstart(3)=1 jexstop(3)=JJ-LM+1 jeystop(3)=JJ-LM jezstop(3)=JJ-LM+1 jhxstop(3)=JJ-LM jhystop(3)=JJ-LM+1 jhzstop(3)=JJ-LM kstart(3)=LM+1 kexstop(3)=KK keystop(3)=KK kezstop(3)=KK khxstop(3)=KK khystop(3)=KK khzstop(3)=KK istart(4)=LM+1 iexstop(4)=II ieystop(4)=II iezstop(4)=II ihxstop(4)=II ihystop(4)=II ihzstop(4)=II jstart(4)=LM+1 jexstop(4)=JJ jeystop(4)=JJ jezstop(4)=JJ jhxstop(4)=JJ jhystop(4)=JJ jhzstop(4)=JJ kstart(4)=1 kexstop(4)=KK-LM+1 keystop(4)=KK-LM+1 kezstop(4)=KK-LM khxstop(4)=KK-LM khystop(4)=KK-LM khzstop(4)=KK-LM+1 istart(5)=1 iexstop(5)=II-LM ieystop(5)=II-LM+1 iezstop(5)=II-LM+1 ihxstop(5)=II-LM+1 ihystop(5)=II-LM ihzstop(5)=II-LM jstart(5)=LM+1 jexstop(5)=JJ jeystop(5)=JJ jezstop(5)=JJ jhxstop(5)=JJ jhystop(5)=JJ jhzstop(5)=JJ kstart(5)=1 kexstop(5)=KK-LM+1 keystop(5)=KK-LM+1 kezstop(5)=KK-LM khxstop(5)=KK-LM khystop(5)=KK-LM khzstop(5)=KK-LM+1 istart(6)=LM+1 iexstop(6)=II ieystop(6)=II iezstop(6)=II ihxstop(6)=II ihystop(6)=II ihzstop(6)=II jstart(6)=1 jexstop(6)=JJ-LM+1 jeystop(6)=JJ-LM jezstop(6)=JJ-LM+1 jhxstop(6)=JJ-LM jhystop(6)=JJ-LM+1 jhzstop(6)=JJ-LM kstart(6)=1 kexstop(6)=KK-LM+1 keystop(6)=KK-LM+1 kezstop(6)=KK-LM khxstop(6)=KK-LM khystop(6)=KK-LM khzstop(6)=KK-LM+1 istart(7)=1 iexstop(7)=II-LM ieystop(7)=II-LM+1 iezstop(7)=II-LM+1 ihxstop(7)=II-LM+1 ihystop(7)=II-LM ihzstop(7)=II-LM jstart(7)=1 jexstop(7)=JJ-LM+1 jeystop(7)=JJ-LM jezstop(7)=JJ-LM+1 jhxstop(7)=JJ-LM jhystop(7)=JJ-LM+1 jhzstop(7)=JJ-LM kstart(7)=1 kexstop(7)=KK-LM+1 keystop(7)=KK-LM+1 kezstop(7)=KK-LM khxstop(7)=KK-LM khystop(7)=KK-LM khzstop(7)=KK-LM+1 return end subroutine setup_pml(Cedx,Cedy,Cedz,Cedxl,Cedyl,Cedzl, & Chdx,Chdy,Chdz,Chdxl,Chdyl,Chdzl,II,JJ,KK,LM,MM,RC0, & EPS0,XMU0,C,DX,DY,DZ,dt, & istepmlx,ispepmlx,jstepmlx,jspepmlx,kstepmlx,kspepmlx, & istepmly,ispepmly,jstepmly,jspepmly,kstepmly,kspepmly, & istepmlz,ispepmlz,jstepmlz,jspepmlz,kstepmlz,kspepmlz, & isthpmlx,isphpmlx,jsthpmlx,jsphpmlx,ksthpmlx,ksphpmlx, & isthpmly,isphpmly,jsthpmly,jsphpmly,ksthpmly,ksphpmly, & isthpmlz,isphpmlz,jsthpmlz,jsphpmlz,ksthpmlz,ksphpmlz, & lpmlst,myrank) integer l integer istepmlx(3,0:7),ispepmlx(3,0:7),jstepmlx(3,0:7), & jspepmlx(3,0:7),kstepmlx(3,0:7),kspepmlx(3,0:7), & istepmly(3,0:7),ispepmly(3,0:7),jstepmly(3,0:7), & jspepmly(3,0:7),kstepmly(3,0:7),kspepmly(3,0:7), & istepmlz(3,0:7),ispepmlz(3,0:7),jstepmlz(3,0:7), & jspepmlz(3,0:7),kstepmlz(3,0:7),kspepmlz(3,0:7), & isthpmlx(3,0:7),isphpmlx(3,0:7),jsthpmlx(3,0:7), & jsphpmlx(3,0:7),ksthpmlx(3,0:7),ksphpmlx(3,0:7), & isthpmly(3,0:7),isphpmly(3,0:7),jsthpmly(3,0:7), & jsphpmly(3,0:7),ksthpmly(3,0:7),ksphpmly(3,0:7), & isthpmlz(3,0:7),isphpmlz(3,0:7),jsthpmlz(3,0:7), & jsphpmlz(3,0:7),ksthpmlz(3,0:7),ksphpmlz(3,0:7), & lpmlst(3) double precision & Cedx(II+1),Cedy(JJ+1),Cedz(KK+1),Cedxl(II+1),Cedyl(JJ+1), & Cedzl(KK+1),Chdx(II+1),Chdy(JJ+1),Chdz(KK+1),Chdxl(II+1), & Chdyl(JJ+1),Chdzl(KK+1),dt,Sglmax,Sgle(LM),Sglh(LM) C GENERATE CONSTATNS FOR FIELD UPDATE EQUATIONS IN PML REGION Sglmax=-(MM+1)*EPS0*C/(2.0*LM*DX)*RC0 do l=1,LM Sgle(l)=Sglmax*((LM-l+1.0)/LM)**MM Sglh(l)=Sglmax*((LM-l+0.5)/LM)**MM end do do l=1,II+1 Cedx(l)=1.0 Chdx(l)=1.0 Cedxl(l)=dt/EPS0/DX Chdxl(l)=dt/XMU0/DX end do do l=1,JJ+1 Cedy(l)=1.0 Chdy(l)=1.0 Cedyl(l)=dt/EPS0/DY Chdyl(l)=dt/XMU0/DY end do do l=1,KK+1 Cedz(l)=1.0 Chdz(l)=1.0 Cedzl(l)=dt/EPS0/DZ Chdzl(l)=dt/XMU0/DZ end do if (myrank.eq.0) then do l=1,LM Cedx(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedy(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedz(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedxl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DX Cedyl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DY Cedzl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DZ Chdx(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdy(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdz(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdxl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DX Chdyl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DY Chdzl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DZ end do end if if (myrank.eq.1) then do l=1,LM Cedx(II-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedy(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedz(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedxl(II-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DX Cedyl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DY Cedzl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DZ Chdx(II-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdy(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdz(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdxl(II-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DX Chdyl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DY Chdzl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DZ end do end if if (myrank.eq.2) then do l=1,LM Cedx(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedy(JJ-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedz(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedxl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DX Cedyl(JJ-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DY Cedzl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DZ Chdx(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdy(JJ-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdz(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdxl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DX Chdyl(JJ-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DY Chdzl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DZ end do end if if (myrank.eq.3) then do l=1,LM Cedx(II-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedy(JJ-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedz(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedxl(II-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DX Cedyl(JJ-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DY Cedzl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DZ Chdx(II-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdy(JJ-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdz(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdxl(II-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DX Chdyl(JJ-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DY Chdzl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DZ end do end if if (myrank.eq.4) then do l=1,LM Cedx(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedy(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedz(KK-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedxl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DX Cedyl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DY Cedzl(KK-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DZ Chdx(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdy(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdz(KK-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdxl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DX Chdyl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DY Chdzl(KK-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DZ end do end if if (myrank.eq.5) then do l=1,LM Cedx(II-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedy(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedz(KK-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedxl(II-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DX Cedyl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DY Cedzl(KK-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DZ Chdx(II-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdy(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdz(KK-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdxl(II-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DX Chdyl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DY Chdzl(KK-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DZ end do end if if (myrank.eq.6) then do l=1,LM Cedx(l)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedy(JJ-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedz(KK-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedxl(l)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DX Cedyl(JJ-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DY Cedzl(KK-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DZ Chdx(l)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdy(JJ-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdz(KK-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdxl(l)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DX Chdyl(JJ-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DY Chdzl(KK-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DZ end do end if if (myrank.eq.7) then do l=1,LM Cedx(II-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedy(JJ-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedz(KK-l+2)=(1.0-Sgle(l)*dt/(2.0*EPS0)) & /(1.0+Sgle(l)*dt/(2.0*EPS0)) Cedxl(II-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DX Cedyl(JJ-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DY Cedzl(KK-l+2)=(dt/EPS0)/(1.0+Sgle(l)*dt/(2.0*EPS0))/DZ Chdx(II-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdy(JJ-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdz(KK-l+2)=(1.0-Sglh(l)*dt/(2.0*EPS0)) & /(1.0+Sglh(l)*dt/(2.0*EPS0)) Chdxl(II-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DX Chdyl(JJ-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DY Chdzl(KK-l+2)=(dt/XMU0)/(1.0+Sglh(l)*dt/(2.0*EPS0))/DZ end do end if c For E-PML c For processor 0 istepmlx(1,0)=1 ispepmlx(1,0)=LM jstepmlx(1,0)=2 jspepmlx(1,0)=JJ kstepmlx(1,0)=2 kspepmlx(1,0)=KK istepmly(1,0)=2 ispepmly(1,0)=LM jstepmly(1,0)=1 jspepmly(1,0)=JJ kstepmly(1,0)=2 kspepmly(1,0)=KK istepmlz(1,0)=2 ispepmlz(1,0)=LM jstepmlz(1,0)=2 jspepmlz(1,0)=JJ kstepmlz(1,0)=1 kspepmlz(1,0)=KK istepmlx(2,0)=1 ispepmlx(2,0)=II jstepmlx(2,0)=2 jspepmlx(2,0)=LM kstepmlx(2,0)=2 kspepmlx(2,0)=KK istepmly(2,0)=2 ispepmly(2,0)=II jstepmly(2,0)=1 jspepmly(2,0)=LM kstepmly(2,0)=2 kspepmly(2,0)=KK istepmlz(2,0)=2 ispepmlz(2,0)=II jstepmlz(2,0)=2 jspepmlz(2,0)=LM kstepmlz(2,0)=1 kspepmlz(2,0)=KK istepmlx(3,0)=1 ispepmlx(3,0)=II jstepmlx(3,0)=2 jspepmlx(3,0)=JJ kstepmlx(3,0)=2 kspepmlx(3,0)=LM istepmly(3,0)=2 ispepmly(3,0)=II jstepmly(3,0)=1 jspepmly(3,0)=JJ kstepmly(3,0)=2 kspepmly(3,0)=LM istepmlz(3,0)=2 ispepmlz(3,0)=II jstepmlz(3,0)=2 jspepmlz(3,0)=JJ kstepmlz(3,0)=1 kspepmlz(3,0)=LM c For processor 1 istepmlx(1,1)=II-LM+1 ispepmlx(1,1)=II jstepmlx(1,1)=2 jspepmlx(1,1)=JJ kstepmlx(1,1)=2 kspepmlx(1,1)=KK istepmly(1,1)=II-LM+2 ispepmly(1,1)=II jstepmly(1,1)=1 jspepmly(1,1)=JJ kstepmly(1,1)=2 kspepmly(1,1)=KK istepmlz(1,1)=II-LM+2 ispepmlz(1,1)=II jstepmlz(1,1)=2 jspepmlz(1,1)=JJ kstepmlz(1,1)=1 kspepmlz(1,1)=KK istepmlx(2,1)=1 ispepmlx(2,1)=II jstepmlx(2,1)=2 jspepmlx(2,1)=LM kstepmlx(2,1)=2 kspepmlx(2,1)=KK istepmly(2,1)=1 ispepmly(2,1)=II jstepmly(2,1)=1 jspepmly(2,1)=LM kstepmly(2,1)=2 kspepmly(2,1)=KK istepmlz(2,1)=1 ispepmlz(2,1)=II jstepmlz(2,1)=2 jspepmlz(2,1)=LM kstepmlz(2,1)=1 kspepmlz(2,1)=KK istepmlx(3,1)=1 ispepmlx(3,1)=II jstepmlx(3,1)=2 jspepmlx(3,1)=JJ kstepmlx(3,1)=2 kspepmlx(3,1)=LM istepmly(3,1)=1 ispepmly(3,1)=II jstepmly(3,1)=1 jspepmly(3,1)=JJ kstepmly(3,1)=2 kspepmly(3,1)=LM istepmlz(3,1)=1 ispepmlz(3,1)=II jstepmlz(3,1)=2 jspepmlz(3,1)=JJ kstepmlz(3,1)=1 kspepmlz(3,1)=LM c For processor 2 istepmlx(1,2)=1 ispepmlx(1,2)=LM jstepmlx(1,2)=1 jspepmlx(1,2)=JJ kstepmlx(1,2)=2 kspepmlx(1,2)=KK istepmly(1,2)=2 ispepmly(1,2)=LM jstepmly(1,2)=1 jspepmly(1,2)=JJ kstepmly(1,2)=2 kspepmly(1,2)=KK istepmlz(1,2)=2 ispepmlz(1,2)=LM jstepmlz(1,2)=1 jspepmlz(1,2)=JJ kstepmlz(1,2)=1 kspepmlz(1,2)=KK istepmlx(2,2)=1 ispepmlx(2,2)=II jstepmlx(2,2)=JJ-LM+2 jspepmlx(2,2)=JJ kstepmlx(2,2)=2 kspepmlx(2,2)=KK istepmly(2,2)=2 ispepmly(2,2)=II jstepmly(2,2)=JJ-LM+1 jspepmly(2,2)=JJ kstepmly(2,2)=2 kspepmly(2,2)=KK istepmlz(2,2)=2 ispepmlz(2,2)=II jstepmlz(2,2)=JJ-LM+2 jspepmlz(2,2)=JJ kstepmlz(2,2)=1 kspepmlz(2,2)=KK istepmlx(3,2)=1 ispepmlx(3,2)=II jstepmlx(3,2)=1 jspepmlx(3,2)=JJ kstepmlx(3,2)=2 kspepmlx(3,2)=LM istepmly(3,2)=2 ispepmly(3,2)=II jstepmly(3,2)=1 jspepmly(3,2)=JJ kstepmly(3,2)=2 kspepmly(3,2)=LM istepmlz(3,2)=2 ispepmlz(3,2)=II jstepmlz(3,2)=1 jspepmlz(3,2)=JJ kstepmlz(3,2)=1 kspepmlz(3,2)=LM c For processor 3 istepmlx(1,3)=II-LM+1 ispepmlx(1,3)=II jstepmlx(1,3)=1 jspepmlx(1,3)=JJ kstepmlx(1,3)=2 kspepmlx(1,3)=KK istepmly(1,3)=II-LM+2 ispepmly(1,3)=II jstepmly(1,3)=1 jspepmly(1,3)=JJ kstepmly(1,3)=2 kspepmly(1,3)=KK istepmlz(1,3)=II-LM+2 ispepmlz(1,3)=II jstepmlz(1,3)=1 jspepmlz(1,3)=JJ kstepmlz(1,3)=1 kspepmlz(1,3)=KK istepmlx(2,3)=1 ispepmlx(2,3)=II jstepmlx(2,3)=JJ-LM+2 jspepmlx(2,3)=JJ kstepmlx(2,3)=2 kspepmlx(2,3)=KK istepmly(2,3)=1 ispepmly(2,3)=II jstepmly(2,3)=JJ-LM+1 jspepmly(2,3)=JJ kstepmly(2,3)=2 kspepmly(2,3)=KK istepmlz(2,3)=1 ispepmlz(2,3)=II jstepmlz(2,3)=JJ-LM+2 jspepmlz(2,3)=JJ kstepmlz(2,3)=1 kspepmlz(2,3)=KK istepmlx(3,3)=1 ispepmlx(3,3)=II jstepmlx(3,3)=1 jspepmlx(3,3)=JJ kstepmlx(3,3)=2 kspepmlx(3,3)=LM istepmly(3,3)=1 ispepmly(3,3)=II jstepmly(3,3)=1 jspepmly(3,3)=JJ kstepmly(3,3)=2 kspepmly(3,3)=LM istepmlz(3,3)=1 ispepmlz(3,3)=II jstepmlz(3,3)=1 jspepmlz(3,3)=JJ kstepmlz(3,3)=1 kspepmlz(3,3)=LM c For processor 4 istepmlx(1,4)=1 ispepmlx(1,4)=LM jstepmlx(1,4)=2 jspepmlx(1,4)=JJ kstepmlx(1,4)=1 kspepmlx(1,4)=KK istepmly(1,4)=2 ispepmly(1,4)=LM jstepmly(1,4)=1 jspepmly(1,4)=JJ kstepmly(1,4)=1 kspepmly(1,4)=KK istepmlz(1,4)=2 ispepmlz(1,4)=LM jstepmlz(1,4)=2 jspepmlz(1,4)=JJ kstepmlz(1,4)=1 kspepmlz(1,4)=KK istepmlx(2,4)=1 ispepmlx(2,4)=II jstepmlx(2,4)=2 jspepmlx(2,4)=LM kstepmlx(2,4)=1 kspepmlx(2,4)=KK istepmly(2,4)=2 ispepmly(2,4)=II jstepmly(2,4)=1 jspepmly(2,4)=LM kstepmly(2,4)=1 kspepmly(2,4)=KK istepmlz(2,4)=2 ispepmlz(2,4)=II jstepmlz(2,4)=2 jspepmlz(2,4)=LM kstepmlz(2,4)=1 kspepmlz(2,4)=KK istepmlx(3,4)=1 ispepmlx(3,4)=II jstepmlx(3,4)=2 jspepmlx(3,4)=JJ kstepmlx(3,4)=KK-LM+2 kspepmlx(3,4)=KK istepmly(3,4)=2 ispepmly(3,4)=II jstepmly(3,4)=1 jspepmly(3,4)=JJ kstepmly(3,4)=KK-LM+2 kspepmly(3,4)=KK istepmlz(3,4)=2 ispepmlz(3,4)=II jstepmlz(3,4)=2 jspepmlz(3,4)=JJ kstepmlz(3,4)=KK-LM+1 kspepmlz(3,4)=KK c For processor 5 istepmlx(1,5)=II-LM+1 ispepmlx(1,5)=II jstepmlx(1,5)=2 jspepmlx(1,5)=JJ kstepmlx(1,5)=1 kspepmlx(1,5)=KK istepmly(1,5)=II-LM+2 ispepmly(1,5)=II jstepmly(1,5)=1 jspepmly(1,5)=JJ kstepmly(1,5)=1 kspepmly(1,5)=KK istepmlz(1,5)=II-LM+2 ispepmlz(1,5)=II jstepmlz(1,5)=2 jspepmlz(1,5)=JJ kstepmlz(1,5)=1 kspepmlz(1,5)=KK istepmlx(2,5)=1 ispepmlx(2,5)=II jstepmlx(2,5)=2 jspepmlx(2,5)=LM kstepmlx(2,5)=1 kspepmlx(2,5)=KK istepmly(2,5)=1 ispepmly(2,5)=II jstepmly(2,5)=1 jspepmly(2,5)=LM kstepmly(2,5)=1 kspepmly(2,5)=KK istepmlz(2,5)=1 ispepmlz(2,5)=II jstepmlz(2,5)=2 jspepmlz(2,5)=LM kstepmlz(2,5)=1 kspepmlz(2,5)=KK istepmlx(3,5)=1 ispepmlx(3,5)=II jstepmlx(3,5)=2 jspepmlx(3,5)=JJ kstepmlx(3,5)=KK-LM+2 kspepmlx(3,5)=KK istepmly(3,5)=1 ispepmly(3,5)=II jstepmly(3,5)=1 jspepmly(3,5)=JJ kstepmly(3,5)=KK-LM+2 kspepmly(3,5)=KK istepmlz(3,5)=1 ispepmlz(3,5)=II jstepmlz(3,5)=2 jspepmlz(3,5)=JJ kstepmlz(3,5)=KK-LM+1 kspepmlz(3,5)=KK c For processor 6 istepmlx(1,6)=1 ispepmlx(1,6)=LM jstepmlx(1,6)=1 jspepmlx(1,6)=JJ kstepmlx(1,6)=1 kspepmlx(1,6)=KK istepmly(1,6)=2 ispepmly(1,6)=LM jstepmly(1,6)=1 jspepmly(1,6)=JJ kstepmly(1,6)=1 kspepmly(1,6)=KK istepmlz(1,6)=2 ispepmlz(1,6)=LM jstepmlz(1,6)=1 jspepmlz(1,6)=JJ kstepmlz(1,6)=1 kspepmlz(1,6)=KK istepmlx(2,6)=1 ispepmlx(2,6)=II jstepmlx(2,6)=JJ-LM+2 jspepmlx(2,6)=JJ kstepmlx(2,6)=1 kspepmlx(2,6)=KK istepmly(2,6)=2 ispepmly(2,6)=II jstepmly(2,6)=JJ-LM+1 jspepmly(2,6)=JJ kstepmly(2,6)=1 kspepmly(2,6)=KK istepmlz(2,6)=2 ispepmlz(2,6)=II jstepmlz(2,6)=JJ-LM+2 jspepmlz(2,6)=JJ kstepmlz(2,6)=1 kspepmlz(2,6)=KK istepmlx(3,6)=1 ispepmlx(3,6)=II jstepmlx(3,6)=1 jspepmlx(3,6)=JJ kstepmlx(3,6)=KK-LM+2 kspepmlx(3,6)=KK istepmly(3,6)=2 ispepmly(3,6)=II jstepmly(3,6)=1 jspepmly(3,6)=JJ kstepmly(3,6)=KK-LM+2 kspepmly(3,6)=KK istepmlz(3,6)=2 ispepmlz(3,6)=II jstepmlz(3,6)=1 jspepmlz(3,6)=JJ kstepmlz(3,6)=KK-LM+1 kspepmlz(3,6)=KK c For processor 7 istepmlx(1,7)=II-LM+1 ispepmlx(1,7)=II jstepmlx(1,7)=1 jspepmlx(1,7)=JJ kstepmlx(1,7)=1 kspepmlx(1,7)=KK istepmly(1,7)=II-LM+2 ispepmly(1,7)=II jstepmly(1,7)=1 jspepmly(1,7)=JJ kstepmly(1,7)=1 kspepmly(1,7)=KK istepmlz(1,7)=II-LM+2 ispepmlz(1,7)=II jstepmlz(1,7)=1 jspepmlz(1,7)=JJ kstepmlz(1,7)=1 kspepmlz(1,7)=KK istepmlx(2,7)=1 ispepmlx(2,7)=II jstepmlx(2,7)=JJ-LM+2 jspepmlx(2,7)=JJ kstepmlx(2,7)=1 kspepmlx(2,7)=KK istepmly(2,7)=1 ispepmly(2,7)=II jstepmly(2,7)=JJ-LM+1 jspepmly(2,7)=JJ kstepmly(2,7)=1 kspepmly(2,7)=KK istepmlz(2,7)=1 ispepmlz(2,7)=II jstepmlz(2,7)=JJ-LM+2 jspepmlz(2,7)=JJ kstepmlz(2,7)=1 kspepmlz(2,7)=KK istepmlx(3,7)=1 ispepmlx(3,7)=II jstepmlx(3,7)=1 jspepmlx(3,7)=JJ kstepmlx(3,7)=KK-LM+2 kspepmlx(3,7)=KK istepmly(3,7)=1 ispepmly(3,7)=II jstepmly(3,7)=1 jspepmly(3,7)=JJ kstepmly(3,7)=KK-LM+2 kspepmly(3,7)=KK istepmlz(3,7)=1 ispepmlz(3,7)=II jstepmlz(3,7)=1 jspepmlz(3,7)=JJ kstepmlz(3,7)=KK-LM+1 kspepmlz(3,7)=KK c For H-PML c For processor 0 isthpmlx(1,0)=2 isphpmlx(1,0)=LM jsthpmlx(1,0)=1 jsphpmlx(1,0)=JJ ksthpmlx(1,0)=1 ksphpmlx(1,0)=KK isthpmly(1,0)=1 isphpmly(1,0)=LM jsthpmly(1,0)=2 jsphpmly(1,0)=JJ ksthpmly(1,0)=1 ksphpmly(1,0)=KK isthpmlz(1,0)=1 isphpmlz(1,0)=LM jsthpmlz(1,0)=1 jsphpmlz(1,0)=JJ ksthpmlz(1,0)=2 ksphpmlz(1,0)=KK isthpmlx(2,0)=2 isphpmlx(2,0)=II jsthpmlx(2,0)=1 jsphpmlx(2,0)=LM ksthpmlx(2,0)=1 ksphpmlx(2,0)=KK isthpmly(2,0)=1 isphpmly(2,0)=II jsthpmly(2,0)=2 jsphpmly(2,0)=LM ksthpmly(2,0)=1 ksphpmly(2,0)=KK isthpmlz(2,0)=1 isphpmlz(2,0)=II jsthpmlz(2,0)=1 jsphpmlz(2,0)=LM ksthpmlz(2,0)=2 ksphpmlz(2,0)=KK isthpmlx(3,0)=2 isphpmlx(3,0)=II jsthpmlx(3,0)=1 jsphpmlx(3,0)=JJ ksthpmlx(3,0)=1 ksphpmlx(3,0)=LM isthpmly(3,0)=1 isphpmly(3,0)=II jsthpmly(3,0)=2 jsphpmly(3,0)=JJ ksthpmly(3,0)=1 ksphpmly(3,0)=LM isthpmlz(3,0)=1 isphpmlz(3,0)=II jsthpmlz(3,0)=1 jsphpmlz(3,0)=JJ ksthpmlz(3,0)=2 ksphpmlz(3,0)=LM c For processor 1 isthpmlx(1,1)=II-LM+2 isphpmlx(1,1)=II jsthpmlx(1,1)=1 jsphpmlx(1,1)=JJ ksthpmlx(1,1)=1 ksphpmlx(1,1)=KK isthpmly(1,1)=II-LM+1 isphpmly(1,1)=II jsthpmly(1,1)=2 jsphpmly(1,1)=JJ ksthpmly(1,1)=1 ksphpmly(1,1)=KK isthpmlz(1,1)=II-LM+1 isphpmlz(1,1)=II jsthpmlz(1,1)=1 jsphpmlz(1,1)=JJ ksthpmlz(1,1)=2 ksphpmlz(1,1)=KK isthpmlx(2,1)=1 isphpmlx(2,1)=II jsthpmlx(2,1)=1 jsphpmlx(2,1)=LM ksthpmlx(2,1)=1 ksphpmlx(2,1)=KK isthpmly(2,1)=1 isphpmly(2,1)=II jsthpmly(2,1)=2 jsphpmly(2,1)=LM ksthpmly(2,1)=1 ksphpmly(2,1)=KK isthpmlz(2,1)=1 isphpmlz(2,1)=II jsthpmlz(2,1)=1 jsphpmlz(2,1)=LM ksthpmlz(2,1)=2 ksphpmlz(2,1)=KK isthpmlx(3,1)=1 isphpmlx(3,1)=II jsthpmlx(3,1)=1 jsphpmlx(3,1)=JJ ksthpmlx(3,1)=1 ksphpmlx(3,1)=LM isthpmly(3,1)=1 isphpmly(3,1)=II jsthpmly(3,1)=2 jsphpmly(3,1)=JJ ksthpmly(3,1)=1 ksphpmly(3,1)=LM isthpmlz(3,1)=1 isphpmlz(3,1)=II jsthpmlz(3,1)=1 jsphpmlz(3,1)=JJ ksthpmlz(3,1)=2 ksphpmlz(3,1)=LM c For processor 2 isthpmlx(1,2)=2 isphpmlx(1,2)=LM jsthpmlx(1,2)=1 jsphpmlx(1,2)=JJ ksthpmlx(1,2)=1 ksphpmlx(1,2)=KK isthpmly(1,2)=1 isphpmly(1,2)=LM jsthpmly(1,2)=1 jsphpmly(1,2)=JJ ksthpmly(1,2)=1 ksphpmly(1,2)=KK isthpmlz(1,2)=1 isphpmlz(1,2)=LM jsthpmlz(1,2)=1 jsphpmlz(1,2)=JJ ksthpmlz(1,2)=2 ksphpmlz(1,2)=KK isthpmlx(2,2)=2 isphpmlx(2,2)=II jsthpmlx(2,2)=JJ-LM+1 jsphpmlx(2,2)=JJ ksthpmlx(2,2)=1 ksphpmlx(2,2)=KK isthpmly(2,2)=1 isphpmly(2,2)=II jsthpmly(2,2)=JJ-LM+2 jsphpmly(2,2)=JJ ksthpmly(2,2)=1 ksphpmly(2,2)=KK isthpmlz(2,2)=1 isphpmlz(2,2)=II jsthpmlz(2,2)=JJ-LM+1 jsphpmlz(2,2)=JJ ksthpmlz(2,2)=2 ksphpmlz(2,2)=KK isthpmlx(3,2)=2 isphpmlx(3,2)=II jsthpmlx(3,2)=1 jsphpmlx(3,2)=JJ ksthpmlx(3,2)=1 ksphpmlx(3,2)=LM isthpmly(3,2)=1 isphpmly(3,2)=II jsthpmly(3,2)=1 jsphpmly(3,2)=JJ ksthpmly(3,2)=1 ksphpmly(3,2)=LM isthpmlz(3,2)=1 isphpmlz(3,2)=II jsthpmlz(3,2)=1 jsphpmlz(3,2)=JJ ksthpmlz(3,2)=2 ksphpmlz(3,2)=LM c For processor 3 isthpmlx(1,3)=II-LM+2 isphpmlx(1,3)=II jsthpmlx(1,3)=1 jsphpmlx(1,3)=JJ ksthpmlx(1,3)=1 ksphpmlx(1,3)=KK isthpmly(1,3)=II-LM+1 isphpmly(1,3)=II jsthpmly(1,3)=1 jsphpmly(1,3)=JJ ksthpmly(1,3)=1 ksphpmly(1,3)=KK isthpmlz(1,3)=II-LM+1 isphpmlz(1,3)=II jsthpmlz(1,3)=1 jsphpmlz(1,3)=JJ ksthpmlz(1,3)=2 ksphpmlz(1,3)=KK isthpmlx(2,3)=1 isphpmlx(2,3)=II jsthpmlx(2,3)=JJ-LM+1 jsphpmlx(2,3)=JJ ksthpmlx(2,3)=1 ksphpmlx(2,3)=KK isthpmly(2,3)=1 isphpmly(2,3)=II jsthpmly(2,3)=JJ-LM+2 jsphpmly(2,3)=JJ ksthpmly(2,3)=1 ksphpmly(2,3)=KK isthpmlz(2,3)=1 isphpmlz(2,3)=II jsthpmlz(2,3)=JJ-LM+1 jsphpmlz(2,3)=JJ ksthpmlz(2,3)=2 ksphpmlz(2,3)=KK isthpmlx(3,3)=1 isphpmlx(3,3)=II jsthpmlx(3,3)=1 jsphpmlx(3,3)=JJ ksthpmlx(3,3)=1 ksphpmlx(3,3)=LM isthpmly(3,3)=1 isphpmly(3,3)=II jsthpmly(3,3)=1 jsphpmly(3,3)=JJ ksthpmly(3,3)=1 ksphpmly(3,3)=LM isthpmlz(3,3)=1 isphpmlz(3,3)=II jsthpmlz(3,3)=1 jsphpmlz(3,3)=JJ ksthpmlz(3,3)=2 ksphpmlz(3,3)=LM c For processor 4 isthpmlx(1,4)=2 isphpmlx(1,4)=LM jsthpmlx(1,4)=1 jsphpmlx(1,4)=JJ ksthpmlx(1,4)=1 ksphpmlx(1,4)=KK isthpmly(1,4)=1 isphpmly(1,4)=LM jsthpmly(1,4)=2 jsphpmly(1,4)=JJ ksthpmly(1,4)=1 ksphpmly(1,4)=KK isthpmlz(1,4)=1 isphpmlz(1,4)=LM jsthpmlz(1,4)=1 jsphpmlz(1,4)=JJ ksthpmlz(1,4)=1 ksphpmlz(1,4)=KK isthpmlx(2,4)=2 isphpmlx(2,4)=II jsthpmlx(2,4)=1 jsphpmlx(2,4)=LM ksthpmlx(2,4)=1 ksphpmlx(2,4)=KK isthpmly(2,4)=1 isphpmly(2,4)=II jsthpmly(2,4)=2 jsphpmly(2,4)=LM ksthpmly(2,4)=1 ksphpmly(2,4)=KK isthpmlz(2,4)=1 isphpmlz(2,4)=II jsthpmlz(2,4)=1 jsphpmlz(2,4)=LM ksthpmlz(2,4)=1 ksphpmlz(2,4)=KK isthpmlx(3,4)=2 isphpmlx(3,4)=II jsthpmlx(3,4)=1 jsphpmlx(3,4)=JJ ksthpmlx(3,4)=KK-LM+1 ksphpmlx(3,4)=KK isthpmly(3,4)=1 isphpmly(3,4)=II jsthpmly(3,4)=2 jsphpmly(3,4)=JJ ksthpmly(3,4)=KK-LM+1 ksphpmly(3,4)=KK isthpmlz(3,4)=1 isphpmlz(3,4)=II jsthpmlz(3,4)=1 jsphpmlz(3,4)=JJ ksthpmlz(3,4)=KK-LM+2 ksphpmlz(3,4)=KK c For processor 5 isthpmlx(1,5)=II-LM+2 isphpmlx(1,5)=II jsthpmlx(1,5)=1 jsphpmlx(1,5)=JJ ksthpmlx(1,5)=1 ksphpmlx(1,5)=KK isthpmly(1,5)=II-LM+1 isphpmly(1,5)=II jsthpmly(1,5)=2 jsphpmly(1,5)=JJ ksthpmly(1,5)=1 ksphpmly(1,5)=KK isthpmlz(1,5)=II-LM+1 isphpmlz(1,5)=II jsthpmlz(1,5)=1 jsphpmlz(1,5)=JJ ksthpmlz(1,5)=1 ksphpmlz(1,5)=KK isthpmlx(2,5)=1 isphpmlx(2,5)=II jsthpmlx(2,5)=1 jsphpmlx(2,5)=LM ksthpmlx(2,5)=1 ksphpmlx(2,5)=KK isthpmly(2,5)=1 isphpmly(2,5)=II jsthpmly(2,5)=2 jsphpmly(2,5)=LM ksthpmly(2,5)=1 ksphpmly(2,5)=KK isthpmlz(2,5)=1 isphpmlz(2,5)=II jsthpmlz(2,5)=1 jsphpmlz(2,5)=LM ksthpmlz(2,5)=1 ksphpmlz(2,5)=KK isthpmlx(3,5)=1 isphpmlx(3,5)=II jsthpmlx(3,5)=1 jsphpmlx(3,5)=JJ ksthpmlx(3,5)=KK-LM+1 ksphpmlx(3,5)=KK isthpmly(3,5)=1 isphpmly(3,5)=II jsthpmly(3,5)=2 jsphpmly(3,5)=JJ ksthpmly(3,5)=KK-LM+1 ksphpmly(3,5)=KK isthpmlz(3,5)=1 isphpmlz(3,5)=II jsthpmlz(3,5)=1 jsphpmlz(3,5)=JJ ksthpmlz(3,5)=KK-LM+2 ksphpmlz(3,5)=KK c For processor 6 isthpmlx(1,6)=2 isphpmlx(1,6)=LM jsthpmlx(1,6)=1 jsphpmlx(1,6)=JJ ksthpmlx(1,6)=1 ksphpmlx(1,6)=KK isthpmly(1,6)=1 isphpmly(1,6)=LM jsthpmly(1,6)=1 jsphpmly(1,6)=JJ ksthpmly(1,6)=1 ksphpmly(1,6)=KK isthpmlz(1,6)=1 isphpmlz(1,6)=LM jsthpmlz(1,6)=1 jsphpmlz(1,6)=JJ ksthpmlz(1,6)=1 ksphpmlz(1,6)=KK isthpmlx(2,6)=2 isphpmlx(2,6)=II jsthpmlx(2,6)=JJ-LM+1 jsphpmlx(2,6)=JJ ksthpmlx(2,6)=1 ksphpmlx(2,6)=KK isthpmly(2,6)=1 isphpmly(2,6)=II jsthpmly(2,6)=JJ-LM+2 jsphpmly(2,6)=JJ ksthpmly(2,6)=1 ksphpmly(2,6)=KK isthpmlz(2,6)=1 isphpmlz(2,6)=II jsthpmlz(2,6)=JJ-LM+1 jsphpmlz(2,6)=JJ ksthpmlz(2,6)=1 ksphpmlz(2,6)=KK isthpmlx(3,6)=2 isphpmlx(3,6)=II jsthpmlx(3,6)=1 jsphpmlx(3,6)=JJ ksthpmlx(3,6)=KK-LM+1 ksphpmlx(3,6)=KK isthpmly(3,6)=1 isphpmly(3,6)=II jsthpmly(3,6)=1 jsphpmly(3,6)=JJ ksthpmly(3,6)=KK-LM+1 ksphpmly(3,6)=KK isthpmlz(3,6)=1 isphpmlz(3,6)=II jsthpmlz(3,6)=1 jsphpmlz(3,6)=JJ ksthpmlz(3,6)=KK-LM+2 ksphpmlz(3,6)=KK c For processor 7 isthpmlx(1,7)=II-LM+2 isphpmlx(1,7)=II jsthpmlx(1,7)=1 jsphpmlx(1,7)=JJ ksthpmlx(1,7)=1 ksphpmlx(1,7)=KK isthpmly(1,7)=II-LM+1 isphpmly(1,7)=II jsthpmly(1,7)=1 jsphpmly(1,7)=JJ ksthpmly(1,7)=1 ksphpmly(1,7)=KK isthpmlz(1,7)=II-LM+1 isphpmlz(1,7)=II jsthpmlz(1,7)=1 jsphpmlz(1,7)=JJ ksthpmlz(1,7)=1 ksphpmlz(1,7)=KK isthpmlx(2,7)=1 isphpmlx(2,7)=II jsthpmlx(2,7)=JJ-LM+1 jsphpmlx(2,7)=JJ ksthpmlx(2,7)=1 ksphpmlx(2,7)=KK isthpmly(2,7)=1 isphpmly(2,7)=II jsthpmly(2,7)=JJ-LM+2 jsphpmly(2,7)=JJ ksthpmly(2,7)=1 ksphpmly(2,7)=KK isthpmlz(2,7)=1 isphpmlz(2,7)=II jsthpmlz(2,7)=JJ-LM+1 jsphpmlz(2,7)=JJ ksthpmlz(2,7)=1 ksphpmlz(2,7)=KK isthpmlx(3,7)=1 isphpmlx(3,7)=II jsthpmlx(3,7)=1 jsphpmlx(3,7)=JJ ksthpmlx(3,7)=KK-LM+1 ksphpmlx(3,7)=KK isthpmly(3,7)=1 isphpmly(3,7)=II jsthpmly(3,7)=1 jsphpmly(3,7)=JJ ksthpmly(3,7)=KK-LM+1 ksphpmly(3,7)=KK isthpmlz(3,7)=1 isphpmlz(3,7)=II jsthpmlz(3,7)=1 jsphpmlz(3,7)=JJ ksthpmlz(3,7)=KK-LM+2 ksphpmlz(3,7)=KK c PML cell number lpmlst(1)=1 lpmlst(2)=lpmlst(1)+(LM+1)*(JJ+1)*(KK+1) lpmlst(3)=lpmlst(2)+(II+1)*(LM+1)*(KK+1) return end c E-field subroutine efld(E,H,Cex,Cey,Cez,Cexl,Ceyl,Cezl,II,JJ,KK,DX,DY,DZ, & istart,iexstop,ieystop,iezstop,jstart,jexstop,jeystop, & jezstop,kstart,kexstop,keystop,kezstop,myrank) integer i,j,k integer istart(0:7),iexstop(0:7),ieystop(0:7),iezstop(0:7), & jstart(0:7),jexstop(0:7),jeystop(0:7),jezstop(0:7), & kstart(0:7),kexstop(0:7),keystop(0:7),kezstop(0:7) double precision & E(II+1,JJ+1,KK+1,3),H(0:II,0:JJ,0:KK,3), & Cex(II+1,JJ+1,KK+1),Cey(II+1,JJ+1,KK+1),Cez(II+1,JJ+1,KK+1), & Cexl(II+1,JJ+1,KK+1),Ceyl(II+1,JJ+1,KK+1), & Cezl(II+1,JJ+1,KK+1) do k=kstart(myrank),kexstop(myrank) do j=jstart(myrank),jexstop(myrank) do i=istart(myrank),iexstop(myrank) E(i,j,k,1)=Cex(i,j,k)*E(i,j,k,1) & +Cexl(i,j,k)/DY*(H(i,j,k,3)-H(i,j-1,k,3)) & -Cexl(i,j,k)/DZ*(H(i,j,k,2)-H(i,j,k-1,2)) end do end do end do do k=kstart(myrank),keystop(myrank) do j=jstart(myrank),jeystop(myrank) do i=istart(myrank),ieystop(myrank) E(i,j,k,2)=Cey(i,j,k)*E(i,j,k,2) & +Ceyl(i,j,k)/DZ*(H(i,j,k,1)-H(i,j,k-1,1)) & -Ceyl(i,j,k)/DX*(H(i,j,k,3)-H(i-1,j,k,3)) end do end do end do do k=kstart(myrank),kezstop(myrank) do j=jstart(myrank),jezstop(myrank) do i=istart(myrank),iezstop(myrank) E(i,j,k,3)=Cez(i,j,k)*E(i,j,k,3) & +Cezl(i,j,k)/DX*(H(i,j,k,2)-H(i-1,j,k,2)) & -Cezl(i,j,k)/DY*(H(i,j,k,1)-H(i,j-1,k,1)) end do end do end do return end c H-field subroutine hfld(E,H,Chx,Chy,Chz,II,JJ,KK, & istart,ihxstop,ihystop,ihzstop,jstart,jhxstop,jhystop, & jhzstop,kstart,khxstop,khystop,khzstop,myrank) integer i,j,k integer istart(0:7),ihxstop(0:7),ihystop(0:7),ihzstop(0:7), & jstart(0:7),jhxstop(0:7),jhystop(0:7),jhzstop(0:7), & kstart(0:7),khxstop(0:7),khystop(0:7),khzstop(0:7) double precision & E(II+1,JJ+1,KK+1,3),H(0:II,0:JJ,0:KK,3),Chx,Chy,Chz do k=kstart(myrank),khxstop(myrank) do j=jstart(myrank),jhxstop(myrank) do i=istart(myrank),ihxstop(myrank) H(i,j,k,1)=H(i,j,k,1) & -Chy*(E(i,j+1,k,3)-E(i,j,k,3)) & +Chz*(E(i,j,k+1,2)-E(i,j,k,2)) end do end do end do do k=kstart(myrank),khystop(myrank) do j=jstart(myrank),jhystop(myrank) do i=istart(myrank),ihystop(myrank) H(i,j,k,2)=H(i,j,k,2) & -Chz*(E(i,j,k+1,1)-E(i,j,k,1)) & +Chx*(E(i+1,j,k,3)-E(i,j,k,3)) end do end do end do do k=kstart(myrank),khzstop(myrank) do j=jstart(myrank),jhzstop(myrank) do i=istart(myrank),ihzstop(myrank) H(i,j,k,3)=H(i,j,k,3) & -Chx*(E(i+1,j,k,2)-E(i,j,k,2)) & +Chy*(E(i,j+1,k,1)-E(i,j,k,1)) end do end do end do return end c E-field pml subroutine epml(E,H,Exy,Exz,Eyx,Eyz,Ezx,Ezy, & Cedx,Cedy,Cedz,Cedxl,Cedyl,Cedzl,II,JJ,KK,NPML, & istepmlx,ispepmlx,jstepmlx,jspepmlx,kstepmlx,kspepmlx, & istepmly,ispepmly,jstepmly,jspepmly,kstepmly,kspepmly, & istepmlz,ispepmlz,jstepmlz,jspepmlz,kstepmlz,kspepmlz, & lpmlst,myrank) integer i,j,k,l,lx,ly,lz integer istepmlx(3,0:7),ispepmlx(3,0:7),jstepmlx(3,0:7), & jspepmlx(3,0:7),kstepmlx(3,0:7),kspepmlx(3,0:7), & istepmly(3,0:7),ispepmly(3,0:7),jstepmly(3,0:7), & jspepmly(3,0:7),kstepmly(3,0:7),kspepmly(3,0:7), & istepmlz(3,0:7),ispepmlz(3,0:7),jstepmlz(3,0:7), & jspepmlz(3,0:7),kstepmlz(3,0:7),kspepmlz(3,0:7),lpmlst(3) double precision E(II+1,JJ+1,KK+1,3),H(0:II,0:JJ,0:KK,3), & Exy(NPML),Exz(NPML),Eyx(NPML),Eyz(NPML),Ezx(NPML),Ezy(NPML), & Cedx(II+1),Cedy(JJ+1),Cedz(KK+1),Cedxl(II+1),Cedyl(JJ+1), & Cedzl(KK+1) do l=1,3 lx=lpmlst(l) do k=kstepmlx(l,myrank),kspepmlx(l,myrank) do j=jstepmlx(l,myrank),jspepmlx(l,myrank) do i=istepmlx(l,myrank),ispepmlx(l,myrank) Exy(lx)=Cedy(j)*Exy(lx) & +Cedyl(j)*(H(i,j,k,3)-H(i,j-1,k,3)) Exz(lx)=Cedz(k)*Exz(lx) & +Cedzl(k)*(H(i,j,k-1,2)-H(i,j,k,2)) E(i,j,k,1)=Exy(lx)+Exz(lx) lx=lx+1 end do end do end do ly=lpmlst(l) do k=kstepmly(l,myrank),kspepmly(l,myrank) do j=jstepmly(l,myrank),jspepmly(l,myrank) do i=istepmly(l,myrank),ispepmly(l,myrank) Eyx(ly)=Cedx(i)*Eyx(ly) & +Cedxl(i)*(H(i-1,j,k,3)-H(i,j,k,3)) Eyz(ly)=Cedz(k)*Eyz(ly) & +Cedzl(k)*(H(i,j,k,1)-H(i,j,k-1,1)) E(i,j,k,2)=Eyx(ly)+Eyz(ly) ly=ly+1 end do end do end do lz=lpmlst(l) do k=kstepmlz(l,myrank),kspepmlz(l,myrank) do j=jstepmlz(l,myrank),jspepmlz(l,myrank) do i=istepmlz(l,myrank),ispepmlz(l,myrank) Ezx(lz)=Cedx(i)*Ezx(lz) & +Cedxl(i)*(H(i,j,k,2)-H(i-1,j,k,2)) Ezy(lz)=Cedy(j)*Ezy(lz) & +Cedyl(j)*(H(i,j-1,k,1)-H(i,j,k,1)) E(i,j,k,3)=Ezx(lz)+Ezy(lz) lz=lz+1 end do end do end do end do return end c H-field pml subroutine hpml(E,H,Hxy,Hxz,Hyx,Hyz,Hzx,Hzy, & Chdx,Chdy,Chdz,Chdxl,Chdyl,Chdzl,II,JJ,KK,NPML, & isthpmlx,isphpmlx,jsthpmlx,jsphpmlx,ksthpmlx,ksphpmlx, & isthpmly,isphpmly,jsthpmly,jsphpmly,ksthpmly,ksphpmly, & isthpmlz,isphpmlz,jsthpmlz,jsphpmlz,ksthpmlz,ksphpmlz, & lpmlst,myrank) integer i,j,k,l,lx,ly,lz integer isthpmlx(3,0:7),isphpmlx(3,0:7),jsthpmlx(3,0:7), & jsphpmlx(3,0:7),ksthpmlx(3,0:7),ksphpmlx(3,0:7), & isthpmly(3,0:7),isphpmly(3,0:7),jsthpmly(3,0:7), & jsphpmly(3,0:7),ksthpmly(3,0:7),ksphpmly(3,0:7), & isthpmlz(3,0:7),isphpmlz(3,0:7),jsthpmlz(3,0:7), & jsphpmlz(3,0:7),ksthpmlz(3,0:7),ksphpmlz(3,0:7),lpmlst(3) double precision E(II+1,JJ+1,KK+1,3),H(0:II,0:JJ,0:KK,3), & Hxy(NPML),Hxz(NPML),Hyx(NPML),Hyz(NPML),Hzx(NPML),Hzy(NPML), & Chdx(II+1),Chdy(JJ+1),Chdz(KK+1),Chdxl(II+1),Chdyl(JJ+1), & Chdzl(KK+1) do l=1,3 lx=lpmlst(l) do k=ksthpmlx(l,myrank),ksphpmlx(l,myrank) do j=jsthpmlx(l,myrank),jsphpmlx(l,myrank) do i=isthpmlx(l,myrank),isphpmlx(l,myrank) Hxy(lx)=Chdy(j)*Hxy(lx) & -Chdyl(j)*(E(i,j+1,k,3)-E(i,j,k,3)) Hxz(lx)=Chdz(k)*Hxz(lx) & +Chdzl(k)*(E(i,j,k+1,2)-E(i,j,k,2)) H(i,j,k,1)=Hxy(lx)+Hxz(lx) lx=lx+1 end do end do end do ly=lpmlst(l) do k=ksthpmly(l,myrank),ksphpmly(l,myrank) do j=jsthpmly(l,myrank),jsphpmly(l,myrank) do i=isthpmly(l,myrank),isphpmly(l,myrank) Hyx(ly)=Chdx(i)*Hyx(ly) & +Chdxl(i)*(E(i+1,j,k,3)-E(i,j,k,3)) Hyz(ly)=Chdz(k)*Hyz(ly) & -Chdzl(k)*(E(i,j,k+1,1)-E(i,j,k,1)) H(i,j,k,2)=Hyx(ly)+Hyz(ly) ly=ly+1 end do end do end do lz=lpmlst(l) do k=ksthpmlz(l,myrank),ksphpmlz(l,myrank) do j=jsthpmlz(l,myrank),jsphpmlz(l,myrank) do i=isthpmlz(l,myrank),isphpmlz(l,myrank) Hzx(lz)=Chdx(i)*Hzx(lz) & -Chdxl(i)*(E(i+1,j,k,2)-E(i,j,k,2)) Hzy(lz)=Chdy(j)*Hzy(lz) & +Chdyl(j)*(E(i,j+1,k,1)-E(i,j,k,1)) H(i,j,k,3)=Hzx(lz)+Hzy(lz) lz=lz+1 end do end do end do end do return end c subroutine for saving temporary voltage and current c save value on own processor subroutine savvi(H,Vgap,Cgap,n,NSTOP, & II,JJ,KK,IA,JA,KGAP,DX,DY) integer n double precision H(0:II,0:JJ,0:KK,3),Vgap(NSTOP),Cgap(NSTOP) Cgap(n)=(-H(IA,JA,KGAP,1)+H(IA,JA-1,KGAP,1))*DX & +(H(IA,JA,KGAP,2)-H(IA-1,JA,KGAP,2))*DY write(02,*) Vgap(n),Cgap(n) return end