[cig-commits] r21959 - seismo/3D/SPECFEM3D/trunk/src/specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Mon Apr 29 05:21:01 PDT 2013
Author: xie.zhinan
Date: 2013-04-29 05:21:01 -0700 (Mon, 29 Apr 2013)
New Revision: 21959
Modified:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90
Log:
clean the code of pml
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90 2013-04-27 22:59:14 UTC (rev 21958)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90 2013-04-29 12:21:01 UTC (rev 21959)
@@ -32,7 +32,8 @@
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner)
+ ispec_is_inner,phase_is_inner,&
+ PML_CONDITIONS,spec_to_CPML,is_CPML)
! returns the updated pressure array: potential_dot_dot_acoustic
@@ -42,7 +43,7 @@
integer :: NSPEC_AB,NGLOB_AB
! displacement and pressure
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic
! global indexing
@@ -59,13 +60,16 @@
logical, dimension(NSPEC_AB) :: ispec_is_inner
logical :: phase_is_inner
+! CPML
+ logical :: PML_CONDITIONS
+ integer :: spec_to_CPML(NSPEC_AB)
+ logical :: is_CPML(NSPEC_AB)
+
! local parameters
real(kind=CUSTOM_REAL) :: displ_x,displ_y,displ_z,displ_n
real(kind=CUSTOM_REAL) :: nx,ny,nz,jacobianw
+ integer :: iface,igll,ispec,iglob,ispec_CPML,i,j,k
- integer :: iface,igll,ispec,iglob
- integer :: i,j,k
-
! loops on all coupling faces
do iface = 1,num_coupling_ac_el_faces
@@ -87,9 +91,21 @@
iglob = ibool(i,j,k,ispec)
! elastic displacement on global point
- displ_x = displ(1,iglob)
- displ_y = displ(2,iglob)
- displ_z = displ(3,iglob)
+ if(PML_CONDITIONS)then
+ if(is_CPML(ispec))then
+ ispec_CPML = spec_to_CPML(ispec)
+ call pml_compute_memory_variables_acoustic_elastic(ispec_CPML,iface,iglob,i,j,k,&
+ displ_x,displ_y,displ_z)
+ else
+ displ_x = displ(1,iglob)
+ displ_y = displ(2,iglob)
+ displ_z = displ(3,iglob)
+ endif
+ else
+ displ_x = displ(1,iglob)
+ displ_y = displ(2,iglob)
+ displ_z = displ(3,iglob)
+ endif
! gets associated normal on GLL point
! (note convention: pointing outwards of acoustic element)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-04-27 22:59:14 UTC (rev 21958)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-04-29 12:21:01 UTC (rev 21959)
@@ -57,6 +57,7 @@
use specfem_par_acoustic
use specfem_par_elastic
use specfem_par_poroelastic
+ use pml_par,only: spec_to_CPML,is_CPML
implicit none
@@ -177,7 +178,8 @@
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner)
+ ispec_is_inner,phase_is_inner,&
+ PML_CONDITIONS,spec_to_CPML,is_CPML)
else
! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
! adjoint definition: \partial_t^2 \bfs^\dagger=-\frac{1}{\rho}\bfnabla\phi^\dagger
@@ -187,7 +189,8 @@
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner)
+ ispec_is_inner,phase_is_inner,&
+ PML_CONDITIONS,spec_to_CPML,is_CPML)
endif
! adjoint/kernel simulations
if( SIMULATION_TYPE == 3 ) &
@@ -197,7 +200,8 @@
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner)
+ ispec_is_inner,phase_is_inner,&
+ PML_CONDITIONS,spec_to_CPML,is_CPML)
else
! on GPU
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90 2013-04-27 22:59:14 UTC (rev 21958)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90 2013-04-29 12:21:01 UTC (rev 21959)
@@ -31,6 +31,8 @@
use pml_par
use specfem_par, only: NSPEC_AB
use constants, only: NDIM,NGLLX,NGLLY,NGLLZ
+ use specfem_par_acoustic, only: ACOUSTIC_SIMULATION,num_coupling_ac_el_faces
+ use specfem_par_elastic, only: ELASTIC_SIMULATION
implicit none
@@ -168,6 +170,12 @@
allocate(potential_dot_dot_acoustic_CPML(NGLLX,NGLLY,NGLLZ,NSPEC_CPML),stat=ier)
if(ier /= 0) stop 'error allocating potential_dot_dot_acoustic_CPML array'
+ ! stores C-PML contribution on elastic/acoustic interface
+ if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
+ allocate(rmemory_coupling_ac_el_displ(3,NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_coupling_ac_el_displ array'
+ endif
+
spec_to_CPML(:) = 0
CPML_type(:) = 0
@@ -236,4 +244,8 @@
potential_dot_dot_acoustic_CPML(:,:,:,:) = 0._CUSTOM_REAL
+ if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
+ rmemory_coupling_ac_el_displ(:,:,:,:,:,:) = 0._CUSTOM_REAL
+ endif
+
end subroutine pml_allocate_arrays
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-04-27 22:59:14 UTC (rev 21958)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-04-29 12:21:01 UTC (rev 21959)
@@ -113,14 +113,6 @@
rmemory_duz_dxl_z(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,2) &
+ PML_duz_dxl_new(i,j,k,ispec_CPML) * coef1_2 + PML_duz_dxl(i,j,k,ispec_CPML) * coef2_2
- duxdxl_x = A6 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,1) + A8 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,2)
- duydxl_y = A6 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,1) + A8 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,2)
- duzdxl_z = A6 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,1) + A8 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,2)
-
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_x(i,j,k,ispec_CPML)
A10 = d_store_x(i,j,k,ispec_CPML)
@@ -149,14 +141,6 @@
+ PML_duz_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dyl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdyl_x = A9 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,1) + A11 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,2)
- duydyl_y = A9 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,1) + A11 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,2)
- duzdyl_z = A9 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,1) + A11 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,2)
-
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_x(i,j,k,ispec_CPML)
A13 = d_store_x(i,j,k,ispec_CPML)
@@ -185,13 +169,6 @@
+ PML_duz_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdzl_x = A12 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,1) + A14 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,2)
- duydzl_y = A12 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,1) + A14 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,2)
- duzdzl_z = A12 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,1) + A14 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A15 and A16 --------------------------
A15 = k_store_x(i,j,k,ispec_CPML)
A16 = d_store_x(i,j,k,ispec_CPML)
@@ -215,11 +192,6 @@
+ PML_duz_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_y = A15 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_y(i,j,k,ispec_CPML,2)
- duzdyl_y = A15 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dyl_y(i,j,k,ispec_CPML,2)
-
rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) &
+ PML_duy_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
@@ -228,12 +200,6 @@
+ PML_duy_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dyl_z(i,j,k,ispec_CPML,2) = 0.d0
- duydzl_z = A15 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dzl_z(i,j,k,ispec_CPML,2)
- duydyl_z = A15 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_z(i,j,k,ispec_CPML,2)
-
-
!---------------------- A17 and A18 --------------------------
A17 = 1.d0
A18 = 0.0
@@ -244,23 +210,12 @@
rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duz_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_x = A17 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_x(i,j,k,ispec_CPML,2)
- duzdxl_x = A17 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dxl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdzl_z = A17 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dzl_z(i,j,k,ispec_CPML,2)
- duxdxl_z = A17 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_z(i,j,k,ispec_CPML,2)
-
-
!---------------------- A19 and A20 --------------------------
A19 = 1.d0
A20 = 0.0
@@ -271,22 +226,12 @@
rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duy_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duydyl_x = A19 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_x(i,j,k,ispec_CPML,2)
- duydxl_x = A19 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dxl_y(i,j,k,ispec_CPML,2) = 0.d0
- duxdyl_y = A19 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
- duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
-
else if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- Y-surface C-PML ---------------------------------
@@ -320,13 +265,6 @@
+ PML_duz_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dxl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdxl_x = A6 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,1) + A8 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,2)
- duydxl_y = A6 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,1) + A8 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,2)
- duzdxl_z = A6 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,1) + A8 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = 1.d0/k_store_y(i,j,k,ispec_CPML)
A10 = 0.d0
@@ -355,13 +293,6 @@
rmemory_duz_dyl_z(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,2) &
+ PML_duz_dyl_new(i,j,k,ispec_CPML) * coef1_2 + PML_duz_dyl(i,j,k,ispec_CPML) * coef2_2
- duxdyl_x = A9 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,1) + A11 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,2)
- duydyl_y = A9 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,1) + A11 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,2)
- duzdyl_z = A9 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,1) + A11 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_y(i,j,k,ispec_CPML)
A13 = d_store_y(i,j,k,ispec_CPML)
@@ -390,13 +321,6 @@
+ PML_duz_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdzl_x = A12 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,1) + A14 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,2)
- duydzl_y = A12 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,1) + A14 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,2)
- duzdzl_z = A12 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,1) + A14 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A15 and A16 --------------------------
A15 = 1.d0
A16 = 0.d0
@@ -407,22 +331,12 @@
rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duz_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_y = A15 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_y(i,j,k,ispec_CPML,2)
- duzdyl_y = A15 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dyl_y(i,j,k,ispec_CPML,2)
-
rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duy_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duy_dyl_z(i,j,k,ispec_CPML,2) = 0.d0
- duydzl_z = A15 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dzl_z(i,j,k,ispec_CPML,2)
- duydyl_z = A15 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A17 and A18 --------------------------
A17 = k_store_y(i,j,k,ispec_CPML)
A18 = d_store_y(i,j,k,ispec_CPML)
@@ -446,11 +360,6 @@
+ PML_duz_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_x = A17 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_x(i,j,k,ispec_CPML,2)
- duzdxl_x = A17 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) &
+ PML_dux_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
@@ -459,12 +368,6 @@
+ PML_dux_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dxl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdzl_z = A17 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dzl_z(i,j,k,ispec_CPML,2)
- duxdxl_z = A17 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_z(i,j,k,ispec_CPML,2)
-
-
!---------------------- A19 and A20--------------------------
A19 = 1.d0
A20 = 0.0
@@ -475,22 +378,12 @@
rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duy_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duydyl_x = A19 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_x(i,j,k,ispec_CPML,2)
- duydxl_x = A19 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dxl_y(i,j,k,ispec_CPML,2) = 0.d0
- duxdyl_y = A19 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
- duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
-
else if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
!------------------------------------------------------------------------------
@@ -525,13 +418,6 @@
+ PML_duz_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dxl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdxl_x = A6 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,1) + A8 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,2)
- duydxl_y = A6 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,1) + A8 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,2)
- duzdxl_z = A6 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,1) + A8 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_z(i,j,k,ispec_CPML)
A10 = d_store_z(i,j,k,ispec_CPML)
@@ -560,13 +446,6 @@
+ PML_duz_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dyl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdyl_x = A9 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,1) + A11 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,2)
- duydyl_y = A9 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,1) + A11 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,2)
- duzdyl_z = A9 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,1) + A11 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = 1.0 / k_store_z(i,j,k,ispec_CPML)
A13 = 0.d0
@@ -595,13 +474,6 @@
rmemory_duz_dzl_z(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2) &
+ PML_duz_dzl_new(i,j,k,ispec_CPML) * coef1_2 + PML_duz_dzl(i,j,k,ispec_CPML) * coef2_2
- duxdzl_x = A12 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,1) + A14 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,2)
- duydzl_y = A12 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,1) + A14 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,2)
- duzdzl_z = A12 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,1) + A14 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A15 and A16 --------------------------
A15 = 1.d0
A16 = 0.d0
@@ -612,22 +484,12 @@
rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duz_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_y = A15 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_y(i,j,k,ispec_CPML,2)
- duzdyl_y = A15 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dyl_y(i,j,k,ispec_CPML,2)
-
rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duy_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duy_dyl_z(i,j,k,ispec_CPML,2) = 0.d0
- duydzl_z = A15 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dzl_z(i,j,k,ispec_CPML,2)
- duydyl_z = A15 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A17 and A18 --------------------------
A17 = 1.d0
A18 = 0.d0
@@ -638,22 +500,12 @@
rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duz_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_x = A17 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_x(i,j,k,ispec_CPML,2)
- duzdxl_x = A17 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dxl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdzl_z = A17 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dzl_z(i,j,k,ispec_CPML,2)
- duxdxl_z = A17 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A19 and A20 --------------------------
A19 = k_store_z(i,j,k,ispec_CPML)
A20 = d_store_z(i,j,k,ispec_CPML)
@@ -677,11 +529,6 @@
+ PML_duy_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duydyl_x = A19 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_x(i,j,k,ispec_CPML,2)
- duydxl_x = A19 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) &
+ PML_dux_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
@@ -690,11 +537,6 @@
+ PML_dux_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dxl_y(i,j,k,ispec_CPML,2) = 0.d0
- duxdyl_y = A19 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
- duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
-
else if( CPML_regions(ispec_CPML) == CPML_XY_ONLY ) then
!------------------------------------------------------------------------------
@@ -730,14 +572,7 @@
rmemory_duz_dxl_z(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,2) &
+ PML_duz_dxl_new(i,j,k,ispec_CPML) * coef1_2 + PML_duz_dxl(i,j,k,ispec_CPML) * coef2_2
- duxdxl_x = A6 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,1) + A8 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,2)
- duydxl_y = A6 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,1) + A8 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,2)
- duzdxl_z = A6 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,1) + A8 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_x(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML)
A10 = 0.d0
@@ -767,13 +602,6 @@
rmemory_duz_dyl_z(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,2) &
+ PML_duz_dyl_new(i,j,k,ispec_CPML) * coef1_2 + PML_duz_dyl(i,j,k,ispec_CPML) * coef2_2
- duxdyl_x = A9 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,1) + A11 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,2)
- duydyl_y = A9 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,1) + A11 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,2)
- duzdyl_z = A9 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,1) + A11 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML)
A13 = d_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
@@ -814,13 +642,6 @@
+ PML_duz_dzl_new(i,j,k,ispec_CPML) * it*deltat * coef1_2 &
+ PML_duz_dzl(i,j,k,ispec_CPML) * it*deltat * coef2_2
- duxdzl_x = A12 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,1) + A14 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,2)
- duydzl_y = A12 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,1) + A14 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,2)
- duzdzl_z = A12 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,1) + A14 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A15 and A16 --------------------------
A15 = k_store_x(i,j,k,ispec_CPML)
A16 = d_store_x(i,j,k,ispec_CPML)
@@ -844,11 +665,6 @@
+ PML_duz_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_y = A15 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_y(i,j,k,ispec_CPML,2)
- duzdyl_y = A15 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dyl_y(i,j,k,ispec_CPML,2)
-
rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) &
+ PML_duy_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
@@ -857,11 +673,6 @@
+ PML_duy_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dyl_z(i,j,k,ispec_CPML,2) = 0.d0
- duydzl_z = A15 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dzl_z(i,j,k,ispec_CPML,2)
- duydyl_z = A15 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A17 and A18 --------------------------
A17 = k_store_y(i,j,k,ispec_CPML)
A18 = d_store_y(i,j,k,ispec_CPML)
@@ -885,11 +696,6 @@
+ PML_duz_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_x = A17 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_x(i,j,k,ispec_CPML,2)
- duzdxl_x = A17 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) &
+ PML_dux_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
@@ -898,11 +704,6 @@
+ PML_dux_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dxl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdzl_z = A17 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dzl_z(i,j,k,ispec_CPML,2)
- duxdxl_z = A17 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A19 and A20--------------------------
A19 = 1.d0
A20 = 0.0
@@ -913,22 +714,12 @@
rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duy_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duydyl_x = A19 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_x(i,j,k,ispec_CPML,2)
- duydxl_x = A19 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dxl_y(i,j,k,ispec_CPML,2) = 0.d0
- duxdyl_y = A19 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
- duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
-
else if( CPML_regions(ispec_CPML) == CPML_XZ_ONLY ) then
!------------------------------------------------------------------------------
@@ -964,13 +755,6 @@
rmemory_duz_dxl_z(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,2) &
+ PML_duz_dxl_new(i,j,k,ispec_CPML) * coef1_2 + PML_duz_dxl(i,j,k,ispec_CPML) * coef2_2
- duxdxl_x = A6 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,1) + A8 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,2)
- duydxl_y = A6 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,1) + A8 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,2)
- duzdxl_z = A6 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,1) + A8 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML)
A10 = d_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) &
@@ -1011,13 +795,6 @@
+ PML_duz_dyl_new(i,j,k,ispec_CPML) * it*deltat * coef1_2 &
+ PML_duz_dyl(i,j,k,ispec_CPML) * it*deltat * coef2_2
- duxdyl_x = A9 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,1) + A11 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,2)
- duydyl_y = A9 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,1) + A11 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,2)
- duzdyl_z = A9 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,1) + A11 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_x(i,j,k,ispec_CPML) / k_store_z(i,j,k,ispec_CPML)
A13 = 0.d0
@@ -1047,13 +824,6 @@
rmemory_duz_dzl_z(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2) &
+ PML_duz_dzl_new(i,j,k,ispec_CPML) * coef1_2 + PML_duz_dzl(i,j,k,ispec_CPML) * coef2_2
- duxdzl_x = A12 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,1) + A14 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,2)
- duydzl_y = A12 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,1) + A14 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,2)
- duzdzl_z = A12 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,1) + A14 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A15 and A16 --------------------------
A15 = k_store_x(i,j,k,ispec_CPML)
A16 = d_store_x(i,j,k,ispec_CPML)
@@ -1077,11 +847,6 @@
+ PML_duz_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_y = A15 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_y(i,j,k,ispec_CPML,2)
- duzdyl_y = A15 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dyl_y(i,j,k,ispec_CPML,2)
-
rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) &
+ PML_duy_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
@@ -1090,11 +855,6 @@
+ PML_duy_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dyl_z(i,j,k,ispec_CPML,2) = 0.d0
- duydzl_z = A15 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dzl_z(i,j,k,ispec_CPML,2)
- duydyl_z = A15 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A17 and A18 --------------------------
A17 = 1.0d0
A18 = 0.d0
@@ -1105,22 +865,12 @@
rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duz_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_x = A17 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_x(i,j,k,ispec_CPML,2)
- duzdxl_x = A17 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_dux_dxl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdzl_z = A17 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dzl_z(i,j,k,ispec_CPML,2)
- duxdxl_z = A17 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A19 and A20 --------------------------
A19 = k_store_z(i,j,k,ispec_CPML)
A20 = d_store_z(i,j,k,ispec_CPML)
@@ -1144,11 +894,6 @@
+ PML_duy_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duydyl_x = A19 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_x(i,j,k,ispec_CPML,2)
- duydxl_x = A19 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) &
+ PML_dux_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
@@ -1157,11 +902,6 @@
+ PML_dux_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dxl_y(i,j,k,ispec_CPML,2) = 0.d0
- duxdyl_y = A19 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
- duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
-
else if( CPML_regions(ispec_CPML) == CPML_YZ_ONLY ) then
!------------------------------------------------------------------------------
@@ -1172,7 +912,7 @@
A6 = k_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML)
A7 = d_store_y(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) &
+ d_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
- + (it+0.0) * deltat * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
+ + it*deltat * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
A8 = - d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
bb = alpha_store(i,j,k,ispec_CPML)
@@ -1208,13 +948,6 @@
+ PML_duz_dxl_new(i,j,k,ispec_CPML) * it*deltat * coef1_2 &
+ PML_duz_dxl(i,j,k,ispec_CPML) * it*deltat * coef2_2
- duxdxl_x = A6 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,1) + A8 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,2)
- duydxl_y = A6 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,1) + A8 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,2)
- duzdxl_z = A6 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,1) + A8 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_z(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML)
A10 = 0.d0
@@ -1244,13 +977,6 @@
rmemory_duz_dyl_z(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,2) &
+ PML_duz_dyl_new(i,j,k,ispec_CPML) * coef1_2 + PML_duz_dyl(i,j,k,ispec_CPML) * coef2_2
- duxdyl_x = A9 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,1) + A11 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,2)
- duydyl_y = A9 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,1) + A11 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,2)
- duzdyl_z = A9 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,1) + A11 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_y(i,j,k,ispec_CPML) / k_store_z(i,j,k,ispec_CPML)
A13 = 0.d0
@@ -1280,13 +1006,6 @@
rmemory_duz_dzl_z(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2) &
+ PML_duz_dzl_new(i,j,k,ispec_CPML) * coef1_2 + PML_duz_dzl(i,j,k,ispec_CPML) * coef2_2
- duxdzl_x = A12 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,1) + A14 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,2)
- duydzl_y = A12 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,1) + A14 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,2)
- duzdzl_z = A12 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,1) + A14 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A15 and A16 --------------------------
A15 = 1.0d0
A16 = 0.0d0
@@ -1297,22 +1016,12 @@
rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duz_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_y = A15 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_y(i,j,k,ispec_CPML,2)
- duzdyl_y = A15 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dyl_y(i,j,k,ispec_CPML,2)
-
rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duy_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) = 0.d0
rmemory_duy_dyl_z(i,j,k,ispec_CPML,2) = 0.d0
- duydzl_z = A15 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dzl_z(i,j,k,ispec_CPML,2)
- duydyl_z = A15 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A17 and A18 --------------------------
A17 = k_store_y(i,j,k,ispec_CPML)
A18 = d_store_y(i,j,k,ispec_CPML)
@@ -1336,11 +1045,6 @@
+ PML_duz_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_x = A17 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_x(i,j,k,ispec_CPML,2)
- duzdxl_x = A17 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) &
+ PML_dux_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
@@ -1349,11 +1053,6 @@
+ PML_dux_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dxl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdzl_z = A17 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dzl_z(i,j,k,ispec_CPML,2)
- duxdxl_z = A17 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A19 and A20--------------------------
A19 = k_store_z(i,j,k,ispec_CPML)
A20 = d_store_z(i,j,k,ispec_CPML)
@@ -1377,11 +1076,6 @@
+ PML_duy_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duydyl_x = A19 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_x(i,j,k,ispec_CPML,2)
- duydxl_x = A19 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) &
+ PML_dux_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
@@ -1390,11 +1084,6 @@
+ PML_dux_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dxl_y(i,j,k,ispec_CPML,2) = 0.d0
- duxdyl_y = A19 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
- duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
-
else if( CPML_regions(ispec_CPML) == CPML_XYZ ) then
!------------------------------------------------------------------------------
@@ -1414,7 +1103,7 @@
A7 = ( d_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) + &
d_store_y(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) ) / &
k_store_x(i,j,k,ispec_CPML) + &
- it*deltat*d_store_y(i,j,k,ispec_CPML)*d_store_z(i,j,k,ispec_CPML)/k_store_x(i,j,k,ispec_CPML)
+ it*deltat * d_store_y(i,j,k,ispec_CPML)*d_store_z(i,j,k,ispec_CPML)/k_store_x(i,j,k,ispec_CPML)
A8 = - d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) / k_store_x(i,j,k,ispec_CPML)
endif
@@ -1466,14 +1155,6 @@
+ PML_duz_dxl(i,j,k,ispec_CPML) * it*deltat * coef2_2
endif
- duxdxl_x = A6 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,1) + A8 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,2)
- duydxl_y = A6 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,1) + A8 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,2)
- duzdxl_z = A6 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,1) + A8 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,2)
-
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML)
if( abs(d_store_y(i,j,k,ispec_CPML)) > 1.d-5 ) then
@@ -1487,7 +1168,7 @@
A10 = ( d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
+ d_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) ) / &
k_store_y(i,j,k,ispec_CPML) + &
- it*deltat*d_store_x(i,j,k,ispec_CPML)*d_store_z(i,j,k,ispec_CPML)/k_store_y(i,j,k,ispec_CPML)
+ it*deltat * d_store_x(i,j,k,ispec_CPML)*d_store_z(i,j,k,ispec_CPML)/k_store_y(i,j,k,ispec_CPML)
A11 = - d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML)
endif
@@ -1539,13 +1220,6 @@
+ PML_duz_dyl(i,j,k,ispec_CPML) * it*deltat * coef2_2
endif
- duxdyl_x = A9 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,1) + A11 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,2)
- duydyl_y = A9 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,1) + A11 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,2)
- duzdyl_z = A9 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,1) + A11 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) / k_store_z(i,j,k,ispec_CPML)
if( abs(d_store_z(i,j,k,ispec_CPML)) > 1.d-5 ) then
@@ -1559,7 +1233,7 @@
A13 = ( d_store_y(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
+ d_store_x(i,j,k,ispec_CPML)*k_store_y(i,j,k,ispec_CPML) ) / &
k_store_z(i,j,k,ispec_CPML) + &
- it*deltat*d_store_x(i,j,k,ispec_CPML)*d_store_y(i,j,k,ispec_CPML)/k_store_z(i,j,k,ispec_CPML)
+ it*deltat * d_store_x(i,j,k,ispec_CPML)*d_store_y(i,j,k,ispec_CPML)/k_store_z(i,j,k,ispec_CPML)
A14 = - d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) / k_store_z(i,j,k,ispec_CPML)
endif
@@ -1612,13 +1286,6 @@
+ PML_duz_dzl(i,j,k,ispec_CPML) * it*deltat * coef2_2
endif
- duxdzl_x = A12 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,1) + A14 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,2)
- duydzl_y = A12 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,1) + A14 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,2)
- duzdzl_z = A12 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,1) + A14 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A15 and A16 --------------------------
A15 = k_store_x(i,j,k,ispec_CPML)
A16 = d_store_x(i,j,k,ispec_CPML)
@@ -1642,11 +1309,6 @@
+ PML_duz_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_y = A15 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_y(i,j,k,ispec_CPML,2)
- duzdyl_y = A15 * PML_duz_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dyl_y(i,j,k,ispec_CPML,2)
-
rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) &
+ PML_duy_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
@@ -1655,11 +1317,6 @@
+ PML_duy_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dyl_z(i,j,k,ispec_CPML,2) = 0.d0
- duydzl_z = A15 * PML_duy_dzl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dzl_z(i,j,k,ispec_CPML,2)
- duydyl_z = A15 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A16 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A17 and A18 --------------------------
A17 = k_store_y(i,j,k,ispec_CPML)
A18 = d_store_y(i,j,k,ispec_CPML)
@@ -1683,11 +1340,6 @@
+ PML_duz_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duz_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_duz_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duzdzl_x = A17 * PML_duz_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_x(i,j,k,ispec_CPML,2)
- duzdxl_x = A17 * PML_duz_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) &
+ PML_dux_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dzl_z(i,j,k,ispec_CPML,2) = 0.d0
@@ -1696,11 +1348,6 @@
+ PML_dux_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dxl_z(i,j,k,ispec_CPML,2) = 0.d0
- duxdzl_z = A17 * PML_dux_dzl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dzl_z(i,j,k,ispec_CPML,2)
- duxdxl_z = A17 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A18 * rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_z(i,j,k,ispec_CPML,2)
-
!---------------------- A19 and A20 --------------------------
A19 = k_store_z(i,j,k,ispec_CPML)
A20 = d_store_z(i,j,k,ispec_CPML)
@@ -1724,11 +1371,6 @@
+ PML_duy_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_duy_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_duy_dxl_x(i,j,k,ispec_CPML,2) = 0.d0
- duydyl_x = A19 * PML_duy_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_x(i,j,k,ispec_CPML,2)
- duydxl_x = A19 * PML_duy_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dxl_x(i,j,k,ispec_CPML,2)
-
rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) &
+ PML_dux_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dyl_y(i,j,k,ispec_CPML,2) = 0.d0
@@ -1737,15 +1379,55 @@
+ PML_dux_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dux_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_dux_dxl_y(i,j,k,ispec_CPML,2) = 0.d0
- duxdyl_y = A19 * PML_dux_dyl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
- duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
- + A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
-
else
stop 'wrong PML flag in PML memory variable calculation routine'
endif
+ duxdxl_x = A6 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ + A7 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,1) + A8 * rmemory_dux_dxl_x(i,j,k,ispec_CPML,2)
+ duxdyl_x = A9 * PML_dux_dyl(i,j,k,ispec_CPML) &
+ + A10 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,1) + A11 * rmemory_dux_dyl_x(i,j,k,ispec_CPML,2)
+ duxdzl_x = A12 * PML_dux_dzl(i,j,k,ispec_CPML) &
+ + A13 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,1) + A14 * rmemory_dux_dzl_x(i,j,k,ispec_CPML,2)
+ duzdzl_x = A17 * PML_duz_dzl(i,j,k,ispec_CPML) &
+ + A18 * rmemory_duz_dzl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_x(i,j,k,ispec_CPML,2)
+ duzdxl_x = A17 * PML_duz_dxl(i,j,k,ispec_CPML) &
+ + A18 * rmemory_duz_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duz_dxl_x(i,j,k,ispec_CPML,2)
+ duydyl_x = A19 * PML_duy_dyl(i,j,k,ispec_CPML) &
+ + A20 * rmemory_duy_dyl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_x(i,j,k,ispec_CPML,2)
+ duydxl_x = A19 * PML_duy_dxl(i,j,k,ispec_CPML) &
+ + A20 * rmemory_duy_dxl_x(i,j,k,ispec_CPML,1) + rmemory_duy_dxl_x(i,j,k,ispec_CPML,2)
+
+ duydxl_y = A6 * PML_duy_dxl(i,j,k,ispec_CPML) &
+ + A7 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,1) + A8 * rmemory_duy_dxl_y(i,j,k,ispec_CPML,2)
+ duydyl_y = A9 * PML_duy_dyl(i,j,k,ispec_CPML) &
+ + A10 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,1) + A11 * rmemory_duy_dyl_y(i,j,k,ispec_CPML,2)
+ duydzl_y = A12 * PML_duy_dzl(i,j,k,ispec_CPML) &
+ + A13 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,1) + A14 * rmemory_duy_dzl_y(i,j,k,ispec_CPML,2)
+ duzdzl_y = A15 * PML_duz_dzl(i,j,k,ispec_CPML) &
+ + A16 * rmemory_duz_dzl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dzl_y(i,j,k,ispec_CPML,2)
+ duzdyl_y = A15 * PML_duz_dyl(i,j,k,ispec_CPML) &
+ + A16 * rmemory_duz_dyl_y(i,j,k,ispec_CPML,1) + rmemory_duz_dyl_y(i,j,k,ispec_CPML,2)
+ duxdyl_y = A19 * PML_dux_dyl(i,j,k,ispec_CPML) &
+ + A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
+ duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ + A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
+
+ duzdxl_z = A6 * PML_duz_dxl(i,j,k,ispec_CPML) &
+ + A7 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,1) + A8 * rmemory_duz_dxl_z(i,j,k,ispec_CPML,2)
+ duzdyl_z = A9 * PML_duz_dyl(i,j,k,ispec_CPML) &
+ + A10 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,1) + A11 * rmemory_duz_dyl_z(i,j,k,ispec_CPML,2)
+ duzdzl_z = A12 * PML_duz_dzl(i,j,k,ispec_CPML) &
+ + A13 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,1) + A14 * rmemory_duz_dzl_z(i,j,k,ispec_CPML,2)
+ duydzl_z = A15 * PML_duy_dzl(i,j,k,ispec_CPML) &
+ + A16 * rmemory_duy_dzl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dzl_z(i,j,k,ispec_CPML,2)
+ duydyl_z = A15 * PML_duy_dyl(i,j,k,ispec_CPML) &
+ + A16 * rmemory_duy_dyl_z(i,j,k,ispec_CPML,1) + rmemory_duy_dyl_z(i,j,k,ispec_CPML,2)
+ duxdzl_z = A17 * PML_dux_dzl(i,j,k,ispec_CPML) &
+ + A18 * rmemory_dux_dzl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dzl_z(i,j,k,ispec_CPML,2)
+ duxdxl_z = A17 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ + A18 * rmemory_dux_dxl_z(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_z(i,j,k,ispec_CPML,2)
+
! compute stress sigma
sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
sigma_yx = mul*duxdyl_x + mul*duydxl_x
@@ -1856,9 +1538,6 @@
rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) &
+ PML_dpotential_dxl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dxl(i,j,k,ispec_CPML) * coef2_2
- dpotentialdxl = A6 * PML_dpotential_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + A8 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_x(i,j,k,ispec_CPML)
A10 = d_store_x(i,j,k,ispec_CPML)
@@ -1879,9 +1558,6 @@
+ PML_dpotential_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dpotential_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) = 0.d0
- dpotentialdyl = A9 * PML_dpotential_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,1) + A11 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_x(i,j,k,ispec_CPML)
A13 = d_store_x(i,j,k,ispec_CPML)
@@ -1902,9 +1578,6 @@
+ PML_dpotential_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dpotential_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) = 0.d0
- dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
-
else if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- Y-surface C-PML ---------------------------------
@@ -1930,9 +1603,6 @@
+ PML_dpotential_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dpotential_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) = 0.d0
- dpotentialdxl = A6 * PML_dpotential_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + A8 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = 1.d0/k_store_y(i,j,k,ispec_CPML)
A10 = 0.d0
@@ -1953,9 +1623,6 @@
rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) &
+ PML_dpotential_dyl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dyl(i,j,k,ispec_CPML) * coef2_2
- dpotentialdyl = A9 * PML_dpotential_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,1) + A11 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_y(i,j,k,ispec_CPML)
A13 = d_store_y(i,j,k,ispec_CPML)
@@ -1976,9 +1643,6 @@
+ PML_dpotential_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dpotential_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) = 0.d0
- dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
-
else if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
!------------------------------------------------------------------------------
@@ -2005,9 +1669,6 @@
+ PML_dpotential_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dpotential_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) = 0.d0
- dpotentialdxl = A6 * PML_dpotential_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + A8 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_z(i,j,k,ispec_CPML)
A10 = d_store_z(i,j,k,ispec_CPML)
@@ -2028,9 +1689,6 @@
+ PML_dpotential_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dpotential_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) = 0.d0
- dpotentialdyl = A9 * PML_dpotential_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,1) + A11 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = 1.0 / k_store_z(i,j,k,ispec_CPML)
A13 = 0.d0
@@ -2051,9 +1709,6 @@
rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) &
+ PML_dpotential_dzl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dzl(i,j,k,ispec_CPML) * coef2_2
- dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
-
else if( CPML_regions(ispec_CPML) == CPML_XY_ONLY ) then
!------------------------------------------------------------------------------
@@ -2081,9 +1736,6 @@
rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) &
+ PML_dpotential_dxl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dxl(i,j,k,ispec_CPML) * coef2_2
- dpotentialdxl = A6 * PML_dpotential_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + A8 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_x(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML)
A10 = 0.d0
@@ -2105,9 +1757,6 @@
rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) &
+ PML_dpotential_dyl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dyl(i,j,k,ispec_CPML) * coef2_2
- dpotentialdyl = A9 * PML_dpotential_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,1) + A11 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML)
A13 = d_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
@@ -2133,12 +1782,9 @@
rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) &
+ PML_dpotential_dzl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dpotential_dzl(i,j,k,ispec_CPML) * coef2_1
rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) &
- + PML_dpotential_dzl_new(i,j,k,ispec_CPML) *it*deltat * coef1_2 &
- + PML_dpotential_dzl(i,j,k,ispec_CPML) *it*deltat * coef2_2
+ + PML_dpotential_dzl_new(i,j,k,ispec_CPML) * it*deltat * coef1_2 &
+ + PML_dpotential_dzl(i,j,k,ispec_CPML) * it*deltat * coef2_2
- dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
-
else if( CPML_regions(ispec_CPML) == CPML_XZ_ONLY ) then
!------------------------------------------------------------------------------
@@ -2166,9 +1812,6 @@
rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) &
+ PML_dpotential_dxl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dxl(i,j,k,ispec_CPML) * coef2_2
- dpotentialdxl = A6 * PML_dpotential_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + A8 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML)
A10 = d_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) &
@@ -2194,12 +1837,9 @@
rmemory_dpotential_dyl(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,1) &
+ PML_dpotential_dyl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dpotential_dyl(i,j,k,ispec_CPML) * coef2_1
rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) &
- + PML_dpotential_dyl_new(i,j,k,ispec_CPML) *it*deltat* coef1_2 &
- + PML_dpotential_dyl(i,j,k,ispec_CPML) *it*deltat* coef2_2
+ + PML_dpotential_dyl_new(i,j,k,ispec_CPML) * it*deltat * coef1_2 &
+ + PML_dpotential_dyl(i,j,k,ispec_CPML) * it*deltat * coef2_2
- dpotentialdyl = A9 * PML_dpotential_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,1) + A11 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_x(i,j,k,ispec_CPML) / k_store_z(i,j,k,ispec_CPML)
A13 = 0.d0
@@ -2221,10 +1861,6 @@
rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) &
+ PML_dpotential_dzl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dzl(i,j,k,ispec_CPML) * coef2_2
- dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
-
-
else if( CPML_regions(ispec_CPML) == CPML_YZ_ONLY ) then
!------------------------------------------------------------------------------
@@ -2256,12 +1892,9 @@
rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) &
+ PML_dpotential_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dpotential_dxl(i,j,k,ispec_CPML) * coef2_1
rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) &
- + PML_dpotential_dxl_new(i,j,k,ispec_CPML) * it*deltat* coef1_2 &
- + PML_dpotential_dxl(i,j,k,ispec_CPML) *it*deltat* coef2_2
+ + PML_dpotential_dxl_new(i,j,k,ispec_CPML) * it*deltat * coef1_2 &
+ + PML_dpotential_dxl(i,j,k,ispec_CPML) * it*deltat * coef2_2
- dpotentialdxl = A6 * PML_dpotential_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + A8 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_z(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML)
A10 = 0.d0
@@ -2283,9 +1916,6 @@
rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) &
+ PML_dpotential_dyl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dyl(i,j,k,ispec_CPML) * coef2_2
- dpotentialdyl = A9 * PML_dpotential_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,1) + A11 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_y(i,j,k,ispec_CPML) / k_store_z(i,j,k,ispec_CPML)
A13 = 0.d0
@@ -2307,9 +1937,6 @@
rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) &
+ PML_dpotential_dzl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dzl(i,j,k,ispec_CPML) * coef2_2
- dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
-
else if( CPML_regions(ispec_CPML) == CPML_XYZ ) then
!------------------------------------------------------------------------------
@@ -2329,7 +1956,7 @@
A7 = (d_store_z(i,j,k,ispec_CPML)*k_store_y(i,j,k,ispec_CPML)+ &
d_store_y(i,j,k,ispec_CPML)*k_store_z(i,j,k,ispec_CPML)) / &
k_store_x(i,j,k,ispec_CPML) + &
- it*deltat*d_store_y(i,j,k,ispec_CPML)*d_store_z(i,j,k,ispec_CPML)/k_store_x(i,j,k,ispec_CPML)
+ it*deltat * d_store_y(i,j,k,ispec_CPML)*d_store_z(i,j,k,ispec_CPML)/k_store_x(i,j,k,ispec_CPML)
A8 = - d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) / k_store_x(i,j,k,ispec_CPML)
endif
@@ -2363,13 +1990,10 @@
+ PML_dpotential_dxl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dxl(i,j,k,ispec_CPML) * coef2_2
else
rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) &
- + PML_dpotential_dxl_new(i,j,k,ispec_CPML) * it*deltat* coef1_2 &
- + PML_dpotential_dxl(i,j,k,ispec_CPML) * it*deltat* coef2_2
+ + PML_dpotential_dxl_new(i,j,k,ispec_CPML) * it*deltat * coef1_2 &
+ + PML_dpotential_dxl(i,j,k,ispec_CPML) * it*deltat * coef2_2
endif
- dpotentialdxl = A6 * PML_dpotential_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + A8 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2)
-
!---------------------- A9, A10 and A11 --------------------------
A9 = k_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML)
if( abs(d_store_y(i,j,k,ispec_CPML)) > 1.d-5 ) then
@@ -2383,7 +2007,7 @@
A10 = (d_store_z(i,j,k,ispec_CPML)*k_store_x(i,j,k,ispec_CPML) &
+d_store_x(i,j,k,ispec_CPML)*k_store_z(i,j,k,ispec_CPML)) / &
k_store_y(i,j,k,ispec_CPML) + &
- it*deltat*d_store_x(i,j,k,ispec_CPML)*d_store_z(i,j,k,ispec_CPML)/k_store_y(i,j,k,ispec_CPML)
+ it*deltat * d_store_x(i,j,k,ispec_CPML)*d_store_z(i,j,k,ispec_CPML)/k_store_y(i,j,k,ispec_CPML)
A11 = - d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML) / k_store_y(i,j,k,ispec_CPML)
endif
@@ -2417,13 +2041,10 @@
+ PML_dpotential_dyl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dyl(i,j,k,ispec_CPML) * coef2_2
else
rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2) &
- + PML_dpotential_dyl_new(i,j,k,ispec_CPML) * it*deltat* coef1_2 &
- + PML_dpotential_dyl(i,j,k,ispec_CPML) * it*deltat* coef2_2
+ + PML_dpotential_dyl_new(i,j,k,ispec_CPML) * it*deltat * coef1_2 &
+ + PML_dpotential_dyl(i,j,k,ispec_CPML) * it*deltat * coef2_2
endif
- dpotentialdyl = A9 * PML_dpotential_dyl(i,j,k,ispec_CPML) &
- + A10 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,1) + A11 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2)
-
!---------------------- A12, A13 and A14 --------------------------
A12 = k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) / k_store_z(i,j,k,ispec_CPML)
if( abs(d_store_z(i,j,k,ispec_CPML)) > 1.d-5 ) then
@@ -2437,7 +2058,7 @@
A13 = (d_store_y(i,j,k,ispec_CPML)*k_store_x(i,j,k,ispec_CPML)&
+d_store_x(i,j,k,ispec_CPML)*k_store_y(i,j,k,ispec_CPML)) / &
k_store_z(i,j,k,ispec_CPML) + &
- it*deltat*d_store_x(i,j,k,ispec_CPML)*d_store_y(i,j,k,ispec_CPML)/k_store_z(i,j,k,ispec_CPML)
+ it*deltat * d_store_x(i,j,k,ispec_CPML)*d_store_y(i,j,k,ispec_CPML)/k_store_z(i,j,k,ispec_CPML)
A14 = - d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML) / k_store_z(i,j,k,ispec_CPML)
endif
@@ -2471,17 +2092,20 @@
+ PML_dpotential_dzl_new(i,j,k,ispec_CPML) * coef1_2 + PML_dpotential_dzl(i,j,k,ispec_CPML) * coef2_2
else
rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) = coef0_2 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2) &
- + PML_dpotential_dzl_new(i,j,k,ispec_CPML) * it*deltat* coef1_2 &
- + PML_dpotential_dzl(i,j,k,ispec_CPML) * it*deltat* coef2_2
+ + PML_dpotential_dzl_new(i,j,k,ispec_CPML) * it*deltat * coef1_2 &
+ + PML_dpotential_dzl(i,j,k,ispec_CPML) * it*deltat * coef2_2
endif
- dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
- + A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
-
else
stop 'wrong PML flag in PML memory variable calculation routine'
endif
+ dpotentialdxl = A6 * PML_dpotential_dxl(i,j,k,ispec_CPML) &
+ + A7 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + A8 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2)
+ dpotentialdyl = A9 * PML_dpotential_dyl(i,j,k,ispec_CPML) &
+ + A10 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,1) + A11 * rmemory_dpotential_dyl(i,j,k,ispec_CPML,2)
+ dpotentialdzl = A12 * PML_dpotential_dzl(i,j,k,ispec_CPML) &
+ + A13 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,1) + A14 * rmemory_dpotential_dzl(i,j,k,ispec_CPML,2)
temp1(i,j,k) = rhoin_jacob_jk * (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
temp2(i,j,k) = rhoin_jacob_ik * (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
temp3(i,j,k) = rhoin_jacob_ij * (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
@@ -2492,3 +2116,582 @@
end subroutine pml_compute_memory_variables_acoustic
+!
+!=====================================================================
+!
+!
+
+subroutine pml_compute_memory_variables_acoustic_elastic(ispec_CPML,iface,iglob,i,j,k,&
+ displ_x,displ_y,displ_z)
+ ! calculates C-PML elastic memory variables and computes stress sigma
+
+ ! second-order accurate convolution term calculation from equation (21) of
+ ! Shumin Wang, Robert Lee, and Fernando L. Teixeira,
+ ! Anisotropic-Medium PML for Vector FETD With Modified Basis Functions,
+ ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006)
+
+ use specfem_par, only: it,deltat
+ use specfem_par_elastic, only: displ,veloc
+ use pml_par
+ use constants, only: CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ,NDIM
+
+ implicit none
+
+ integer, intent(in) :: ispec_CPML,iface,iglob
+
+ ! local parameters
+ integer :: i,j,k
+ real(kind=CUSTOM_REAL) :: bb,coef0_1,coef1_1,coef2_1,coef0_2,coef1_2,coef2_2
+ real(kind=CUSTOM_REAL) :: A6,A7,A8,A9,A10,A11,A12,A13,A14
+ real(kind=CUSTOM_REAL) :: displ_x,displ_y,displ_z
+
+
+ if( CPML_regions(ispec_CPML) == CPML_X_ONLY ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- X-surface C-PML ---------------------------------
+ !------------------------------------------------------------------------------
+
+ ! displ_x
+ A6 = 1.d0
+ A7 = 0.d0
+ A8 = 0.d0
+
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) = 0.d0
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = 0.d0
+
+ displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) + &
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+
+ ! displ_y
+ A9 = k_store_x(i,j,k,ispec_CPML)
+ A10 = d_store_x(i,j,k,ispec_CPML)
+ A11 = 0.d0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * coef1_1 + (displ(2,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = 0.d0
+
+ displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) + &
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+
+ ! displ_z
+ A12 = k_store_x(i,j,k,ispec_CPML)
+ A13 = d_store_x(i,j,k,ispec_CPML)
+ A14 = 0.d0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) &
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * coef1_1 + (displ(3,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,2) = 0.d0
+
+ displ_z = A12 * displ(3,iglob) + A13 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) + &
+ + A14 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
+ !------------------------------------------------------------------------------
+ !---------------------------- Y-surface C-PML ---------------------------------
+ !------------------------------------------------------------------------------
+
+ ! displ_x
+ A6 = k_store_y(i,j,k,ispec_CPML)
+ A7 = d_store_y(i,j,k,ispec_CPML)
+ A8 = 0.0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_1 + (displ(1,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = 0.d0
+
+ displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) + &
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+
+ ! displ_y
+ A9 = 1.d0
+ A10 = 0.d0
+ A11 = 0.d0
+
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) = 0.d0
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = 0.d0
+
+ displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) + &
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+
+ ! displ_z
+ A12 = k_store_y(i,j,k,ispec_CPML)
+ A13 = d_store_y(i,j,k,ispec_CPML)
+ A14 = 0.0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) &
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * coef1_1 + (displ(3,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,2) = 0.d0
+
+ displ_z = A12 * displ(3,iglob) + A13 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) + &
+ + A14 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- Z-surface C-PML ---------------------------------
+ !------------------------------------------------------------------------------
+
+ ! displ_x
+ A6 = k_store_z(i,j,k,ispec_CPML)
+ A7 = d_store_z(i,j,k,ispec_CPML)
+ A8 = 0.0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_1 + (displ(1,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = 0.d0
+
+ displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) + &
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+
+ ! displ_y
+ A9 = k_store_z(i,j,k,ispec_CPML)
+ A10 = d_store_z(i,j,k,ispec_CPML)
+ A11 = 0.0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * coef1_1 + (displ(2,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = 0.d0
+
+ displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) + &
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+
+ ! displ_z
+ A12 = 1.d0
+ A13 = 0.d0
+ A14 = 0.d0
+
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) = 0.d0
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,2) = 0.d0
+
+ displ_z = A12 * displ(3,iglob) + A13 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) + &
+ + A14 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_XY_ONLY ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- XY-edge C-PML -----------------------------------
+ !------------------------------------------------------------------------------
+ ! displ_x
+ A6 = k_store_y(i,j,k,ispec_CPML)
+ A7 = d_store_y(i,j,k,ispec_CPML)
+ A8 = 0.0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_1 + (displ(1,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = 0.d0
+
+ displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) + &
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+
+ ! displ_y
+ A9 = k_store_x(i,j,k,ispec_CPML)
+ A10 = d_store_x(i,j,k,ispec_CPML)
+ A11 = 0.d0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * coef1_1 + (displ(2,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = 0.d0
+
+ displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) + &
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+
+ ! displ_z
+ A12 = k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML)
+ A13 = d_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
+ + d_store_y(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
+ + it*deltat * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
+ A14 = - d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ coef0_2 = coef0_1
+ coef1_2 = coef1_1
+ coef2_2 = coef2_1
+
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) &
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * coef1_1 + (displ(3,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,2) = coef0_2 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2) &
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * it*deltat * coef1_2 &
+ + (displ(3,iglob)) * it*deltat * coef2_2
+
+ displ_z = A12 * displ(3,iglob) + A13 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) + &
+ + A14 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_XZ_ONLY ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- XZ-edge C-PML -----------------------------------
+ !------------------------------------------------------------------------------
+
+ ! displ_x
+ A6 = k_store_z(i,j,k,ispec_CPML)
+ A7 = d_store_z(i,j,k,ispec_CPML)
+ A8 = 0.0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_1 + (displ(1,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = 0.d0
+
+ displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) + &
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+
+ ! displ_y
+ A9 = k_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML)
+ A10 = d_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) &
+ + d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
+ + it*deltat * d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
+ A11 = - d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ coef0_2 = coef0_1
+ coef1_2 = coef1_1
+ coef2_2 = coef2_1
+
+
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * coef1_1 + (displ(2,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = coef0_2 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) &
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * it*deltat * coef1_2 &
+ + (displ(2,iglob)) * it*deltat * coef2_2
+
+ displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) + &
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+
+ ! displ_z
+ A12 = k_store_x(i,j,k,ispec_CPML)
+ A13 = d_store_x(i,j,k,ispec_CPML)
+ A14 = 0.0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) &
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * coef1_1 + (displ(3,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,2) = 0.d0
+
+ displ_z = A12 * displ(3,iglob) + A13 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) + &
+ + A14 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_YZ_ONLY ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- YZ-edge C-PML -----------------------------------
+ !------------------------------------------------------------------------------
+
+ ! displ_x
+ A6 = k_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML)
+ A7 = d_store_y(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) &
+ + d_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
+ + it*deltat * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
+ A8 = - d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ coef0_2 = coef0_1
+ coef1_2 = coef1_1
+ coef2_2 = coef2_1
+
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_1 + (displ(1,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = coef0_2 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) &
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_2 + (displ(1,iglob)) * coef2_2
+
+ displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) + &
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+
+ ! displ_y
+ A9 = k_store_z(i,j,k,ispec_CPML)
+ A10 = d_store_z(i,j,k,ispec_CPML)
+ A11 = 0.d0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * coef1_1 + (displ(2,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = 0.d0
+
+ displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) + &
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+
+ ! displ_z
+ A12 = k_store_y(i,j,k,ispec_CPML)
+ A13 = d_store_y(i,j,k,ispec_CPML)
+ A14 = 0.0
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) &
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * coef1_1 + (displ(3,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,2) = 0.d0
+
+ displ_z = A12 * displ(3,iglob) + A13 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) + &
+ + A14 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2)
+
+
+ else if( CPML_regions(ispec_CPML) == CPML_XYZ ) then
+
+ !------------------------------------------------------------------------------
+ !---------------------------- XYZ-corner C-PML --------------------------------
+ !------------------------------------------------------------------------------
+ ! displ_x
+ A6 = k_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML)
+ A7 = d_store_y(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) &
+ + d_store_z(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
+ + it*deltat * d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
+ A8 = - d_store_y(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ coef0_2 = coef0_1
+ coef1_2 = coef1_1
+ coef2_2 = coef2_1
+
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_1 + (displ(1,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = coef0_2 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) &
+ + (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_2 + (displ(1,iglob)) * coef2_2
+
+ displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) + &
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+
+ ! displ_y
+ A9 = k_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML)
+ A10 = d_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML) &
+ + d_store_z(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
+ + it*deltat * d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
+ A11 = - d_store_x(i,j,k,ispec_CPML) * d_store_z(i,j,k,ispec_CPML)
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ coef0_2 = coef0_1
+ coef1_2 = coef1_1
+ coef2_2 = coef2_1
+
+
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * coef1_1 + (displ(2,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = coef0_2 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) &
+ + (displ(2,iglob) + deltat * veloc(2,iglob)) * it*deltat * coef1_2 &
+ + (displ(2,iglob)) * it*deltat * coef2_2
+
+ displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) + &
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+
+ ! displ_z
+ A12 = k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML)
+ A13 = d_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML) &
+ + d_store_y(i,j,k,ispec_CPML) * k_store_x(i,j,k,ispec_CPML) &
+ + it*deltat * d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
+ A14 = - d_store_x(i,j,k,ispec_CPML) * d_store_y(i,j,k,ispec_CPML)
+
+ bb = alpha_store(i,j,k,ispec_CPML)
+ coef0_1 = exp(-bb * deltat)
+
+ if( abs(bb) > 1.d-5 ) then
+ coef1_1 = (1.d0 - exp(-bb * deltat/2.d0)) / bb
+ coef2_1 = (1.d0 - exp(-bb * deltat/2.d0)) * exp(-bb * deltat/2.d0) / bb
+ else
+ coef1_1 = deltat/2.0d0
+ coef2_1 = deltat/2.0d0
+ endif
+
+ coef0_2 = coef0_1
+ coef1_2 = coef1_1
+ coef2_2 = coef2_1
+
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) = coef0_1 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) &
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * coef1_1 + (displ(3,iglob)) * coef2_1
+ rmemory_coupling_ac_el_displ(3,i,j,k,iface,2) = coef0_2 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2) &
+ + (displ(3,iglob) + deltat * veloc(3,iglob)) * it*deltat * coef1_2 &
+ + (displ(3,iglob)) * it*deltat * coef2_2
+
+ displ_z = A12 * displ(3,iglob) + A13 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) + &
+ + A14 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2)
+ else
+ stop 'wrong PML flag in PML memory variable calculation routine'
+ endif
+
+end subroutine pml_compute_memory_variables_acoustic_elastic
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-04-27 22:59:14 UTC (rev 21958)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-04-29 12:21:01 UTC (rev 21959)
@@ -53,7 +53,6 @@
logical, dimension(:), allocatable :: is_CPML
! thickness of C-PML layers
-!ZN real(CUSTOM_REAL) :: CPML_width,CPML_width_x,CPML_width_y,CPML_width_z
real(CUSTOM_REAL) :: CPML_width_x,CPML_width_y,CPML_width_z
! C-PML damping profile arrays
@@ -106,4 +105,7 @@
! C-PML contribution to update the second derivative of the potential to the global mesh
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: potential_dot_dot_acoustic_CPML
+ ! C-PML contribution to update displacement on elastic/acoustic interface
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:,:), allocatable :: rmemory_coupling_ac_el_displ
+
end module pml_par
More information about the CIG-COMMITS
mailing list