[cig-commits] r21227 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Mon Jan 14 04:52:32 PST 2013
Author: xie.zhinan
Date: 2013-01-14 04:52:31 -0800 (Mon, 14 Jan 2013)
New Revision: 21227
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
Log:
solve data type inconsistency when set CUSTOM_REAL=real in elastic/acoustic simulation with PML
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2013-01-14 10:04:45 UTC (rev 21226)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2013-01-14 12:52:31 UTC (rev 21227)
@@ -141,9 +141,10 @@
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: potential_dot_dot_acoustic_PML
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) ::PML_dux_dxl,PML_dux_dzl,&
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: PML_dux_dxl,PML_dux_dzl,&
PML_dux_dxl_new,PML_dux_dzl_new
- real(kind=CUSTOM_REAL) :: coef0, coef1, coef2,bb,deltat
+ real(kind=CUSTOM_REAL) :: coef0, coef1, coef2,bb
+ double precision :: deltat
real(kind=CUSTOM_REAL) :: A0, A1, A2, A3, A4, A5, A6, A7, A8
logical :: PML_BOUNDARY_CONDITIONS
@@ -166,16 +167,19 @@
PML_dux_dzl_new = 0._CUSTOM_REAL
endif
-! loop over spectral elements
+! loop over spectral elementsbb
do ispec = ifirstelem,ilastelem
!---
!--- acoustic spectral element
!---
if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rhol = sngl(density(1,kmato(ispec)))
+ else
+ rhol = density(1,kmato(ispec))
+ endif
- rhol = density(1,kmato(ispec))
-
! first double loop over GLL points to compute and store gradients
do j = 1,NGLLZ
do i = 1,NGLLX
@@ -201,7 +205,7 @@
dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
- ! derivative along x and along z
+ ! derivative along x and along zbb
if(PML_BOUNDARY_CONDITIONS .and. is_PML(ispec))then
ispec_PML=spec_to_PML(ispec)
@@ -224,7 +228,6 @@
xizl = xiz(i,j,ispec)
gammaxl = gammax(i,j,ispec)
gammazl = gammaz(i,j,ispec)
-
! derivatives of potential
dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
@@ -247,12 +250,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(-bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
- coef1 = (1.d0 - exp(-bb * deltat / 2.d0)) / bb
- coef2 = (1.d0 - exp(-bb* deltat / 2.d0)) * exp(-bb * deltat / 2.d0) / bb
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
+ coef1 = (1._CUSTOM_REAL - exp(-bb * deltat / 2._CUSTOM_REAL)) / bb
+ coef2 = (1._CUSTOM_REAL - exp(-bb* deltat / 2._CUSTOM_REAL)) * exp(-bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2.0_CUSTOM_REAL
+ coef2 = deltat / 2.0_CUSTOM_REAL
end if
rmemory_acoustic_dux_dx(i,j,ispec_PML) = coef0*rmemory_acoustic_dux_dx(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) * coef1 + PML_dux_dxl(i,j,ispec_PML) * coef2
@@ -276,12 +279,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
- if ( abs( bb ) > 1.d-3) then
- coef1 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) / bb
- coef2 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) / bb
+ if ( abs( bb ) > 0.001_CUSTOM_REAL) then
+ coef1 = (1.0_CUSTOM_REAL - exp(- bb * deltat / 2.0_CUSTOM_REAL) ) / bb
+ coef2 = (1.0_CUSTOM_REAL - exp(- bb * deltat / 2.0_CUSTOM_REAL) ) * exp(- bb * deltat / 2.0_CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2.0_CUSTOM_REAL
+ coef2 = deltat / 2.0_CUSTOM_REAL
end if
rmemory_acoustic_dux_dz(i,j,ispec_PML) = coef0 * rmemory_acoustic_dux_dz(i,j,ispec_PML) &
+ PML_dux_dzl_new(i,j,ispec_PML) *coef1 + PML_dux_dzl(i,j,ispec_PML) * coef2
@@ -313,12 +316,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
- coef1 = ( 1.d0 - exp(- bb * deltat / 2.d0) ) / bb
- coef2 = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * exp(- bb * deltat / 2.d0) / bb
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
+ coef1 = ( 1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) / bb
+ coef2 = ( 1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) * exp(- bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2.0_CUSTOM_REAL
+ coef2 = deltat / 2.0_CUSTOM_REAL
end if
rmemory_acoustic_dux_dx(i,j,ispec_PML) = coef0*rmemory_acoustic_dux_dx(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) * coef1 + PML_dux_dxl(i,j,ispec_PML) * coef2
@@ -344,12 +347,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
- coef1 = ( 1.d0 - exp(- bb * deltat / 2.d0) ) / bb
- coef2 = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * exp(- bb * deltat / 2.d0) / bb
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
+ coef1 = ( 1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) / bb
+ coef2 = ( 1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) * exp(- bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2.0_CUSTOM_REAL
+ coef2 = deltat / 2.0_CUSTOM_REAL
end if
rmemory_acoustic_dux_dz(i,j,ispec_PML) = coef0 * rmemory_acoustic_dux_dz(i,j,ispec_PML) &
@@ -379,12 +382,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
- if ( abs( bb ) > 1.d-3) then
- coef1 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) / bb
- coef2 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) / bb
+ if ( abs( bb ) > 0.001_CUSTOM_REAL) then
+ coef1 = (1.0_CUSTOM_REAL - exp(- bb * deltat / 2.0_CUSTOM_REAL) ) / bb
+ coef2 = (1.0_CUSTOM_REAL - exp(- bb * deltat / 2.0_CUSTOM_REAL) ) * exp(- bb * deltat / 2.0_CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2.0_CUSTOM_REAL
+ coef2 = deltat / 2.0_CUSTOM_REAL
end if
rmemory_acoustic_dux_dx(i,j,ispec_PML) = coef0*rmemory_acoustic_dux_dx(i,j,ispec_PML) &
@@ -407,12 +410,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(-bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
- coef1 = (1.d0 - exp(-bb * deltat / 2.d0)) / bb
- coef2 = (1.d0 - exp(-bb* deltat / 2.d0)) * exp(-bb * deltat / 2.d0) / bb
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
+ coef1 = (1._CUSTOM_REAL - exp(-bb * deltat / 2._CUSTOM_REAL)) / bb
+ coef2 = (1._CUSTOM_REAL - exp(-bb* deltat / 2._CUSTOM_REAL)) * exp(-bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2.0_CUSTOM_REAL
+ coef2 = deltat / 2.0_CUSTOM_REAL
end if
rmemory_acoustic_dux_dz(i,j,ispec_PML) = coef0 * rmemory_acoustic_dux_dz(i,j,ispec_PML) &
@@ -435,7 +438,13 @@
jacobianl = jacobian(i,j,ispec)
! if external density model
- if(assign_external_model) rhol = rhoext(i,j,ispec)
+ if(assign_external_model)then
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rhol = sngl(rhoext(i,j,ispec))
+ else
+ rhol = rhoext(i,j,ispec)
+ endif
+ endif
! for acoustic medium
! also add GLL integration weights
tempx1(i,j) = wzgll(j)*jacobianl*(xixl*dux_dxl + xizl*dux_dzl) / rhol
@@ -448,11 +457,19 @@
do j = 1,NGLLZ
do i = 1,NGLLX
if(assign_external_model) then
- rhol = rhoext(i,j,ispec)
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rhol = sngl(rhoext(i,j,ispec))
+ else
+ rhol = rhoext(i,j,ispec)
+ endif
else
-! rhol = density(1,kmato(ispec))
- lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
- mul_relaxed = poroelastcoef(2,1,kmato(ispec))
+ if(CUSTOM_REAL == SIZE_REAL) then
+ lambdal_relaxed = sngl(poroelastcoef(1,1,kmato(ispec)))
+ mul_relaxed = sngl(poroelastcoef(2,1,kmato(ispec)))
+ else
+ lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
+ mul_relaxed = poroelastcoef(2,1,kmato(ispec))
+ endif
kappal = lambdal_relaxed + TWO*mul_relaxed/3._CUSTOM_REAL
rhol = density(1,kmato(ispec))
endif
@@ -470,12 +487,12 @@
bb = alpha_x_store(i,j,ispec_PML)
coef0 = exp(- bb * deltat)
- if ( abs( bb ) > 1.d-3) then
- coef1 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) / bb
- coef2 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) / bb
+ if ( abs( bb ) > 0.001_CUSTOM_REAL) then
+ coef1 = (1._CUSTOM_REAL - exp(- bb * deltat / 2.0_CUSTOM_REAL) ) / bb
+ coef2 = (1._CUSTOM_REAL - exp(- bb * deltat / 2.0_CUSTOM_REAL) ) * exp(- bb * deltat / 2.0_CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2.0_CUSTOM_REAL
+ coef2 = deltat / 2.0_CUSTOM_REAL
end if
rmemory_potential_acoustic(1,i,j,ispec_PML)=coef0 * rmemory_potential_acoustic(1,i,j,ispec_PML) &
@@ -495,12 +512,12 @@
bb = alpha_x_store(i,j,ispec_PML)
coef0 = exp(- bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
- coef1 = ( 1 - exp(- bb * deltat / 2) ) / bb
- coef2 = ( 1 - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) / bb
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
+ coef1 = ( 1._CUSTOM_REAL - exp(- bb * deltat / 2) ) / bb
+ coef2 = ( 1._CUSTOM_REAL - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2.0_CUSTOM_REAL
+ coef2 = deltat / 2.0_CUSTOM_REAL
end if
rmemory_potential_acoustic(1,i,j,ispec_PML)=&
@@ -512,12 +529,12 @@
bb = alpha_z_store(i,j,ispec_PML)
coef0 = exp(- bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
- coef1 = ( 1 - exp(- bb * deltat / 2) ) / bb
- coef2 = ( 1 - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) / bb
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
+ coef1 = ( 1.0_CUSTOM_REAL - exp(- bb * deltat / 2) ) / bb
+ coef2 = ( 1.0_CUSTOM_REAL - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2.0_CUSTOM_REAL
+ coef2 = deltat / 2.0_CUSTOM_REAL
end if
rmemory_potential_acoustic(2,i,j,ispec_PML)=&
@@ -535,16 +552,16 @@
bb = alpha_z_store(i,j,ispec_PML)
coef0 = exp(- bb * deltat)
- if ( abs( bb ) > 1.d-3) then
- coef1 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) / bb
- coef2 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) / bb
+ if ( abs( bb ) > 0.001_CUSTOM_REAL) then
+ coef1 = (1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) / bb
+ coef2 = (1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) * exp(- bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
end if
- rmemory_potential_acoustic(1,i,j,ispec_PML)=0.d0
+ rmemory_potential_acoustic(1,i,j,ispec_PML)=0._CUSTOM_REAL
rmemory_potential_acoustic(2,i,j,ispec_PML)=coef0 * rmemory_potential_acoustic(2,i,j,ispec_PML) &
+ (potential_acoustic(iglob)+deltat*potential_dot_acoustic(iglob)) * coef1 &
+ potential_acoustic(iglob) * coef2
@@ -561,7 +578,7 @@
A1 = d_x_store(i,j,ispec_PML)
A2 = k_x_store(i,j,ispec_PML)
A3 = d_x_store(i,j,ispec_PML) * alpha_x_store(i,j,ispec_PML) ** 2
- A4 = 0.d0
+ A4 = 0._CUSTOM_REAL
if(stage_time_scheme == 6) then
bb = alpha_x_store(i,j,ispec_PML)
@@ -572,7 +589,7 @@
rmemory_potential_acoustic(1,i,j,ispec_PML) = rmemory_potential_acoustic(1,i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_potential_acoust_LDDRK(1,i,j,ispec_PML)
- rmemory_potential_acoustic(2,i,j,ispec_PML) =0.d0
+ rmemory_potential_acoustic(2,i,j,ispec_PML) =0._CUSTOM_REAL
end if
potential_dot_dot_acoustic_PML(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
@@ -596,7 +613,7 @@
A3 = alpha_x_store(i,j,ispec_PML) ** 2*(d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML)+ &
d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)) &
- -2.d0 * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)+ &
+ -2._CUSTOM_REAL * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)+ &
(it+0.5)*deltat*alpha_x_store(i,j,ispec_PML)**2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
A4 = -alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
@@ -604,7 +621,7 @@
if(stage_time_scheme == 6) then
A3 = alpha_x_store(i,j,ispec_PML) ** 2*(d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML)+ &
d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)) &
- -2.d0 * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
+ -2._CUSTOM_REAL * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
A4 = alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
end if
@@ -647,13 +664,13 @@
A0 = - alpha_z_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML)
A1 = d_z_store(i,j,ispec_PML)
A2 = k_z_store(i,j,ispec_PML)
- A3 = 0.d0
+ A3 = 0._CUSTOM_REAL
A4 = d_z_store(i,j,ispec_PML) * alpha_z_store(i,j,ispec_PML) ** 2
if(stage_time_scheme == 6) then
bb = alpha_z_store(i,j,ispec_PML)
- rmemory_potential_acoustic(1,i,j,ispec_PML) =0.d0
+ rmemory_potential_acoustic(1,i,j,ispec_PML) =0._CUSTOM_REAL
rmemory_potential_acoust_LDDRK(2,i,j,ispec_PML) = &
alpha_LDDRK(i_stage) * rmemory_potential_acoust_LDDRK(2,i,j,ispec_PML) &
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-01-14 10:04:45 UTC (rev 21226)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-01-14 12:52:31 UTC (rev 21227)
@@ -425,12 +425,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(-bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
- coef1 = (1.d0 - exp(-bb * deltat / 2.d0)) / bb
- coef2 = (1.d0 - exp(-bb* deltat / 2.d0)) * exp(-bb * deltat / 2.d0)/ bb
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
+ coef1 = (1._CUSTOM_REAL - exp(-bb * deltat / 2._CUSTOM_REAL)) / bb
+ coef2 = (1._CUSTOM_REAL - exp(-bb* deltat / 2._CUSTOM_REAL)) * exp(-bb * deltat / 2._CUSTOM_REAL)/ bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
end if
if(ROTATE_PML_ACTIVATE)then
@@ -495,12 +495,12 @@
coef0 = exp(- bb * deltat)
- if ( abs( bb ) > 1.d-3) then
- coef1 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) / bb
- coef2 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) / bb
+ if ( abs( bb ) > 0.001_CUSTOM_REAL) then
+ coef1 = (1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) / bb
+ coef2 = (1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) * exp(- bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
end if
if(ROTATE_PML_ACTIVATE)then
@@ -573,12 +573,12 @@
coef0 = exp(- bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
- coef1 = ( 1.d0 - exp(- bb * deltat / 2.d0) ) / bb
- coef2 = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * exp(- bb * deltat / 2.d0) / bb
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
+ coef1 = ( 1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) / bb
+ coef2 = ( 1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) * exp(- bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
end if
if(ROTATE_PML_ACTIVATE)then
@@ -644,12 +644,12 @@
coef0 = exp(- bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
- coef1 = ( 1.d0 - exp(- bb * deltat / 2.d0) ) / bb
- coef2 = ( 1.d0 - exp(- bb * deltat / 2.d0) ) * exp(- bb * deltat / 2.d0) / bb
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
+ coef1 = ( 1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) / bb
+ coef2 = ( 1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) * exp(- bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
end if
if(ROTATE_PML_ACTIVATE)then
@@ -718,12 +718,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
- if ( abs( bb ) > 1.d-3) then
- coef1 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) / bb
- coef2 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) / bb
+ if ( abs( bb ) > 0.001_CUSTOM_REAL) then
+ coef1 = (1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) / bb
+ coef2 = (1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) * exp(- bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
end if
if(ROTATE_PML_ACTIVATE)then
@@ -785,12 +785,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(-bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
- coef1 = (1.d0 - exp(-bb * deltat / 2.d0)) / bb
- coef2 = (1.d0 - exp(-bb* deltat / 2.d0)) * exp(-bb * deltat / 2.d0) / bb
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
+ coef1 = (1._CUSTOM_REAL - exp(-bb * deltat / 2._CUSTOM_REAL)) / bb
+ coef2 = (1._CUSTOM_REAL - exp(-bb* deltat / 2._CUSTOM_REAL)) * exp(-bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
end if
if(ROTATE_PML_ACTIVATE)then
@@ -935,7 +935,7 @@
if(PML_BOUNDARY_CONDITIONS .and. is_PML(ispec)) then
ispec_PML=spec_to_PML(ispec)
if(ROTATE_PML_ACTIVATE)then
- theta = -ROTATE_PML_ANGLE/180.d0*Pi
+ theta = -ROTATE_PML_ANGLE/180._CUSTOM_REAL*Pi
if(it==1)write(*,*)theta,ROTATE_PML_ACTIVATE,cos(theta),sin(theta)
ct=cos(theta)
st=sin(theta)
@@ -1087,12 +1087,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
- if ( abs( bb ) > 1.d-3) then
- coef1 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) / bb
- coef2 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) / bb
+ if ( abs( bb ) > 0.001_CUSTOM_REAL) then
+ coef1 = (1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) / bb
+ coef2 = (1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) * exp(- bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
end if
rmemory_displ_elastic(1,1,i,j,ispec_PML)=coef0 * rmemory_displ_elastic(1,1,i,j,ispec_PML) &
@@ -1116,12 +1116,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
coef1 = ( 1 - exp(- bb * deltat / 2) ) / bb
coef2 = ( 1 - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
end if
rmemory_displ_elastic(1,1,i,j,ispec_PML)= &
@@ -1138,12 +1138,12 @@
if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
- if ( abs(bb) > 1.d-3 ) then
+ if ( abs(bb) > 0.001_CUSTOM_REAL ) then
coef1 = ( 1 - exp(- bb * deltat / 2) ) / bb
coef2 = ( 1 - exp(- bb * deltat / 2) ) * exp(- bb * deltat / 2) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
end if
rmemory_displ_elastic(2,1,i,j,ispec_PML)=coef0 * rmemory_displ_elastic(2,1,i,j,ispec_PML) &
@@ -1165,19 +1165,19 @@
if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
- if ( abs( bb ) > 1.d-3) then
- coef1 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) / bb
- coef2 = (1.0d0 - exp(- bb * deltat / 2.0d0) ) * exp(- bb * deltat / 2.0d0) / bb
+ if ( abs( bb ) > 0.001_CUSTOM_REAL) then
+ coef1 = (1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) / bb
+ coef2 = (1._CUSTOM_REAL - exp(- bb * deltat / 2._CUSTOM_REAL) ) * exp(- bb * deltat / 2._CUSTOM_REAL) / bb
else
- coef1 = deltat / 2.0d0
- coef2 = deltat / 2.0d0
+ coef1 = deltat / 2._CUSTOM_REAL
+ coef2 = deltat / 2._CUSTOM_REAL
end if
- rmemory_displ_elastic(1,1,i,j,ispec_PML)=0.d0
+ rmemory_displ_elastic(1,1,i,j,ispec_PML)=0._CUSTOM_REAL
rmemory_displ_elastic(2,1,i,j,ispec_PML)=coef0 * rmemory_displ_elastic(2,1,i,j,ispec_PML) &
+ displ_elastic_new(1,iglob) * coef1 + displ_elastic(1,iglob) * coef2
- rmemory_displ_elastic(1,3,i,j,ispec_PML)=0.d0
+ rmemory_displ_elastic(1,3,i,j,ispec_PML)=0._CUSTOM_REAL
rmemory_displ_elastic(2,3,i,j,ispec_PML)=coef0 * rmemory_displ_elastic(2,3,i,j,ispec_PML) &
+ displ_elastic_new(3,iglob) * coef1 + displ_elastic(3,iglob) * coef2
end if
@@ -1191,7 +1191,7 @@
A1 = d_x_store(i,j,ispec_PML)
A2 = k_x_store(i,j,ispec_PML)
A3 = d_x_store(i,j,ispec_PML) * alpha_x_store(i,j,ispec_PML) ** 2
- A4 = 0.d0
+ A4 = 0._CUSTOM_REAL
if(stage_time_scheme == 6) then
@@ -1202,7 +1202,7 @@
+ deltat * (-bb * rmemory_displ_elastic(1,1,i,j,ispec_PML) + displ_elastic(1,iglob))
rmemory_displ_elastic(1,1,i,j,ispec_PML) = rmemory_displ_elastic(1,1,i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,1,i,j,ispec_PML)
- rmemory_displ_elastic(2,1,i,j,ispec_PML) =0.d0
+ rmemory_displ_elastic(2,1,i,j,ispec_PML) =0._CUSTOM_REAL
rmemory_displ_elastic_LDDRK(1,3,i,j,ispec_PML) = &
alpha_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,3,i,j,ispec_PML) &
@@ -1210,7 +1210,7 @@
rmemory_displ_elastic(1,3,i,j,ispec_PML) = rmemory_displ_elastic(1,3,i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,3,i,j,ispec_PML)
- rmemory_displ_elastic(2,3,i,j,ispec_PML) =0.d0
+ rmemory_displ_elastic(2,3,i,j,ispec_PML) =0._CUSTOM_REAL
end if
@@ -1243,14 +1243,14 @@
A3 = alpha_x_store(i,j,ispec_PML) ** 2*(d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML)+ &
d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)) &
- -2.d0 * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)+ &
+ -2._CUSTOM_REAL * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)+ &
(it+0.5)*deltat*alpha_x_store(i,j,ispec_PML)**2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
A4 = -alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
if(stage_time_scheme == 6) then
A3 = alpha_x_store(i,j,ispec_PML) ** 2*(d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML)+ &
d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)) &
- -2.d0 * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
+ -2._CUSTOM_REAL * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
A4 = alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
end if
@@ -1307,21 +1307,21 @@
A0 = - alpha_z_store(i,j,ispec_PML) * d_z_store(i,j,ispec_PML)
A1 = d_z_store(i,j,ispec_PML)
A2 = k_z_store(i,j,ispec_PML)
- A3 = 0.d0
+ A3 = 0._CUSTOM_REAL
A4 = d_z_store(i,j,ispec_PML) * alpha_z_store(i,j,ispec_PML) ** 2
if(stage_time_scheme == 6) then
bb = alpha_z_store(i,j,ispec_PML)
- rmemory_displ_elastic(1,1,i,j,ispec_PML) =0.d0
+ rmemory_displ_elastic(1,1,i,j,ispec_PML) =0._CUSTOM_REAL
rmemory_displ_elastic_LDDRK(2,1,i,j,ispec_PML) = &
alpha_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(2,1,i,j,ispec_PML) &
+ deltat * (-bb * rmemory_displ_elastic(2,1,i,j,ispec_PML) + displ_elastic(1,iglob))
rmemory_displ_elastic(2,1,i,j,ispec_PML) = rmemory_displ_elastic(2,1,i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(2,1,i,j,ispec_PML)
- rmemory_displ_elastic(1,3,i,j,ispec_PML) =0.d0
+ rmemory_displ_elastic(1,3,i,j,ispec_PML) =0._CUSTOM_REAL
rmemory_displ_elastic_LDDRK(2,3,i,j,ispec_PML) = &
alpha_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(2,3,i,j,ispec_PML) &
+ deltat * (-bb * rmemory_displ_elastic(2,3,i,j,ispec_PML) + displ_elastic(3,iglob))
@@ -1483,8 +1483,8 @@
ty = rho_vs*vy
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- displtx=0.0d0
- displtz=0.0d0
+ displtx=0._CUSTOM_REAL
+ displtz=0._CUSTOM_REAL
if(ADD_SPRING_TO_STACEY)then
@@ -1598,8 +1598,8 @@
ty = rho_vs*vy
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- displtx=0.0d0
- displtz=0.0d0
+ displtx=0._CUSTOM_REAL
+ displtz=0._CUSTOM_REAL
if(ADD_SPRING_TO_STACEY)then
@@ -1728,8 +1728,8 @@
tz = 0
endif
- displtx=0.0d0
- displtz=0.0d0
+ displtx=0._CUSTOM_REAL
+ displtz=0._CUSTOM_REAL
if(ADD_SPRING_TO_STACEY)then
@@ -1858,8 +1858,8 @@
tz = 0
endif
- displtx=0.0d0
- displtz=0.0d0
+ displtx=0._CUSTOM_REAL
+ displtz=0._CUSTOM_REAL
if(ADD_SPRING_TO_STACEY)then
@@ -2120,9 +2120,9 @@
e1_force_RK(i,j,ispec,i_sls,i_stage) = deltat * (theta_n * temp_time_scheme - Un * tauinv)
if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5d0
- if(i_stage == 2)weight_rk = 0.5d0
- if(i_stage == 3)weight_rk = 1.0d0
+ if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
if(i_stage==1)then
e1_initial_rk(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls)
@@ -2133,9 +2133,9 @@
elseif(i_stage==4)then
- e1(i,j,ispec,i_sls) = e1_initial_rk(i,j,ispec,i_sls) + 1.0d0 / 6.0d0 * &
- (e1_force_RK(i,j,ispec,i_sls,1) + 2.0d0 * e1_force_RK(i,j,ispec,i_sls,2) + &
- 2.0d0 * e1_force_RK(i,j,ispec,i_sls,3) + e1_force_RK(i,j,ispec,i_sls,4))
+ e1(i,j,ispec,i_sls) = e1_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
+ (e1_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e1_force_RK(i,j,ispec,i_sls,2) + &
+ 2._CUSTOM_REAL * e1_force_RK(i,j,ispec,i_sls,3) + e1_force_RK(i,j,ispec,i_sls,4))
endif
endif
@@ -2173,18 +2173,18 @@
e11(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5d0
- if(i_stage == 2)weight_rk = 0.5d0
- if(i_stage == 3)weight_rk = 1.0d0
+ if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
if(i_stage==1)then
e11_initial_rk(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)
endif
e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) &
+ weight_rk * e11_force_RK(i,j,ispec,i_sls,i_stage)
elseif(i_stage==4)then
- e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + 1.0d0 / 6.0d0 * &
- (e11_force_RK(i,j,ispec,i_sls,1) + 2.0d0 * e11_force_RK(i,j,ispec,i_sls,2) + &
- 2.0d0 * e11_force_RK(i,j,ispec,i_sls,3) + e11_force_RK(i,j,ispec,i_sls,4))
+ e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
+ (e11_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e11_force_RK(i,j,ispec,i_sls,2) + &
+ 2._CUSTOM_REAL * e11_force_RK(i,j,ispec,i_sls,3) + e11_force_RK(i,j,ispec,i_sls,4))
endif
endif
@@ -2219,18 +2219,18 @@
e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * (temp_time_scheme*temper_time_scheme- &
e13(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5d0
- if(i_stage == 2)weight_rk = 0.5d0
- if(i_stage == 3)weight_rk = 1.0d0
+ if(i_stage == 1)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 2)weight_rk = 0.5_CUSTOM_REAL
+ if(i_stage == 3)weight_rk = 1._CUSTOM_REAL
if(i_stage==1)then
e13_initial_rk(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)
endif
e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) &
+ weight_rk * e13_force_RK(i,j,ispec,i_sls,i_stage)
elseif(i_stage==4)then
- e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + 1.0d0 / 6.0d0 * &
- (e13_force_RK(i,j,ispec,i_sls,1) + 2.0d0 * e13_force_RK(i,j,ispec,i_sls,2) + &
- 2.0d0 * e13_force_RK(i,j,ispec,i_sls,3) + e13_force_RK(i,j,ispec,i_sls,4))
+ e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + 1._CUSTOM_REAL / 6._CUSTOM_REAL * &
+ (e13_force_RK(i,j,ispec,i_sls,1) + 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,2) + &
+ 2._CUSTOM_REAL * e13_force_RK(i,j,ispec,i_sls,3) + e13_force_RK(i,j,ispec,i_sls,4))
endif
endif
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2013-01-14 10:04:45 UTC (rev 21226)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90 2013-01-14 12:52:31 UTC (rev 21227)
@@ -90,7 +90,7 @@
! inverse mass matrices
integer :: nglob_elastic
- real(kind=CUSTOM_REAL), dimension(nglob_elastic) :: rmass_inverse_elastic_one,&
+ real(kind=CUSTOM_REAL), dimension(nglob_elastic) :: rmass_inverse_elastic_one,&
rmass_inverse_elastic_three
integer :: nglob_acoustic
@@ -192,15 +192,11 @@
if (this_element_has_PML) then
ispec_PML=spec_to_PML(ispec)
-! if ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
-! .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec))) then
if (region_CPML(ispec) == CPML_LEFT .or. region_CPML(ispec) == CPML_RIGHT) then
rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
+ wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML)&
+ d_x_store(i,j,ispec_PML) * deltat / 2.d0)
rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
-! elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
-! (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
elseif (region_CPML(ispec) == CPML_TOP_LEFT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_BOTTOM_RIGHT) then
rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob) &
@@ -245,14 +241,10 @@
if (this_element_has_PML) then
ispec_PML=spec_to_PML(ispec)
-! if ((which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
-! .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
if (region_CPML(ispec) == CPML_LEFT .or. region_CPML(ispec) == CPML_RIGHT) then
rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
+ wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(i,j,ispec_PML)&
+ d_x_store(i,j,ispec_PML) * deltat / 2.d0)
-! elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
-! (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
elseif (region_CPML(ispec) == CPML_TOP_LEFT .or. region_CPML(ispec) == CPML_TOP_RIGHT .or. &
region_CPML(ispec) == CPML_BOTTOM_LEFT .or. region_CPML(ispec) == CPML_BOTTOM_RIGHT) then
rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2013-01-14 10:04:45 UTC (rev 21226)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2013-01-14 12:52:31 UTC (rev 21227)
@@ -487,12 +487,12 @@
write(IOUT,*) ' d0_left =', d0_x_left
! endif
- d_x = ZERO
- d_z = ZERO
- K_x = ZERO
- K_z = ZERO
- alpha_x = ZERO
- alpha_z = ZERO
+ d_x = 0.d0
+ d_z = 0.d0
+ K_x = 0.d0
+ K_z = 0.d0
+ alpha_x = 0.d0
+ alpha_z = 0.d0
! origin of the PML layer (position of right edge minus thickness, in meters)
xoriginleft = thickness_PML_x_left+xmin
@@ -506,12 +506,12 @@
if (is_PML(ispec)) then
do j=1,NGLLZ
do i=1,NGLLX
- K_x_store(i,j,ispec_PML) = 0.0
- K_z_store(i,j,ispec_PML) = 0.0
- d_x_store(i,j,ispec_PML) = 0.0
- d_z_store(i,j,ispec_PML) = 0.0
- alpha_x_store(i,j,ispec_PML) = 0.0
- alpha_z_store(i,j,ispec_PML) = 0.0
+ K_x_store(i,j,ispec_PML) = 0._CUSTOM_REAL
+ K_z_store(i,j,ispec_PML) = 0._CUSTOM_REAL
+ d_x_store(i,j,ispec_PML) = 0._CUSTOM_REAL
+ d_z_store(i,j,ispec_PML) = 0._CUSTOM_REAL
+ alpha_x_store(i,j,ispec_PML) = 0._CUSTOM_REAL
+ alpha_z_store(i,j,ispec_PML) = 0._CUSTOM_REAL
iglob=ibool(i,j,ispec)
! abscissa of current grid point along the damping profile
@@ -616,12 +616,21 @@
endif
endif
+ if(CUSTOM_REAL == SIZE_REAL) then
+ K_x_store(i,j,ispec_PML) = sngl(K_x)
+ K_z_store(i,j,ispec_PML) = sngl(K_z)
+ d_x_store(i,j,ispec_PML) = sngl(d_x)
+ d_z_store(i,j,ispec_PML) = sngl(d_z)
+ alpha_x_store(i,j,ispec_PML) = sngl(alpha_x)
+ alpha_z_store(i,j,ispec_PML) = sngl(alpha_z)
+ else
K_x_store(i,j,ispec_PML) = K_x
K_z_store(i,j,ispec_PML) = K_z
d_x_store(i,j,ispec_PML) = d_x
d_z_store(i,j,ispec_PML) = d_z
alpha_x_store(i,j,ispec_PML) = alpha_x
alpha_z_store(i,j,ispec_PML) = alpha_z
+ endif
enddo
enddo
More information about the CIG-COMMITS
mailing list