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

!*********************************************************************** 3Ett9fBd  
!     3-D FDTD code with CPML absorbing boundary conditions EwOi` g  
!*********************************************************************** >iWw i'T=  
! /A%31WE&1  
!     Program author: Jamesina J. Simpson, Ph.D. Graduate Student PZCOJK
 
!                     National Science Foundation Fellow { Em fw9L  
!                     NU Computational Electromagnetics Laboratory /q6 ^.>b  
!                     Director: Allen Taflove tt]ZGn*  
!                     Department of Electrical Engineering and Y:K1v:Knw  
!                                                      Computer Science |z.Z='`  
!                     Northwestern University f`,isy[  
!                     2145 Sheridan Road P%kJq^&  
!                     Evanston, IL 60208 l Gy`{E|  
!                     [url=mailto:j-simpson@northwestern.edu]j-simpson@northwestern.edu[/url] 1}I%yOi)  
! P&8QKX3 j^  
!     Copyright 2005 7@9R^,M4:  
! ]1i1_AR'`  
!     This FORTRAN 90 code implements the finite-difference time-domain gn4Sz")  
!     solution of Maxwell's curl equations over a three-dimensional w 8oIq*  
!     Cartesian space lattice.  The grid is terminated by CPML absorbing @:G#[>nKe  
!     boundary conditions.  However, in a straightforward manner it can |gJI}"T  
!     be altered to have PEC planes at any of the outer grid boundaries.  $N17GqoC  
!     Also, the number of grid cells as well as the thickness of EN@
Pr `R  
!     the PML in each Cartesian direction can be varied independently. +Zi@+|"BCN  
! sQr |3}I(  
!     For illustrative purposes, the code supplied here models the PML .s>PDzM$  
!     numerical experiment described in Section 7.11.2.  Three output eTY""EWU  
!     files are generated having the following data: /m+\oZ ]d  
! A>Oi9%OY:  
!       (1) Ez at the source for each time step; ZHOh(  
!       (2) Ey at the probe point for each time step;  and >=3ay^(Y2D  
!       (3) The plane of Ez values 1 mm above the source in the k- UhF+},gU  
!           direction recorded at the time step, "record_grid".  This sT%^W  
!           data can be viewed in Matlab using the following commands: )|i]"8I  
! :~dI2e\:  
!             >>  load view_grid.dat P{h$> 6c  
!             >>  for j = 1:Jmax W*NK-F[  
!             >>      image(1:Imax,j) = view_grid(1+(j-1)*Imax:j*Imax); t(.xEl;Ma  
!             >>  end ]$ iqJL  
!             >>  pcolor(image'); shading flat Ni)/L( &  
! R"k}wRnxY  
!     The relative error (equation 7.135) graphed in Fig. 7.6 on 4l"oq"uc  
!     page 318 can be reproduced by comparing the output file Q.])En >i  
!     "probe.dat" with that of a much larger benchmark grid having ?Y#x`DMh  
!     Imax, Jmax, and Kmax increased to the values mentioned in the text. C\C*'l6d  
! Qo \;)  
!     This code has been tested using the Intel Fortran Compiler 8.0 for I'T@}{h  
!     Linux.  The executable can be created by typing
%oN5jt  
!     "ifort fdtd3D_CPML.f90" at the command prompt. LuUfdzH  
! q.g0Oz@z  
!     This program is not guaranteed to be free of defects or bugs. 9..! g:  
!     Please report any bugs that may exist to: 9~C$C  
!                                          [url=mailto:j-simpson@northwestern.edu]j-simpson@northwestern.edu[/url] [%yCnt  
! t"L-9kCM  
! A@k`$xevVj  
!*********************************************************************** Nh/B8:035  
~OQ/ |ws  
   PROGRAM fdtd3D_CPML o+.LG($+U  
   IMPLICIT NONE n*GsM6Y&  
!  .................................. m%G:|`f7  
!  Input Fundamental Constants (MKS units) IA?v[xu  
   REAL, PARAMETER ::                            & K trR+:  
      pi = 3.14159265358979, C = 2.99792458E8, & j5qrM_Chg  
      muO = 4.0 * pi * 1.0E-7, epsO = 1.0/(C*C*muO) -H(\[{3{V  
!  .................................. ]W>kbHImz  
!  Specify Material Relative Permittivity and Conductivity p$OkWSi~
 
   REAL, PARAMETER::                      & >`!Lh`n7_  
      epsR = 1.0, sigM1 = 0.0   ! free space lQnqPQY  
!  .................................. 2K&5Kt/  
!  Specify Number of Time Steps and Grid Size Parameters Dwp-*QK^G  
   INTEGER, PARAMETER ::                                     & j@+QwZL|  
      nmax = 2100, &  ! total number of time steps eZ[O:Wvk:  
      Imax = 51, Jmax = 126, Kmax = 26  ! grid size corresponding to the mtg3}etA  
                                        ! number of Ez field components @ wJ|vW_.  
!  .................................. #99fFs`w  
!  Specify Grid Cell Size in Each Direction and Calculate the 8<Yqpb  
!  Resulting Courant-Stable Time Step ~^pV>>LX|  
   REAL, PARAMETER ::                                        & q+/7v9  
      dx = 1.0E-3, dy = 1.0E-3, dz = 1.0E-3,  &  ! cell size in each direction auV<=1<zJ  
      dt = 0.99 / (C*(1.0/dx**2+1.0/dy**2+1.0/dz**2)**0.5) ;/]vmgl2  
                                                 ! time step increment g4 G?hv`R  
!  .................................. \z>L,U  
!  Specify the Impulsive Source (See Equation 7.134) f)/5%W7n}  
   REAL, PARAMETER ::                                        & n%Gk {h5  
      tw = 53.0E-12, tO = 4.0*tw  exQU  
!  .................................. 5t&;>-A'?'  
!  Specify the Time Step at which the Grid is Recorded for Visualization !XJS"owr  
   INTEGER, PARAMETER ::                                        & c$[cDf~  
      record_grid = 300 zw,=mpf3_  
!  ..................................  W.t`  
!  Specify the PEC Plate Boundaries and the Source/Recording Points 3BWYSJ|  
   INTEGER, PARAMETER ::                                    & XrD@q  
      istart = (Imax-1)/2-11, iend = istart+24, jstart = istart,   & )"Vd8*e  
      jend = jstart + 99, kstart = Kmax/2, kend = kstart,      & DU>#eR0G  
      isource = istart, jsource = jstart, ksource = kend+1,     & XAf,k&f3  
      irecv1 = iend, jrecv1 = jend+1, krecv1 = kend  ! Ey at probe point \9GJa"xA`  
!  .................................. BVk&TGa;[$  
!  Specify the CPML Thickness in Each Direction (Value of Zero lYt|C^  
!  Corresponds to No PML, and the Grid is Terminated with a PEC) )cd5iE:FO  
   INTEGER, PARAMETER ::                        & %RF9R"t$  
      ! PML thickness in each direction 1:r
8p6  
      nxPML_1 = 11, nxPML_2 = nxPML_1, nyPML_1 = nxPML_1,      & eL7\})!W  
      nyPML_2 = nxPML_1, nzPML_1 = nxPML_1, nzPML_2 = nxPML_2 Q(J6;s#b  
!  .................................. %qf ?_2v  
!  Specify the CPML Order and Other Parameters Y &"rf   
   INTEGER, PARAMETER ::                        & l;C00ZBOc  
      m = 3, ma = 1  !xz0zT.  
   REAL, PARAMETER  ::                     & ic_q<Y}  
      sig_x_max = 0.75 * (0.8*(m+1)/(dx*(muO/epsO*epsR)**0.5)),   & ndU<,{r  
      sig_y_max = 0.75 * (0.8*(m+1)/(dy*(muO/epsO*epsR)**0.5)),   & bF X0UE>  
      sig_z_max = 0.75 * (0.8*(m+1)/(dz*(muO/epsO*epsR)**0.5)),   & ^dF?MQA<@  
      alpha_x_max = 0.24,   & 20XN5dTFT  
      alpha_y_max = alpha_x_max, alpha_z_max = alpha_x_max, & P5^<c\Mr,Y  
      kappa_x_max = 15.0, & ;"77?)  
      kappa_y_max = kappa_x_max, kappa_z_max = kappa_x_max W~e/3#R\=  
   INTEGER ::                                                & OL
S.0UEc  
 i,j,ii,jj,k,kk,n p Y[dJxB  
   REAL  ::                                                   & y6*9, CF  
      source, P1, P2 _
:0  
!     TM components `swf~  
   REAL,DIMENSION(Imax, Jmax, Kmax)  ::                      & b^~"4fU  
      Ez, CA, CB, sig, eps N4JJA+  
   REAL,DIMENSION(Imax-1, Jmax, Kmax)  ::                      & c0%.GcF0{  
      Hy Fs,#d%4@%  
   REAL,DIMENSION(Imax,Jmax-1, Kmax)  ::                      & ev_4!+ko  
      Hx Dy[_Ix/Y,  
   REAL  ::                        & $ 3/G)
/A  
      DA, DB +jPs0?}s  
!     TE components vdLBf+Zi  
   REAL,DIMENSION(Imax-1, Jmax-1, Kmax-1)  ::                      & IyfhVk?  
      Hz `x0GT\O2-  
   REAL,DIMENSION(Imax-1, Jmax, Kmax-1)  ::                      & U94Tp A6  
      Ex w`")^KXi  
   REAL,DIMENSION(Imax,Jmax-1, Kmax-1)  ::                      & #M5R>&?Jqz  
      Ey E 9LKVs}  
!  PML 
Nhnw'9  
   REAL ,DIMENSION(nxPML_1,Jmax,Kmax) ::                       & Y"RjMyQh  
      psi_Ezx_1 nA#FGfZ{Ge  
   REAL ,DIMENSION(nxPML_2,Jmax,Kmax) ::                       & ^t'mfG|DV  
      psi_Ezx_2 mDT"%I"4j  
   REAL ,DIMENSION(nxPML_1-1,Jmax,Kmax) ::                       & }c$@0x;YQ  
      psi_Hyx_1 oju}0h'1  
   REAL ,DIMENSION(nxPML_2-1,Jmax,Kmax) ::                       & ~ojH$=K>d  
      psi_Hyx_2 E0Y>2HOuL  
   REAL ,DIMENSION(Imax,nyPML_1,Kmax) ::                       & '4""Gz  
      psi_Ezy_1                               "`jey)&H*M  
   REAL ,DIMENSION(Imax,nyPML_2,Kmax) ::                       & B]o5HA<k  
      psi_Ezy_2 (` N@4w=  
   REAL ,DIMENSION(Imax,nyPML_1-1,Kmax) ::                       & +hJ@w-u,G  
      psi_Hxy_1                               muc>4!Q  
   REAL ,DIMENSION(Imax,nyPML_2-1,Kmax) ::                       & l6wN&JHTh  
      psi_Hxy_2 _RHB ^y;-  
   REAL ,DIMENSION(Imax,Jmax-1,nzPML_1-1) ::                       & cn~M:LW23  
      psi_Hxz_1 ?'dsiA[  
   REAL ,DIMENSION(Imax,Jmax-1,nzPML_2-1) ::                       & Vk>m/"  
      psi_Hxz_2 Hmi]qK[F  
   REAL ,DIMENSION(Imax-1,Jmax,nzPML_1-1) ::                       & dfss_}R  
      psi_Hyz_1 cy6lsJ"?  
   REAL ,DIMENSION(Imax-1,Jmax,nzPML_2-1) ::                       & K~ 6[zJ4  
      psi_Hyz_2 y@'m D*z  
   REAL ,DIMENSION(Imax-1,Jmax,nzPML_1) ::                       & _4]GP3`  
      psi_Exz_1 };z[x2l^  
   REAL ,DIMENSION(Imax-1,Jmax,nzPML_2) ::                       & GV1SKa  
      psi_Exz_2 {xzs{)9|Y4  
   REAL ,DIMENSION(Imax,Jmax-1,nzPML_1) ::                       & CE'd`_;HLn  
      psi_Eyz_1 $ MN1:ih  
   REAL ,DIMENSION(Imax,Jmax-1,nzPML_2) ::                       & #:yh2y7a%  
      psi_Eyz_2 +l "z  
   REAL ,DIMENSION(nxPML_1-1,Jmax-1,Kmax-1) ::                       & dP0%<Q|  
      psi_Hzx_1 kDz!v?Z2+B  
   REAL ,DIMENSION(nxPML_1,Jmax-1,Kmax-1) ::                       & }?0At<(d  
      psi_Eyx_1 ,'E+f%  
   REAL ,DIMENSION(nxPML_2-1,Jmax-1,Kmax-1) ::                       & l^B PTg)X@  
      psi_Hzx_2 sKvz<7pag  
   REAL ,DIMENSION(nxPML_2,Jmax-1,Kmax-1) ::                       & 2tq~NA\#t  
      psi_Eyx_2 j6NK7Li  
   REAL ,DIMENSION(Imax-1,nyPML_1-1,Kmax-1) ::                       & 7vRFF@eq}  
      psi_Hzy_1                               d0Jaa1b~O  
   REAL ,DIMENSION(Imax-1,nyPML_1,Kmax-1) ::                       & ;c!}'2>vM  
      psi_Exy_1                               =68
CR[H  
   REAL ,DIMENSION(Imax-1,nyPML_2-1,Kmax-1) ::                       & rk$$gXg9/  
      psi_Hzy_2 <]jKpJ{3N  
   REAL ,DIMENSION(Imax-1,nyPML_2,Kmax-1) ::                       & |NsrO8H   
      psi_Exy_2 {b<;?Dus^  
   REAL ,DIMENSION(nxPML_1) ::                       & q (?%$u.  
      be_x_1, ce_x_1, alphae_x_PML_1, sige_x_PML_1, kappae_x_PML_1 W.fsW<{4j  
   REAL ,DIMENSION(nxPML_1-1) ::                       & ^HYrJr$y  
      bh_x_1, ch_x_1, alphah_x_PML_1, sigh_x_PML_1, kappah_x_PML_1 h*2NFL~#  
   REAL ,DIMENSION(nxPML_2) ::                       & gg-4ce/  
      be_x_2, ce_x_2, alphae_x_PML_2, sige_x_PML_2, kappae_x_PML_2 +~6Nq(kV  
   REAL ,DIMENSION(nxPML_2-1) ::                       & BO3#*J5S\  
      bh_x_2, ch_x_2, alphah_x_PML_2, sigh_x_PML_2, kappah_x_PML_2 <u'q._m  
   REAL ,DIMENSION(nyPML_1) ::                       & eB$S d  
      be_y_1, ce_y_1, alphae_y_PML_1, sige_y_PML_1, kappae_y_PML_1 V@QWJZ"  
   REAL ,DIMENSION(nyPML_1-1) ::                       & Aw38Tw  
      bh_y_1, ch_y_1, alphah_y_PML_1, sigh_y_PML_1, kappah_y_PML_1 bP4<q?FKcN  
   REAL ,DIMENSION(nyPML_2) ::                       & xzqgem`[\  
      be_y_2, ce_y_2, alphae_y_PML_2, sige_y_PML_2, kappae_y_PML_2 O$2= Z  
   REAL ,DIMENSION(nyPML_2-1) ::                       & %SHjJCS3  
      bh_y_2, ch_y_2, alphah_y_PML_2, sigh_y_PML_2, kappah_y_PML_2 )9.i'{{ 0  
   REAL ,DIMENSION(nzPML_1) ::                       & <)vjoRv  
      be_z_1, ce_z_1, alphae_z_PML_1, sige_z_PML_1, kappae_z_PML_1 jI-\~  
   REAL ,DIMENSION(nzPML_1-1) ::                       & ?G48GxJ  
      bh_z_1, ch_z_1, alphah_z_PML_1, sigh_z_PML_1, kappah_z_PML_1 23F<f+2S  
   REAL ,DIMENSION(nzPML_2) ::                       & u0)~Im,X  
      be_z_2, ce_z_2, alphae_z_PML_2, sige_z_PML_2, kappae_z_PML_2 6WN1DW  
   REAL ,DIMENSION(nzPML_2-1) ::                       & `pHlGbrW  
      bh_z_2, ch_z_2, alphah_z_PML_2, sigh_z_PML_2, kappah_z_PML_2 ] X9
e|  
!     denominators for the update equations &K1\"  
   REAL,DIMENSION(Imax-1)  ::                      & mkR1iY  
      den_ex, den_hx QL<uQ`>(  
   REAL,DIMENSION(Jmax-1)  ::                      & 4!%TY4bJ  
      den_ey, den_hy kFJ sB,2-  
   REAL,DIMENSION(Kmax-1)  ::                      & u;[*Z  
      den_ez, den_hz -W^jmwM   
zPC&p{S>  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ : " ([i"  
!  OPEN OUTPUT FILES Vh.9/$xQ  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Guc~]
B  
   OPEN (UNIT = 30, FILE = "source.dat") %(
c5T
)B9  
   OPEN (UNIT = 31, FILE = "probe.dat") )N8[@  
   OPEN (UNIT = 33, FILE = "view_grid.dat") [%8t~zg  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ LV
^V`m0#  
!  INITIALIZE VARIABLES xS4
B"/  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ L7_qs+  
   P1 = 0.0 vbDSNm#Yv  
   P2 = 0.0 yK$.wd2,  
   Ez(:,:,:) = 0.0 px!
TRbf  
   Hz(:,:,:) = 0.0 :|GC~JElo5  
   Ex(:,:,:) = 0.0 ~F</s.  
   Hx(:,:,:) = 0.0 D2U")g}U  
   Ey(:,:,:) = 0.0 `eF&|3!IYQ  
   Hy(:,:,:) = 0.0 TTJFF\$?  
   sig(:,:,:) = sigM1 y?t2@f]!XK  
   eps(:,:,:) = epsR*epsO VUo7Evc:.P  
   psi_Exy_1(:,:,:) = 0.0 RY>BP[h  
   psi_Exy_2(:,:,:) = 0.0 %|6t\[gn  
   psi_Exz_1(:,:,:) = 0.0 \Kx@?,  
   psi_Exz_2(:,:,:) = 0.0 YkPc&&#  
   psi_Eyx_1(:,:,:) = 0.0 tP?pN]Q$,  
   psi_Eyx_2(:,:,:) = 0.0 !zhg3B#p  
   psi_Eyz_1(:,:,:) = 0.0 E5#Dn.!~  
   psi_Eyz_2(:,:,:) = 0.0 ^* xhbM;  
   psi_Ezy_1(:,:,:) = 0.0 lHV&8fny  
   psi_Ezy_2(:,:,:) = 0.0 "|2|Vju%  
   psi_Ezx_1(:,:,:) = 0.0 xHA6  
   psi_Ezx_2(:,:,:) = 0.0 Zd~l_V f  
   psi_Hxy_1(:,:,:) = 0.0 ;W/K7}  
   psi_Hxy_2(:,:,:) = 0.0 j2UiZLuV  
   psi_Hxz_1(:,:,:) = 0.0 [>B`"nyNQ  
   psi_Hxz_2(:,:,:) = 0.0 1@ e22\  
   psi_Hyx_1(:,:,:) = 0.0 [a\
U8 w  
   psi_Hyx_2(:,:,:) = 0.0 q,ry3Nr4n  
   psi_Hyz_1(:,:,:) = 0.0 ^`W8>czi  
   psi_Hyz_2(:,:,:) = 0.0 BD)5br].  
   psi_Hzy_1(:,:,:) = 0.0 Q: H`TSR]  
   psi_Hzy_2(:,:,:) = 0.0 4QK~qAi  
   psi_Hzx_1(:,:,:) = 0.0 w3,1ImrXp  
   psi_Hzx_2(:,:,:) = 0.0 "Y9PS_u(~  
   write(*,*)"Imax: ", Imax [<HU~PP  
   write(*,*)"Jmax: ", Jmax 4`x.d  
   write(*,*)"Kmax: ", Kmax P,$|.pd'  
   write(*,*)"dt: ", dt ^!={=No]  
   write(*,*)"nmax: ", nmax S(K}.C1x  
   write(*,*)"max time: ", nmax*dt 'u
\my
 
   write(*,*)"record grid after ", record_grid, "dt" *1}UK9X;  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ^yO+-A2zC  
!  SET CPML PARAMETERS IN EACH DIRECTION kl0!*j  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ zL1*w@6  
   DO i = 1,nxPML_1 -P+@n)?T6  
      sige_x_PML_1(i) = sig_x_max * ( (nxPML_1 - i) / (nxPML_1 - 1.0) )**m qdlz#-B  
      alphae_x_PML_1(i) = alpha_x_max*((i-1.0)/(nxPML_1-1.0))**ma .,)C^hs@  
      kappae_x_PML_1(i) = 1.0+(kappa_x_max-1.0)*   & sXD.*D  
                                 ((nxPML_1 - i) / (nxPML_1 - 1.0))**m n/ \{}9  
      be_x_1(i) = EXP(-(sige_x_PML_1(i) / kappae_x_PML_1(i) +   & $\W|{u`  
                                 alphae_x_PML_1(i))*dt/epsO) O4Wn+$AN  
      if ((sige_x_PML_1(i) == 0.0) .and.        & sHk,#EsKH  
         (alphae_x_PML_1(i) == 0.0) .and. (i == nxPML_1)) then ewo
1^&#>  
         ce_x_1(i) = 0.0 _K4Igq  
      else ?0v(_ v  
         ce_x_1(i) = sige_x_PML_1(i)*(be_x_1(i)-1.0)/       & CXb-{|I}d  
            (sige_x_PML_1(i)+kappae_x_PML_1(i)*alphae_x_PML_1(i)) & '7+e!>"  
            / kappae_x_PML_1(i) nLA8Hy"8z  
      endif qjK'sge/  
   ENDDO  d"E@e21  
   DO i = 1,nxPML_1-1 52*9q!  
      sigh_x_PML_1(i) = sig_x_max * ( (nxPML_1 - i - 0.5)/(nxPML_1-1.0))**m J 8/]&Ow  
      alphah_x_PML_1(i) = alpha_x_max*((i-0.5)/(nxPML_1-1.0))**ma n9Mi?#xIp  
      kappah_x_PML_1(i) = 1.0+(kappa_x_max-1.0)*   & E3pnu.;U:_  
                            ((nxPML_1 - i - 0.5) / (nxPML_1 - 1.0))**m F,t ,Ja  
      bh_x_1(i) = EXP(-(sigh_x_PML_1(i) / kappah_x_PML_1(i) +   & X+'z@xpj  
                                 alphah_x_PML_1(i))*dt/epsO) ]kJinXHW  
      ch_x_1(i) = sigh_x_PML_1(i)*(bh_x_1(i)-1.0)/      & Dr;-2$Kt/&  
                  (sigh_x_PML_1(i)+kappah_x_PML_1(i)*alphah_x_PML_1(i)) & >)/,5VSE  
                  / kappah_x_PML_1(i) !+&Rn\e%7  
   ENDDO *5SOXrvhu6  
   DO i = 1,nxPML_2 To x{Sk3L  
      sige_x_PML_2(i) = sig_x_max * ( (nxPML_2 - i) / (nxPML_2 - 1.0) )**m B.wRZDEvc  
      alphae_x_PML_2(i) = alpha_x_max*((i-1.0)/(nxPML_2-1.0))**ma j}0W|*  
      kappae_x_PML_2(i) = 1.0+(kappa_x_max-1.0)*   & ZyJdz+L{@V  
                                 ((nxPML_2 - i) / (nxPML_2 - 1.0))**m i.eu$~F  
      be_x_2(i) = EXP(-(sige_x_PML_2(i) / kappae_x_PML_2(i) +   & bNUb  
                                 alphae_x_PML_2(i))*dt/epsO) *,/ADtL  
      if ((sige_x_PML_2(i) == 0.0) .and.        & SDdefB  
         (alphae_x_PML_2(i) == 0.0) .and. (i == nxPML_2)) then fM[Qn*.  
         ce_x_2(i) = 0.0 Ueq*R(9>  
      else eXHk6[%[  
         ce_x_2(i) = sige_x_PML_2(i)*(be_x_2(i)-1.0)/       & +=XDNSw  
            (sige_x_PML_2(i)+kappae_x_PML_2(i)*alphae_x_PML_2(i)) & QAp+LSm  
            / kappae_x_PML_2(i) ZT UaF4k j  
      endif EcBJ-j6d  
   ENDDO j@JhxCe1+R  
   DO i = 1,nxPML_2-1 On[:]#  
      sigh_x_PML_2(i) = sig_x_max * ( (nxPML_2 - i - 0.5)/(nxPML_2-1.0))**m ]Iku(<*Ya  
      alphah_x_PML_2(i) = alpha_x_max*((i-0.5)/(nxPML_2-1.0))**ma 3?Ml]=u  
      kappah_x_PML_2(i) = 1.0+(kappa_x_max-1.0)*   & &P+7U
m(  
                            ((nxPML_2 - i - 0.5) / (nxPML_2 - 1.0))**m Q(P'4XCm  
      bh_x_2(i) = EXP(-(sigh_x_PML_2(i) / kappah_x_PML_2(i) +   & 24ojjxz+  
                                 alphah_x_PML_2(i))*dt/epsO) Kq`C5  
      ch_x_2(i) = sigh_x_PML_2(i)*(bh_x_2(i)-1.0)/      & $1QQidB  
                  (sigh_x_PML_2(i)+kappah_x_PML_2(i)*alphah_x_PML_2(i)) & s0}OsHAj  
                  / kappah_x_PML_2(i) 5t`:=@u  
   ENDDO d)'am 3Q  
   DO j = 1,nyPML_1 x%acWeV5  
      sige_y_PML_1(j) = sig_y_max * ( (nyPML_1 - j ) / (nyPML_1 - 1.0) )**m :
Hk:Goo2  
      alphae_y_PML_1(j) = alpha_y_max*((j-1)/(nyPML_1-1.0))**ma  2EG`  
      kappae_y_PML_1(j) = 1.0+(kappa_y_max-1.0)*   & bCM&Fe0GM  
                                 ((nyPML_1 - j) / (nyPML_1 - 1.0))**m 8hx4s(1!  
      be_y_1(j) = EXP(-(sige_y_PML_1(j) / kappae_y_PML_1(j) +   & c:  /Wk  
                                 alphae_y_PML_1(j))*dt/epsO) B{\cV-X$0  
      if ((sige_y_PML_1(j) == 0.0) .and.        & =")}wl=s  
         (alphae_y_PML_1(j) == 0.0) .and. (j == nyPML_1)) then ohK_~  
         ce_y_1(j) = 0.0 NVEjUt/  
      else 0KW@j>=jK  
         ce_y_1(j) = sige_y_PML_1(j)*(be_y_1(j)-1.0)/       & J3oEN'
8S  
            (sige_y_PML_1(j)+kappae_y_PML_1(j)*alphae_y_PML_1(j)) &  E*[dc  
            / kappae_y_PML_1(j) CNC3">Dk~9  
      endif _JlbVe[<  
   ENDDO MYQZqlV  
   DO j = 1,nyPML_1-1 '@Zau\xC  
      sigh_y_PML_1(j) = sig_y_max * ( (nyPML_1 - j - 0.5)/(nyPML_1-1.0))**m fS08q9,S/  
      alphah_y_PML_1(j) = alpha_y_max*((j-0.5)/(nyPML_1-1.0))**ma 7j9:s>D  
      kappah_y_PML_1(j) = 1.0+(kappa_y_max-1.0)*   & vqq7IV)|  
                            ((nyPML_1 - j - 0.5) / (nyPML_1 - 1.0))**m I~LN)hqdo  
      bh_y_1(j) = EXP(-(sigh_y_PML_1(j) / kappah_y_PML_1(j) +   & RC_w 1:h  
                                 alphah_y_PML_1(j))*dt/epsO) jOJ$QT  
      ch_y_1(j) = sigh_y_PML_1(j)*(bh_y_1(j)-1.0)/      & KO`ftz3 +  
                  (sigh_y_PML_1(j)+kappah_y_PML_1(j)*alphah_y_PML_1(j)) & }GIwYh/  
                  / kappah_y_PML_1(j) XcoV27  
   ENDDO zEG6T*  
   DO j = 1,nyPML_2 D&Xh|}2A  
      sige_y_PML_2(j) = sig_y_max * ( (nyPML_2 - j ) / (nyPML_2 - 1.0) )**m \a:#e%]qz9  
      alphae_y_PML_2(j) = alpha_y_max*((j-1)/(nyPML_2-1.0))**ma H%,
jB<-.A  
      kappae_y_PML_2(j) = 1.0+(kappa_y_max-1.0)*   & (=v :@\r  
                                 ((nyPML_2 - j) / (nyPML_2 - 1.0))**m ^/|agQ7D2  
      be_y_2(j) = EXP(-(sige_y_PML_2(j) / kappae_y_PML_2(j) +   & <a)B5B>  
                                 alphae_y_PML_2(j))*dt/epsO) `
[hc{ynO|  
      if ((sige_y_PML_2(j) == 0.0) .and.        & QXXB>gOY5  
         (alphae_y_PML_2(j) == 0.0) .and. (j == nyPML_2)) then 4)L(41h  
         ce_y_2(j) = 0.0 hdL/zW7]  
      else 9(;5!q,Gsg  
         ce_y_2(j) = sige_y_PML_2(j)*(be_y_2(j)-1.0)/       & T#;W5<"  
            (sige_y_PML_2(j)+kappae_y_PML_2(j)*alphae_y_PML_2(j)) & TO&ohATp  
            / kappae_y_PML_2(j) /S32)=(  
      endif RlRkw+%m  
   ENDDO Z`5jX;Z!  
   DO j = 1,nyPML_2-1 #;9I3,@/Y  
      sigh_y_PML_2(j) = sig_y_max * ( (nyPML_2 - j - 0.5)/(nyPML_2-1.0))**m oF1{/ERS  
      alphah_y_PML_2(j) = alpha_y_max*((j-0.5)/(nyPML_2-1.0))**ma axJuJ`+Y  
      kappah_y_PML_2(j) = 1.0+(kappa_y_max-1.0)*   & X<\E 'v`~  
                            ((nyPML_2 - j - 0.5) / (nyPML_2 - 1.0))**m /}G+PUk7  
      bh_y_2(j) = EXP(-(sigh_y_PML_2(j) / kappah_y_PML_2(j) +   & t')47k\  
                                 alphah_y_PML_2(j))*dt/epsO) 9FB[`}  
      ch_y_2(j) = sigh_y_PML_2(j)*(bh_y_2(j)-1.0)/      & V$y6=Q<c  
                  (sigh_y_PML_2(j)+kappah_y_PML_2(j)*alphah_y_PML_2(j)) & q=NI}k  
                  / kappah_y_PML_2(j) Lu.D,oP  
   ENDDO en"]u,!  
   DO k = 1,nzPML_1 >s 4"2X  
      sige_z_PML_1(k) = sig_z_max * ( (nzPML_1 - k ) / (nzPML_1 - 1.0) )**m (@t O1g  
      alphae_z_PML_1(k) = alpha_z_max*((k-1)/(nzPML_1-1.0))**ma ~U] "
dbQ  
      kappae_z_PML_1(k) = 1.0+(kappa_z_max-1.0)*   & XCCh*qym  
                                 ((nzPML_1 - k) / (nzPML_1 - 1.0))**m 2TH13k$  
      be_z_1(k) = EXP(-(sige_z_PML_1(k) / kappae_z_PML_1(k) +   & k[TVu5R  
                                 alphae_z_PML_1(k))*dt/epsO) hu?Q,[+o  
      if ((sige_z_PML_1(k) == 0.0) .and.        & tUS)1*{_  
         (alphae_z_PML_1(k) == 0.0) .and. (k == nzPML_1)) then XILB>o.^3  
         ce_z_1(k) = 0.0 qPh @Bl3  
      else }ll&qb
 
         ce_z_1(k) = sige_z_PML_1(k)*(be_z_1(k)-1.0)/       & 3AKT>Wy =  
            (sige_z_PML_1(k)+kappae_z_PML_1(k)*alphae_z_PML_1(k)) & ]Gk;n/! B  
            / kappae_z_PML_1(k) D?+\"lI  
      endif uvNLm]*  
   ENDDO )e#KL$B)v  
   DO k = 1,nzPML_1-1 cT2&nZ  
      sigh_z_PML_1(k) = sig_z_max * ( (nzPML_1 - k - 0.5)/(nzPML_1-1.0))**m i!H)@4jX  
      alphah_z_PML_1(k) = alpha_z_max*((k-0.5)/(nzPML_1-1.0))**ma K U 2LJ_~Y  
      kappah_z_PML_1(k) = 1.0+(kappa_z_max-1.0)*   & {WBe(dc_%  
                            ((nzPML_1 - k - 0.5) / (nzPML_1 - 1.0))**m nz{ ;]U1  
      bh_z_1(k) = EXP(-(sigh_z_PML_1(k) / kappah_z_PML_1(k) +   & "r!>p\.0O  
                                 alphah_z_PML_1(k))*dt/epsO) bW"bkA80  
      ch_z_1(k) = sigh_z_PML_1(k)*(bh_z_1(k)-1.0)/      & j7(sYo@x7  
                  (sigh_z_PML_1(k)+kappah_z_PML_1(k)*alphah_z_PML_1(k)) & 8Ld`$_E  
                  / kappah_z_PML_1(k) jZjWz1+  
   ENDDO `pMI@"m  
   DO k = 1,nzPML_2 UL9
]LEGG  
      sige_z_PML_2(k) = sig_z_max * ( (nzPML_2 - k ) / (nzPML_2 - 1.0) )**m Rm@#GP`  
      alphae_z_PML_2(k) = alpha_z_max*((k-1)/(nzPML_2-1.0))**ma [v@3|@  
      kappae_z_PML_2(k) = 1.0+(kappa_z_max-1.0)*   & L<t>o":o  
                                 ((nzPML_2 - k) / (nzPML_2 - 1.0))**m oRf.34  
      be_z_2(k) = EXP(-(sige_z_PML_2(k) / kappae_z_PML_2(k) +   & fL0dy[Ch@  
                                 alphae_z_PML_2(k))*dt/epsO) U0%T<6*H  
      if ((sige_z_PML_2(k) == 0.0) .and.        & t

cDWx:Q  
         (alphae_z_PML_2(k) == 0.0) .and. (k == nzPML_2)) then Ck:RlF[6C  
         ce_z_2(k) = 0.0 dTV4 Q`Z  
      else |PGF g0li  
         ce_z_2(k) = sige_z_PML_2(k)*(be_z_2(k)-1.0)/       & " ^v/Y  
            (sige_z_PML_2(k)+kappae_z_PML_2(k)*alphae_z_PML_2(k)) &
a ~s:f5S>  
            / kappae_z_PML_2(k) ^
Rr!YnEN  
      endif <WXGDCj  
   ENDDO JD`IPQb~E  
   DO k = 1,nzPML_2-1 ;rggO0Y  
      sigh_z_PML_2(k) = sig_z_max * ( (nzPML_2 - k - 0.5)/(nzPML_2-1.0))**m s$%t*T2J>  
      alphah_z_PML_2(k) = alpha_z_max*((k-0.5)/(nzPML_2-1.0))**ma a0wSXd  
      kappah_z_PML_2(k) = 1.0+(kappa_z_max-1.0)*   & AnF"+<  
                            ((nzPML_2 - k - 0.5) / (nzPML_2 - 1.0))**m pUcN-WA  
      bh_z_2(k) = EXP(-(sigh_z_PML_2(k) / kappah_z_PML_2(k) +   & ^T?zR7r  
                                 alphah_z_PML_2(k))*dt/epsO) X;e=d+pw  
      ch_z_2(k) = sigh_z_PML_2(k)*(bh_z_2(k)-1.0)/      & |`T$Iq  
                  (sigh_z_PML_2(k)+kappah_z_PML_2(k)*alphah_z_PML_2(k)) & 4*AkUkP:T  
                  / kappah_z_PML_2(k) JfOBZQ  
   ENDDO kM:Z(Z7$  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ x;^DlyyYU  
!  FILL IN UPDATING COEFFICIENTS /M,C%.-  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 7 XNZEi9o  
   DA = 1.0 L3/SIoqd  
   DB = (dt/(muO)) O/XG}G.x|  
   DO i = 1,Imax d}RU-uiW  
      DO j = 1,Jmax r\PO?1  
         DO k = 1,Kmax (DLk+N4UHA  
            CA(i,j,k) = (1.0 - sig(i,j,k)*dt / (2.0*eps(i,j,k))) /        & ~U;M1>  
               (1.0 + sig(i,j,k) * dt / (2.0*eps(i,j,k))) g~7x+cu0  
            CB(i,j,k) = (dt/(eps(i,j,k))) /                            & <?2g\+{s9  
               (1.0 + sig(i,j,k)*dt / (2.0*eps(i,j,k))) }WowgY  
            ENDDO xK*G'3Ge  
 ENDDO c-`'`L^J  
   ENDDO o+],L_Ab  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ;\'d9C  
!  FILL IN DENOMINATORS FOR FIELD UPDATES  ff;9P5X  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 2_/H,  
   ii = nxPML_2-1 >='y+68  
   DO i = 1,Imax-1 T'ko
=k  
      if (i <= nxPML_1-1) then mm dQ\\  
         den_hx(i) = 1.0/(kappah_x_PML_1(i)*dx) ]$7yB3S,B  
      elseif (i >= Imax+1-nxPML_2) then .4[\%r\i  
         den_hx(i) = 1.0/(kappah_x_PML_2(ii)*dx) @u>:(9bp  
         ii = ii-1 =x xN3Ay  
      else Y5ZBP?P  
         den_hx(i) = 1.0/dx *'-[J2  
      endif [PIh^DhK  
   ENDDO kC_Kb&Q0  
   jj = nyPML_2-1 ~D0e\Q(A  
   DO j = 1,Jmax-1 EXF|;@-"  
      if (j <= nyPML_1-1) then $~ >/_<~  
         den_hy(j) = 1.0/(kappah_y_PML_1(j)*dy) rq#\x{l  
      elseif (j >= Jmax+1-nyPML_2) then tLE7s_^  
         den_hy(j) = 1.0/(kappah_y_PML_2(jj)*dy) g`Kh&|GU  
         jj = jj-1 <"hq}B  
      else ip+?k<]z  
         den_hy(j) = 1.0/dy .)g7s? K  
      endif I^GZ9@UE  
   ENDDO ^KF%Z2:$  
   kk = nzPML_2-1 @$7'{*  
   DO k = 2,Kmax-1 , 3&DA  
      if (k <= nzPML_1) then _'mK=`>u  
         den_hz(k) = 1.0/(kappah_z_PML_1(k-1)*dz) >2s31 {  
      elseif (k >= Kmax+1-nzPML_2) then oypF0?!m  
         den_hz(k) = 1.0/(kappah_z_PML_2(kk)*dz) ,k9xI<i  
         kk = kk - 1 f-BPT2U+  
      else
8C=8Wjm  
         den_hz(k) = 1.0/dz 3j6Am{9  
      endif Z[G[.\0  
   ENDDO D_Zt:tzO  
   ii = nxPML_2 FyhLMW3  
   DO i = 1,Imax-1 oazY?E]}3  
      if (i <= nxPML_1) then 5M&<tj/[a0  
         den_ex(i) = 1.0/(kappae_x_PML_1(i)*dx) o)L)|  
      elseif (i >= Imax+1-nxPML_2) then ?APCDZ^  
         den_ex(i) = 1.0/(kappae_x_PML_2(ii)*dx) @hF$qevX  
         ii = ii-1 >^~W'etX|  
      else "Ng%"N
z  
         den_ex(i) = 1.0/dx 4n0Iw  I  
      endif v.6K;TY.  
   ENDDO Vv&GyqoO]  
   jj = nyPML_2 "oT
&KW   
   DO j = 1,Jmax-1 nIqNhJ+  
      if (j <= nyPML_1) then vnx+1T  
         den_ey(j) = 1.0/(kappae_y_PML_1(j)*dy) S}3?  
      elseif (j >= Jmax+1-nyPML_2) then ._[uSBR'  
         den_ey(j) = 1.0/(kappae_y_PML_2(jj)*dy) Ew
)1O9f  
         jj = jj-1 >sGiDK @  
      else z7s}-w,  
         den_ey(j) = 1.0/dy 8\S$iGd  
      endif ,UNnz&H+f
 
   ENDDO zi_0*znw  
   kk = nzPML_2 a|-ozBFR  
   DO k = 1,Kmax-1 O%1uBc  
      if (k <= nzPML_1) then Tl L,dPM  
         den_ez(k) = 1.0/(kappae_z_PML_1(k)*dz) pS~=T}o  
      elseif (k >= Kmax-nzPML_2) then A~@x8  
         den_ez(k) = 1.0/(kappae_z_PML_2(kk)*dz) bMB@${i}  
         kk = kk - 1 bo-lT-I  
      else ]1s6=  
         den_ez(k) = 1.0/dz Kut@z>SK  
      endif cS%;JV>C  
   ENDDO MmQk@~  
x"P@[T  
!.:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. HP,{/ $i:  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 4C }#lW9  
!  BEGIN TIME STEP sdBB(  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !PfdY&.)  
  write(*,*)"begin time-stepping" W+i&!'  
  DO n = 1,nmax 9*S9~  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ gLwrYG7@  
!  UPDATE Hx b"y4-KV  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ eyuQ}
R  
   DO k = 2,Kmax-1 h,ipQ>  
      DO i = 1,Imax-1 ,0@QBr5P  
    DO j = 1,Jmax-1 CsJ&,(s(  
       Hx(i,j,k) = DA * Hx(i,j,k) + DB *       & K2gF;(  
   ( (Ez(i,j,k) - Ez(i,j+1,k))*den_hy(j)  +    & t+#vcg,G  
     (Ey(i,j,k) - Ey(i,j,k-1))*den_hz(k) ) r<[G~n  
    ENDDO mq+x=  
 ENDDO :{ Lihe~\  
      DO i = 1,Imax-1 =@P]eK/  
!..................................................................... lvH} 8lJ  
!  PML for bottom Hx, j-direction =J`M}BBx  
!..................................................................... W.l#@p  
         DO j = 1,nyPML_1-1 r8>Qs RnU%  
         psi_Hxy_1(i,j,k) = bh_y_1(j)*psi_Hxy_1(i,j,k)                 & q4SEvP}fLx  
     + ch_y_1(j) *(Ez(i,j,k) - Ez(i,j+1,k))/dy m#^;V  
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxy_1(i,j,k) pB0 SCS*  
         ENDDO hm0MO,i"  
!..................................................................... ,rX|_4n*  
!  PML for top Hx, j-direction @!*I mN
MI  
!..................................................................... #'97
mg  
         jj = nyPML_2-1 H`4KhdqR
 
         DO j = Jmax+1-nyPML_2,Jmax-1 S])*LUi  
         psi_Hxy_2(i,jj,k) = bh_y_2(jj)*psi_Hxy_2(i,jj,k)       & G%>{Z?!B  
     + ch_y_2(jj) *(Ez(i,j,k) -       & t;}`~B  
                                Ez(i,(j+1),k))/dy )T@?.J`  
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxy_2(i,jj,k) t4UL|fI  
            jj = jj-1 \s6VOR/  
         ENDDO T%yGSk  
 ENDDO CQs,G8\/  
   ENDDO "%w E>E  
   DO i = 1,Imax-1 JvF0s}#4  
      DO j = 1,Jmax-1 QsBC[7<jd-  
!..................................................................... RBpv40n0  
!  PML for bottom Hx, k-direction <l<y R?  
!..................................................................... O f]/tdPp  
         DO k = 2,nzPML_1 \F14]`i  
         psi_Hxz_1(i,j,k-1) = bh_z_1(k-1)*psi_Hxz_1(i,j,k-1)            & }J6 y NoXu  
     + ch_z_1(k-1) *(Ey(i,j,k) - Ey(i,j,k-1))/dz /H$:Q|T}  
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxz_1(i,j,k-1) UlPhW~F)  
         ENDDO N j:W6? A  
!..................................................................... X2`n&JE  
!  PML for top Hx, k-direction ~En]sj  
!..................................................................... H#3Ma1z  
         kk = nzPML_2-1 23X-h#w  
         DO k = Kmax+1-nzPML_2,Kmax-1 ft$!u-`  
         psi_Hxz_2(i,j,kk) = bh_z_2(kk)*psi_Hxz_2(i,j,kk)             & @+sYwlA~  
     + ch_z_2(kk) *(Ey(i,j,k) -       & !`dMTW  
                                Ey(i,j,k-1))/dz Fr#QM0--B  
            Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxz_2(i,j,kk) p: u@? k  
           kk = kk-1 z8j7K'vV1  
         ENDDO ]f6,4[  
 ENDDO W^i[7 r  
   ENDDO "(iQ-g Mm  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ H-3*},9  
!  UPDATE Hy 6('CB|ga  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ xmejoOF  
   DO k = 2,Kmax-1 h*?]A  
      DO i = 1,Imax-1 zh5ovA%  
    DO j = 1,Jmax-1 ~&D5RfK5f  
            Hy(i,j,k) = DA * Hy(i,j,k) + DB *             & @xr}(.  
   ( (Ez(i+1,j,k) - Ez(i,j,k))*den_hx(i) +      & 4veXg/l  
     (Ex(i,j,k-1) - Ex(i,j,k))*den_hz(k) ) 4Mox
P  
    ENDDO ywj'O e41  
      ENDDO e 3TKg  
      DO j = 1,Jmax-1 K@?K4o   
!..................................................................... @2sr/gX^  
!  PML for bottom Hy, i-direction \=H+m%  
!..................................................................... nnTiu,2R  
         DO i = 1,nxPML_1-1 Bqw/\Lxwlf  
       psi_Hyx_1(i,j,k) = bh_x_1(i)*psi_Hyx_1(i,j,k)                   & S<g~VK!Tt  
    + ch_x_1(i)*(Ez(i+1,j,k) - Ez(i,j,k))/dx %>1C($^  
       Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyx_1(i,j,k) D2f~*!vEnA  
         ENDDO @%8$k[
  
!..................................................................... X$=/H 6R5Z  
!  PML for top Hy, i-direction nq\~`vH|Gd  
!..................................................................... e\}'i-  
         ii = nxPML_2-1 M %!O)r#Pn  
         DO i = Imax+1-nxPML_2,Imax-1 6 )lWuY]e  
       psi_Hyx_2(ii,j,k) = bh_x_2(ii)*psi_Hyx_2(ii,j,k)     & &X,
6v  
    + ch_x_2(ii)*(Ez(i+1,j,k) -    & X.AE>fx*h  
                                Ez(i,j,k))/dx f.:0T&%G  
       Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyx_2(ii,j,k) #e[igxwi  
            ii = ii-1 K/$5SN1  
    ENDDO TVK*l*  
     ENDDO ])
#?rRw  
   ENDDO ^xq)Q?[{  
   DO i = 1,Imax-1 fr&K^je\  
      DO j = 1,Jmax-1 "Mh}n-oju  
!..................................................................... u6 4{w,  
!  PML for bottom Hy, k-direction S~B{G T\M  
!..................................................................... 26MoYO!k  
         DO k = 2,nzPML_1 $PMD$c  
       psi_Hyz_1(i,j,k-1) = bh_z_1(k-1)*psi_Hyz_1(i,j,k-1)          & zANsv9R~  
    + ch_z_1(k-1)*(Ex(i,j,k-1) - Ex(i,j,k))/dz 8Rnq &8A  
       Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyz_1(i,j,k-1) G1:"Gxja  
         ENDDO wq]vcY9^  
!..................................................................... 8-9<r  
!  PML for top Hy, k-direction T/tCX[}  
!..................................................................... 'JJ :  
         kk = nzPML_2-1 juMHc$d17  
         DO k = Kmax+1-nzPML_2,Kmax-1 WL;2&S/{@  
     psi_Hyz_2(i,j,kk) = bh_z_2(kk)*psi_Hyz_2(i,j,kk)               & < Q6  
    + ch_z_2(kk)*(Ex(i,j,k-1) -    & vb$i00?  
                                Ex(i,j,k))/dz Q 822 #  
     Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyz_2(i,j,kk) ",]A.,  
            kk = kk-1 SKo*8r   
         ENDDO %Lom#:L'  
     ENDDO Ci?
RuZ"  
   ENDDO 8aVj@x$'  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ AIw~@*T  
!  UPDATE Hz WyUa3$[gO  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ :`S\p[5  
   DO k = 1,Kmax-1 rEjEz+wu  
      DO i = 1,Imax-1 '\P+Bu]6&  
        DO j = 1,Jmax-1 oTveY  
            Hz(i,j,k) = DA * Hz(i,j,k) + DB       & n66b(6"mO2  
                  * ((Ey(i,j,k) - Ey(i+1,j,k))*den_hx(i) +        & X(Lz&fkd  
       (Ex(i,j+1,k) - Ex(i,j,k))*den_hy(j)) "sed{?  
    ENDDO )Mh5q&ow  
      ENDDO +<)H2  
      DO j = 1,Jmax-1 ootkf=  
!..................................................................... c@$W]o"A  
!  PML for bottom Hz, x-direction #yOeL3|b'  
!..................................................................... BX2}ar  
         DO i = 1,nxPML_1-1 cUwR6I9  
          psi_Hzx_1(i,j,k) = bh_x_1(i)*psi_Hzx_1(i,j,k)                 & v
%S$5  
     + ch_x_1(i) *(Ey(i,j,k) - Ey(i+1,j,k))/dx ?}No'E1!I  
       Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzx_1(i,j,k) nbBox,zW  
         ENDDO } A}Vd:#  
!..................................................................... Q7_#k66gb7  
!  PML for top Hz, x-direction IeB^BD+j  
!..................................................................... A s8IjGNs{  
         ii = nxPML_2-1 @1V?94T1  
         DO i = Imax+1-nxPML_2,Imax-1 x";4)u=  
          psi_Hzx_2(ii,j,k) = bh_x_2(ii)*psi_Hzx_2(ii,j,k)            & 6\x/Z=}L  
     + ch_x_2(ii) *(Ey(i,j,k) -       & [pz1f!Wn  
                                Ey(i+1,j,k))/dx q5.5%W  
       Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzx_2(ii,j,k) rGzGbI=  
            ii = ii-1 h4c4!S  
    ENDDO ki8Jl}dr  
      ENDDO }sM_^&e4X  
      DO i = 1,Imax-1 64SRW8AH  
!..................................................................... :!fU+2$`^(  
!  PML for bottom Hz, y-direction c22L]Sxo  
!..................................................................... VVc-Dx  
         DO j = 1,nyPML_1-1 y7LM}dH#m  
            psi_Hzy_1(i,j,k) = bh_y_1(j)*psi_Hzy_1(i,j,k)                   & j:0< tjE  
    + ch_y_1(j)*(Ex(i,j+1,k) - Ex(i,j,k))/dy ;JX2ebx  
       Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzy_1(i,j,k) o:@A%*jg  
         ENDDO "MDy0Tj8EN  
!..................................................................... DS7Pioa86  
!  PML for top Hz, y-direction Se[>z(  
!..................................................................... 6\m'MV`R!  
         jj = nyPML_2-1 p e$WSS J  
         DO j = Jmax+1-nyPML_2,Jmax-1 4TcW%  
            psi_Hzy_2(i,jj,k) = bh_y_2(jj)*psi_Hzy_2(i,jj,k)               & Kk>va->R  
    + ch_y_2(jj)*(Ex(i,j+1,k) -    & jtPHk*>^wu  
                                Ex(i,j,k))/dy >t4<2|!(M  
       Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzy_2(i,jj,k) UM. Se(kS  
            jj = jj-1 ;eMb$px  
    ENDDO hmv*IF.  
      ENDDO 
=9A!5  
   ENDDO QS<)*  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Xliw(B'\a4  
!  UPDATE Ex ! ;R}=  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Q'rX]kk_  
   DO k = 1,Kmax-1 "SpsSQ  
      DO i = 1,Imax-1 qcfg 55]'c  
    DO j = 2,Jmax-1 }K
V)F,`  
              IF (i >= istart-1 .and. i <= iend .and. j >= jstart .and.  & *!,k`=.([#  
                   j <= jend .and. k >= kstart .and. k <= kend) THEN )Hqn  
                 Ex(i,j,k) = 0.0 FwDEYG  
              ELSE on f
7V  
          Ex(i,j,k) = CA(i,j,k) * Ex(i,j,k) + CB(i,j,k) *       & l(~NpT{=V  
     ( (Hz(i,j,k) - Hz(i,j-1,k))*den_ey(j)  +    & Olr'n% }  
     (Hy(i,j,k) - Hy(i,j,k+1))*den_ez(k) ) &izk$~  
              ENDIF MYUL y2)  
    ENDDO Zc=#Y  
 ENDDO Wbi12{C  
      DO i = 1,Imax-1 $ 3R5p  
!..................................................................... /re0"!0y  
!  PML for bottom Ex, j-direction $~T|v7Y%  
!..................................................................... gO bP  
         DO j = 2,nyPML_1 6W)#FO`  
         psi_Exy_1(i,j,k) = be_y_1(j)*psi_Exy_1(i,j,k)                 & ?d)FYB  
     + ce_y_1(j) *(Hz(i,j,k) - Hz(i,j-1,k))/dy G4"[ynlWV  
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exy_1(i,j,k) T>m|C}yy  
         ENDDO 'L2[^iF9  
!..................................................................... vcSb:('  
!  PML for top Ex, j-direction 8%;]]{(B  
!..................................................................... LHq*E`  
         jj = nyPML_2 ":L d}~>  
         DO j = Jmax+1-nyPML_2,Jmax-1 9iGp0_J  
         psi_Exy_2(i,jj,k) = be_y_2(jj)*psi_Exy_2(i,jj,k)       & /[|A(,N}{  
     + ce_y_2(jj) *(Hz(i,j,k) -       & n&FRjq9y  
                                Hz(i,(j-1),k))/dy sl*&.F,v=  
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exy_2(i,jj,k) ,772$
7x  
            jj = jj-1 8VMD304  
         ENDDO df9$k0Fx  
 ENDDO !-8y;,P  
   ENDDO  bRx}ih  
   DO i = 1,Imax-1 _x<7^^VT  
      DO j = 2,Jmax-1 5,du2  
!..................................................................... /+f3jy:d  
!  PML for bottom Ex, k-direction D#11 N^-K  
!..................................................................... B{NGrC`5)  
         DO k = 1,nzPML_1 =]-D_$S~  
         psi_Exz_1(i,j,k) = be_z_1(k)*psi_Exz_1(i,j,k)                 & qhNYQ/uS  
     + ce_z_1(k) *(Hy(i,j,k) - Hy(i,j,k+1))/dz U.J/ "}5`T  
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exz_1(i,j,k) ;uqi  
         ENDDO .7n`]S/  
!..................................................................... cB7'>L  
!  PML for top Ex, k-direction uD4$<rSHb  
!..................................................................... ^/7Y3n!|3  
         kk = nzPML_2 eit%U  
         DO k = Kmax-nzPML_2,Kmax-1 }M${ _D  
         psi_Exz_2(i,j,kk) = be_z_2(kk)*psi_Exz_2(i,j,kk)             & (?z"_\^n/  
     + ce_z_2(kk) *(Hy(i,j,k) -       & ?YnB:z*eV  
                                Hy(i,j,k+1))/dz @*JS[w$1  
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exz_2(i,j,kk) 3NrWt2?  
            kk = kk-1 gm}zF%B"  
         ENDDO -aIB_  
 ENDDO t,qz%J&a  
   ENDDO zl$z>z)  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ vVH*\&H\T  
!  UPDATE Ey 9j$J}=y  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ z?(QM:  
   DO k = 1,Kmax-1 OyyE0  
      DO i = 2,Imax-1 x-SYfvYY  
    DO j = 1,Jmax-1 {y
|j**NZ  
              IF (i >= istart .and. i <= iend .and. j >= jstart-1 .and. & D=~3N  
                   j <= jend .and. k >= kstart .and. k <= kend) THEN tZA%^Y  
                 Ey(i,j,k) = 0.0 't3nh  
              ELSE 7niI65  
                 Ey(i,j,k) = CA(i,j,k) * Ey(i,j,k) + CB(i,j,k) *    & 8-q4'@(  
   ( (Hz(i-1,j,k) - Hz(i,j,k))*den_ex(i) +         & :JzJ(q/  
     (Hx(i,j,k+1) - Hx(i,j,k))*den_ez(k) ) zG~nRt{4  
              ENDIF qa5 T(
:8  
    ENDDO nPjN\Es6  
      ENDDO 9w$+Qc  
      DO j = 1,Jmax-1 M;E$ ]Z9  
!..................................................................... sT>l ?L  
!  PML for bottom Ey, i-direction jn>RE   
!..................................................................... j.UO>1{7  
         DO i = 2,nxPML_1 '];=1loD  
       psi_Eyx_1(i,j,k) = be_x_1(i)*psi_Eyx_1(i,j,k)                   & &(A'uX.>pr  
    + ce_x_1(i)*(Hz(i-1,j,k) - Hz(i,j,k))/dx RGLi#:0_.x  
       Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyx_1(i,j,k) }E

\u2]  
         ENDDO .Yxf0y?uv  
!..................................................................... Med0O~T%  
!  PML for top Ey, i-direction jNP%BNd1f  
!..................................................................... oY7 eVuz  
         ii = nxPML_2 *vu  
         DO i = Imax+1-nxPML_2,Imax-1 oqy}?<SQ  
       psi_Eyx_2(ii,j,k) = be_x_2(ii)*psi_Eyx_2(ii,j,k)              & >2K:O\&  
    + ce_x_2(ii)*(Hz(i-1,j,k) -    & xBAASy  
                               Hz(i,j,k))/dx ]PZ\N~T
 
       Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyx_2(ii,j,k) r'& 6P-Vm  
            ii = ii-1 =kn-F T  
    ENDDO IHj9n>c)[  
     ENDDO VY "i>Ae  
   ENDDO >FReGiK$T  
   DO i = 2,Imax-1 0/zgjT|fe  
      DO j = 1,Jmax-1 p}h9>R  
!..................................................................... w(nQ:;oC  
!  PML for bottom Ey, k-direction tCr?!Y~  
!..................................................................... 9C\@10D  
         DO k = 1,nzPML_1 \rmge4`4  
       psi_Eyz_1(i,j,k) = be_z_1(k)*psi_Eyz_1(i,j,k)        .. [2nP
r^  
{iVmae  

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

不错，可以很好的参考

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

下载看看吧，已经编译通过。 @"s\eL,r  
Code只能学习，不要拿来用（为了提高自己的能力）。 J1<fE(X  
实际上，只有自己编写Code,才会灵活运用，发现PML的问题和缺陷，并创新，然后发论文。 JXeqVKF  
我自己编写了许多PML的code，以此发表了许多SCI检索的论文。 1V`]sfRK  
71tMX[x  
请问楼主，从哪搞到的这个code？

这么厉害。写PML的code，能发很多SCI的文章。 CQ( @7  
0KQ8;&a|  
介绍一下经验吧。 =hs@W)-O  

5nj~RUK  

引用第3楼ljx690529于2008-12-08 20:06发表的 : GC\/B0!  
下载看看吧，已经编译通过。 ^w12k2a  
Code只能学习，不要拿来用（为了提高自己的能力）。 m~U{ V9;*  
实际上，只有自己编写Code,才会灵活运用，发现PML的问题和缺陷，并创新，然后发论文。 Nz}Q"6L  
我自己编写了许多PML的code，以此发表了许多SCI检索的论文。
.......

引用第3楼ljx690529于2008-12-08 20:06发表的 : ^"l$p,P+  
下载看看吧，已经编译通过。 tQ~B!j]  
Code只能学习，不要拿来用（为了提高自己的能力）。 0\#Q;Z2  
实际上，只有自己编写Code,才会灵活运用，发现PML的问题和缺陷，并创新，然后发论文。 g7P1]CZ}  
我自己编写了许多PML的code，以此发表了许多SCI检索的论文。
.......

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

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

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

谢谢分享！

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