displacement.f90

1D有限差分波动方程模拟

!=======================================================================

SUBROUTINE  DISPLACEMENT (D, DM)

  USE PRECISION   , ONLY:                                              &
                          WP     
  USE CONTROL_DATA, ONLY:                                              &
                          KTTO, MZ
  USE GRID_MEDIUM , ONLY:                                              &
                          DEN
  USE NONREF_BOUND, ONLY:                                              &
                          AZ, SZ_D, SZ_DP
  USE AUXIL       , ONLY:                                              &
                          A, B
  USE WAVEFIELD   , ONLY:                                              &
                          T

!-----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER :: L
  
  REAL(WP), DIMENSION(0:MZ), INTENT(INOUT) :: D, DM

!-----------------------------------------------------------------------


  IF      ( KTTO < 0 ) THEN   ! top boundary is free surface: use AFDA

    DM(0) = 2._WP*D (0) - DM(0) + DEN (0) *                            &
               (   35._WP/ 8._WP * T (0)                               &
                 - 35._WP/24._WP * T (1)                               &
                 + 21._WP/40._WP * T (2)                               &
                 -  5._WP/56._WP * T (3) )

    DM(1) = 2._WP*D (1) - DM(1) + DEN (1) *                            &
               ( - 31._WP/ 24._WP * T (0)                              &
                 + 29._WP/ 24._WP * T (1)                              &
                 -  3._WP/ 40._WP * T (2)                              &
                 +  1._WP/168._WP * T (3) )

  ELSE IF ( KTTO == 0 ) THEN  ! top boundary is plane of symmetry: use symmetric imaging

    DM(0) = DM(2)

    DM(1) = 2._WP*D (1) - DM(1) + DEN (1) *                            &
          ( A*( T (  2) - T (  0) )                                    &
          + B*( T (  1) - T (  0) ) )

  ELSE                        ! top boundary is 'non-reflecting': use 2nd order FDA
  
    L = 1
    DM(L) = 2._WP*D (L) - DM(L) + DEN (L) *                            &
          (     T (L  ) - T (L-1)   )

  END IF
      
  DO  L = 2, MZ-2             ! interior FD scheme

    DM(L) = 2._WP*D (L) - DM(L) + DEN (L) *                            &
          ( A*( T (L+1) - T (L-2) )                                    &
          + B*( T (L  ) - T (L-1) ) )

  END DO

  IF ( KTTO > 0 ) THEN        ! top boundary is 'non-reflecting': use FDA for NR boundary  
  
    DM(0) =                      AZ%ATOP01 *   DM(1)                   &
                               + AZ%ATOP02 *   DM(2)                   &
        + AZ%ATOP10 *  SZ_D%T0 + AZ%ATOP11 *  SZ_D%T1                  &
                               + AZ%ATOP12 *  SZ_D%T2                  &
        + AZ%ATOP20 * SZ_DP%T0 + AZ%ATOP21 * SZ_DP%T1                  &
                               + AZ%ATOP22 * SZ_DP%T2
                               
  END IF
  
                              ! bottom boundary is 'non-reflecting': use 2nd order FDA
  L = MZ-1

    DM(L) = 2._WP*D (L) - DM(L) + DEN (L) *                            &
          (     T (L  ) - T (L-1)   )
          

  L = MZ                      ! bottom boundary is 'non-reflecting': use FDA for NR boundary
  
    DM(MZ) =                       AZ%ABOT01 *   DM(MZ-1)              &
                                 + AZ%ABOT02 *   DM(MZ-2)              &
        + AZ%ABOT10 *  SZ_D%TMZ0 + AZ%ABOT11 *  SZ_D%TMZ1              &
                                 + AZ%ABOT12 *  SZ_D%TMZ2              &
        + AZ%ABOT20 * SZ_DP%TMZ0 + AZ%ABOT21 * SZ_DP%TMZ1              &
                                 + AZ%ABOT22 * SZ_DP%TMZ2
  


END  SUBROUTINE  DISPLACEMENT