pml_update.f90

Sfdtd Simple finite-difference time-domain

            Hyx_1 = b%top_y%H(ix-1, iy, iz, 4)
          ENDIF
          IF(iy .EQ. g%nygh) THEN
            IF (iz < g%nzl) THEN
              Hxy_1 = b%bas_z%H(ix, iy-1, iz, 1)
              Hxz_1 = b%bas_z%H(ix, iy-1, iz, 2)
            ELSEIF(iz > g%nzyh) THEN
              Hxy_1 = b%top_z%H(ix, iy-1, iz, 1)
              Hxz_1 = b%top_z%H(ix, iy-1, iz, 2)
            ELSE
              Hxy_1 = g%H(ix, iy-1, iz, 1)
              Hxz_1 = 0.0d0
            ENDIF
          ELSE
            Hxy_1 = b%top_y%H(ix, iy-1, iz, 1)
            Hxz_1 = b%top_y%H(ix, iy-1, iz, 2)
          ENDIF

          c1 = 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), c1, c2)

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

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

          b%top_y%E(ix, iy, iz, 6) =     c3 *   b%top_y%E(ix, iy, iz, 6)                &
                                       - c4 * ( b%top_y%H(ix, iy, iz, 1)                &
                                              - Hxy_1                                   &
                                              + b%top_y%H(ix, iy, iz, 2)                &
                                              - Hxz_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 (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_z%H(ix-1, iy, iz, 3)
          Hyx_1 = b%bas_z%H(ix-1, iy, iz, 4)
        ENDIF
        IF (iy .EQ. g%nyl) THEN
           Hxy_1 = b%bas_y%H(ix, iy-1, iz, 1)
           Hxz_1 = b%bas_y%H(ix, iy-1, iz, 2)
        ELSE
          Hxy_1 = b%bas_z%H(ix, iy-1, iz, 1)
          Hxz_1 = b%bas_z%H(ix, iy-1, iz, 2)
        ENDIF

        c1 = 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), c1, c2)

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

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

        b%bas_z%E(ix, iy, iz, 6) =     c3 *   b%bas_z%E(ix, iy, iz, 6)                &
                                     - c4 * ( b%bas_z%H(ix, iy, iz, 1)                &
                                            - Hxy_1                                   &
                                            + b%bas_z%H(ix, iy, iz, 2)                &
                                            - Hxz_1 )
      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, 1 

        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_z%H(ix-1, iy, iz, 3)
          Hyx_1 = b%top_z%H(ix-1, iy, iz, 4)
        ENDIF
        IF (iy .EQ. g%nyl) THEN
          Hxy_1 = b%bas_y%H(ix, iy-1, iz, 1)
          Hxz_1 = b%bas_y%H(ix, iy-1, iz, 2)
        ELSE
          Hxy_1 = b%top_z%H(ix, iy-1, iz, 1)
          Hxz_1 = b%top_z%H(ix, iy-1, iz, 2)
        ENDIF

        c1 = 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), c1, c2)

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

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

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

END SUBROUTINE PML_update_Ez

SUBROUTINE PML_update_Hx(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                                              :: Ezx_1, Ezy_1
  DOUBLE PRECISION                                              :: Eyz_1, Eyx_1

  ! bas_x
  DO ix = b%min_xo, g%nxl-1, 1
    DO iy = b%min_yo, b%max_yo-1, 1
      DO iz = b%min_zo, b%max_zo-1, 1 
         
        c1 = EXP(- b%bas_x%sigma(ix, iy, iz, 5) * g%dt/g%eps)
        CALL compute_H_factor(g, b%bas_x%sigma(ix, iy, iz, 5), c1, c2)

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

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

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

  ! top_x
  DO ix = g%nxgh, b%max_xo-1, 1
    DO iy = b%min_yo, b%max_yo-1, 1
      DO iz = b%min_zo, b%max_zo-1, 1

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

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

        b%top_x%H(ix, iy, iz, 1) =     c1 *   b%top_x%H(ix, iy, iz, 1)                     &
                                     - c2 * ( b%top_x%E(ix, iy+1, iz, 5)                   &
                                            - b%top_x%E(ix, iy, iz, 5)                     &
                                            + b%top_x%E(ix, iy+1, iz, 6)                   &
                                            - b%top_x%E(ix, iy, iz, 6))

        b%top_x%H(ix, iy, iz, 2) =     c3 *   b%top_x%H(ix, iy, iz, 2)                     &
                                     + c4 * ( b%top_x%E(ix, iy, iz+1, 3)                   &
                                            - b%top_x%E(ix, iy, iz, 3)                     &
                                            + b%top_x%E(ix, iy, iz+1, 4)                   &
                                            - b%top_x%E(ix, iy, iz, 4))
      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. (g%nyl-1)) THEN
          IF (iz < g%nzl) THEN
            Ezx_1 = b%bas_z%E(ix, iy+1, iz, 5)
            Ezy_1 = b%bas_z%E(ix, iy+1, iz, 6)
          ELSEIF(iz > g%nzyh) THEN
            Ezx_1 = b%top_z%E(ix, iy+1, iz, 5)
            Ezy_1 = b%top_z%E(ix, iy+1, iz, 6)
          ELSE
            Ezx_1 = g%E(ix, iy+1, iz, 3)
            Ezy_1 = 0.0d0
          ENDIF
        ELSE
          Ezx_1 = b%bas_y%E(ix, iy+1, iz, 5)
          Ezy_1 = b%bas_y%E(ix, iy+1, iz, 6)
        ENDIF

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

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

        b%bas_y%H(ix, iy, iz, 1) =     c1 *   b%bas_y%H(ix, iy, iz, 1)                     &
                                     - c2 * ( Ezx_1                                        &
                                            - b%bas_y%E(ix, iy, iz, 5)                     &
                                            + Ezy_1                                        &
                                            - b%bas_y%E(ix, iy, iz, 6))

        b%bas_y%H(ix, iy, iz, 2) =     c3 *   b%bas_y%H(ix, iy, iz, 2)                     &
                                     + c4 * ( b%bas_y%E(ix, iy, iz+1, 3)                   &
                                            - b%bas_y%E(ix, iy, iz, 3)                     &
                                            + b%bas_y%E(ix, iy, iz+1, 4)                   &
                                            - b%bas_y%E(ix, iy, iz, 4))
      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, 1

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

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

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

        b%top_y%H(ix, iy, iz, 2) =     c3 *   b%top_y%H(ix, iy, iz, 2)                     &
                                     + c4 * ( b%top_y%E(ix, iy, iz+1, 3)                   &
                                            - b%top_y%E(ix, iy, iz, 3)                     &
                                            + b%top_y%E(ix, iy, iz+1, 4)                   &
                                            - b%top_y%E(ix, iy, iz, 4))
      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(iy .EQ. g%nyyh) THEN
          Ezx_1 = b%top_y%E(ix, iy+1, iz, 5)
          Ezy_1 = b%top_y%E(ix, iy+1, iz, 6)
        ELSE
          Ezx_1 = b%bas_z%E(ix, iy+1, iz, 5)
          Ezy_1 = b%bas_z%E(ix, iy+1, iz, 6)
        ENDIF
        IF(iz .EQ. (g%nzl-1)) THEN
          Eyz_1 = g%E(ix, iy, iz+1, 2)
          Eyx_1 = 0.0d0
        ELSE
          Eyz_1 = b%bas_z%E(ix, iy, iz+1, 3)
          Eyx_1 = b%bas_z%E(ix, iy, iz+1, 4)
        ENDIF

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

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

        b%bas_z%H(ix, iy, iz, 1) =     c1 *   b%bas_z%H(ix, iy, iz, 1)                     &
                                     - c2 * ( Ezx_1                                        &
                                            - b%bas_z%E(ix, iy, iz, 5)                     &
                                            + Ezy_1                                        &
                                            - b%bas_z%E(ix, iy, iz, 6))

        b%bas_z%H(ix, iy, iz, 2) =     c3 *   b%bas_z%H(ix, iy, iz, 2)                     &
                                     + c4 * ( Eyz_1                                        &
                                            - b%bas_z%E(ix, iy, iz, 3)                     &
                                            + Eyx_1                                        &
                                            - b%bas_z%E(ix, iy, iz, 4))
      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, 1 

       IF(iy .EQ. g%nyyh) THEN
          Ezx_1 = b%top_y%E(ix, iy+1, iz, 5)
          Ezy_1 = b%top_y%E(ix, iy+1, iz, 6)
        ELSE
          Ezx_1 = b%top_z%E(ix, iy+1, iz, 5)
          Ezy_1 = b%top_z%E(ix, iy+1, iz, 6)
        ENDIF

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

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