fdtd3d_cpml.f90

FDTD的fortran程序

!.....................................................................
         DO i = 1,nxPML_1-1
	      psi_Hyx_1(i,j,k) = bh_x_1(i)*psi_Hyx_1(i,j,k)                   &
				+ ch_x_1(i)*(Ez(i+1,j,k) - Ez(i,j,k))/dx
	      Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyx_1(i,j,k)
         ENDDO
!.....................................................................
!  PML for top Hy, i-direction
!.....................................................................
         ii = nxPML_2-1
         DO i = Imax+1-nxPML_2,Imax-1
	      psi_Hyx_2(ii,j,k) = bh_x_2(ii)*psi_Hyx_2(ii,j,k)     &
				+ ch_x_2(ii)*(Ez(i+1,j,k) -    &
                                Ez(i,j,k))/dx
	      Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyx_2(ii,j,k)
            ii = ii-1
	   ENDDO
     ENDDO
   ENDDO
   DO i = 1,Imax-1
      DO j = 1,Jmax-1
!.....................................................................
!  PML for bottom Hy, k-direction
!.....................................................................
         DO k = 2,nzPML_1
	      psi_Hyz_1(i,j,k-1) = bh_z_1(k-1)*psi_Hyz_1(i,j,k-1)          &
				+ ch_z_1(k-1)*(Ex(i,j,k-1) - Ex(i,j,k))/dz
	      Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyz_1(i,j,k-1)
         ENDDO
!.....................................................................
!  PML for top Hy, k-direction
!.....................................................................
         kk = nzPML_2-1
         DO k = Kmax+1-nzPML_2,Kmax-1
	    psi_Hyz_2(i,j,kk) = bh_z_2(kk)*psi_Hyz_2(i,j,kk)               &
				+ ch_z_2(kk)*(Ex(i,j,k-1) -    &
                                Ex(i,j,k))/dz
	    Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyz_2(i,j,kk)
            kk = kk-1
         ENDDO
     ENDDO
   ENDDO

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  UPDATE Hz
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DO k = 1,Kmax-1
      DO i = 1,Imax-1
        DO j = 1,Jmax-1
            Hz(i,j,k) = DA * Hz(i,j,k) + DB       &
                  * ((Ey(i,j,k) - Ey(i+1,j,k))*den_hx(i) +        &
			    (Ex(i,j+1,k) - Ex(i,j,k))*den_hy(j))
	   ENDDO
      ENDDO
      DO j = 1,Jmax-1
!.....................................................................
!  PML for bottom Hz, x-direction
!.....................................................................
         DO i = 1,nxPML_1-1
   	      psi_Hzx_1(i,j,k) = bh_x_1(i)*psi_Hzx_1(i,j,k)                 &
	 			+ ch_x_1(i) *(Ey(i,j,k) - Ey(i+1,j,k))/dx
	      Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzx_1(i,j,k)
         ENDDO
!.....................................................................
!  PML for top Hz, x-direction
!.....................................................................
         ii = nxPML_2-1
         DO i = Imax+1-nxPML_2,Imax-1
   	      psi_Hzx_2(ii,j,k) = bh_x_2(ii)*psi_Hzx_2(ii,j,k)            &
	 			+ ch_x_2(ii) *(Ey(i,j,k) -       &
                                Ey(i+1,j,k))/dx
	      Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzx_2(ii,j,k)
            ii = ii-1
	   ENDDO
      ENDDO
      DO i = 1,Imax-1
!.....................................................................
!  PML for bottom Hz, y-direction
!.....................................................................
         DO j = 1,nyPML_1-1
            psi_Hzy_1(i,j,k) = bh_y_1(j)*psi_Hzy_1(i,j,k)                   &
				+ ch_y_1(j)*(Ex(i,j+1,k) - Ex(i,j,k))/dy
	      Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzy_1(i,j,k)
         ENDDO
!.....................................................................
!  PML for top Hz, y-direction
!.....................................................................
         jj = nyPML_2-1
         DO j = Jmax+1-nyPML_2,Jmax-1
            psi_Hzy_2(i,jj,k) = bh_y_2(jj)*psi_Hzy_2(i,jj,k)               &
				+ ch_y_2(jj)*(Ex(i,j+1,k) -    &
                                Ex(i,j,k))/dy
	      Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzy_2(i,jj,k)
            jj = jj-1
	   ENDDO
      ENDDO
   ENDDO

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  UPDATE Ex
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DO k = 1,Kmax-1
      DO i = 1,Imax-1
	   DO j = 2,Jmax-1
              IF (i >= istart-1 .and. i <= iend .and. j >= jstart .and.  &
                   j <= jend .and. k >= kstart .and. k <= kend) THEN
                 Ex(i,j,k) = 0.0
              ELSE
	         Ex(i,j,k) = CA(i,j,k) * Ex(i,j,k) + CB(i,j,k) *       &
  			( (Hz(i,j,k) - Hz(i,j-1,k))*den_ey(j)  +    &
			  (Hy(i,j,k) - Hy(i,j,k+1))*den_ez(k) )
              ENDIF
	   ENDDO
	ENDDO 
      DO i = 1,Imax-1
!.....................................................................
!  PML for bottom Ex, j-direction
!.....................................................................
         DO j = 2,nyPML_1
  	      psi_Exy_1(i,j,k) = be_y_1(j)*psi_Exy_1(i,j,k)                 &
	 			+ ce_y_1(j) *(Hz(i,j,k) - Hz(i,j-1,k))/dy
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exy_1(i,j,k)
         ENDDO
!.....................................................................
!  PML for top Ex, j-direction
!.....................................................................
         jj = nyPML_2
         DO j = Jmax+1-nyPML_2,Jmax-1
  	      psi_Exy_2(i,jj,k) = be_y_2(jj)*psi_Exy_2(i,jj,k)       &
	 			+ ce_y_2(jj) *(Hz(i,j,k) -       &
                                Hz(i,(j-1),k))/dy
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exy_2(i,jj,k)
            jj = jj-1
         ENDDO
	ENDDO
   ENDDO
   DO i = 1,Imax-1
      DO j = 2,Jmax-1
!.....................................................................
!  PML for bottom Ex, k-direction
!.....................................................................
         DO k = 1,nzPML_1
  	      psi_Exz_1(i,j,k) = be_z_1(k)*psi_Exz_1(i,j,k)                 &
	 			+ ce_z_1(k) *(Hy(i,j,k) - Hy(i,j,k+1))/dz
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exz_1(i,j,k)
         ENDDO
!.....................................................................
!  PML for top Ex, k-direction
!.....................................................................
         kk = nzPML_2
         DO k = Kmax-nzPML_2,Kmax-1
  	      psi_Exz_2(i,j,kk) = be_z_2(kk)*psi_Exz_2(i,j,kk)             &
	 			+ ce_z_2(kk) *(Hy(i,j,k) -       &
                                Hy(i,j,k+1))/dz
            Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exz_2(i,j,kk)
            kk = kk-1
         ENDDO
	ENDDO
   ENDDO

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  UPDATE Ey
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DO k = 1,Kmax-1
      DO i = 2,Imax-1
	   DO j = 1,Jmax-1
              IF (i >= istart .and. i <= iend .and. j >= jstart-1 .and. &
                   j <= jend .and. k >= kstart .and. k <= kend) THEN
                 Ey(i,j,k) = 0.0
              ELSE
                 Ey(i,j,k) = CA(i,j,k) * Ey(i,j,k) + CB(i,j,k) *    &
			( (Hz(i-1,j,k) - Hz(i,j,k))*den_ex(i) +         &
			  (Hx(i,j,k+1) - Hx(i,j,k))*den_ez(k) )
              ENDIF
	   ENDDO 
      ENDDO
      DO j = 1,Jmax-1
!.....................................................................
!  PML for bottom Ey, i-direction
!.....................................................................
         DO i = 2,nxPML_1
	      psi_Eyx_1(i,j,k) = be_x_1(i)*psi_Eyx_1(i,j,k)                   &
				+ ce_x_1(i)*(Hz(i-1,j,k) - Hz(i,j,k))/dx
	      Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyx_1(i,j,k)
         ENDDO
!.....................................................................
!  PML for top Ey, i-direction
!.....................................................................
         ii = nxPML_2
         DO i = Imax+1-nxPML_2,Imax-1
	      psi_Eyx_2(ii,j,k) = be_x_2(ii)*psi_Eyx_2(ii,j,k)              &
				+ ce_x_2(ii)*(Hz(i-1,j,k) -    &
                               Hz(i,j,k))/dx
	      Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyx_2(ii,j,k)
            ii = ii-1
	   ENDDO
     ENDDO
   ENDDO
   DO i = 2,Imax-1
      DO j = 1,Jmax-1
!.....................................................................
!  PML for bottom Ey, k-direction
!.....................................................................
         DO k = 1,nzPML_1
	      psi_Eyz_1(i,j,k) = be_z_1(k)*psi_Eyz_1(i,j,k)                   &
				+ ce_z_1(k)*(Hx(i,j,k+1) - Hx(i,j,k))/dz
	      Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyz_1(i,j,k)
         ENDDO
!.....................................................................
!  PML for top Ey, k-direction
!.....................................................................
         kk = nzPML_2
         DO k = Kmax-nzPML_2,Kmax-1
	    psi_Eyz_2(i,j,kk) = be_z_2(kk)*psi_Eyz_2(i,j,kk)               &
				+ ce_z_2(kk)*(Hx(i,j,k+1) -    &
                                Hx(i,j,k))/dz
	    Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyz_2(i,j,kk)
            kk = kk-1
         ENDDO
     ENDDO
   ENDDO
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  UPDATE Ez
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DO k = 2,Kmax-1
      DO i = 2,Imax-1
         DO j = 2,Jmax-1
            Ez(i,j,k) = CA(i,j,k) * Ez(i,j,k) + CB(i,j,k)       &
                  * ((Hy(i,j,k) - Hy(i-1,j,k))*den_ex(i) +        &
			    (Hx(i,j-1,k) - Hx(i,j,k))*den_ey(j))
	   ENDDO
      ENDDO
      DO j = 2,Jmax-1
!.....................................................................
!  PML for bottom Ez, x-direction
!.....................................................................
         DO i = 2,nxPML_1
   	      psi_Ezx_1(i,j,k) = be_x_1(i)*psi_Ezx_1(i,j,k)             &
	 			+ ce_x_1(i) *(Hy(i,j,k) - Hy(i-1,j,k))/dx
	      Ez(i,j,k) = Ez(i,j,k) + CB(i,j,k)*psi_Ezx_1(i,j,k)
         ENDDO
!.....................................................................
!  PML for top Ez, x-direction
!.....................................................................
         ii = nxPML_2
         DO i = Imax+1-nxPML_2,Imax-1
   	      psi_Ezx_2(ii,j,k) = be_x_2(ii)*psi_Ezx_2(ii,j,k)       &
	 			+ ce_x_2(ii) *(Hy(i,j,k) -       &
                                Hy(i-1,j,k))/dx
	      Ez(i,j,k) = Ez(i,j,k) + CB(i,j,k)*psi_Ezx_2(ii,j,k)
            ii = ii-1
	   ENDDO
      ENDDO
      DO i = 2,Imax-1
!.....................................................................
!  PML for bottom Ez, y-direction
!.....................................................................
         DO j = 2,nyPML_1
            psi_Ezy_1(i,j,k) = be_y_1(j)*psi_Ezy_1(i,j,k)            &
				+ ce_y_1(j)*(Hx(i,j-1,k) - Hx(i,j,k))/dy
	      Ez(i,j,k) = Ez(i,j,k) + CB(i,j,k)*psi_Ezy_1(i,j,k)
         ENDDO
!.....................................................................
!  PML for top Ez, y-direction
!.....................................................................
         jj = nyPML_2
         DO j = Jmax+1-nyPML_2,Jmax-1
            psi_Ezy_2(i,jj,k) = be_y_2(jj)*psi_Ezy_2(i,jj,k)         &
				+ ce_y_2(jj)*(Hx(i,j-1,k) -    &
                                Hx(i,j,k))/dy
	      Ez(i,j,k) = Ez(i,j,k) + CB(i,j,k)*psi_Ezy_2(i,jj,k)
            jj = jj-1
	   ENDDO
      ENDDO
   ENDDO

!-----------------------------------------------------------------------
!   SOURCE
!-----------------------------------------------------------------------
   i = isource
   j = jsource
   k = ksource
   source = -2.0*((n*dt-tO)/tw) * exp(-((n*dt-tO)/tw)**2.0)
   Ez(i,j,k) = Ez(i,j,k) - CB(i,j,k)*source

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  RECORD GRID FOR VISUALIZATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   IF (n == record_grid) then
      DO j = 1,Jmax
         DO i = 1,Imax      
            write(33,*)Ez(i,j,ksource+1)
	   ENDDO
	ENDDO
   CLOSE(UNIT = 33)
   endif   
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  WRITE TO OUTPUT FILES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
   P1 = Ez(isource,jsource,ksource)
   write(30,*)P1      
   P2 = Ey(irecv1,jrecv1,krecv1)
   write(31,*)P2      

   IF (mod(n,10) == 0) then
    WRITE(*,*)n, " of ", nmax
   ENDIF

   ENDDO
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  END TIME STEP
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!.:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:.
    WRITE(*,*)"done time-stepping"

!-----------------------------------------------------------------------
! CLOSE OUTPUT FILES
!-----------------------------------------------------------------------   
   CLOSE(UNIT = 30)
   CLOSE(UNIT = 31)

   END PROGRAM fdtd3D_CPML
!cccccccc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5cccgtgccc6ccccccccc7cc