total2.f90

fortran source code for fdtd parall computing

  end do
end do

!*********************************连接边界引入入射波Ez分量*********************

do I=Itmin+1,Itmax-1
  Ez(I,Jtmax)= Ez(I,Jtmax)-FE*HCB(I,3)
  Ez(I,Jtmin)= Ez(I,Jtmin)+FE*HCB(I,1)
end do

do J=Jtmin+1,Jtmax-1
  Ez(Itmax,J)= Ez(Itmax,J)+FE*HCB(J,2)
  Ez(Itmin,J)= Ez(Itmin,J)-FE*HCB(J,4)
end do


!*********************************连接边界角点引入入射波************************
T1=HCB(Jtmax,2)-HCB(Itmax,3)
Ez(Itmax,Jtmax)= Ez(Itmax,Jtmax)+FE*T1
 
T1=-HCB(Jtmin,4)+HCB(Itmin,1)
Ez(Itmin,Jtmin)= Ez(Itmin,Jtmin)+FE*T1

T1=HCB(Jtmin,2)+HCB(Itmax,1)
Ez(Itmax,Jtmin)= Ez(Itmax,Jtmin)+FE*T1

T1=-HCB(Itmin,3)-HCB(Jtmax,4)
Ez(Itmin,Jtmax)= Ez(Itmin,Jtmax)+FE*T1

!*********************************吸收边界条件Ez*******************************
do I=Imin+1,Imax-1
   
   T1=Hy(I,Jmin)-Hy(I-1,Jmin)
   T1=T1+Hy(I,Jmin+1)-Hy(I-1,Jmin+1)
   T1=T1*TEMFlag
   T2=Ez(I,Jmin+1)-Ez(I,Jmin)
   Ez(I,Jmin)=EAB(I,1)-.33333333*T2
   Ez(I,Jmin)=Ez(I,Jmin)+.16666667*T1
   EAB(I,1)=Ez(I,Jmin+1)  

   T1=Hy(I,Jmax)-Hy(I-1,Jmax)
   T1=T1+Hy(I,Jmax-1)-Hy(I-1,Jmax-1)
   T1=T1*TEMFlag
   T2=Ez(I,Jmax-1)-Ez(I,Jmax)
   Ez(I,Jmax)=EAB(I,3)-.33333333*T2
   Ez(I,Jmax)=Ez(I,Jmax)+.16666667*T1
   EAB(I,3)=Ez(I,Jmax-1)

end do

do J=Jmin+1,Jmax-1

  T1=Hx(Imin,J)-Hx(Imin,J-1)
  T1=T1+Hx(Imin+1,J)-Hx(Imin+1,J-1)
  T1=T1*TEMFlag
  T2=Ez(Imin+1,J)-Ez(Imin,J)
  Ez(Imin,J)=EAB(J,4)-.33333333*T2
  Ez(Imin,J)=Ez(Imin,J)-.16666667*T1
  EAB(J,4)=Ez(Imin+1,J)

  
  T1=Hx(Imax,J)-Hx(Imax,J-1)
  T1=T1+Hx(Imax-1,J)-Hx(Imax-1,J-1)
  T1=T1*TEMFlag
  T2=Ez(Imax-1,J)-Ez(Imax,J)
  Ez(Imax,J)=EAB(J,2)-.33333333*T2
  Ez(Imax,J)=Ez(Imax,J)-.16666667*T1
  EAB(J,2)=Ez(Imax-1,J)
end do

! *********************吸收边界角点Ez************

Iflag=1-Iflag         !交替存储标志

Ez(Imin,Jmin)=.2928932*EAC(1,1,Iflag)+.7071068*EAC(1,2,Iflag)
EAC(1,1,Iflag)=Ez(Imin,Jmin)
EAC(1,2,Iflag)=Ez(Imin+1,Jmin+1)

Ez(Imax,Jmin)=.2928932*EAC(2,1,Iflag)+.7071068*EAC(2,2,Iflag)
EAC(2,1,Iflag)=Ez(Imax,Jmin)
EAC(2,2,Iflag)=Ez(Imax-1,Jmin+1)

Ez(Imin,Jmax)=.2928932*EAC(3,1,Iflag)+.7071068*EAC(3,2,Iflag)
EAC(3,1,Iflag)=Ez(Imin,Jmax)
EAC(3,2,Iflag)=Ez(Imin+1,Jmax-1)

Ez(Imax,Jmax)=.2928932*EAC(4,1,Iflag)+.7071068*EAC(4,2,Iflag)
EAC(4,1,Iflag)=Ez(Imax,Jmax)
EAC(4,2,Iflag)=Ez(Imax-1,Jmax-1)
! **************Hx,Hy的FDTD迭代********************

do I=Imin,Imax
   do J=Jmin,Jmax-1
      T1=Ez(I,J+1)-Ez(I,J)
      Hx(I,J)=FH1(ob(i,j),ob(i,j+1))*Hx(I,J)-FH2(ob(i,j),ob(i,j+1))*T1
      end do
end do

do I=Imin,Imax-1
  do J=Jmin,Jmax
     T1=(Ez(I+1,J)-Ez(I,J))
     Hy(I,J)=FH1(ob(i,j),ob(i+1,j))*Hy(I,J)+FH2(ob(i,j),ob(i+1,j))*T1
   end do
end do

! *******************连接边界引入入射波Hx,Hy分量*****************

do I=Itmin,Itmax
  Hx(I,Jtmin-1)=Hx(I,Jtmin-1)+FH*ECB(I,1)
  Hx(I,Jtmax)=Hx(I,Jtmax)-FH*ECB(I,3)
end do

do J=Jtmin,Jtmax
  Hy(Itmin-1,J)=Hy(Itmin-1,J)-FH*ECB(J,4)
  Hy(Itmax,J)=Hy(Itmax,J)+FH*ECB(J,2)
end do



end do                       ! 时间步主循环结束

! **********************输出结果*****************************

if(OUTflag.EQ.'I') then            !如果输出标志为‘I’,则输出场量为虚部
  write(*,'(1H+,10X,32h Now outputing IMAGE part ......)')
  OPEN(1,FILE='FDTD.EI',ACCESS='DIRECT',RECL=4)
  OPEN(4,FILE='FDTD.HxI',ACCESS='DIRECT',RECL=4)
  OPEN(5,FILE='FDTD.HyI',ACCESS='DIRECT',RECL=4)
  OPEN(2,FILE='FDTD.BI')
  EI0=Ein(0)
else
  write(*,'(1H+,10X,31H Now outputing REAL part ......)')
  OPEN(1,FILE='FDTD.ER',ACCESS='DIRECT',RECL=4)
  OPEN(4,FILE='FDTD.HxR',ACCESS='DIRECT',RECL=4)
  OPEN(5,FILE='FDTD.HyR',ACCESS='DIRECT',RECL=4)
  OPEN(2,FILE='FDTD.BR')
  ER0=Ein(0)
end if

NN=0
do I=Imin,Imax                          !输出整个计算域的场量的实部或虚部到临时文件
  do J=Jmin,Jmax
    NN=NN+1
    write(1,REC=NN)Ez(I,J)
    write(4,REC=NN)Hx(I,J)/Z
    write(5,REC=NN)Hy(I,J)/Z
  end do
end do
CLOSE(1)
CLOSE(4)
CLOSE(5)


do I=Iomin,Iomax                         !输出输出边界上场量的实部或虚部到临时文件中 
   T1=.5*(Ez(I,Jomin)+Ez(I,Jomin-1))     !下边界
   write(2,*)T1,Hx(I,Jomin-1)/Z
   T1=.5*(Ez(I,Jomax)+Ez(I,Jomax+1))     !上边界
   write(2,*)T1,Hx(I,Jomax)/Z
end do
do J=Jomin,Jomax                          
   T1=.5*(Ez(Iomax,J)+Ez(Iomax+1,J))     !右边界
   write(2,*)T1,Hy(Iomax,J)/Z
   T1=.5*(Ez(Iomin,J)+Ez(Iomin-1,J))     !上边界
   write(2,*)T1,Hy(Iomin-1,J)/Z
end do
CLOSE(2)

if(OUTflag.EQ.'I')then
  OUTflag='R'
  TimeStop=TimeStop+WL/2                 !延迟WL/2时间步以输出实部
  GOTO 999
end if

! ***************************将输出结果组成复数*********************************************

write(*,'(1H+,10X,40h Make.PIC files ......                   )')

CT1 = (1.,0.)/CMPLX(ER0,EI0)                              !入射波电场在原点处的值的倒数
CT2 = CEXP(CMPLX(0,-.5*OMIGA))*CT1                        !CT1比CT2倒数滞后半个时间步
                                                          !CT1,CT2将用来分别对电场和磁场归一
OPEN(1,FILE='FDTD.EI',ACCESS = 'DIRECT',RECL = 4)
OPEN(2,FILE='FDTD.ER',ACCESS = 'DIRECT',RECL = 4)
FileName=FileName(1:FileNameLength)//'AZ.PIC'
OPEN(4,FILE=FileName, ACCESS = 'DIRECT',RECL = 4)
FileName=FileName(1:FileNameLength)//'PZ.PIC'
OPEN(5,FILE=FileName, ACCESS = 'DIRECT',RECL = 4)

NN = 0
do I = Imin,Imax                                     !输出整个计算域的Ez(幅值和相位,用CT1归一)
  do J = Jmin,Jmax
  NN = NN + 1
  read(1,REC=NN) T1
  read(2,REC=NN) T2
  CT3 = CMPLX(T2,T1)*CT1
  write(4,REC = NN) ABS(CT3)
  T1 = ATAN2(IMAG(CT3),real(CT3))
  write(5,REC = NN) T1
 end do 
end do
CLOSE(1)
CLOSE(2)
CLOSE(4)
CLOSE(5)


OPEN(6,FILE='FDTD.HxI',ACCESS = 'DIRECT',RECL=4)
OPEN(7,FILE='FDTD.HxR',ACCESS = 'DIRECT',RECL=4)
OPEN(8,FILE='FDTD.HyI',ACCESS = 'DIRECT',RECL=4)
OPEN(9,FILE='FDTD.HyR',ACCESS = 'DIRECT',RECL=4)
FileName=FileName(1:FileNameLength)//'AX.PIC'
OPEN(10,FILE=FileName,ACCESS='DIRECT',RECL=4)
FileName=FileName(1:FileNameLength)//'PX.PIC'
OPEN(11,FILE=FileName,ACCESS='DIRECT',RECL=4)
FileName=FileName(1:FileNameLength)//'AY.PIC'
OPEN(12,FILE=FileName,ACCESS='DIRECT',RECL=4)
FileName=FileName(1:FileNameLength)//'PY.PIC'
OPEN(13,FILE=FileName,ACCESS='DIRECT',RECL=4)


NN=0
do I=Imin,Imax                                     !输出整个计算域的Hx,Hy(幅值和相位,用CT2归一)
  do J=Jmin,Jmax
    NN=NN+1
    read(6,REC=NN) T1
    read(7,REC=NN) T2
    CT3=CMPLX(T2,T1)*CT2
    write(10,REC=NN)ABS(CT3)
    T1=ATAN2(IMAG(CT3),real(CT3))
    write(11,REC=NN) T1

    read(8,REC=NN) T1
    read(9,REC=NN) T2
    CT3=CMPLX(T2,T1)*CT2
    write(12,REC=NN)ABS(CT3)
    T1=ATAN2(IMAG(CT3),real(CT3))
    write(13,REC=NN) T1

  end do
end do
CLOSE(6)
CLOSE(7)
CLOSE(8)
CLOSE(9)
CLOSE(10)
CLOSE(11)
CLOSE(12)
CLOSE(13)

OPEN(1,FILE='FDTD.BI')
OPEN(2,FILE='FDTD.BR')
FileName=FileName(1:FileNameLength)//'.BND'

    OPEN(3,FILE=filename,access='direct',recl=8)


    write(3,rec=1)WaveLength                   !输出输出边界的场值(幅值和相位，用CT1或CT2归一)
    write(3,rec=2)WL
    write(3,rec=3)TEMFlag
    write(3,rec=4)Phi*180/Pi
    write(3,rec=5)Imin
    write(3,rec=6)Imax
    write(3,rec=7)Jmin
    write(3,rec=8)Jmax
    write(3,rec=9)Itmin
    write(3,rec=10)Itmax
    write(3,rec=11)Jtmin
    write(3,rec=12)Jtmax
    write(3,rec=13)Iomin
    write(3,rec=14)Iomax
    write(3,rec=15)Jomin
    write(3,rec=16)Jomax
    write(3,rec=17)TimeStop-WL/2


    NN=17
    do iii=1,10000
      read(1,*,end = 1111)T1,T2
      read(2,*)T3,T4
      NN=NN+1
      CT3=CMPLX(T3,T1)*CT1
      write(3,rec=NN)CT3
      NN=NN+1
      CT4=CMPLX(T4,T2)*CT2
      write(3,rec=NN)CT4
    end do
1111 continue
     CLOSE(1)
     CLOSE(2)
     CLOSE(3)
     
                                  ! 输出.REM(与.PIC对应)文件以便用Photo.exe显示计算区域场值
     Step=WaveLength/WL
     MaxX=Imax*Step
     MixY=Imin*Step
     MaxY=Jmax*Step
     MixY=Jmin*Step
     XNum=Imax-Imin+1
     YNum=Jmax-Jmin+1


FileName=FileName(1:FileNameLength)//'AZ.REM'
open(1,file=FileName)

    write(1,*)'''Title'''''''''''''''
    write(1,*)''' '''
    write(rem,101)MaxX
101 format(1x,''' MaxX:',f9.4,'m''')
    write(1,*)rem
    write(rem,102)MinX
102 format(1x,'''MinX:',f9.4,'m''')
    write(1,*)rem
    write(rem,103)MaxY
103 format(1x,'''MaxY:',f9.4,'m''')
    write(1,*)rem
    write(rem,104)MinY
104 format(1x,'''MinY:',f9.4,'m''')
    write(1,*)rem
    write(rem,105)Step
105 format(1x,'''Step:',f9.4,'m''')
    write(1,*)rem
    write(rem,106)XNum
106 format(1x,'''XNum:',I4,'(grid)''')
    write(1,*)rem
    write(rem,107)YNum
107 format(1x,'''YNum:',I4,'(grid)''')
    write(1,*)rem
    write(1,*)''' '''   

    write(1,*)'''WaveMode:'//WaveMode,' '''

    write(rem,108)Wavelength
108 format(1x,'''WL =',f9.4,'m''')
    write(1,*)rem
    write(rem,109)WL
109 format(1x,'''=',I4,'grid''')
    write(1,*)rem
    write(1,*)'''Incident Angle '''

    write(rem,110)Phi*180/Pi
110 format(1x,''' ',f7.2,'degree''')
    write(1,*)rem

    write(rem,111)Timestop-W1/2
111 format(1x,'''TimeStep =',I5,' ''')
    write(1,*)rem

    write(1,*)''' '''
    write(1,*)''' '''
    write(1,*)MaxX,MinX,MaxY,MinY
    write(1,*)XNum,YNum,Step
    close(1)


    FileName = FileName(1:FileNameLength)//'PZ.REM'
    open(1,file = FileName)
    write(1,*)'''Title'''''''''''''''
    write(1,*)''' '''
    write(rem,101)MaxX
    write(1,*)rem
    write(rem,102)MinX
    write(1,*)rem
    write(rem,103)MaxY
    write(1,*)rem
    write(rem,104)MinY
    write(1,*)rem
    write(rem,105)Step
    write(1,*)rem
    write(rem,106)XNum
    write(1,*)rem
    write(rem,107)YNum
    write(1,*)rem
    write(1,*)''' '''

    write(1,*)'''WaveMode:'//WaveMode,' '''


    write(rem,108)Wavelength
    write(1,*)rem
    write(rem,109)WL
    write(1,*)rem
    write(1,*)'''Incident Angle '''
    write(rem,110)Phi*180/Pi
    write(1,*)rem
    write(rem,111)Timestop-WL/2
    write(1,*)rem
    write(1,*)''' '''
    write(1,*)''' '''
    write(1,*)MaxX,MinX,MaxY,MinY
    write(1,*)XNum,YNum,Step
    close(1)


    open(1,file = 'FDTD,EI')       !删除临时文件 
    close(1,status = 'delete')
    open(2,file = 'FDTD.ER')
    close(2,status = 'delete')
    open(3,file = 'FDTD.HXI')
    close(3,status = 'delete')
    open(4,file = 'FDTD.HXR')
    close(4,status = 'delete')
    open(5,file = 'FDTD.HYI')
    close(5,status = 'delete')
    open(6,file = 'FDTD.HYR')
    close(6,status = 'delete')
    open(7,file = 'FDTD.BI')
    close(7,status = 'delete')
    open(8,file = 'FDTD.BR')
    close(8,status = 'delete')


    end