pml_update.f90

Sfdtd Simple finite-difference time-domain

        ELSE
          IF(ix .EQ. g%nxl) THEN
            Hzx_1 = b%bas_x%H(ix-1, iy, iz, 5)
            Hzy_1 = b%bas_x%H(ix-1, iy, iz, 6)
          ELSE
            Hzx_1 = b%bas_y%H(ix-1, iy, iz, 5)
            Hzy_1 = b%bas_y%H(ix-1, iy, iz, 6)
          ENDIF

          c1 = EXP(- b%bas_y%sigma(ix, iy, iz, 3)*g%dt/g%eps)
          CALL compute_E_factor(g, b%bas_y%sigma(ix, iy, iz, 3), c1, c2)

          c3 = EXP(- b%bas_y%sigma(ix, iy, iz, 1)*g%dt/g%eps)
          CALL compute_E_factor(g, b%bas_y%sigma(ix, iy, iz, 1), c3, c4)

          b%bas_y%E(ix, iy, iz, 3) =     c1 *   b%bas_y%E(ix, iy, iz, 3)                &
                                       + c2 * ( b%bas_y%H(ix, iy, iz, 1)                &
                                              - b%bas_y%H(ix, iy, iz-1, 1)              &
                                              + b%bas_y%H(ix, iy, iz, 2)                &
                                              - b%bas_y%H(ix, iy, iz-1, 2))

          b%bas_y%E(ix, iy, iz, 4) =     c3 *   b%bas_y%E(ix, iy, iz, 4)                &
                                       - c4 * ( b%bas_y%H(ix, iy, iz, 5)                &
                                              - Hzx_1                                   &
                                              + b%bas_y%H(ix, iy, iz, 6)                &
                                              - Hzy_1)
        ENDIF
      ENDDO
    ENDDO
  ENDDO

  ! top_y
  DO ix = g%nxl, g%nxyh, 1
    DO iy = g%nygh, b%max_yo-1, 1
      DO iz = b%min_zo, b%max_zo, 1
    
        IF (iz .EQ. b%min_zo .OR. iz .EQ. b%max_zo) THEN
          b%top_y%E(ix, iy, iz, 3) = 0.0d0
          b%top_y%E(ix, iy, iz, 4) = 0.0d0
        ELSE 
          IF(ix .EQ. g%nxl) THEN
            Hzx_1 = b%bas_x%H(ix-1, iy, iz, 5)
            Hzy_1 = b%bas_x%H(ix-1, iy, iz, 6)
          ELSE
            Hzx_1 = b%top_y%H(ix-1, iy, iz, 5)
            Hzy_1 = b%top_y%H(ix-1, iy, iz, 6)
          ENDIF

          c1 = EXP(- b%top_y%sigma(ix, iy, iz, 3)*g%dt/g%eps)
          CALL compute_E_factor(g, b%top_y%sigma(ix, iy, iz, 3), c1, c2)

          c3 = EXP(- b%top_y%sigma(ix, iy, iz, 1)*g%dt/g%eps)
          CALL compute_E_factor(g, b%top_y%sigma(ix, iy, iz, 1), c3, c4)

          b%top_y%E(ix, iy, iz, 3) =     c1 *   b%top_y%E(ix, iy, iz, 3)                &
                                       + c2 * ( b%top_y%H(ix, iy, iz, 1)                &
                                              - b%top_y%H(ix, iy, iz-1, 1)              &
                                              + b%top_y%H(ix, iy, iz, 2)                &
                                              - b%top_y%H(ix, iy, iz-1, 2))

          b%top_y%E(ix, iy, iz, 4) =     c3 *   b%top_y%E(ix, iy, iz, 4)                &
                                       - c4 * ( b%top_y%H(ix, iy, iz, 5)                &
                                              - Hzx_1                                   &
                                              + b%top_y%H(ix, iy, iz, 6)                &
                                              - Hzy_1)
        ENDIF
      ENDDO
    ENDDO
  ENDDO

  ! bas_z
  DO ix = g%nxl, g%nxyh, 1
    DO iy = g%nyl, g%nyyh, 1
      DO iz = b%min_zo, g%nzl-1, 1 

        IF (iz .EQ. b%min_zo) THEN
          b%bas_z%E(ix, iy, iz, 3) = 0.0d0
          b%bas_z%E(ix, iy, iz, 4) = 0.0d0
        ELSE
          IF (ix .EQ. g%nxl) THEN
            Hzx_1 = b%bas_x%H(ix-1, iy, iz, 5)
            Hzy_1 = b%bas_x%H(ix-1, iy, iz, 6)
          ELSE
            Hzx_1 = b%bas_z%H(ix-1, iy, iz, 5)
            Hzy_1 = b%bas_z%H(ix-1, iy, iz, 6)
          ENDIF

          c1 = EXP(- b%bas_z%sigma(ix, iy, iz, 3)*g%dt/g%eps)
          CALL compute_E_factor(g, b%bas_z%sigma(ix, iy, iz, 3), c1, c2)

          c3 = EXP(- b%bas_z%sigma(ix, iy, iz, 1)*g%dt/g%eps)
          CALL compute_E_factor(g, b%bas_z%sigma(ix, iy, iz, 1), c3, c4)

          b%bas_z%E(ix, iy, iz, 3) =     c1 *   b%bas_z%E(ix, iy, iz, 3)                &
                                       + c2 * ( b%bas_z%H(ix, iy, iz, 1)                &
                                              - b%bas_z%H(ix, iy, iz-1, 1)              &
                                              + b%bas_z%H(ix, iy, iz, 2)                &
                                              - b%bas_z%H(ix, iy, iz-1, 2))

          b%bas_z%E(ix, iy, iz, 4) =     c3 *   b%bas_z%E(ix, iy, iz, 4)                &
                                       - c4 * ( b%bas_z%H(ix, iy, iz, 5)                &
                                              - Hzx_1                                   &
                                              + b%bas_z%H(ix, iy, iz, 6)                &
                                              - Hzy_1)
        ENDIF
      ENDDO
    ENDDO
  ENDDO

  ! top_z
  DO ix = g%nxl, g%nxyh, 1
    DO iy = g%nyl, g%nyyh, 1
      DO iz = g%nzgh, b%max_zo, 1 

        IF (iz .EQ. b%max_zo) THEN
          b%top_z%E(ix, iy, iz, 3) = 0.0d0
          b%top_z%E(ix, iy, iz, 4) = 0.0d0
        ELSE
          IF (ix .EQ. g%nxl) THEN
            Hzx_1 = b%bas_x%H(ix-1, iy, iz, 5)
            Hzy_1 = b%bas_x%H(ix-1, iy, iz, 6)
          ELSE
            Hzx_1 = b%top_z%H(ix-1, iy, iz, 5)
            Hzy_1 = b%top_z%H(ix-1, iy, iz, 6)
          ENDIF
          IF(iz .EQ. g%nzgh) THEN
            Hxy_1 = 0.0d0
            Hxz_1 = g%H(ix, iy, iz-1, 1)
          ELSE
            Hxy_1 = b%top_z%H(ix, iy, iz-1, 1)
            Hxz_1 = b%top_z%H(ix, iy, iz-1, 2)
          ENDIF

          c1 = EXP(- b%top_z%sigma(ix, iy, iz, 3)*g%dt/g%eps)
          CALL compute_E_factor(g, b%top_z%sigma(ix, iy, iz, 3), c1, c2)

          c3 = EXP(- b%top_z%sigma(ix, iy, iz, 1)*g%dt/g%eps)
          CALL compute_E_factor(g, b%top_z%sigma(ix, iy, iz, 1), c3, c4)

          b%top_z%E(ix, iy, iz, 3) =     c1 *   b%top_z%E(ix, iy, iz, 3)                &
                                       + c2 * ( b%top_z%H(ix, iy, iz, 1)                &
                                              - Hxy_1                                   &
                                              + b%top_z%H(ix, iy, iz, 2)                &
                                              - Hxz_1)

          b%top_z%E(ix, iy, iz, 4) =     c3 *   b%top_z%E(ix, iy, iz, 4)                &
                                       - c4 * ( b%top_z%H(ix, iy, iz, 5)                &
                                              - Hzx_1                                   &
                                              + b%top_z%H(ix, iy, iz, 6)                &
                                              - Hzy_1)
        ENDIF
      ENDDO
    ENDDO
  ENDDO

END SUBROUTINE PML_update_Ey

SUBROUTINE PML_update_Ez(g, b)

  TYPE(gitter), INTENT(IN)                                      :: g
  TYPE(pml_boundary), INTENT(INOUT)                             :: b
  INTEGER                                                       :: ix, iy, iz
  DOUBLE PRECISION                                              :: c1, c2, c3, c4
  DOUBLE PRECISION                                              :: Hyz_1, Hyx_1
  DOUBLE PRECISION                                              :: Hxy_1, Hxz_1
  
  ! bas_x
  DO ix = b%min_xo, g%nxl-1, 1
    DO iy = b%min_yo, b%max_yo, 1
      DO iz = b%min_zo, b%max_zo-1, 1 

        IF(ix .EQ. b%min_xo) THEN
          b%bas_x%E(ix, iy, iz, 5) = 0.0d0
          b%bas_x%E(ix, iy, iz, 6) = 0.0d0
        ELSEIF (iy .EQ. b%min_yo  .OR. iy .EQ. b%max_yo) THEN
          b%bas_x%E(ix, iy, iz, 5) = 0.0d0
          b%bas_x%E(ix, iy, iz, 6) = 0.0d0
        ELSE

          c1 = EXP(- b%bas_x%sigma(ix, iy, iz, 1)*g%dt/g%eps)
          CALL compute_E_factor(g, b%bas_x%sigma(ix, iy, iz, 1), c1, c2)

          c3 = EXP(- b%bas_x%sigma(ix, iy, iz, 2)*g%dt/g%eps)
          CALL compute_E_factor(g, b%bas_x%sigma(ix, iy, iz, 2), c3, c4)

          b%bas_x%E(ix, iy, iz, 5) =     c1 *   b%bas_x%E(ix, iy, iz, 5)                &
                                       + c2 * ( b%bas_x%H(ix, iy, iz, 3)                &
                                              - b%bas_x%H(ix-1, iy, iz, 3)              &
                                              + b%bas_x%H(ix, iy, iz, 4)                &
                                              - b%bas_x%H(ix-1, iy, iz, 4))

          b%bas_x%E(ix, iy, iz, 6) =     c3 *   b%bas_x%E(ix, iy, iz, 6)                &
                                       - c4 * ( b%bas_x%H(ix, iy, iz, 1)                &
                                              - b%bas_x%H(ix, iy-1, iz, 1)              &
                                              + b%bas_x%H(ix, iy, iz, 2)                &
                                              - b%bas_x%H(ix, iy-1, iz, 2))
        ENDIF
      ENDDO
    ENDDO
  ENDDO

  ! top_x
  DO ix = g%nxgh, b%max_xo, 1
    DO iy = b%min_yo, b%max_yo, 1
      DO iz = b%min_zo, b%max_zo-1, 1 
      
        IF (ix .EQ. b%max_xo) THEN
          b%top_x%E(ix, iy, iz, 5) = 0.0d0
          b%top_x%E(ix, iy, iz, 6) = 0.0d0
        ELSEIF (iy .EQ. b%min_yo .OR. iy .EQ. b%max_yo) THEN
          b%top_x%E(ix, iy, iz, 5) = 0.0d0
          b%top_x%E(ix, iy, iz, 6) = 0.0d0
        ELSE
          IF(ix .EQ. g%nxgh) THEN
            IF (iy < g%nyl) THEN
              Hyz_1 = b%bas_y%H(ix-1, iy, iz, 3)
              Hyx_1 = b%bas_y%H(ix-1, iy, iz, 4)
            ELSEIF (iy > g%nyyh) THEN
              Hyz_1 = b%top_y%H(ix-1, iy, iz, 3)
              Hyx_1 = b%top_y%H(ix-1, iy, iz, 4)
            ELSE
              IF (iz < g%nzl) THEN
                Hyz_1 = b%bas_z%H(ix-1, iy, iz, 3)
                Hyx_1 = b%bas_z%H(ix-1, iy, iz, 4)
              ELSEIF (iz > g%nzyh) THEN
                Hyz_1 = b%top_z%H(ix-1, iy, iz, 3)
                Hyx_1 = b%top_z%H(ix-1, iy, iz, 4) 
              ELSE
                Hyz_1 = g%H(ix-1, iy, iz, 2)
                Hyx_1 = 0.0d0
              ENDIF
            ENDIF
          ELSE
            Hyz_1 = b%top_x%H(ix-1, iy, iz, 3)
            Hyx_1 = b%top_x%H(ix-1, iy, iz, 4)
          ENDIF

          c1 = EXP(- b%top_x%sigma(ix, iy, iz, 1)*g%dt/g%eps)
          CALL compute_E_factor(g, b%top_x%sigma(ix, iy, iz, 1), c1, c2)

          c3 = EXP(- b%top_x%sigma(ix, iy, iz, 2)*g%dt/g%eps)
          CALL compute_E_factor(g, b%top_x%sigma(ix, iy, iz, 2), c3, c4)

          b%top_x%E(ix, iy, iz, 5) =     c1 *   b%top_x%E(ix, iy, iz, 5)                &
                                       + c2 * ( b%top_x%H(ix, iy, iz, 3)                &
                                              - Hyz_1                                   &
                                              + b%top_x%H(ix, iy, iz, 4)                &
                                              - Hyx_1 )
  
          b%top_x%E(ix, iy, iz, 6) =     c3 *   b%top_x%E(ix, iy, iz, 6)                &
                                       - c4 * ( b%top_x%H(ix, iy, iz, 1)                &
                                              - b%top_x%H(ix, iy-1, iz, 1)              &
                                              + b%top_x%H(ix, iy, iz, 2)                &
                                              - b%top_x%H(ix, iy-1, iz, 2))
        ENDIF
      ENDDO
    ENDDO
  ENDDO

  ! bas_y
  DO ix = g%nxl, g%nxyh, 1
    DO iy = b%min_yo, g%nyl-1, 1
      DO iz = b%min_zo, b%max_zo-1, 1 

        IF (iy .EQ. b%min_yo) THEN
          b%bas_y%E(ix, iy, iz, 5) = 0.0d0
          b%bas_y%E(ix, iy, iz, 6) = 0.0d0
        ELSE
          IF(ix .EQ. g%nxl) THEN
            Hyz_1 = b%bas_x%H(ix-1, iy, iz, 3)
            Hyx_1 = b%bas_x%H(ix-1, iy, iz, 4)
          ELSE
            Hyz_1 = b%bas_y%H(ix-1, iy, iz, 3)
            Hyx_1 = b%bas_y%H(ix-1, iy, iz, 4)
          ENDIF

          c1 = EXP(- b%bas_y%sigma(ix, iy, iz, 1)*g%dt/g%eps)
          CALL compute_E_factor(g, b%bas_y%sigma(ix, iy, iz, 1), c1, c2)

          c3 = EXP(- b%bas_y%sigma(ix, iy, iz, 2)*g%dt/g%eps)
          CALL compute_E_factor(g, b%bas_y%sigma(ix, iy, iz, 2), c3, c4)

          b%bas_y%E(ix, iy, iz, 5) =     c1 *   b%bas_y%E(ix, iy, iz, 5)                &
                                       + c2 * ( b%bas_y%H(ix, iy, iz, 3)                &
                                              - Hyz_1                                   &
                                              + b%bas_y%H(ix, iy, iz, 4)                &
                                              - Hyx_1 )

          b%bas_y%E(ix, iy, iz, 6) =     c3 *   b%bas_y%E(ix, iy, iz, 6)                &
                                       - c4 * ( b%bas_y%H(ix, iy, iz, 1)                &
                                              - b%bas_y%H(ix, iy-1, iz, 1)              &
                                              + b%bas_y%H(ix, iy, iz, 2)                &
                                              - b%bas_y%H(ix, iy-1, iz, 2))
        ENDIF
      ENDDO
    ENDDO
  ENDDO

  ! top_y
  DO ix = g%nxl, g%nxyh, 1
    DO iy = g%nygh, b%max_yo, 1
      DO iz = b%min_zo, b%max_zo-1, 1 

        IF (iy .EQ. b%max_yo) THEN
          b%top_y%E(ix, iy, iz, 5) = 0.0d0
          b%top_y%E(ix, iy, iz, 6) = 0.0d0
        ELSE
          IF(ix .EQ. g%nxl) THEN
            Hyz_1 = b%bas_x%H(ix-1, iy, iz, 3)
            Hyx_1 = b%bas_x%H(ix-1, iy, iz, 4)
          ELSE
            Hyz_1 = b%top_y%H(ix-1, iy, iz, 3)