[资料共享]fortran编写的3D的PML程序铁出来，大家都看看

!*********************************************************************** @?kJ).  
!     3-D FDTD code with CPML absorbing boundary conditions 8)j@aiF`  
!*********************************************************************** h4)Bs\==mT  
! 7TX2&kMoc  
!     Program author: Jamesina J. Simpson, Ph.D. Graduate Student CvDxq:x  
!                     National Science Foundation Fellow <RXwM6G2  
!                     NU Computational Electromagnetics Laboratory `@Oa lg  
!                     Director: Allen Taflove `yO'[2  
!                     Department of Electrical Engineering and `<frgXu64  
!                                                      Computer Science vScjq5"p  
!                     Northwestern University [k<1`z3  
!                     2145 Sheridan Road d`y!cu2}  
!                     Evanston, IL 60208 B^v8,;jZT  
!                     [url=mailto:j-simpson@northwestern.edu]j-simpson@northwestern.edu[/url] @'6"7g  
! "_1)CDqP  
!     Copyright 2005 *O~e T  
! `]i []|  
!     This FORTRAN 90 code implements the finite-difference time-domain Bc5+ss  
!     solution of Maxwell's curl equations over a three-dimensional WsGths+[  
!     Cartesian space lattice.  The grid is terminated by CPML absorbing h6:#!Rg  
!     boundary conditions.  However, in a straightforward manner it can iE].&>w  
!     be altered to have PEC planes at any of the outer grid boundaries.  `3cCH  
!     Also, the number of grid cells as well as the thickness of !t#F/C  
!     the PML in each Cartesian direction can be varied independently. E4z)Mr#  
! Fhf<T`  
!     For illustrative purposes, the code supplied here models the PML 6%Mt  
!     numerical experiment described in Section 7.11.2.  Three output 'b&yrBFD  
!     files are generated having the following data: A8r^)QJP{  
! |nUl\WRd\  
!       (1) Ez at the source for each time step; K.~q+IYP[  
!       (2) Ey at the probe point for each time step;  and j5$GFi\kB  
!       (3) The plane of Ez values 1 mm above the source in the k- !l@zT}i??  
!           direction recorded at the time step, "record_grid".  This GVGlVAo|@  
!           data can be viewed in Matlab using the following commands:
==N{1gO]  
! ;C.S3}  
!             >>  load view_grid.dat $IQw=w7p  
!             >>  for j = 1:Jmax 0#_'o ,  
!             >>      image(1:Imax,j) = view_grid(1+(j-1)*Imax:j*Imax); pzZk\-0R  
!             >>  end j<<3Pr  
!             >>  pcolor(image'); shading flat C[sh,  
! }UwO<#  
!     The relative error (equation 7.135) graphed in Fig. 7.6 on EgjR^A1W2  
!     page 318 can be reproduced by comparing the output file 4rg2y]  
!     "probe.dat" with that of a much larger benchmark grid having LeOP;#  
!     Imax, Jmax, and Kmax increased to the values mentioned in the text. Z4b||  
! \ 0W!4
D  
!     This code has been tested using the Intel Fortran Compiler 8.0 for 8' DW#%  
!     Linux.  The executable can be created by typing \M<3}t  
!     "ifort fdtd3D_CPML.f90" at the command prompt. p4!:]0c  
! WKfkKk;G  
!     This program is not guaranteed to be free of defects or bugs. aB~S?.l  
!     Please report any bugs that may exist to: V)ig)(CT  
!                                          [url=mailto:j-simpson@northwestern.edu]j-simpson@northwestern.edu[/url] (V:E2WR  
! 8^/I>0EZ  
! Ifc]K?  
!*********************************************************************** u\yVR$pQ  
uVE.,)xz  
   PROGRAM fdtd3D_CPML W~1~k{A  
   IMPLICIT NONE c2nKPEX&5  
!  .................................. M5357Q  
!  Input Fundamental Constants (MKS units)
=FP0\cQ.  
   REAL, PARAMETER ::                            & Z @ dC+0[=  
      pi = 3.14159265358979, C = 2.99792458E8, & FS6`6M.K  
      muO = 4.0 * pi * 1.0E-7, epsO = 1.0/(C*C*muO) ').}Nz  
!  .................................. Yr.sm!xA  
!  Specify Material Relative Permittivity and Conductivity 4<dcB@v  
   REAL, PARAMETER::                      & w obgu  
      epsR = 1.0, sigM1 = 0.0   ! free space 'a6<ixgo0  
!  .................................. (V`Md\NL`  
!  Specify Number of Time Steps and Grid Size Parameters FzQTDu9  
   INTEGER, PARAMETER ::                                     & hj%ye~|~  
      nmax = 2100, &  ! total number of time steps :Qt  
      Imax = 51, Jmax = 126, Kmax = 26  ! grid size corresponding to the W !w,f;  
                                        ! number of Ez field components DyJ.BQdk)  
!  .................................. b;sVls  
!  Specify Grid Cell Size in Each Direction and Calculate the fxaJZz$o  
!  Resulting Courant-Stable Time Step {76c%<`WaP  
   REAL, PARAMETER ::                                        & /P|fB]p  
      dx = 1.0E-3, dy = 1.0E-3, dz = 1.0E-3,  &  ! cell size in each direction u$#Wv2|mk  
      dt = 0.99 / (C*(1.0/dx**2+1.0/dy**2+1.0/dz**2)**0.5) Yb3mP!3q8Z  
                                                 ! time step increment CC=d I  
!  .................................. RGKYW>$0RR  
!  Specify the Impulsive Source (See Equation 7.134) CAq/K?:8  
   REAL, PARAMETER ::                                        & H,3\0BKk  
      tw = 53.0E-12, tO = 4.0*tw  nfSbM3D]h  
!  .................................. PXzsj
.  
!  Specify the Time Step at which the Grid is Recorded for Visualization x+8_4>,>Y7  
   INTEGER, PARAMETER ::                                        & .IBp\7W!?E  
      record_grid = 300 c/N@zum,{  
!  .................................. IYn]U4P.  
!  Specify the PEC Plate Boundaries and the Source/Recording Points }o4N<%/+  
   INTEGER, PARAMETER ::                                    & D"(L5jR8m@  
      istart = (Imax-1)/2-11, iend = istart+24, jstart = istart,   & EP'h@zdz  
      jend = jstart + 99, kstart = Kmax/2, kend = kstart,      & JxE53ev  
      isource = istart, jsource = jstart, ksource = kend+1,     & R0ID2:i]F  
      irecv1 = iend, jrecv1 = jend+1, krecv1 = kend  ! Ey at probe point ~c^>54  
!  .................................. g- AHdYJ  
!  Specify the CPML Thickness in Each Direction (Value of Zero W&*f#
E  
!  Corresponds to No PML, and the Grid is Terminated with a PEC) TxjYrzC  
   INTEGER, PARAMETER ::                        & zi R5:d3   
      ! PML thickness in each direction 9vUO
*D  
      nxPML_1 = 11, nxPML_2 = nxPML_1, nyPML_1 = nxPML_1,      & fZ04!R  
      nyPML_2 = nxPML_1, nzPML_1 = nxPML_1, nzPML_2 = nxPML_2 ^z1&8k"[^  
!  .................................. Qd}m`YW-f$  
!  Specify the CPML Order and Other Parameters f"7O "6  
   INTEGER, PARAMETER ::                        & @AHm!9?o  
      m = 3, ma = 1 9EDfd NN  
   REAL, PARAMETER  ::                     & y&}E~5O  
      sig_x_max = 0.75 * (0.8*(m+1)/(dx*(muO/epsO*epsR)**0.5)),   & voP7"Dl[  
      sig_y_max = 0.75 * (0.8*(m+1)/(dy*(muO/epsO*epsR)**0.5)),   & A1'hlAGF  
      sig_z_max = 0.75 * (0.8*(m+1)/(dz*(muO/epsO*epsR)**0.5)),   & &,A64y  
      alpha_x_max = 0.24,   & 3vhnwDcK  
      alpha_y_max = alpha_x_max, alpha_z_max = alpha_x_max, & [[PEa-992  
      kappa_x_max = 15.0, & "Ve.cP,7(  
      kappa_y_max = kappa_x_max, kappa_z_max = kappa_x_max AjcX  N  
   INTEGER ::                                                & 0"2=n.##  
 i,j,ii,jj,k,kk,n U*Ge<(v$  
   REAL  ::                                                   & @k[R/,#'[t  
      source, P1, P2 F<>!kK/c  
!     TM components SDu#Yt&mhh  
   REAL,DIMENSION(Imax, Jmax, Kmax)  ::                      & .+7;)K
  
      Ez, CA, CB, sig, eps {!6/x9>  
   REAL,DIMENSION(Imax-1, Jmax, Kmax)  ::                      & fkf1m:Ckh  
      Hy 5;0g!&-t#  
   REAL,DIMENSION(Imax,Jmax-1, Kmax)  ::                      & Bjk]ZU0T  
      Hx ]8q3>  
   REAL  ::                        & t!*+8Q!e  
      DA, DB x~xa6  
!     TE components /|#&px)G  
   REAL,DIMENSION(Imax-1, Jmax-1, Kmax-1)  ::                      & @1*ohdHH  
      Hz 5` Te\H  
   REAL,DIMENSION(Imax-1, Jmax, Kmax-1)  ::                      & Ee4&g<X.  
      Ex l^ZI* z7N  
   REAL,DIMENSION(Imax,Jmax-1, Kmax-1)  ::                      & O"Ku1t!  
      Ey yjfat&$  
!  PML .CL^BiD.D  
   REAL ,DIMENSION(nxPML_1,Jmax,Kmax) ::                       &
.ObZ\.I  
      psi_Ezx_1 2gO2jJlv  
   REAL ,DIMENSION(nxPML_2,Jmax,Kmax) ::                       & 0/S_e)U  
      psi_Ezx_2 S
*j6OwZ
 
   REAL ,DIMENSION(nxPML_1-1,Jmax,Kmax) ::                       & hX`}Q4(k  
      psi_Hyx_1 0#w?HCx=  
   REAL ,DIMENSION(nxPML_2-1,Jmax,Kmax) ::                       & U2uF&6v  
      psi_Hyx_2 8} =JKR^cK  
   REAL ,DIMENSION(Imax,nyPML_1,Kmax) ::                       & >e\9Bf_  
      psi_Ezy_1                               @d Jr/6Yx  
   REAL ,DIMENSION(Imax,nyPML_2,Kmax) ::                       & DXz}YIEC  
      psi_Ezy_2 wP[t0/dl  
   REAL ,DIMENSION(Imax,nyPML_1-1,Kmax) ::                       & f$NMM >z  
      psi_Hxy_1                               zSo)k~&[3  
   REAL ,DIMENSION(Imax,nyPML_2-1,Kmax) ::                       & I%- " |]$  
      psi_Hxy_2 64hk2a8  
   REAL ,DIMENSION(Imax,Jmax-1,nzPML_1-1) ::                       & j'9"cE5_  
      psi_Hxz_1 X[f=h=|  
   REAL ,DIMENSION(Imax,Jmax-1,nzPML_2-1) ::                       & 6-?66gmT  
      psi_Hxz_2 (Qa/EkE^*w  
   REAL ,DIMENSION(Imax-1,Jmax,nzPML_1-1) ::                       & 311LC cRp  
      psi_Hyz_1 ylt`*|$  
   REAL ,DIMENSION(Imax-1,Jmax,nzPML_2-1) ::                       & $O9^SB  
      psi_Hyz_2 fS~;>n%R  
   REAL ,DIMENSION(Imax-1,Jmax,nzPML_1) ::                       & sR/b$j>i3  
      psi_Exz_1 626Z5Afg  
   REAL ,DIMENSION(Imax-1,Jmax,nzPML_2) ::                       & N<QXmgqx  
      psi_Exz_2 ^i:`ZfA#
 
   REAL ,DIMENSION(Imax,Jmax-1,nzPML_1) ::                       & =-_B:d;  
      psi_Eyz_1 Cvn#=6V3  
   REAL ,DIMENSION(Imax,Jmax-1,nzPML_2) ::                       & {?'fyEeg  
      psi_Eyz_2 B}PIRk@a1  
   REAL ,DIMENSION(nxPML_1-1,Jmax-1,Kmax-1) ::                       & TMYd47  
      psi_Hzx_1 \[Z?&  
   REAL ,DIMENSION(nxPML_1,Jmax-1,Kmax-1) ::                       & | WvUq  
      psi_Eyx_1 3p0LN'q]A  
   REAL ,DIMENSION(nxPML_2-1,Jmax-1,Kmax-1) ::                       & @V03a )6,h  
      psi_Hzx_2 Eb=}FuV  
   REAL ,DIMENSION(nxPML_2,Jmax-1,Kmax-1) ::                       & =[F<7pvE  
      psi_Eyx_2 cA%%IL$R  
   REAL ,DIMENSION(Imax-1,nyPML_1-1,Kmax-1) ::                       & G:zua`u[  
      psi_Hzy_1                               
s kg*  
   REAL ,DIMENSION(Imax-1,nyPML_1,Kmax-1) ::                       & 2WU@*%sk"  
      psi_Exy_1                               H$I=W>;  
   REAL ,DIMENSION(Imax-1,nyPML_2-1,Kmax-1) ::                       & r*Yi1j/  
      psi_Hzy_2 [Q=dCX9%  
   REAL ,DIMENSION(Imax-1,nyPML_2,Kmax-1) ::                       & 76u&EG%  
      psi_Exy_2 4:U?u  
   REAL ,DIMENSION(nxPML_1) ::                       & 5nsq[Q`  
      be_x_1, ce_x_1, alphae_x_PML_1, sige_x_PML_1, kappae_x_PML_1 ]('D^Ro  
   REAL ,DIMENSION(nxPML_1-1) ::                       & ?O>V%@  
      bh_x_1, ch_x_1, alphah_x_PML_1, sigh_x_PML_1, kappah_x_PML_1 rETRTp0HT  
   REAL ,DIMENSION(nxPML_2) ::                       & X7s `U5'l  
      be_x_2, ce_x_2, alphae_x_PML_2, sige_x_PML_2, kappae_x_PML_2 FJasS8  
   REAL ,DIMENSION(nxPML_2-1) ::                       & lP<:tR~K  
      bh_x_2, ch_x_2, alphah_x_PML_2, sigh_x_PML_2, kappah_x_PML_2 P2bZ65>3y  
   REAL ,DIMENSION(nyPML_1) ::                       & xO-+i\ ZV  
      be_y_1, ce_y_1, alphae_y_PML_1, sige_y_PML_1, kappae_y_PML_1 27;ci:5  
   REAL ,DIMENSION(nyPML_1-1) ::                       & 7)s^8+  
      bh_y_1, ch_y_1, alphah_y_PML_1, sigh_y_PML_1, kappah_y_PML_1 =$]uoA  
   REAL ,DIMENSION(nyPML_2) ::                       & &^W|iXi
#  
      be_y_2, ce_y_2, alphae_y_PML_2, sige_y_PML_2, kappae_y_PML_2 w8n|B?Sr  
   REAL ,DIMENSION(nyPML_2-1) ::                       & l
sJnI|  
      bh_y_2, ch_y_2, alphah_y_PML_2, sigh_y_PML_2, kappah_y_PML_2 cReB~wk  
   REAL ,DIMENSION(nzPML_1) ::                       & !Q#{o^{Y~  
      be_z_1, ce_z_1, alphae_z_PML_1, sige_z_PML_1, kappae_z_PML_1 {mAU3x  
   REAL ,DIMENSION(nzPML_1-1) ::                       & ,?L2wl[  
      bh_z_1, ch_z_1, alphah_z_PML_1, sigh_z_PML_1, kappah_z_PML_1 1Tu *79A  
   REAL ,DIMENSION(nzPML_2) ::                       & 8MSC.0  
      be_z_2, ce_z_2, alphae_z_PML_2, sige_z_PML_2, kappae_z_PML_2 ~8lwe*lNV  
   REAL ,DIMENSION(nzPML_2-1) ::                       & XLH0 ;+CL{  
      bh_z_2, ch_z_2, alphah_z_PML_2, sigh_z_PML_2, kappah_z_PML_2 s8}@=]aA  
!     denominators for the update equations OF\rgz
 
   REAL,DIMENSION(Imax-1)  ::                      & dPhQ :sd>  
      den_ex, den_hx 9uGrk^<t  
   REAL,DIMENSION(Jmax-1)  ::                      & j?MAED  
      den_ey, den_hy Q[nEsYP  
   REAL,DIMENSION(Kmax-1)  ::                      & U"
Zmv  
      den_ez, den_hz Uw47LP  
[7[0^ad  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ V O3x~E  
!  OPEN OUTPUT FILES qHxqQ'ks;  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ TA{\PKA)  
   OPEN (UNIT = 30, FILE = "source.dat") ^Vpq$'!  
   OPEN (UNIT = 31, FILE = "probe.dat") xD6@Qk  
   OPEN (UNIT = 33, FILE = "view_grid.dat") t+eVR8  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ WL'P)lI5  
!  INITIALIZE VARIABLES >Q(3*d >  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ )kP5u`v  
   P1 = 0.0 )b:7-}d  
   P2 = 0.0 Y'
ow  
   Ez(:,:,:) = 0.0 KQ~i<1&j  
   Hz(:,:,:) = 0.0 ~$~5qwl  
   Ex(:,:,:) = 0.0 /n>qCuw  
   Hx(:,:,:) = 0.0 z_|/5$T>U  
   Ey(:,:,:) = 0.0
dH-s2r%s  
   Hy(:,:,:) = 0.0 #FNcF>3>  
   sig(:,:,:) = sigM1 @;'o2   
   eps(:,:,:) = epsR*epsO iRrl^\qn  
   psi_Exy_1(:,:,:) = 0.0 &uO%_6J  
   psi_Exy_2(:,:,:) = 0.0 Y./2Ely  
   psi_Exz_1(:,:,:) = 0.0 \_lod kf  
   psi_Exz_2(:,:,:) = 0.0 d+'p@!W_  
   psi_Eyx_1(:,:,:) = 0.0 HaP}Y:p  
   psi_Eyx_2(:,:,:) = 0.0 1;{Rhu7* k  
   psi_Eyz_1(:,:,:) = 0.0  -4
cXRv]  
   psi_Eyz_2(:,:,:) = 0.0 NUp<e%zB  
   psi_Ezy_1(:,:,:) = 0.0 yxH ( c  
   psi_Ezy_2(:,:,:) = 0.0 m%76i;uP  
   psi_Ezx_1(:,:,:) = 0.0 zRyuq1Zyc,  
   psi_Ezx_2(:,:,:) = 0.0 OVV]x{  
   psi_Hxy_1(:,:,:) = 0.0 RO.k]x6  
   psi_Hxy_2(:,:,:) = 0.0 <WcR,d  
   psi_Hxz_1(:,:,:) = 0.0 ^Y'HaneoM  
   psi_Hxz_2(:,:,:) = 0.0 :53)Nv  
   psi_Hyx_1(:,:,:) = 0.0 oZAB_A)[-  
   psi_Hyx_2(:,:,:) = 0.0 Ry'= ke  
   psi_Hyz_1(:,:,:) = 0.0
[Nk3|u`h  
   psi_Hyz_2(:,:,:) = 0.0 V>b\[(=s  
   psi_Hzy_1(:,:,:) = 0.0 wYmM"60  
   psi_Hzy_2(:,:,:) = 0.0 ;1a~pF S  
   psi_Hzx_1(:,:,:) = 0.0 ~P@Q7T*  
   psi_Hzx_2(:,:,:) = 0.0 Cvf[/C+  
   write(*,*)"Imax: ", Imax 7R`mf   
   write(*,*)"Jmax: ", Jmax CO`_^7o9(  
   write(*,*)"Kmax: ", Kmax mFi&YpHu3  
   write(*,*)"dt: ", dt % ClHCoyA  
   write(*,*)"nmax: ", nmax 49vcoHlf  
   write(*,*)"max time: ", nmax*dt '#&os`mQ  
   write(*,*)"record grid after ", record_grid, "dt" o|$AyS{1  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ } q$ WvY/  
!  SET CPML PARAMETERS IN EACH DIRECTION e:5bzk!~  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ~>Aj  
   DO i = 1,nxPML_1 `|/<\  
      sige_x_PML_1(i) = sig_x_max * ( (nxPML_1 - i) / (nxPML_1 - 1.0) )**m 9o+)?1\  
      alphae_x_PML_1(i) = alpha_x_max*((i-1.0)/(nxPML_1-1.0))**ma ]bK=FIK2  
      kappae_x_PML_1(i) = 1.0+(kappa_x_max-1.0)*   & ~x|F)~:0=  
                                 ((nxPML_1 - i) / (nxPML_1 - 1.0))**m #*KNPh  
      be_x_1(i) = EXP(-(sige_x_PML_1(i) / kappae_x_PML_1(i) +   & AT%u%cE-  
                                 alphae_x_PML_1(i))*dt/epsO) f`;j:O  
      if ((sige_x_PML_1(i) == 0.0) .and.        & D ^x-^6^  
         (alphae_x_PML_1(i) == 0.0) .and. (i == nxPML_1)) then {QwHc5Bf  
         ce_x_1(i) = 0.0 D}mjN=Y  
      else Ba!J"b]  
         ce_x_1(i) = sige_x_PML_1(i)*(be_x_1(i)-1.0)/       & ;LMJd@  
            (sige_x_PML_1(i)+kappae_x_PML_1(i)*alphae_x_PML_1(i)) & PBp^|t]E>  
            / kappae_x_PML_1(i) 3nhXZOO1  
      endif 7:2WgLo  
   ENDDO lMBLIB
]i  
   DO i = 1,nxPML_1-1 43@{JK9G  
      sigh_x_PML_1(i) = sig_x_max * ( (nxPML_1 - i - 0.5)/(nxPML_1-1.0))**m K2<9mDn&  
      alphah_x_PML_1(i) = alpha_x_max*((i-0.5)/(nxPML_1-1.0))**ma lqm1!5dt  
      kappah_x_PML_1(i) = 1.0+(kappa_x_max-1.0)*   & m8o(J\]  
                            ((nxPML_1 - i - 0.5) / (nxPML_1 - 1.0))**m 6E^.7%3  
      bh_x_1(i) = EXP(-(sigh_x_PML_1(i) / kappah_x_PML_1(i) +   & [Lz
w#XE  
                                 alphah_x_PML_1(i))*dt/epsO) rsy'q(N[  
      ch_x_1(i) = sigh_x_PML_1(i)*(bh_x_1(i)-1.0)/      & smnSDS  
                  (sigh_x_PML_1(i)+kappah_x_PML_1(i)*alphah_x_PML_1(i)) & |b)Y#)C;  
                  / kappah_x_PML_1(i) u4/kR  
   ENDDO ]4pkcV P  
   DO i = 1,nxPML_2 {$wjO7Glp  
      sige_x_PML_2(i) = sig_x_max * ( (nxPML_2 - i) / (nxPML_2 - 1.0) )**m fe9LEM8j  
      alphae_x_PML_2(i) = alpha_x_max*((i-1.0)/(nxPML_2-1.0))**ma !Y[lQXv  
      kappae_x_PML_2(i) = 1.0+(kappa_x_max-1.0)*   & W`u[h0\c  
                                 ((nxPML_2 - i) / (nxPML_2 - 1.0))**m fyByz=pl  
      be_x_2(i) = EXP(-(sige_x_PML_2(i) / kappae_x_PML_2(i) +   & z[3L2U~6  
                                 alphae_x_PML_2(i))*dt/epsO) c\Q7"!e  
      if ((sige_x_PML_2(i) == 0.0) .and.        & hX=A)73(  
         (alphae_x_PML_2(i) == 0.0) .and. (i == nxPML_2)) then BYMi6wts  
         ce_x_2(i) = 0.0  {PVWD7  
      else #n#@fAY  
         ce_x_2(i) = sige_x_PML_2(i)*(be_x_2(i)-1.0)/       & ekWePL;rR2  
            (sige_x_PML_2(i)+kappae_x_PML_2(i)*alphae_x_PML_2(i)) & K[!&b0O  
            / kappae_x_PML_2(i) <x<"n t  
      endif 3yB!M  
   ENDDO vuoQz\  
   DO i = 1,nxPML_2-1 =_:et0  
      sigh_x_PML_2(i) = sig_x_max * ( (nxPML_2 - i - 0.5)/(nxPML_2-1.0))**m egq67S  
      alphah_x_PML_2(i) = alpha_x_max*((i-0.5)/(nxPML_2-1.0))**ma o*o/q],C9-  
      kappah_x_PML_2(i) = 1.0+(kappa_x_max-1.0)*   & ;<(W% _  
                            ((nxPML_2 - i - 0.5) / (nxPML_2 - 1.0))**m L{8xlx`  
      bh_x_2(i) = EXP(-(sigh_x_PML_2(i) / kappah_x_PML_2(i) +   & Y4v|ko`l%  
                                 alphah_x_PML_2(i))*dt/epsO) !K'}K>iT  
      ch_x_2(i) = sigh_x_PML_2(i)*(bh_x_2(i)-1.0)/      & JW3B'_0  
                  (sigh_x_PML_2(i)+kappah_x_PML_2(i)*alphah_x_PML_2(i)) & FSA1gAW6g  
                  / kappah_x_PML_2(i) <=>=.kmGt  
   ENDDO ;JayoJ  
   DO j = 1,nyPML_1 P7{gfiB  
      sige_y_PML_1(j) = sig_y_max * ( (nyPML_1 - j ) / (nyPML_1 - 1.0) )**m Qx+%"YO  
      alphae_y_PML_1(j) = alpha_y_max*((j-1)/(nyPML_1-1.0))**ma %'X[^W  
      kappae_y_PML_1(j) = 1.0+(kappa_y_max-1.0)*   & }Nf%n@  
                                 ((nyPML_1 - j) / (nyPML_1 - 1.0))**m AD|2qM))  
      be_y_1(j) = EXP(-(sige_y_PML_1(j) / kappae_y_PML_1(j) +   & :>Ay^{vf=  
                                 alphae_y_PML_1(j))*dt/epsO) bp~g;h*E2  
      if ((sige_y_PML_1(j) == 0.0) .and.        & R~S
;sJ& c  
         (alphae_y_PML_1(j) == 0.0) .and. (j == nyPML_1)) then ~,i-8jl,  
         ce_y_1(j) = 0.0 }n4
V|f-  
      else `.i!NBA'6  
         ce_y_1(j) = sige_y_PML_1(j)*(be_y_1(j)-1.0)/       & /N&CaH\;^$  
            (sige_y_PML_1(j)+kappae_y_PML_1(j)*alphae_y_PML_1(j)) & C,N
Jb+J  
            / kappae_y_PML_1(j) BS:+~|3w  
      endif _Nqt21sL  
   ENDDO j4^97  
   DO j = 1,nyPML_1-1 pP*a  
      sigh_y_PML_1(j) = sig_y_max * ( (nyPML_1 - j - 0.5)/(nyPML_1-1.0))**m *tK\R&4,4s  
      alphah_y_PML_1(j) = alpha_y_max*((j-0.5)/(nyPML_1-1.0))**ma '2[albxSc  
      kappah_y_PML_1(j) = 1.0+(kappa_y_max-1.0)*   & %>g3~yl  
                            ((nyPML_1 - j - 0.5) / (nyPML_1 - 1.0))**m +Zaj,oEE  
      bh_y_1(j) = EXP(-(sigh_y_PML_1(j) / kappah_y_PML_1(j) +   & 2fWTY0  
                                 alphah_y_PML_1(j))*dt/epsO) R5X.^u  
      ch_y_1(j) = sigh_y_PML_1(j)*(bh_y_1(j)-1.0)/      & Ndi9FD3im  
                  (sigh_y_PML_1(j)+kappah_y_PML_1(j)*alphah_y_PML_1(j)) & Yi$vg  
                  / kappah_y_PML_1(j) 61)-cVC  
   ENDDO x3o]U)^  
   DO j = 1,nyPML_2 (adyZ/j  
      sige_y_PML_2(j) = sig_y_max * ( (nyPML_2 - j ) / (nyPML_2 - 1.0) )**m
HBu
[gh;b  
      alphae_y_PML_2(j) = alpha_y_max*((j-1)/(nyPML_2-1.0))**ma SUU !7Yd|  
      kappae_y_PML_2(j) = 1.0+(kappa_y_max-1.0)*   & TzNn^ir=HX  
                                 ((nyPML_2 - j) / (nyPML_2 - 1.0))**m u[DfzH  
      be_y_2(j) = EXP(-(sige_y_PML_2(j) / kappae_y_PML_2(j) +   & E.J
kf\  
                                 alphae_y_PML_2(j))*dt/epsO) d'G0m9u2  
      if ((sige_y_PML_2(j) == 0.0) .and.        & ~wkj&yVT  
         (alphae_y_PML_2(j) == 0.0) .and. (j == nyPML_2)) then =*r])Vg^  
         ce_y_2(j) = 0.0 yrC7F`.  
      else C<m{*C-`a  
         ce_y_2(j) = sige_y_PML_2(j)*(be_y_2(j)-1.0)/       & Y07ZB'K  
            (sige_y_PML_2(j)+kappae_y_PML_2(j)*alphae_y_PML_2(j)) & i?
 K|TC`  
            / kappae_y_PML_2(j) fmb} 2h  
      endif x7eQ2h6O  
   ENDDO V0P>YQq9s  
   DO j = 1,nyPML_2-1 P_Gw-`L5T  
      sigh_y_PML_2(j) = sig_y_max * ( (nyPML_2 - j - 0.5)/(nyPML_2-1.0))**m @Bf%s(Uj+  
      alphah_y_PML_2(j) = alpha_y_max*((j-0.5)/(nyPML_2-1.0))**ma lXjXqk\  
      kappah_y_PML_2(j) = 1.0+(kappa_y_max-1.0)*   & .O0+H+  
                            ((nyPML_2 - j - 0.5) / (nyPML_2 - 1.0))**m -)RH5WGS  
      bh_y_2(j) = EXP(-(sigh_y_PML_2(j) / kappah_y_PML_2(j) +   & MP4z-4Y  
                                 alphah_y_PML_2(j))*dt/epsO) #0y)U;dA+w  
      ch_y_2(j) = sigh_y_PML_2(j)*(bh_y_2(j)-1.0)/      & /#m=*&!CB  
                  (sigh_y_PML_2(j)+kappah_y_PML_2(j)*alphah_y_PML_2(j)) & >8qQK r\"  
                  / kappah_y_PML_2(j) )|DM~%$QM  
   ENDDO 4
}eepJOn  
   DO k = 1,nzPML_1 /Sc l#4bW  
      sige_z_PML_1(k) = sig_z_max * ( (nzPML_1 - k ) / (nzPML_1 - 1.0) )**m
,#L=v]  
      alphae_z_PML_1(k) = alpha_z_max*((k-1)/(nzPML_1-1.0))**ma g.zEn/SM  
      kappae_z_PML_1(k) = 1.0+(kappa_z_max-1.0)*   & F7O(Cy"1  
                                 ((nzPML_1 - k) / (nzPML_1 - 1.0))**m {l
/-LZ.  
      be_z_1(k) = EXP(-(sige_z_PML_1(k) / kappae_z_PML_1(k) +   & Jc|6&  
                                 alphae_z_PML_1(k))*dt/epsO) JNJ=e,O,  
      if ((sige_z_PML_1(k) == 0.0) .and.        & ".<DAs j  
         (alphae_z_PML_1(k) == 0.0) .and. (k == nzPML_1)) then }wHW7SJ  
         ce_z_1(k) = 0.0 <U!`J[n%  
      else Na?!;1]_  
         ce_z_1(k) = sige_z_PML_1(k)*(be_z_1(k)-1.0)/       & *fn*h[pV&  
            (sige_z_PML_1(k)+kappae_z_PML_1(k)*alphae_z_PML_1(k)) & 5*,f Fib  
            / kappae_z_PML_1(k) IHcD*zQ  
      endif )~HUo9K9  
   ENDDO Wb?8j M  
   DO k = 1,nzPML_1-1 6 1F(
<!  
      sigh_z_PML_1(k) = sig_z_max * ( (nzPML_1 - k - 0.5)/(nzPML_1-1.0))**m 93` AWg/T  
      alphah_z_PML_1(k) = alpha_z_max*((k-0.5)/(nzPML_1-1.0))**ma 3tZ]4ms}  
      kappah_z_PML_1(k) = 1.0+(kappa_z_max-1.0)*   & tavpq.0O  
                            ((nzPML_1 - k - 0.5) / (nzPML_1 - 1.0))**m 3< 2}V  
      bh_z_1(k) = EXP(-(sigh_z_PML_1(k) / kappah_z_PML_1(k) +   & \kU &^Hi  
                                 alphah_z_PML_1(k))*dt/epsO) VrxQc qPr`  
      ch_z_1(k) = sigh_z_PML_1(k)*(bh_z_1(k)-1.0)/      & :[hgxJu+  
                  (sigh_z_PML_1(k)+kappah_z_PML_1(k)*alphah_z_PML_1(k)) & {6_|/KE9_  
                  / kappah_z_PML_1(k) C'o64+W^  
   ENDDO 4'!c*@Y  
   DO k = 1,nzPML_2 Te<}*qvD  
      sige_z_PML_2(k) = sig_z_max * ( (nzPML_2 - k ) / (nzPML_2 - 1.0) )**m q@nP}Pv&5  
      alphae_z_PML_2(k) = alpha_z_max*((k-1)/(nzPML_2-1.0))**ma }Du}c3  
      kappae_z_PML_2(k) = 1.0+(kappa_z_max-1.0)*   & ?6fnpGX@a  
                                 ((nzPML_2 - k) / (nzPML_2 - 1.0))**m 8>WC5%f*  
      be_z_2(k) = EXP(-(sige_z_PML_2(k) / kappae_z_PML_2(k) +   & w[IE
 
                                 alphae_z_PML_2(k))*dt/epsO) c0;t4( &8  
      if ((sige_z_PML_2(k) == 0.0) .and.        & ZI
3Nq  
         (alphae_z_PML_2(k) == 0.0) .and. (k == nzPML_2)) then \4-"L>  
         ce_z_2(k) = 0.0 =Q+= f  
      else %\H|B0  
         ce_z_2(k) = sige_z_PML_2(k)*(be_z_2(k)-1.0)/       & (}EB2V9Hh  
            (sige_z_PML_2(k)+kappae_z_PML_2(k)*alphae_z_PML_2(k)) & #py[  
            / kappae_z_PML_2(k) XbD4:i%  
      endif t\%HX.8[;%  
   ENDDO 'Pn3%&O$  
   DO k = 1,nzPML_2-1 R}4So1  
      sigh_z_PML_2(k) = sig_z_max * ( (nzPML_2 - k - 0.5)/(nzPML_2-1.0))**m uFPF!Ern  
      alphah_z_PML_2(k) = alpha_z_max*((k-0.5)/(nzPML_2-1.0))**ma q gLaa  
      kappah_z_PML_2(k) = 1.0+(kappa_z_max-1.0)*   & DW_1,:,?7l  
                            ((nzPML_2 - k - 0.5) / (nzPML_2 - 1.0))**m p%*%n3bw  
      bh_z_2(k) = EXP(-(sigh_z_PML_2(k) / kappah_z_PML_2(k) +   & JAQb{KefdO  
                                 alphah_z_PML_2(k))*dt/epsO) jO1r)hw N>  
      ch_z_2(k) = sigh_z_PML_2(k)*(bh_z_2(k)-1.0)/      & qT`k*i?  
                  (sigh_z_PML_2(k)+kappah_z_PML_2(k)*alphah_z_PML_2(k)) & nysUZB  
                  / kappah_z_PML_2(k) n=c 2Kc  
   ENDDO O"c;|zCc>  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .MJofE;Jn  
!  FILL IN UPDATING COEFFICIENTS 06N}k<10O  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ a6WI170^1  
   DA = 1.0 1>E<8&2[L  
   DB = (dt/(muO)) v,4pp@8rv  
   DO i = 1,Imax /$'AjIg4:&  
      DO j = 1,Jmax GJtZ&H  
         DO k = 1,Kmax t 0|
!(3  
            CA(i,j,k) = (1.0 - sig(i,j,k)*dt / (2.0*eps(i,j,k))) /        & R)RG[F#  
               (1.0 + sig(i,j,k) * dt / (2.0*eps(i,j,k))) s1h/}  
            CB(i,j,k) = (dt/(eps(i,j,k))) /                            & n3D;"
a3  
               (1.0 + sig(i,j,k)*dt / (2.0*eps(i,j,k))) [ji#U s:h  
            ENDDO b{]z wpf  
 ENDDO IbP#_Vt  
   ENDDO @]-jl}:]  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ F=a
<~EpZ  
!  FILL IN DENOMINATORS FOR FIELD UPDATES .QN>z-YA6:  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \0vr>C  
   ii = nxPML_2-1 fDvl/|62{  
   DO i = 1,Imax-1 )Ccq4i  
      if (i <= nxPML_1-1) then ),ma_{$N  
         den_hx(i) = 1.0/(kappah_x_PML_1(i)*dx) w &(|e <  
      elseif (i >= Imax+1-nxPML_2) then N=QfP  
         den_hx(i) = 1.0/(kappah_x_PML_2(ii)*dx) Y!gCMLL  
         ii = ii-1 qpE&go=k'  
      else 0^RXGN  
         den_hx(i) = 1.0/dx fH>NJK;  
      endif A#95&kJpy  
   ENDDO {yB&xj[z  
   jj = nyPML_2-1 <\}KT*Xp  
   DO j = 1,Jmax-1 #R.-KUW:  
      if (j <= nyPML_1-1) then % <qw  
         den_hy(j) = 1.0/(kappah_y_PML_1(j)*dy) HLZ;8/|48m  
      elseif (j >= Jmax+1-nyPML_2) then +L9Eqll  
         den_hy(j) = 1.0/(kappah_y_PML_2(jj)*dy) 7U2J xE  
         jj = jj-1 Gq }U|Z  
      else */|9= $54  
         den_hy(j) = 1.0/dy .b6VQCS~9  
      endif #EsNeBu  
   ENDDO }`,t$NV`  
   kk = nzPML_2-1 !=:$lzS^  
   DO k = 2,Kmax-1 kAC&S!n  
      if (k <= nzPML_1) then >w^YO25q  
         den_hz(k) = 1.0/(kappah_z_PML_1(k-1)*dz) Ndcg/d  
      elseif (k >= Kmax+1-nzPML_2) then B3 5E8/  
         den_hz(k) = 1.0/(kappah_z_PML_2(kk)*dz) ^]a#7/]o  
         kk = kk - 1 %]8qAtV^3j  
      else [3--(#R\}?  
         den_hz(k) = 1.0/dz jhg;%+KB  
      endif R]btAu;Z  
   ENDDO Wu_kx2h  
   ii = nxPML_2 3YFU*f,  
   DO i = 1,Imax-1 R5 9S@MsuD  
      if (i <= nxPML_1) then !qN||mCH  
         den_ex(i) = 1.0/(kappae_x_PML_1(i)*dx) 2VJR$Pao  
      elseif (i >= Imax+1-nxPML_2) then iMP]W_  
         den_ex(i) = 1.0/(kappae_x_PML_2(ii)*dx) J1:1B,^y  
         ii = ii-1 ktdW`R\+  
      else <c,u3cp  
         den_ex(i) = 1.0/dx ~ y;6W0x  
      endif vbD{N3p)?n  
   ENDDO ~;&m*2 |V  
   jj = nyPML_2 /O^RF}  
   DO j = 1,Jmax-1 \DD0s8  
      if (j <= nyPML_1) then (C@~3!AVa  
         den_ey(j) = 1.0/(kappae_y_PML_1(j)*dy) qycI(5S,  
      elseif (j >= Jmax+1-nyPML_2) then i2&ed_h<?  
         den_ey(j) = 1.0/(kappae_y_PML_2(jj)*dy) 5_z33,q2  
         jj = jj-1 MvWaB  
      else pD01,5/  
         den_ey(j) = 1.0/dy }=8B*  
      endif hijgF@  
   ENDDO p-/|mL  
   kk = nzPML_2 ?}a;}Q6  
   DO k = 1,Kmax-1 frT]5?{  
      if (k <= nzPML_1) then `/ q|@B7  
         den_ez(k) = 1.0/(kappae_z_PML_1(k)*dz) mHBnC&-/  
      elseif (k >= Kmax-nzPML_2) then K?zH35f$  
         den_ez(k) = 1.0/(kappae_z_PML_2(kk)*dz) K/3)g9Z&io  
         kk = kk - 1 y1bbILWej  
      else yzLpK;  
         den_ez(k) = 1.0/dz :
9_L6  
      endif j"|=C$Kn/  
   ENDDO /i-J&*6_  
HsHB!mQV  
!.:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. t7?Zxq  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ HOE_S!N  
!  BEGIN TIME STEP H+ 0$tHi  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Y>EzTV  
  write(*,*)"begin time-stepping" 0 rXxRQ  
  DO n = 1,nmax |<.lW  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ J~1r{5V4{  
!  UPDATE Hx ZOQTINf  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ F>-B3x  
   DO k = 2,Kmax-1 QZ:8+[oy  
      DO i = 1,Imax-1 NX/;+{  
    DO j = 1,Jmax-1 yLdVd P  
       Hx(i,j,k) = DA * Hx(i,j,k) + DB *       & r;Gi+Ca5  
   ( (Ez(i,j,k) - Ez(i,j+1,k))*den_hy(j)  +    & f/Hm{<BY  
     (Ey(i,j,k) - Ey(i,j,k-1))*den_hz(k) ) %JHGiCv|  
    ENDDO sh`s/JRf  
 ENDDO d\61;C  
      DO i = 1,Imax-1 X'`~s}vGO  
!..................................................................... 6u`)QUmItg  
!  PML for bottom Hx, j-direction :SGF45>B@  
!..................................................................... qLT>Mz)$%  
         DO j = 1,nyPML_1-1 "Za>ZRR  
         psi_Hxy_1(i,j,k) = bh_y_1(j)*psi_Hxy_1(i,j,k)                 & {jho&Ai  
     + ch_y_1(j) *(Ez(i,j,k) - Ez(i,j+1,k))/dy !>e5z|1   
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxy_1(i,j,k) t$xY #:  
         ENDDO G1"zElug  
!..................................................................... 0\@dYPa&C  
!  PML for top Hx, j-direction i]|Yg$  
!..................................................................... weE/TW\e  
         jj = nyPML_2-1 r]l!WRn  
         DO j = Jmax+1-nyPML_2,Jmax-1 Ar:*oiU  
         psi_Hxy_2(i,jj,k) = bh_y_2(jj)*psi_Hxy_2(i,jj,k)       & #&m0WI1  
     + ch_y_2(jj) *(Ez(i,j,k) -       & Rx);7j/5  
                                Ez(i,(j+1),k))/dy $n+w$CI)  
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxy_2(i,jj,k) 2GWMlI  
            jj = jj-1 \CXQo4P  
         ENDDO !\"5rNy  
 ENDDO QzzW x2  
   ENDDO {FIXc^m'  
   DO i = 1,Imax-1 %QKRFPYhS  
      DO j = 1,Jmax-1 k-HCeZ  
!..................................................................... 2cqI[t@0  
!  PML for bottom Hx, k-direction { Q`QX`#  
!..................................................................... nM-h
&na{s  
         DO k = 2,nzPML_1 rZfN+S,g  
         psi_Hxz_1(i,j,k-1) = bh_z_1(k-1)*psi_Hxz_1(i,j,k-1)            & L2pp6bW  
     + ch_z_1(k-1) *(Ey(i,j,k) - Ey(i,j,k-1))/dz j|/]#@
Yr  
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxz_1(i,j,k-1) <X7FMNr[  
         ENDDO }2uI?i8  
!..................................................................... S`"LV $8  
!  PML for top Hx, k-direction %9M~f*  
!..................................................................... @*Wh  
         kk = nzPML_2-1 T:IW%?M  
         DO k = Kmax+1-nzPML_2,Kmax-1 H*Kj3NgY  
         psi_Hxz_2(i,j,kk) = bh_z_2(kk)*psi_Hxz_2(i,j,kk)             & YR>B_,Gl  
     + ch_z_2(kk) *(Ey(i,j,k) -       & E !a|Xp  
                                Ey(i,j,k-1))/dz \yd s5g!:  
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxz_2(i,j,kk) f~dd3m('  
           kk = kk-1 )YX 'N<[  
         ENDDO YE=q:Bv  
 ENDDO ZP;j9T!  
   ENDDO %ix)8+Eb  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (A?H1 9  
!  UPDATE Hy }*ZHgf]~#  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ /u0' 6V  
   DO k = 2,Kmax-1 3v mjCm  
      DO i = 1,Imax-1 =7C%P%yt  
    DO j = 1,Jmax-1 E7CH^]x  
            Hy(i,j,k) = DA * Hy(i,j,k) + DB *             & :bWUuXVtJ  
   ( (Ez(i+1,j,k) - Ez(i,j,k))*den_hx(i) +      & q@@T]V6  
     (Ex(i,j,k-1) - Ex(i,j,k))*den_hz(k) ) (q4),y<:[  
    ENDDO VGceD$<  
      ENDDO nD)K}4  
      DO j = 1,Jmax-1 8g:;)u4$P  
!..................................................................... ltgc:&=|@  
!  PML for bottom Hy, i-direction 8 XICF  
!..................................................................... dd\n8f  
         DO i = 1,nxPML_1-1 xZQg'IT  
       psi_Hyx_1(i,j,k) = bh_x_1(i)*psi_Hyx_1(i,j,k)                   & 1 29q`u;  
    + ch_x_1(i)*(Ez(i+1,j,k) - Ez(i,j,k))/dx YQpSlCCo 3  
       Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyx_1(i,j,k) cu%
C"  
         ENDDO <%Ostqj  
!..................................................................... m^H21P"z  
!  PML for top Hy, i-direction tX@_fYb  
!..................................................................... j8fpj{hp  
         ii = nxPML_2-1 d8m6B6 CW  
         DO i = Imax+1-nxPML_2,Imax-1 ) :\xHR4  
       psi_Hyx_2(ii,j,k) = bh_x_2(ii)*psi_Hyx_2(ii,j,k)     & (d<4"!  
    + ch_x_2(ii)*(Ez(i+1,j,k) -    & "bm  
                                Ez(i,j,k))/dx j~+>o[c  
       Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyx_2(ii,j,k) #ssN027  
            ii = ii-1 UQ+?\wi*  
    ENDDO ;Y;
qg  
     ENDDO >j'ZPwj^  
   ENDDO 8:#rA*Y  
   DO i = 1,Imax-1 :k3Nt5t!  
      DO j = 1,Jmax-1 qOUqs'7/]  
!..................................................................... .9"Y_/0   
!  PML for bottom Hy, k-direction g+ MdHn[  
!..................................................................... ]]&M@FM2z  
         DO k = 2,nzPML_1 ,WG<hgg-U)  
       psi_Hyz_1(i,j,k-1) = bh_z_1(k-1)*psi_Hyz_1(i,j,k-1)          & ,daKC  
    + ch_z_1(k-1)*(Ex(i,j,k-1) - Ex(i,j,k))/dz @EDs~ lPv  
       Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyz_1(i,j,k-1) raCi 8  
         ENDDO >zhO7,=,  
!..................................................................... _x?
uU  
!  PML for top Hy, k-direction 66'?&Xx'  
!..................................................................... `8I&7c  
         kk = nzPML_2-1 BWM YpZom  
         DO k = Kmax+1-nzPML_2,Kmax-1 *0|IXGr  
     psi_Hyz_2(i,j,kk) = bh_z_2(kk)*psi_Hyz_2(i,j,kk)               & :sP!p`dl  
    + ch_z_2(kk)*(Ex(i,j,k-1) -    & Uv?|G%cD-  
                                Ex(i,j,k))/dz :LQ5u[g$\  
     Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyz_2(i,j,kk) K,e w>U  
            kk = kk-1 Wa1, p  
         ENDDO e+j)~RBnu3  
     ENDDO ecgtUb8K  
   ENDDO lr{?"tl_  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ _^
'I  
!  UPDATE Hz ]\xy\\b/`  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '2laTl]`  
   DO k = 1,Kmax-1 S3@|Q\*r  
      DO i = 1,Imax-1 ph_4q@  
        DO j = 1,Jmax-1 gCN$}  
            Hz(i,j,k) = DA * Hz(i,j,k) + DB       & KrXdnY8  
                  * ((Ey(i,j,k) - Ey(i+1,j,k))*den_hx(i) +        & Cz W:L&t  
       (Ex(i,j+1,k) - Ex(i,j,k))*den_hy(j)) ]ZoD'-,  
    ENDDO w1Ec_y{  
      ENDDO VB?Ohk]<  
      DO j = 1,Jmax-1 DF g,Xa#  
!..................................................................... X5`AGyX  
!  PML for bottom Hz, x-direction l%`F&8K  
!..................................................................... !Yw3 d   
         DO i = 1,nxPML_1-1 +Ag!?T  
          psi_Hzx_1(i,j,k) = bh_x_1(i)*psi_Hzx_1(i,j,k)                 & 9&2Vm;F_  
     + ch_x_1(i) *(Ey(i,j,k) - Ey(i+1,j,k))/dx Udc=,yo3Qm  
       Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzx_1(i,j,k) r)1'eP
I"  
         ENDDO z9p05NFH  
!..................................................................... PmR~c,  
!  PML for top Hz, x-direction X5khCLHi  
!..................................................................... c= x,ijY "  
         ii = nxPML_2-1 oh KCdT~  
         DO i = Imax+1-nxPML_2,Imax-1 |0\0a&tkPl  
          psi_Hzx_2(ii,j,k) = bh_x_2(ii)*psi_Hzx_2(ii,j,k)            & (lGaPMEU}  
     + ch_x_2(ii) *(Ey(i,j,k) -       & :Mcu  
                                Ey(i+1,j,k))/dx {XWZ<OjG  
       Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzx_2(ii,j,k) aj"M>zd*}  
            ii = ii-1 !j\" w p  
    ENDDO E^rbcGJ  
      ENDDO js8{]04y  
      DO i = 1,Imax-1 <N`J`J-[  
!..................................................................... ^<:sdv>Y5  
!  PML for bottom Hz, y-direction 0z[dlHi  
!.....................................................................
XgxX.`H7  
         DO j = 1,nyPML_1-1 3"kdjOB  
            psi_Hzy_1(i,j,k) = bh_y_1(j)*psi_Hzy_1(i,j,k)                   & ]/!<PF  
    + ch_y_1(j)*(Ex(i,j+1,k) - Ex(i,j,k))/dy Pf$pt  
       Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzy_1(i,j,k)
1$2D O  
         ENDDO sU"}-de  
!.....................................................................
=;3fq-  
!  PML for top Hz, y-direction ncX/L[L  
!..................................................................... 5hbJOo0BZ  
         jj = nyPML_2-1 JGFt0He]  
         DO j = Jmax+1-nyPML_2,Jmax-1 
Z1h]  
            psi_Hzy_2(i,jj,k) = bh_y_2(jj)*psi_Hzy_2(i,jj,k)               & uu/+.9  
    + ch_y_2(jj)*(Ex(i,j+1,k) -    & +MB!B9M@  
                                Ex(i,j,k))/dy #!9S}b$  
       Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzy_2(i,jj,k) Kv@eI$t5  
            jj = jj-1 Bcjx>#3?L  
    ENDDO 7-~Q5Kr.  
      ENDDO dDD<E?TjD  
   ENDDO #9
m$ N  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ N~?{UOZd  
!  UPDATE Ex %',F  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Qf|U0  
   DO k = 1,Kmax-1 cDoo
*  
      DO i = 1,Imax-1 eVR5Xar  
    DO j = 2,Jmax-1 Maqf[ Vky  
              IF (i >= istart-1 .and. i <= iend .and. j >= jstart .and.  & k`{@pt.  
                   j <= jend .and. k >= kstart .and. k <= kend) THEN 7nVRn9Hn  
                 Ex(i,j,k) = 0.0 0}:2Q#  
              ELSE {66fG53x  
          Ex(i,j,k) = CA(i,j,k) * Ex(i,j,k) + CB(i,j,k) *       & ~+H"
-+  
     ( (Hz(i,j,k) - Hz(i,j-1,k))*den_ey(j)  +    & N,J9Wu ZJ\  
     (Hy(i,j,k) - Hy(i,j,k+1))*den_ez(k) ) "iM~Hy  
              ENDIF QeDQo  
    ENDDO a2f^x@0k  
 ENDDO ~fa(=.h  
      DO i = 1,Imax-1 p:OPwD+  
!..................................................................... P1U*g!  
!  PML for bottom Ex, j-direction mrJQ#  
!..................................................................... i`0v#P  
         DO j = 2,nyPML_1 ooj~&fu  
         psi_Exy_1(i,j,k) = be_y_1(j)*psi_Exy_1(i,j,k)                 & y?)}8T^  
     + ce_y_1(j) *(Hz(i,j,k) - Hz(i,j-1,k))/dy 4#(ZNP  
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exy_1(i,j,k) H Y ynMP  
         ENDDO O Lc}_  
!..................................................................... JI/_ce  
!  PML for top Ex, j-direction yE\wj  
!..................................................................... #$1og=  
         jj = nyPML_2 e2+BWKaU  
         DO j = Jmax+1-nyPML_2,Jmax-1 $x#FgD(iI  
         psi_Exy_2(i,jj,k) = be_y_2(jj)*psi_Exy_2(i,jj,k)       & 20TCG0%x  
     + ce_y_2(jj) *(Hz(i,j,k) -       & @<@SMK)  
                                Hz(i,(j-1),k))/dy ZO7&vF}  
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exy_2(i,jj,k) }%-t+Tf,  
            jj = jj-1 NJ
MJ  
         ENDDO ,=~z6[  
 ENDDO tj*y)28
-  
   ENDDO 9ktEm|F3  
   DO i = 1,Imax-1 Z Dhx5SL&  
      DO j = 2,Jmax-1 8*SP~q  
!..................................................................... LrCk*@  
!  PML for bottom Ex, k-direction m3^/:<  
!.....................................................................

Gs*G<P"  
         DO k = 1,nzPML_1 ;c
-3g]  
         psi_Exz_1(i,j,k) = be_z_1(k)*psi_Exz_1(i,j,k)                 & m))<!3  
     + ce_z_1(k) *(Hy(i,j,k) - Hy(i,j,k+1))/dz $2MAZGJV
 
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exz_1(i,j,k) 4<X!<]3]  
         ENDDO ZL< MC~  
!..................................................................... 2T)sXBu  
!  PML for top Ex, k-direction ^#):c`  
!..................................................................... zD)pF1,7:8  
         kk = nzPML_2 K?4FT$9G  
         DO k = Kmax-nzPML_2,Kmax-1 e/8z+H^H  
         psi_Exz_2(i,j,kk) = be_z_2(kk)*psi_Exz_2(i,j,kk)             & OI0B:()  
     + ce_z_2(kk) *(Hy(i,j,k) -       & W><dYy=z5
 
                                Hy(i,j,k+1))/dz 5y} v{Ijt  
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exz_2(i,j,kk) Pz_NDI  
            kk = kk-1 :+%Zh@u\  
         ENDDO }Z`(aDH  
 ENDDO D%7
kBfCb  
   ENDDO @cq`:_.[  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ;XC@=RpX  
!  UPDATE Ey UzKFf&-:;K  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ D\~e&0*  
   DO k = 1,Kmax-1 Ao7`G':  
      DO i = 2,Imax-1 o+?rI p  
    DO j = 1,Jmax-1 vU*x2fVb}  
              IF (i >= istart .and. i <= iend .and. j >= jstart-1 .and. & A K/z6XGy  
                   j <= jend .and. k >= kstart .and. k <= kend) THEN 70B)|<$  
                 Ey(i,j,k) = 0.0 k]rLjcB  
              ELSE Qpt&3_   
                 Ey(i,j,k) = CA(i,j,k) * Ey(i,j,k) + CB(i,j,k) *    & ;U`X 6d  
   ( (Hz(i-1,j,k) - Hz(i,j,k))*den_ex(i) +         & e#/kNHl  
     (Hx(i,j,k+1) - Hx(i,j,k))*den_ez(k) ) J.RAmU<  
              ENDIF l+wc'=]  
    ENDDO [`Cq\mI-W  
      ENDDO aj~bt-cE  
      DO j = 1,Jmax-1 ue8qIZH  
!..................................................................... %gcc y|  
!  PML for bottom Ey, i-direction ~/IexQB&  
!..................................................................... p
\bFdxv#  
         DO i = 2,nxPML_1 33~qgK1>  
       psi_Eyx_1(i,j,k) = be_x_1(i)*psi_Eyx_1(i,j,k)                   & .1QgK  
    + ce_x_1(i)*(Hz(i-1,j,k) - Hz(i,j,k))/dx 1j\aH&)GH  
       Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyx_1(i,j,k) [<wbbvXR  
         ENDDO 6<5Jq\-h  
!..................................................................... R[1BfZ6s  
!  PML for top Ey, i-direction ~"!a9GZ  
!..................................................................... ~SjZk|  
         ii = nxPML_2 eX2<}'W<  
         DO i = Imax+1-nxPML_2,Imax-1 \;]kYO}  
       psi_Eyx_2(ii,j,k) = be_x_2(ii)*psi_Eyx_2(ii,j,k)              & 15zrrU~D  
    + ce_x_2(ii)*(Hz(i-1,j,k) -    & }Uf<ZXW  
                               Hz(i,j,k))/dx ]$M<]w,IJ2  
       Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyx_2(ii,j,k) C2<CWPn<  
            ii = ii-1 *OdX u&5  
    ENDDO 0+b0<  
     ENDDO R:aYL~  
   ENDDO }wSy  
   DO i = 2,Imax-1 0m+8P$)C%  
      DO j = 1,Jmax-1 q9dLHi<1  
!..................................................................... n~Szf  
!  PML for bottom Ey, k-direction |z^pL1Z]5  
!..................................................................... \h3e-)  
         DO k = 1,nzPML_1 K`nJVc  
       psi_Eyz_1(i,j,k) = be_z_1(k)*psi_Eyz_1(i,j,k)        .. ACV ek  
P`'$  

未注册仅能浏览部分内容，查看全部内容及附件请先 登录 或 注册

不错，可以很好的参考

谢谢，正在学习PML层的东西

下载看看吧，已经编译通过。 ZXh~79  
Code只能学习，不要拿来用（为了提高自己的能力）。 bSY;[{Kl  
实际上，只有自己编写Code,才会灵活运用，发现PML的问题和缺陷，并创新，然后发论文。 n'ehB%"  
我自己编写了许多PML的code，以此发表了许多SCI检索的论文。 [h^f%  
5pB^Y MP  
请问楼主，从哪搞到的这个code？

这么厉害。写PML的code，能发很多SCI的文章。 g&{gD^9)4  
.Ln;m8  
介绍一下经验吧。 NQDLI 1o  
kX@bv"i  

引用第3楼ljx690529于2008-12-08 20:06发表的 : \d `dV0X  
下载看看吧，已经编译通过。 \_lG#p|  
Code只能学习，不要拿来用（为了提高自己的能力）。 l.XknF  
实际上，只有自己编写Code,才会灵活运用，发现PML的问题和缺陷，并创新，然后发论文。 "}0)YRz%  
我自己编写了许多PML的code，以此发表了许多SCI检索的论文。
.......

引用第3楼ljx690529于2008-12-08 20:06发表的 : ybD{4&ZE  
下载看看吧，已经编译通过。 l4iuu  
Code只能学习，不要拿来用（为了提高自己的能力）。 WpXODkQL  
实际上，只有自己编写Code,才会灵活运用，发现PML的问题和缺陷，并创新，然后发论文。 -s9Y(>  
我自己编写了许多PML的code，以此发表了许多SCI检索的论文。
.......

jf$JaY  
呵呵，别人给传的，我想拿出来大家共享一下，能否加你qq好友，共同交流一下

请公布qq号码，我也要和你交流

我只是把这个程序转贴而已,我用M写程序的

谢谢分享！

REAL,DIMENSION(Imax-1, Jmax, Kmax-1)  ::   & 7E75s)KH  
      Ex
.Qi1I  
这句话怎么这么写呢？我觉得应该是REAL,DIMENSION（Imax-1，Jmax，Kmax）啊！能否解释一下？