! pml_update.f90! ! Original Berenger's Perfectly Matched layer: update rules !!    Copyright (C) 2007  Paul Panserrieu, < peutetre@cs.tu-berlin.de >!!    This program is free software: you can redistribute it and/or modify!    it under the terms of the GNU General Public License as published by!    the Free Software Foundation, either version 3 of the License.! ! last modified: 19-10-2007 12:08:51 PM CESTMODULE pml_updateUSE fdtd_gitterUSE pmlIMPLICIT NONECONTAINSSUBROUTINE ex_inner_update(g, b)  TYPE(gitter), INTENT(INOUT)      :: g  TYPE(pml_boundary), INTENT(IN)   :: b  INTEGER                          :: ix, iy, iz  DOUBLE PRECISION                 :: Hy_1, Hz_1  DO iz=g%nzl, g%nzyh, 1    DO iy=g%nyl, g%nyyh, 1      DO ix=g%nxl, g%nxyh, 1                IF (iz .EQ. g%nzl) THEN          Hy_1 = b%bas_z%H(ix, iy, iz-1, 3) + b%bas_z%H(ix, iy, iz-1, 4)        ELSE          Hy_1 = g%H(ix, iy, iz-1, 2)        ENDIF        IF (iy .EQ. g%nyl) THEN          Hz_1 = b%bas_y%H(ix, iy-1, iz, 5) + b%bas_y%H(ix, iy-1, iz, 6)        ELSE          Hz_1 = g%H(ix, iy-1, iz, 3)        ENDIF        g%E(ix, iy, iz, 1) = g%E(ix, iy, iz, 1)                             &                             + g%MAT_EPS                                    &                               * ( g%H(ix, iy, iz, 3)                       &                                 - g%H(ix, iy, iz, 2)                       &                                 - Hz_1                                     &                                 + Hy_1                                     &                                 )      ENDDO    ENDDO  ENDDOEND SUBROUTINE ex_inner_updateSUBROUTINE ey_inner_update(g, b)  TYPE(gitter), INTENT(INOUT)      :: g  TYPE(pml_boundary), INTENT(IN)   :: b  INTEGER                          :: ix, iy, iz  DOUBLE PRECISION                 :: Hx_1, Hz_1  DO iz=g%nzl, g%nzyh, 1    DO iy=g%nyl, g%nyyh, 1      DO ix=g%nxl, g%nxyh, 1        IF (iz .EQ. g%nzl) THEN          Hx_1 = b%bas_z%H(ix, iy, iz-1, 1) + b%bas_z%H(ix, iy, iz-1, 2)        ELSE          Hx_1 = g%H(ix, iy, iz-1, 1)        ENDIF        IF (ix .EQ. g%nxl) THEN          Hz_1 = b%bas_x%H(ix-1, iy, iz, 5) + b%bas_x%H(ix-1, iy, iz, 6)        ELSE          Hz_1 = g%H(ix-1, iy, iz, 3)        ENDIF        g%E(ix, iy, iz, 2)= g%E(ix, iy, iz, 2)                              &                             + g%MAT_EPS                                    &                               * ( g%H(ix, iy, iz, 1)                       &                                 - g%H(ix, iy, iz, 3)                       &                                 - Hx_1                                     &                                 + Hz_1                                     &                                 )      ENDDO    ENDDO  ENDDOEND SUBROUTINE ey_inner_updateSUBROUTINE ez_inner_update(g, b)  TYPE(gitter), INTENT(INOUT)      :: g  TYPE(pml_boundary), INTENT(IN)   :: b  INTEGER                          :: ix, iy, iz  DOUBLE PRECISION                 :: Hx_1, Hy_1  DO iz=g%nzl, g%nzyh, 1    DO iy=g%nyl, g%nyyh, 1      DO ix=g%nxl, g%nxyh, 1        IF (iy .EQ. g%nyl) THEN          Hx_1 = b%bas_y%H(ix, iy-1, iz, 1) + b%bas_y%H(ix, iy-1, iz, 2)        ELSE          Hx_1 = g%H(ix, iy-1, iz, 1)        ENDIF        IF (ix .EQ. g%nxl) THEN          Hy_1 = b%bas_x%H(ix-1, iy, iz, 3) + b%bas_x%H(ix-1, iy, iz, 4)        ELSE          Hy_1 = g%H(ix-1, iy, iz, 2)        ENDIF        g%E(ix, iy, iz, 3)= g%E(ix, iy, iz, 3)                              &                             + g%MAT_EPS                                    &                               * ( g%H(ix, iy, iz, 2)                       &                                 - g%H(ix, iy, iz, 1)                       &                                 - Hy_1                                     &                                 + Hx_1                                     &                                 )      ENDDO    ENDDO  ENDDOEND SUBROUTINE ez_inner_updateSUBROUTINE hx_inner_update(g, b)  TYPE(gitter), INTENT(INOUT)      :: g  TYPE(pml_boundary), INTENT(IN)   :: b  INTEGER                          :: ix, iy, iz  DOUBLE PRECISION                 :: Ey, Ez  DO iz=g%nzl, g%nzyh, 1    DO iy=g%nyl, g%nyyh, 1      DO ix=g%nxl, g%nxyh, 1                  IF (iy < g%nyyh) THEN          Ez = g%E(ix, iy+1, iz  , 3)        ELSE          Ez = b%top_y%E(ix, iy+1, iz, 5) + b%top_y%E(ix, iy+1, iz, 6)        ENDIF        IF (iz < g%nzyh) THEN          Ey = g%E(ix, iy, iz+1, 2)        ELSE          Ey = b%top_z%E(ix, iy, iz+1, 3) + b%top_z%E(ix, iy, iz+1, 4)        ENDIF        g%H(ix, iy, iz, 1) = g%H(ix, iy, iz, 1)                           &                            - g%MAT_MUE                                   &                              * ( g%E(ix, iy  , iz  , 2)                  &                                + Ez                                      &                                - Ey                                      &                                - g%E(ix, iy  , iz  , 3)                  &                                )      ENDDO    ENDDO  ENDDOEND SUBROUTINE hx_inner_updateSUBROUTINE hy_inner_update(g, b)  TYPE(gitter), INTENT(INOUT)      :: g  TYPE(pml_boundary), INTENT(IN)   :: b  INTEGER                          :: ix, iy, iz  DOUBLE PRECISION                 :: Ex, Ez  DO iz=g%nzl, g%nzyh, 1    DO iy=g%nyl, g%nyyh, 1      DO ix=g%nxl, g%nxyh, 1        IF (iz< g%nzyh) THEN          Ex = g%E(ix, iy, iz+1, 1)        ELSE          Ex = b%top_z%E(ix, iy, iz+1, 1) + b%top_z%E(ix, iy, iz+1, 2)        ENDIF        IF (ix < g%nxyh) THEN          Ez = g%E(ix+1, iy, iz, 3)        ELSE          Ez = b%top_x%E(ix+1, iy, iz, 5) + b%top_x%E(ix+1, iy, iz, 6)        ENDIF        g%H(ix, iy, iz, 2) = g%H(ix, iy, iz, 2)                             &                              - g%MAT_MUE                                   &                                * ( g%E(ix, iy, iz, 3)                      &                                  + Ex                                      &                                  - Ez                                      &                                  - g%E(ix, iy, iz, 1)                      &                                  )      ENDDO    ENDDO  ENDDOEND SUBROUTINE hy_inner_updateSUBROUTINE hz_inner_update(g, b)  TYPE(gitter), INTENT(INOUT)      :: g  TYPE(pml_boundary), INTENT(IN)   :: b  INTEGER                          :: ix, iy, iz  DOUBLE PRECISION                 :: Ex, Ey  DO iz=g%nzl, g%nzyh, 1    DO iy=g%nyl, g%nyyh, 1      DO ix=g%nxl, g%nxyh, 1        IF (ix < g%nxyh) THEN          Ey = g%E(ix+1, iy, iz, 2)        ELSE          Ey = b%top_x%E(ix+1, iy, iz, 3) + b%top_x%E(ix+1, iy, iz, 4)        ENDIF        IF (iy < g%nyyh) THEN          Ex = g%E(ix, iy+1, iz, 1)        ELSE          Ex = b%top_y%E(ix, iy+1, iz, 1) +  b%top_y%E(ix, iy+1, iz, 2)        ENDIF                  g%H(ix, iy, iz, 3) = g%H(ix, iy, iz, 3)                             &                              - g%MAT_MUE                                   &                                * ( g%E(ix, iy, iz, 1)                      &                                  + Ey                                      &                                  - Ex                                      &                                  - g%E(ix, iy, iz, 2)                      &                                  )      ENDDO    ENDDO  ENDDOEND SUBROUTINE hz_inner_update! avoid NaN when sigma = 0 with the right limit of expressionSUBROUTINE compute_E_factor(g, sigma, past_factor, new_factor)  TYPE(gitter), INTENT(IN)                                      :: g  DOUBLE PRECISION, INTENT(IN)                                  :: sigma  DOUBLE PRECISION, INTENT(IN)                                  :: past_factor  DOUBLE PRECISION, INTENT(INOUT)                               :: new_factor  IF (sigma .NE. 0.0d0) THEN    new_factor = (1.0d0 - past_factor) / (sigma * g%dx)  ELSE    new_factor = g%dt / (g%eps * g%dx)  ENDIFEND SUBROUTINE compute_E_factor SUBROUTINE compute_H_factor(g, sigma, past_factor, new_factor)  TYPE(gitter), INTENT(IN)                                      :: g  DOUBLE PRECISION, INTENT(IN)                                  :: sigma  DOUBLE PRECISION, INTENT(IN)                                  :: past_factor  DOUBLE PRECISION, INTENT(INOUT)                               :: new_factor  IF (sigma .NE. 0.0d0) THEN    new_factor = (1.0d0 - past_factor) / (sigma * (g%mue/g%eps) * g%dx)  ELSE    new_factor = g%dt / (g%mue * g%dx)  ENDIFEND SUBROUTINE compute_H_factor! PML Algorithm with PEC at last layerSUBROUTINE PML_update_Ex(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                                              :: Hzx_1, Hzy_1  DOUBLE PRECISION                                              :: Hyz_1, Hyx_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         IF (iy .EQ. b%min_yo .OR. iy .EQ. b%max_yo) THEN          b%bas_x%E(ix, iy, iz, 1) = 0.0d0          b%bas_x%E(ix, iy, iz, 2) = 0.0d0        ELSEIF (iz .EQ. b%min_zo .OR. iz .EQ. b%max_zo) THEN          b%bas_x%E(ix, iy, iz, 1) = 0.0d0          b%bas_x%E(ix, iy, iz, 2) = 0.0d0        ELSE          c1 = 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), c1, c2)          c3 = EXP(- b%bas_x%sigma(ix, iy, iz, 3)*g%dt/g%eps)          CALL compute_E_factor(g, b%bas_x%sigma(ix, iy, iz, 3), c3, c4)          b%bas_x%E(ix, iy, iz, 1) =    c1 *    b%bas_x%E(ix, iy, iz, 1)                &                                      + c2 *  ( b%bas_x%H(ix, iy, iz, 5)                &                                              - b%bas_x%H(ix, iy-1, iz, 5)              &                                              + b%bas_x%H(ix, iy, iz, 6)                &                                              - b%bas_x%H(ix, iy-1, iz, 6))          b%bas_x%E(ix, iy, iz, 2) =    c3 *   b%bas_x%E(ix, iy, iz, 2)                 &                                      - c4 * ( b%bas_x%H(ix, iy, iz, 3)                 &                                              - b%bas_x%H(ix, iy, iz-1, 3)              &                                              + b%bas_x%H(ix, iy, iz, 4)                &                                              - b%bas_x%H(ix, iy, iz-1, 4))        ENDIF      ENDDO    ENDDO  ENDDO  ! top_x  DO ix = g%nxgh, b%max_xo-1, 1    DO iy = b%min_yo, b%max_yo, 1      DO iz = b%min_zo, b%max_zo, 1