[cig-commits] r21998 - in seismo/3D/SPECFEM3D/trunk/src: generate_databases specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Tue May 7 12:56:47 PDT 2013
Author: xie.zhinan
Date: 2013-05-07 12:56:46 -0700 (Tue, 07 May 2013)
New Revision: 21998
Modified:
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
Log:
commit the first stable edition of CPML for solid-fluid simulation
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -92,9 +92,36 @@
allocate(rmass_acoustic(nglob),stat=ier); if(ier /= 0) stop 'error allocating array rmass_acoustic'
rmass_acoustic(:) = 0._CUSTOM_REAL
+ allocate(rmass_acoustic_interface(nglob),stat=ier); if(ier /= 0) stop 'error allocating array rmass_acoustic'
+ rmass_acoustic_interface(:) = 0._CUSTOM_REAL
+
! returns acoustic mass matrix
if( PML_CONDITIONS ) then
call create_mass_matrices_pml_acoustic(nspec,ibool)
+
+ do ispec=1,nspec
+ if( ispec_is_acoustic(ispec) ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+ jacobianl = jacobianstore(i,j,k,ispec)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_acoustic_interface(iglob) = rmass_acoustic_interface(iglob) + &
+ sngl( dble(jacobianl) * weight / dble(kappastore(i,j,k,ispec)) )
+ else
+ rmass_acoustic_interface(iglob) = rmass_acoustic_interface(iglob) + &
+ jacobianl * weight / kappastore(i,j,k,ispec)
+ endif
+ enddo
+ enddo
+ enddo
+ endif
+ enddo
else
do ispec=1,nspec
if( ispec_is_acoustic(ispec) ) then
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -300,6 +300,7 @@
endif
if( ACOUSTIC_SIMULATION) then
deallocate(rmass_acoustic)
+ deallocate(rmass_acoustic_interface) !ZN
deallocate(rmassz_acoustic)
endif
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -203,6 +203,9 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass,rmass_acoustic,&
rmass_solid_poroelastic,rmass_fluid_poroelastic
+! mass matrix for interface !ZN
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_acoustic_interface !ZN
+
! mass matrix contributions
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx,rmassy,rmassz
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassz_acoustic
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -33,17 +33,16 @@
use generate_databases_par, only: ibool,NGLOB_AB,d_store_x,d_store_y,d_store_z, &
K_store_x,K_store_y,K_store_z,alpha_store,CPML_to_spec, &
CPML_width_x,CPML_width_y,CPML_width_z,NPOWER,K_MAX_PML, &
- CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,nspec_cpml,PML_INSTEAD_OF_FREE_SURFACE, &
+ CUSTOM_REAL,SIZE_REAL,NGLLX,NGLLY,NGLLZ,nspec_cpml,PML_INSTEAD_OF_FREE_SURFACE, &
IMAIN,FOUR_THIRDS,CPML_REGIONS,f0_FOR_PML,PI, &
CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ
- use create_regions_mesh_ext_par, only: kappastore,mustore,rhostore,rho_vp,ispec_is_acoustic,ispec_is_elastic, &
+ use create_regions_mesh_ext_par, only: rhostore,rho_vp,ispec_is_acoustic,ispec_is_elastic, &
ELASTIC_SIMULATION, ACOUSTIC_SIMULATION
implicit none
integer, intent(in) :: myrank
-
real(kind=CUSTOM_REAL), dimension(NGLOB_AB), intent(in) :: xstore,ystore,zstore
! local parameters
@@ -68,7 +67,8 @@
CPML_z_top,CPML_z_bottom,&
CPML_width_x_left_max_all, CPML_width_x_right_max_all,&
CPML_width_y_front_max_all,CPML_width_y_back_max_all,&
- CPML_width_z_top_max_all,CPML_width_z_bottom_max_all
+ CPML_width_z_top_max_all,CPML_width_z_bottom_max_all,&
+ vp_elastic,vp_acoustic,vp_max,vp_max_all
! stores damping profiles
allocate(d_store_x(NGLLX,NGLLY,NGLLZ,nspec_cpml),stat=ier)
@@ -101,6 +101,8 @@
! from Festa and Vilotte (2005)
ALPHA_MAX_PML = PI*f0_FOR_PML
+! Assuming the computational domain is convex and can be approximatly seen as a box
+! Calculation of origin of whole computational domain
x_min = minval(xstore(:))
x_max = maxval(xstore(:))
y_min = minval(ystore(:))
@@ -108,13 +110,24 @@
z_min = minval(zstore(:))
z_max = maxval(zstore(:))
- x_min_all = 10.e30
- x_max_all = -10.e30
- y_min_all = 10.e30
- y_max_all = -10.e30
- z_min_all = 10.e30
- z_max_all = -10.e30
+ if(CUSTOM_REAL == SIZE_REAL) then
+ x_min_all = 10.e30
+ y_min_all = 10.e30
+ z_min_all = 10.e30
+ x_max_all = -10.e30
+ y_max_all = -10.e30
+ z_max_all = -10.e30
+ else
+ x_min_all = 10.d30
+ y_min_all = 10.d30
+ z_min_all = 10.d30
+
+ x_max_all = -10.d30
+ y_max_all = -10.d30
+ z_max_all = -10.d30
+ endif
+
call min_all_all_cr(x_min,x_min_all)
call min_all_all_cr(y_min,y_min_all)
call min_all_all_cr(z_min,z_min_all)
@@ -123,20 +136,19 @@
call max_all_all_cr(y_max,y_max_all)
call max_all_all_cr(z_max,z_max_all)
- x_origin = (x_min_all + x_max_all)/2.0
- y_origin = (y_min_all + y_max_all)/2.0
- z_origin = (z_max_all + z_min_all)/2.0
+ x_origin = (x_min_all + x_max_all)/2._CUSTOM_REAL
+ y_origin = (y_min_all + y_max_all)/2._CUSTOM_REAL
+ z_origin = (z_max_all + z_min_all)/2._CUSTOM_REAL
-!Calculation of CPML_width_x,CPML_width_y,CPML_width_Z
-!we assume CPML_width_x,CPML_width_y,CPML_width_Z are constants inside PML layer
+! Assuming CPML_width_x,CPML_width_y,CPML_width_Z are constants inside PML layer
+! Calculation of width of PML along x, y and z direction, such as CPML_width_x,CPML_width_y,CPML_width_Z
+ CPML_width_x_left = 0._CUSTOM_REAL
+ CPML_width_x_right = 0._CUSTOM_REAL
+ CPML_width_y_front = 0._CUSTOM_REAL
+ CPML_width_y_back = 0._CUSTOM_REAL
+ CPML_width_z_top = 0._CUSTOM_REAL
+ CPML_width_z_bottom = 0._CUSTOM_REAL
- CPML_width_x_left = 0.0
- CPML_width_x_right = 0.0
- CPML_width_y_front = 0.0
- CPML_width_y_back = 0.0
- CPML_width_z_top = 0.0
- CPML_width_z_bottom = 0.0
-
CPML_x_right = x_max_all
CPML_x_left = x_min_all
CPML_y_front = y_max_all
@@ -151,7 +163,7 @@
do i=1,NGLLX
iglob = ibool(i,j,k,ispec)
if( CPML_regions(ispec_CPML) == CPML_X_ONLY ) then
- if(xstore(iglob) - x_origin > 0.d0)then
+ if(xstore(iglob) - x_origin > 0._CUSTOM_REAL)then
if(xstore(iglob) - x_origin <= CPML_x_right - x_origin )then
CPML_x_right = xstore(iglob)
endif
@@ -163,7 +175,7 @@
endif
if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
- if(ystore(iglob) - y_origin > 0.d0)then
+ if(ystore(iglob) - y_origin > 0._CUSTOM_REAL)then
if(ystore(iglob) - y_origin <= CPML_y_front - y_origin )then
CPML_y_front = ystore(iglob)
endif
@@ -175,7 +187,7 @@
endif
if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
- if(zstore(iglob) - z_origin > 0.d0)then
+ if(zstore(iglob) - z_origin > 0._CUSTOM_REAL)then
if(zstore(iglob) - z_origin <= CPML_z_top - z_origin )then
CPML_z_top = zstore(iglob)
endif
@@ -219,6 +231,31 @@
zorigintop = z_max_all - CPML_width_z_top_max_all
endif
+! Calculation of maximum p velocity inside PML
+ vp_max = 0._CUSTOM_REAL
+ do ispec_CPML=1,nspec_cpml
+ ispec = CPML_to_spec(ispec_CPML)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ vp_elastic = rho_vp(i,j,k,ispec)/rhostore(i,j,k,ispec)
+ vp_acoustic = rho_vp(i,j,k,ispec)/rhostore(i,j,k,ispec)
+
+ if(vp_acoustic .ge. vp_max)then
+ vp_max = vp_acoustic
+ endif
+ if(vp_elastic .ge. vp_max)then
+ vp_max = vp_acoustic
+ endif
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ call max_all_all_cr(vp_max,vp_max_all)
+
+
! user output
if( myrank == 0 ) then
write(IMAIN,*)
@@ -250,9 +287,15 @@
do i=1,NGLLX
! calculates P-velocity
if( ispec_is_acoustic(ispec) ) then
- vp = sqrt( kappastore(i,j,k,ispec)/rhostore(i,j,k,ispec) )
+! vp = rho_vp(i,j,k,ispec)/rhostore(i,j,k,ispec)
+! For convenience only, when computing the damping profile inside PML,
+! we set the required variable "vp" to be constant and equal to "vp_max_all"
+ vp = vp_max_all
else if( ispec_is_elastic(ispec) ) then
- vp = (FOUR_THIRDS * mustore(i,j,k,ispec) + kappastore(i,j,k,ispec)) / rho_vp(i,j,k,ispec)
+! vp = rho_vp(i,j,k,ispec)/rhostore(i,j,k,ispec)
+! For convenience only, when computing the damping profile inside PML,
+! we set the required variable "vp" to be constant and equal to "vp_max_all"
+ vp = vp_max_all
else
print*,'element index',ispec
print*,'C-PML element index ',ispec_CPML
@@ -266,7 +309,7 @@
!---------------------------- X-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
- if( xstore(iglob) - x_origin > 0.d0 ) then
+ if( xstore(iglob) - x_origin > 0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
@@ -276,20 +319,21 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ ! avoid d_x to be less than zero due to
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- else if( xstore(iglob) - x_origin < 0.d0 ) then
+ else if( xstore(iglob) - x_origin < 0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
@@ -298,17 +342,17 @@
! gets damping profile at the C-PML grid point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
else
@@ -320,11 +364,11 @@
K_store_x(i,j,k,ispec_CPML) = K_x
d_store_x(i,j,k,ispec_CPML) = d_x
- K_store_y(i,j,k,ispec_CPML) = 1.d0
- d_store_y(i,j,k,ispec_CPML) = 0.d0
+ K_store_y(i,j,k,ispec_CPML) = 1._CUSTOM_REAL
+ d_store_y(i,j,k,ispec_CPML) = 0._CUSTOM_REAL
- K_store_z(i,j,k,ispec_CPML) = 1.d0
- d_store_z(i,j,k,ispec_CPML) = 0.d0
+ K_store_z(i,j,k,ispec_CPML) = 1._CUSTOM_REAL
+ d_store_z(i,j,k,ispec_CPML) = 0._CUSTOM_REAL
alpha_store(i,j,k,ispec_CPML) = alpha_x
@@ -333,7 +377,7 @@
!---------------------------- Y-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
- if( ystore(iglob) - y_origin > 0.d0 ) then
+ if( ystore(iglob) - y_origin > 0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
@@ -342,20 +386,20 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- else if( ystore(iglob) - y_origin < 0.d0 ) then
+ else if( ystore(iglob) - y_origin < 0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
@@ -364,17 +408,17 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
else
@@ -384,14 +428,14 @@
!! DK DK define an alias for y and z variable names (which are the same)
! stores damping profiles and auxiliary coefficients at the C-PML element's GLL points
- K_store_x(i,j,k,ispec_CPML) = 1.d0
- d_store_x(i,j,k,ispec_CPML) = 0.d0
+ K_store_x(i,j,k,ispec_CPML) = 1._CUSTOM_REAL
+ d_store_x(i,j,k,ispec_CPML) = 0._CUSTOM_REAL
K_store_y(i,j,k,ispec_CPML) = K_y
d_store_y(i,j,k,ispec_CPML) = d_y
- K_store_z(i,j,k,ispec_CPML) = 1.d0
- d_store_z(i,j,k,ispec_CPML) = 0.d0
+ K_store_z(i,j,k,ispec_CPML) = 1._CUSTOM_REAL
+ d_store_z(i,j,k,ispec_CPML) = 0._CUSTOM_REAL
alpha_store(i,j,k,ispec_CPML) = alpha_y
@@ -400,7 +444,7 @@
!---------------------------- Z-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
- if( zstore(iglob) - z_origin > 0.d0 ) then
+ if( zstore(iglob) - z_origin > 0._CUSTOM_REAL ) then
if( PML_INSTEAD_OF_FREE_SURFACE ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
@@ -410,20 +454,20 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
endif
- else if( zstore(iglob) - z_origin < 0.d0 ) then
+ else if( zstore(iglob) - z_origin < 0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
@@ -432,17 +476,17 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
else
@@ -451,11 +495,11 @@
!! DK DK define an alias for y and z variable names (which are the same)
! stores damping profiles and auxiliary coefficients at the C-PML element's GLL points
- K_store_x(i,j,k,ispec_CPML) = 1.d0
- d_store_x(i,j,k,ispec_CPML) = 0.d0
+ K_store_x(i,j,k,ispec_CPML) = 1._CUSTOM_REAL
+ d_store_x(i,j,k,ispec_CPML) = 0._CUSTOM_REAL
- K_store_y(i,j,k,ispec_CPML) = 1.d0
- d_store_y(i,j,k,ispec_CPML) = 0.d0
+ K_store_y(i,j,k,ispec_CPML) = 1._CUSTOM_REAL
+ d_store_y(i,j,k,ispec_CPML) = 0._CUSTOM_REAL
K_store_z(i,j,k,ispec_CPML) = K_z
d_store_z(i,j,k,ispec_CPML) = d_z
@@ -467,7 +511,7 @@
!---------------------------- XY-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
- if( xstore(iglob) - x_origin>0.d0 .and. ystore(iglob) - y_origin>0.d0 ) then
+ if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. ystore(iglob) - y_origin>0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
@@ -476,8 +520,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
@@ -487,30 +531,30 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- else if( xstore(iglob) - x_origin>0.d0 .and. ystore(iglob) - y_origin<0.d0 ) then
+ else if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. ystore(iglob) - y_origin<0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
@@ -519,8 +563,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
@@ -530,30 +574,30 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- else if( xstore(iglob) - x_origin<0.d0 .and. ystore(iglob) - y_origin>0.d0 ) then
+ else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. ystore(iglob) - y_origin>0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
@@ -562,8 +606,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
@@ -573,30 +617,30 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- else if( xstore(iglob) - x_origin <0.d0 .and. ystore(iglob) - y_origin <0.d0 ) then
+ else if( xstore(iglob) - x_origin <0._CUSTOM_REAL .and. ystore(iglob) - y_origin <0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
@@ -605,8 +649,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
@@ -616,27 +660,27 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
else
@@ -652,17 +696,17 @@
d_store_y(i,j,k,ispec_CPML) = d_y
- K_store_z(i,j,k,ispec_CPML) = 1.d0
- d_store_z(i,j,k,ispec_CPML) = 0.d0
+ K_store_z(i,j,k,ispec_CPML) = 1._CUSTOM_REAL
+ d_store_z(i,j,k,ispec_CPML) = 0._CUSTOM_REAL
- alpha_store(i,j,k,ispec_CPML) = ALPHA_MAX_PML / 2.d0
+ alpha_store(i,j,k,ispec_CPML) = ALPHA_MAX_PML / 2._CUSTOM_REAL
else if( CPML_regions(ispec_CPML) == CPML_XZ_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- XZ-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
- if( xstore(iglob) - x_origin>0.d0 .and. zstore(iglob) - z_origin>0.d0 ) then
+ if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
if( PML_INSTEAD_OF_FREE_SURFACE ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
@@ -672,8 +716,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
@@ -683,32 +727,32 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
endif
- else if( xstore(iglob) - x_origin>0.d0 .and. zstore(iglob) - z_origin<0.d0 ) then
+ else if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. zstore(iglob) - z_origin<0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
@@ -717,8 +761,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
@@ -728,30 +772,30 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- else if( xstore(iglob) - x_origin<0.d0 .and. zstore(iglob) - z_origin>0.d0 ) then
+ else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
if( PML_INSTEAD_OF_FREE_SURFACE ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
@@ -761,8 +805,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
@@ -772,32 +816,32 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
endif
- else if( xstore(iglob) - x_origin<0.d0 .and. zstore(iglob) - z_origin<0.d0 ) then
+ else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. zstore(iglob) - z_origin<0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
@@ -806,8 +850,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
@@ -817,27 +861,27 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
else
@@ -849,20 +893,20 @@
K_store_x(i,j,k,ispec_CPML) = K_x
d_store_x(i,j,k,ispec_CPML) = d_x
- K_store_y(i,j,k,ispec_CPML) = 1.d0
- d_store_y(i,j,k,ispec_CPML) = 0.d0
+ K_store_y(i,j,k,ispec_CPML) = 1._CUSTOM_REAL
+ d_store_y(i,j,k,ispec_CPML) = 0._CUSTOM_REAL
K_store_z(i,j,k,ispec_CPML) = K_z
d_store_z(i,j,k,ispec_CPML) = d_z
- alpha_store(i,j,k,ispec_CPML) = ALPHA_MAX_PML / 2.d0
+ alpha_store(i,j,k,ispec_CPML) = ALPHA_MAX_PML / 2._CUSTOM_REAL
else if( CPML_regions(ispec_CPML) == CPML_YZ_ONLY ) then
!------------------------------------------------------------------------------
!---------------------------- YZ-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
- if( ystore(iglob) - y_origin>0.d0 .and. zstore(iglob) - z_origin>0.d0 ) then
+ if( ystore(iglob) - y_origin>0._CUSTOM_REAL .and. zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
if( PML_INSTEAD_OF_FREE_SURFACE ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
@@ -872,8 +916,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
@@ -883,31 +927,31 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
endif
- else if( ystore(iglob) - y_origin>0.d0 .and. zstore(iglob) - z_origin<0.d0 ) then
+ else if( ystore(iglob) - y_origin>0._CUSTOM_REAL .and. zstore(iglob) - z_origin<0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
@@ -916,8 +960,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
@@ -927,30 +971,30 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- else if( ystore(iglob) - y_origin<0.d0 .and. zstore(iglob) - z_origin>0.d0 ) then
+ else if( ystore(iglob) - y_origin<0._CUSTOM_REAL .and. zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
if( PML_INSTEAD_OF_FREE_SURFACE ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
@@ -960,8 +1004,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
@@ -971,32 +1015,32 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
endif
- else if( ystore(iglob) - y_origin<0.d0 .and. zstore(iglob) - z_origin<0.d0 ) then
+ else if( ystore(iglob) - y_origin<0._CUSTOM_REAL .and. zstore(iglob) - z_origin<0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
@@ -1005,8 +1049,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
@@ -1016,27 +1060,27 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
else
@@ -1044,8 +1088,8 @@
endif
!! DK DK define an alias for y and z variable names (which are the same)
- K_store_x(i,j,k,ispec_CPML) = 1.d0
- d_store_x(i,j,k,ispec_CPML) = 0.d0
+ K_store_x(i,j,k,ispec_CPML) = 1._CUSTOM_REAL
+ d_store_x(i,j,k,ispec_CPML) = 0._CUSTOM_REAL
K_store_y(i,j,k,ispec_CPML) = K_y
d_store_y(i,j,k,ispec_CPML) = d_y
@@ -1053,16 +1097,16 @@
K_store_z(i,j,k,ispec_CPML) = K_z
d_store_z(i,j,k,ispec_CPML) = d_z
- alpha_store(i,j,k,ispec_CPML) = ALPHA_MAX_PML / 2.d0
+ alpha_store(i,j,k,ispec_CPML) = ALPHA_MAX_PML / 2._CUSTOM_REAL
else if( CPML_regions(ispec_CPML) == CPML_XYZ ) then
!------------------------------------------------------------------------------
!---------------------------- XYZ-corner C-PML --------------------------------
!------------------------------------------------------------------------------
- if( xstore(iglob) - x_origin>0.d0 .and. &
- ystore(iglob) - y_origin>0.d0 .and. &
- zstore(iglob) - z_origin>0.d0 ) then
+ if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. &
+ ystore(iglob) - y_origin>0._CUSTOM_REAL .and. &
+ zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
if( PML_INSTEAD_OF_FREE_SURFACE ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
@@ -1072,8 +1116,8 @@
! gets damping profile at the C-PML grid point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
@@ -1083,8 +1127,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
@@ -1094,43 +1138,43 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
endif
- else if( xstore(iglob) - x_origin>0.d0 .and. &
- ystore(iglob) - y_origin>0.d0 .and. &
- zstore(iglob) - z_origin<0.d0 ) then
+ else if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. &
+ ystore(iglob) - y_origin>0._CUSTOM_REAL .and. &
+ zstore(iglob) - z_origin<0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
@@ -1139,8 +1183,8 @@
! gets damping profile at the C-PML grid point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
@@ -1150,8 +1194,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
@@ -1161,43 +1205,43 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- else if( xstore(iglob) - x_origin>0.d0 .and. &
- ystore(iglob) - y_origin<0.d0 .and. &
- zstore(iglob) - z_origin>0.d0 ) then
+ else if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. &
+ ystore(iglob) - y_origin<0._CUSTOM_REAL .and. &
+ zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
if( PML_INSTEAD_OF_FREE_SURFACE ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
@@ -1207,8 +1251,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
@@ -1218,8 +1262,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
@@ -1230,44 +1274,44 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
endif
- else if( xstore(iglob) - x_origin>0.d0 .and. &
- ystore(iglob) - y_origin<0.d0 .and. &
- zstore(iglob) - z_origin < 0.d0 ) then
+ else if( xstore(iglob) - x_origin>0._CUSTOM_REAL .and. &
+ ystore(iglob) - y_origin<0._CUSTOM_REAL .and. &
+ zstore(iglob) - z_origin < 0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xstore(iglob) - xoriginright
@@ -1276,8 +1320,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
@@ -1288,8 +1332,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
@@ -1300,42 +1344,42 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- else if( xstore(iglob) - x_origin<0.d0 .and. &
- ystore(iglob) - y_origin>0.d0 .and. &
- zstore(iglob) - z_origin>0.d0 ) then
+ else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. &
+ ystore(iglob) - y_origin>0._CUSTOM_REAL .and. &
+ zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
if( PML_INSTEAD_OF_FREE_SURFACE ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
@@ -1345,8 +1389,8 @@
! gets damping profile at the C-PML grid point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
@@ -1356,8 +1400,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
@@ -1367,44 +1411,44 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
endif
- else if( xstore(iglob) - x_origin<0.d0 .and. &
- ystore(iglob) - y_origin>0.d0 .and. &
- zstore(iglob) - z_origin<0.d0 ) then
+ else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. &
+ ystore(iglob) - y_origin>0._CUSTOM_REAL .and. &
+ zstore(iglob) - z_origin<0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
@@ -1413,8 +1457,8 @@
! gets damping profile at the C-PML grid point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = ystore(iglob) - yoriginfront
@@ -1424,8 +1468,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
@@ -1435,42 +1479,42 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- else if( xstore(iglob) - x_origin<0.d0 .and. &
- ystore(iglob) - y_origin<0.d0 .and. &
- zstore(iglob) - z_origin>0.d0 ) then
+ else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. &
+ ystore(iglob) - y_origin<0._CUSTOM_REAL .and. &
+ zstore(iglob) - z_origin>0._CUSTOM_REAL ) then
if( PML_INSTEAD_OF_FREE_SURFACE ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
@@ -1480,8 +1524,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
@@ -1491,8 +1535,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zstore(iglob) - zorigintop
@@ -1502,44 +1546,44 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
endif
- else if( xstore(iglob) - x_origin<0.d0 .and. &
- ystore(iglob) - y_origin<0.d0 .and. &
- zstore(iglob) - z_origin < 0.d0 ) then
+ else if( xstore(iglob) - x_origin<0._CUSTOM_REAL .and. &
+ ystore(iglob) - y_origin<0._CUSTOM_REAL .and. &
+ zstore(iglob) - z_origin < 0._CUSTOM_REAL ) then
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_x = xoriginleft - xstore(iglob)
@@ -1548,8 +1592,8 @@
! gets damping profile at the C-PML element's GLL point
d_x = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_x)
- alpha_x = ALPHA_MAX_PML / 2.d0
- K_x = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_x = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_y = yoriginback - ystore(iglob)
@@ -1559,8 +1603,8 @@
! gets damping profile at the C-PML element's GLL point
d_y = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_y)
- alpha_y = ALPHA_MAX_PML / 2.d0
- K_y = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_y = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
! gets abscissa of current grid point along the damping profile
abscissa_in_PML_z = zoriginbottom - zstore(iglob)
@@ -1570,37 +1614,37 @@
! gets damping profile at the C-PML element's GLL point
d_z = pml_damping_profile_l(myrank,iglob,dist,vp,CPML_width_z)
- alpha_z = ALPHA_MAX_PML / 2.d0
- K_z = 1.d0 + (K_MAX_PML - 1.d0) * dist**NPOWER
+ alpha_z = ALPHA_MAX_PML / 2._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL + (K_MAX_PML - 1._CUSTOM_REAL) * dist**NPOWER
- if( d_x < 0.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( d_x < 0._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( K_x < 1.d0 ) then
- d_x = 0.d0
- K_x = 1.d0
+ if( K_x < 1._CUSTOM_REAL ) then
+ d_x = 0._CUSTOM_REAL
+ K_x = 1._CUSTOM_REAL
endif
- if( d_y < 0.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( d_y < 0._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( K_y < 1.d0 ) then
- d_y = 0.d0
- K_y = 1.d0
+ if( K_y < 1._CUSTOM_REAL ) then
+ d_y = 0._CUSTOM_REAL
+ K_y = 1._CUSTOM_REAL
endif
- if( d_z < 0.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( d_z < 0._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
- if( K_z < 1.d0 ) then
- d_z = 0.d0
- K_z = 1.d0
+ if( K_z < 1._CUSTOM_REAL ) then
+ d_z = 0._CUSTOM_REAL
+ K_z = 1._CUSTOM_REAL
endif
else
@@ -1617,7 +1661,7 @@
K_store_z(i,j,k,ispec_CPML) = K_z
d_store_z(i,j,k,ispec_CPML) = d_z
- alpha_store(i,j,k,ispec_CPML) = ALPHA_MAX_PML / 2.d0
+ alpha_store(i,j,k,ispec_CPML) = ALPHA_MAX_PML / 2._CUSTOM_REAL
endif
enddo
@@ -1651,9 +1695,10 @@
! gets damping profile
if( NPOWER >= 1 ) then
! In INRIA research report section 6.1: http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
- ! pml_damping_profile_l = - ((NPOWER + 1) * vp * log(CPML_Rcoef) / (2.d0 * delta)) * dist**(NPOWER)
+ ! pml_damping_profile_l = - ((NPOWER + 1) * vp * log(CPML_Rcoef) / (2._CUSTOM_REAL * delta)) * dist**(NPOWER)
! due to tests it is more accurate to use following definition in case NPOWER = 1 defined in constants.h.in
- pml_damping_profile_l = - ((NPOWER + 1) * vp * log(CPML_Rcoef) / (2.d0 * delta)) * dist**(1.2*NPOWER)
+ pml_damping_profile_l = - ((NPOWER + 1) * vp * log(CPML_Rcoef) / (2._CUSTOM_REAL * delta)) &
+ * dist**(1.2_CUSTOM_REAL*NPOWER)
else
call exit_mpi(myrank,'C-PML error: NPOWER must be greater than or equal to 1')
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -98,6 +98,7 @@
! acoustic
if( ACOUSTIC_SIMULATION ) then
write(IOUT) rmass_acoustic
+ write(IOUT) rmass_acoustic_interface !ZN
endif
! this array is needed for acoustic simulations but also for elastic simulations with CPML,
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -33,7 +33,7 @@
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
- PML_CONDITIONS,spec_to_CPML,is_CPML)
+ PML_CONDITIONS,spec_to_CPML,is_CPML,potential_dot_dot_acoustic_interface)
! returns the updated pressure array: potential_dot_dot_acoustic
@@ -44,7 +44,8 @@
! displacement and pressure
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic,&
+ potential_dot_dot_acoustic_interface
! global indexing
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
@@ -130,6 +131,8 @@
! it also means you have to calculate and update this here first before
! calculating the coupling on the elastic side for the acceleration...
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + jacobianw*displ_n
+ potential_dot_dot_acoustic_interface(iglob) = potential_dot_dot_acoustic_interface(iglob) &
+ + jacobianw*displ_n
enddo ! igll
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -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,is_CPML,potential_dot_dot_acoustic_interface)
! returns the updated acceleration array: accel
@@ -43,7 +44,8 @@
! displacement and pressure
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: accel
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic,&
+ potential_dot_dot_acoustic_interface
! global indexing
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
@@ -66,6 +68,10 @@
integer :: iface,igll,ispec,iglob
integer :: i,j,k
+! CPML
+ logical :: PML_CONDITIONS
+ logical :: is_CPML(NSPEC_AB)
+
! loops on all coupling faces
do iface = 1,num_coupling_ac_el_faces
@@ -87,7 +93,15 @@
iglob = ibool(i,j,k,ispec)
! acoustic pressure on global point
- pressure = - potential_dot_dot_acoustic(iglob)
+ if(PML_CONDITIONS)then
+ if(is_CPML(ispec))then
+ pressure = - potential_dot_dot_acoustic_interface(iglob)
+ else
+ pressure = - potential_dot_dot_acoustic(iglob)
+ endif
+ else
+ pressure = - potential_dot_dot_acoustic(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-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -62,7 +62,7 @@
implicit none
! local parameters
- integer:: iphase
+ integer:: iphase,iface,ispec,iglob,igll,i,j,k
logical:: phase_is_inner
! enforces free surface (zeroes potentials at free surface)
@@ -116,7 +116,7 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
rhostore,jacobian,ibool,deltat, &
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
- phase_ispec_inner_acoustic)
+ phase_ispec_inner_acoustic,ELASTIC_SIMULATION,potential_dot_dot_acoustic_interface)
endif
! adjoint simulations
@@ -140,7 +140,7 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
rhostore,jacobian,ibool,deltat, &
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
- phase_ispec_inner_acoustic)
+ phase_ispec_inner_acoustic,ELASTIC_SIMULATION,potential_dot_dot_acoustic_interface)
endif
endif
@@ -177,7 +177,7 @@
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
- PML_CONDITIONS,spec_to_CPML,is_CPML)
+ PML_CONDITIONS,spec_to_CPML,is_CPML,potential_dot_dot_acoustic_interface)
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
@@ -188,7 +188,7 @@
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
- PML_CONDITIONS,spec_to_CPML,is_CPML)
+ PML_CONDITIONS,spec_to_CPML,is_CPML,potential_dot_dot_acoustic_interface)
endif
! adjoint/kernel simulations
if( SIMULATION_TYPE == 3 ) &
@@ -199,7 +199,7 @@
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner,&
- PML_CONDITIONS,spec_to_CPML,is_CPML)
+ PML_CONDITIONS,spec_to_CPML,is_CPML,potential_dot_dot_acoustic_interface)
else
! on GPU
@@ -346,6 +346,13 @@
if(.NOT. GPU_MODE) then
! divides pressure with mass matrix
potential_dot_dot_acoustic(:) = potential_dot_dot_acoustic(:) * rmass_acoustic(:)
+
+ if(PML_CONDITIONS)then
+ if(ELASTIC_SIMULATION ) then
+ potential_dot_dot_acoustic_interface(:) = potential_dot_dot_acoustic_interface(:) * &
+ rmass_acoustic_interface(:)
+ endif
+ endif
! adjoint simulations
if (SIMULATION_TYPE == 3) &
@@ -355,6 +362,44 @@
call kernel_3_a_acoustic_cuda(Mesh_pointer,NGLOB_AB)
endif
+! The outer boundary condition to use for PML elements in fluid layers is Neumann for the potential
+! because we need Dirichlet conditions for the displacement vector, which means Neumann for the potential.
+! Thus, there is nothing to enforce explicitly here.
+! There is something to enforce explicitly only in the case of elastic elements, for which a Dirichlet
+! condition is needed for the displacement vector, which is the vectorial unknown for these elements.
+
+! However, enforcing explicitly potential_dot_dot_acoustic, potential_dot_acoustic, potential_acoustic
+! to be zero on outer boundary of PML help to improve the accuracy of absorbing low-frequency wave components
+! in case of long-time simulation
+
+ ! C-PML boundary
+ if(PML_CONDITIONS)then
+ do iface=1,num_abs_boundary_faces
+ ispec = abs_boundary_ispec(iface)
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+ if( ispec_is_acoustic(ispec) .and. is_CPML(ispec) ) then
+ ! reference gll points on boundary face
+ do igll = 1,NGLLSQUARE
+
+ ! gets local indices for GLL point
+ i = abs_boundary_ijk(1,igll,iface)
+ j = abs_boundary_ijk(2,igll,iface)
+ k = abs_boundary_ijk(3,igll,iface)
+
+ iglob=ibool(i,j,k,ispec)
+
+ potential_dot_dot_acoustic(iglob) = 0.0
+ potential_dot_acoustic(iglob) = 0.0
+ potential_acoustic(iglob) = 0.0
+ if(ELASTIC_SIMULATION ) then
+ potential_dot_dot_acoustic_interface(iglob) = 0.0
+ endif
+ enddo
+ endif ! ispec_is_acoustic
+ endif
+ enddo
+ endif
+
! update velocity
! note: Newmark finite-difference time scheme with acoustic domains:
! (see e.g. Hughes, 1987; Chaljub et al., 2003)
@@ -365,7 +410,6 @@
!
! where
! chi, chi_dot, chi_dot_dot are acoustic (fluid) potentials ( dotted with respect to time)
-! u, v, a are displacement,velocity & acceleration
! M is mass matrix, K stiffness matrix and B boundary term
! f denotes a source term
!
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -34,7 +34,7 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
rhostore,jacobian,ibool,deltat, &
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
- phase_ispec_inner_acoustic)
+ phase_ispec_inner_acoustic,ELASTIC_SIMULATION,potential_dot_dot_acoustic_interface)
! computes forces for acoustic elements
!
@@ -50,7 +50,7 @@
! acoustic potentials
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: &
- potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic
+ potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic
! time step
real(kind=CUSTOM_REAL) :: deltat
@@ -75,6 +75,10 @@
integer :: num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic
integer, dimension(num_phase_ispec_acoustic,2) :: phase_ispec_inner_acoustic
+!CPML fluid-solid interface
+ logical :: ELASTIC_SIMULATION
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic_interface
+
! local variables
real(kind=CUSTOM_REAL) :: hp1,hp2,hp3
@@ -93,7 +97,6 @@
! local C-PML absorbing boundary conditions parameters
integer :: ispec_CPML
-
if( iphase == 1 ) then
num_elements = nspec_outer_acoustic
else
@@ -252,6 +255,9 @@
! do not merge this second line with the first using an ".and." statement
! because array is_CPML() is unallocated when PML_CONDITIONS is false
if(is_CPML(ispec)) then
+ if(ELASTIC_SIMULATION)then
+ potential_dot_dot_acoustic_interface(iglob) = potential_dot_dot_acoustic(iglob)
+ endif
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - &
potential_dot_dot_acoustic_CPML(i,j,k,ispec_CPML)
endif
@@ -263,11 +269,6 @@
enddo ! end of loop over all spectral elements
-! The outer boundary condition to use for PML elements in fluid layers is Neumann for the potential
-! because we need Dirichlet conditions for the displacement vector, which means Neumann for the potential.
-! Thus, there is nothing to enforce explicitly here.
-! There is something to enforce explicitly only in the case of elastic elements, for which a Dirichlet
-! condition is needed for the displacement vector, which is the vectorial unknown for these elements.
end subroutine compute_forces_acoustic_noDev
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -32,6 +32,7 @@
use specfem_par_acoustic
use specfem_par_elastic
use specfem_par_poroelastic
+ use pml_par
use fault_solver_dynamic, only : bc_dynflt_set3d_all,SIMULATION_TYPE_DYN
use fault_solver_kinematic, only : bc_kinflt_set_all,SIMULATION_TYPE_KIN
@@ -39,6 +40,7 @@
integer:: iphase
logical:: phase_is_inner
+ integer:: iface,ispec,igll,i,j,k,iglob
! distinguishes two runs: for points on MPI interfaces, and points within the partitions
do iphase=1,2
@@ -88,8 +90,7 @@
dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic,ispec_is_elastic,&
- abs_boundary_ijk,num_abs_boundary_faces,phase_is_inner,ispec_is_inner,abs_boundary_ispec)
+ phase_ispec_inner_elastic)
! adjoint simulations: backward/reconstructed wavefield
if( SIMULATION_TYPE == 3 ) &
@@ -119,8 +120,7 @@
b_dsdx_top,b_dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic,ispec_is_elastic,&
- abs_boundary_ijk,num_abs_boundary_faces,phase_is_inner,ispec_is_inner,abs_boundary_ispec)
+ phase_ispec_inner_elastic)
endif
@@ -184,7 +184,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,is_CPML,potential_dot_dot_acoustic_interface)
else
! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
! adoint definition: pressure^\dagger=potential^\dagger
@@ -194,7 +195,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,is_CPML,potential_dot_dot_acoustic_interface)
endif
! adjoint simulations
@@ -205,7 +207,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,is_CPML,potential_dot_dot_acoustic_interface)
else
! on GPU
@@ -367,6 +370,32 @@
endif
endif
+ ! C-PML boundary
+ if(PML_CONDITIONS)then
+ do iface=1,num_abs_boundary_faces
+ ispec = abs_boundary_ispec(iface)
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+ if( ispec_is_elastic(ispec) .and. is_CPML(ispec)) then
+ ! reference gll points on boundary face
+ do igll = 1,NGLLSQUARE
+
+ ! gets local indices for GLL point
+ i = abs_boundary_ijk(1,igll,iface)
+ j = abs_boundary_ijk(2,igll,iface)
+ k = abs_boundary_ijk(3,igll,iface)
+
+ iglob=ibool(i,j,k,ispec)
+
+ accel(:,iglob) = 0.0
+ veloc(:,iglob) = 0.0
+ displ(:,iglob) = 0.0
+
+ enddo
+ endif ! ispec_is_elastic
+ endif
+ enddo
+ endif
+
! updates velocities
! Newmark finite-difference time scheme with elastic domains:
! (see e.g. Hughes, 1987; Chaljub et al., 2003)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -50,10 +50,9 @@
dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic,ispec_is_elastic,&
- abs_boundary_ijk,num_abs_boundary_faces,phase_is_inner,ispec_is_inner,abs_boundary_ispec)
+ phase_ispec_inner_elastic)
- use constants, only: NGLLX,NGLLY,NGLLZ,NDIM,N_SLS,SAVE_MOHO_MESH,ONE_THIRD,FOUR_THIRDS,IOUT,NGLLSQUARE
+ use constants, only: NGLLX,NGLLY,NGLLZ,NDIM,N_SLS,SAVE_MOHO_MESH,ONE_THIRD,FOUR_THIRDS
use pml_par
use fault_solver_dynamic, only : Kelvin_Voigt_eta
use specfem_par, only : FULL_ATTENUATION_SOLID
@@ -128,8 +127,6 @@
! C-PML absorbing boundary conditions
logical :: PML_CONDITIONS
- logical, dimension(NSPEC_AB) :: ispec_is_elastic
-
! local parameters
integer :: i_SLS,imodulo_N_SLS
integer :: ispec,iglob,ispec_p,num_elements
@@ -181,16 +178,6 @@
! local C-PML absorbing boundary conditions parameters
integer :: ispec_CPML
-! communication overlap
- logical, dimension(NSPEC_AB) :: ispec_is_inner
- logical :: phase_is_inner
-
-! outer boundary of CPML
- integer :: num_abs_boundary_faces
- integer :: abs_boundary_ijk(3,NGLLSQUARE,num_abs_boundary_faces)
- integer :: abs_boundary_ispec(num_abs_boundary_faces)
- integer :: igll,iface
-
imodulo_N_SLS = mod(N_SLS,3)
! choses inner/outer elements
@@ -869,40 +856,6 @@
enddo ! spectral element loop
-!! DK DK to Jo, to debug CPML, 22 March 2013:
-!! DK DK I think that there is an error in the loops below, you should also check if ispec is a CPML element,
-!! DK DK and also if ispec is an elastic or viscoelastic element (and NOT for instance an acoustic element)
-!! DK DK
-!! DK DK thus test something like: if (is_CPML(ispec) .and. elastic(ispec)) then
-!! DK DK or something like that
-
- ! C-PML boundary
- if(PML_CONDITIONS)then
- do iface=1,num_abs_boundary_faces
- ispec = abs_boundary_ispec(iface)
- if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
- if( ispec_is_elastic(ispec) .and. is_CPML(ispec)) then
- ! reference gll points on boundary face
- do igll = 1,NGLLSQUARE
-
- ! gets local indices for GLL point
- i = abs_boundary_ijk(1,igll,iface)
- j = abs_boundary_ijk(2,igll,iface)
- k = abs_boundary_ijk(3,igll,iface)
-
- iglob=ibool(i,j,k,ispec)
-
- accel(:,iglob) = 0.0
- veloc(:,iglob) = 0.0
- displ(:,iglob) = 0.0
-
- enddo
- endif ! ispec_is_elastic
- endif
- enddo
- endif
-
-
end subroutine compute_forces_viscoelastic_noDev
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -277,9 +277,6 @@
if( .not. GPU_MODE .and. GRAVITY ) &
stop 'GRAVITY only supported in GPU mode'
- if(STACEY_ABSORBING_CONDITIONS .and. PML_CONDITIONS) &
- stop 'STACEY_ABSORBING_CONDITIONS and PML_CONDITIONS are mutually exclusive but are both set to .true.'
-
! absorbing surfaces
if( STACEY_ABSORBING_CONDITIONS ) then
! for arbitrary orientation of elements, which face belongs to xmin,xmax,etc... -
@@ -288,23 +285,26 @@
! just to be sure for now..
if( NGLLX /= NGLLY .and. NGLLY /= NGLLZ ) &
stop 'STACEY_ABSORBING_CONDITIONS must have NGLLX = NGLLY = NGLLZ'
- if( PML_INSTEAD_OF_FREE_SURFACE ) then
- print*, 'please modify Par_file'
+ if( PML_CONDITIONS ) then
+ print*, 'please modify Par_file and recompile solver'
+ stop 'STACEY_ABSORBING_CONDITIONS and PML_CONDITIONS are both set to .true.'
+ else if( PML_INSTEAD_OF_FREE_SURFACE ) then
+ print*, 'please modify Par_file and recompile solver'
stop 'PML_INSTEAD_OF_FREE_SURFACE = .true. is incompatible with STACEY_ABSORBING_CONDITIONS = .true.'
endif
else
if( STACEY_INSTEAD_OF_FREE_SURFACE ) then
- print*, 'please modify Par_file'
+ print*, 'please modify Par_file and recompile solver'
stop 'STACEY_ABSORBING_CONDITIONS must be activated when STACEY_INSTEAD_OF_FREE_SURFACE is set to .true.'
endif
endif
if( PML_CONDITIONS ) then
if( STACEY_INSTEAD_OF_FREE_SURFACE ) then
- print*, 'please modify Par_file'
+ print*, 'please modify Par_file and recompile solver'
stop 'STACEY_INSTEAD_OF_FREE_SURFACE = .true. is incompatible with PML_CONDITIONS = .true.'
else if( .not. SUPPRESS_UTM_PROJECTION ) then
- print*, 'please modify Par_file'
+ print*, 'please modify Par_file and recompile solver'
stop 'SUPPRESS_UTM_PROJECTION must be activated when PML_CONDITIONS is set to .true.'
endif
else
@@ -313,13 +313,13 @@
endif
if( STACEY_INSTEAD_OF_FREE_SURFACE .and. PML_INSTEAD_OF_FREE_SURFACE ) then
- print*, 'please modify Par_file'
+ print*, 'please modify Par_file and recompile solver'
stop 'error: STACEY_INSTEAD_OF_FREE_SURFACE and PML_INSTEAD_OF_FREE_SURFACE are both set to .true.'
endif
! checks the MOVIE_TYPE parameter
if( MOVIE_TYPE /= 1 .and. MOVIE_TYPE /= 2 ) then
- stop 'error: MOVIE_TYPE must be either 1 or 2! Please modify Par_file'
+ stop 'error: MOVIE_TYPE must be either 1 or 2! Please modify Par_file and recompile solver'
endif
! check that the code has been compiled with the right values
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -243,14 +243,22 @@
deallocate(rmassz_acoustic)
endif
- call assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,rmass_acoustic, &
+ call assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,rmass_acoustic,&
num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
my_neighbours_ext_mesh)
+ call assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,rmass_acoustic_interface, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
+ my_neighbours_ext_mesh)
+
! fill mass matrix with fictitious non-zero values to make sure it can be inverted globally
where(rmass_acoustic <= 0._CUSTOM_REAL) rmass_acoustic = 1._CUSTOM_REAL
rmass_acoustic(:) = 1._CUSTOM_REAL / rmass_acoustic(:)
+
+ where(rmass_acoustic_interface <= 0._CUSTOM_REAL) rmass_acoustic_interface = 1._CUSTOM_REAL
+ rmass_acoustic_interface(:) = 1._CUSTOM_REAL / rmass_acoustic_interface(:)
endif
if(ELASTIC_SIMULATION) then
@@ -462,7 +470,7 @@
double precision :: MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD
real(kind=CUSTOM_REAL):: scale_factorl
integer :: i,j,k,ispec,ier
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: scale_factor,scale_factor_kappa !ZN
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: scale_factor,scale_factor_kappa
! if attenuation is on, shift shear moduli to center frequency of absorption period band, i.e.
! rescale mu to average (central) frequency for attenuation
@@ -476,13 +484,12 @@
if( ier /= 0 ) call exit_mpi(myrank,'error allocation scale_factor')
scale_factor(:,:,:,:) = 1._CUSTOM_REAL
- one_minus_sum_beta_kappa(:,:,:,:) = 1._CUSTOM_REAL !ZN
- factor_common_kappa(:,:,:,:,:) = 1._CUSTOM_REAL !ZN
- allocate( scale_factor_kappa(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa),stat=ier) !ZN
- if( ier /= 0 ) call exit_mpi(myrank,'error allocation scale_factor_kappa') !ZN
- scale_factor_kappa(:,:,:,:) = 1._CUSTOM_REAL !ZN
+ one_minus_sum_beta_kappa(:,:,:,:) = 1._CUSTOM_REAL
+ factor_common_kappa(:,:,:,:,:) = 1._CUSTOM_REAL
+ allocate( scale_factor_kappa(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_kappa),stat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error allocation scale_factor_kappa')
+ scale_factor_kappa(:,:,:,:) = 1._CUSTOM_REAL
-
! reads in attenuation arrays
open(unit=27, file=prname(1:len_trim(prname))//'attenuation.bin', &
status='old',action='read',form='unformatted',iostat=ier)
@@ -500,11 +507,11 @@
read(27) factor_common
read(27) scale_factor
- if(FULL_ATTENUATION_SOLID)then !ZN
- read(27) one_minus_sum_beta_kappa !ZN
- read(27) factor_common_kappa !ZN
- read(27) scale_factor_kappa !ZN
- endif !ZN
+ if(FULL_ATTENUATION_SOLID)then
+ read(27) one_minus_sum_beta_kappa
+ read(27) factor_common_kappa
+ read(27) scale_factor_kappa
+ endif
close(27)
@@ -538,11 +545,11 @@
scale_factorl = scale_factor(i,j,k,ispec)
mustore(i,j,k,ispec) = mustore(i,j,k,ispec) * scale_factorl
- if(FULL_ATTENUATION_SOLID)then !ZN
+ if(FULL_ATTENUATION_SOLID)then
! scales kappa moduli
scale_factorl = scale_factor_kappa(i,j,k,ispec)
kappastore(i,j,k,ispec) = kappastore(i,j,k,ispec) * scale_factorl
- endif !ZN
+ endif
enddo
enddo
@@ -550,7 +557,7 @@
enddo
deallocate(scale_factor)
- deallocate(scale_factor_kappa) !ZN
+ deallocate(scale_factor_kappa)
! statistics
! user output
@@ -567,14 +574,14 @@
! clear memory variables if attenuation
! initialize memory variables for attenuation
- epsilondev_trace(:,:,:,:) = 0._CUSTOM_REAL !ZN
+ epsilondev_trace(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_xx(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_yy(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_xy(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_xz(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_yz(:,:,:,:) = 0._CUSTOM_REAL
- R_trace(:,:,:,:,:) = 0._CUSTOM_REAL !ZN
+ R_trace(:,:,:,:,:) = 0._CUSTOM_REAL
R_xx(:,:,:,:,:) = 0._CUSTOM_REAL
R_yy(:,:,:,:,:) = 0._CUSTOM_REAL
R_xy(:,:,:,:,:) = 0._CUSTOM_REAL
@@ -582,7 +589,7 @@
R_yz(:,:,:,:,:) = 0._CUSTOM_REAL
if(FIX_UNDERFLOW_PROBLEM) then
- R_trace(:,:,:,:,:) = VERYSMALLVAL !ZN
+ R_trace(:,:,:,:,:) = VERYSMALLVAL
R_xx(:,:,:,:,:) = VERYSMALLVAL
R_yy(:,:,:,:,:) = VERYSMALLVAL
R_xy(:,:,:,:,:) = VERYSMALLVAL
@@ -841,13 +848,13 @@
! memory variables if attenuation
if( ATTENUATION ) then
- b_R_trace = 0._CUSTOM_REAL !ZN
+ b_R_trace = 0._CUSTOM_REAL
b_R_xx = 0._CUSTOM_REAL
b_R_yy = 0._CUSTOM_REAL
b_R_xy = 0._CUSTOM_REAL
b_R_xz = 0._CUSTOM_REAL
b_R_yz = 0._CUSTOM_REAL
- b_epsilondev_trace = 0._CUSTOM_REAL !ZN
+ b_epsilondev_trace = 0._CUSTOM_REAL
b_epsilondev_xx = 0._CUSTOM_REAL
b_epsilondev_yy = 0._CUSTOM_REAL
b_epsilondev_xy = 0._CUSTOM_REAL
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -92,6 +92,8 @@
if( ier /= 0 ) stop 'error allocating array potential_dot_acoustic'
allocate(potential_dot_dot_acoustic(NGLOB_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array potential_dot_dot_acoustic'
+ allocate(potential_dot_dot_acoustic_interface(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array potential_dot_dot_acoustic_interface'
if( SIMULATION_TYPE /= 1 ) then
allocate(potential_acoustic_adj_coupling(NGLOB_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array potential_acoustic_adj_coupling'
@@ -99,7 +101,10 @@
! mass matrix, density
allocate(rmass_acoustic(NGLOB_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array rmass_acoustic'
+ allocate(rmass_acoustic_interface(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rmass_acoustic_interface'
read(27) rmass_acoustic
+ read(27) rmass_acoustic_interface
! initializes mass matrix contribution
allocate(rmassz_acoustic(NGLOB_AB),stat=ier)
@@ -741,13 +746,11 @@
if( ier /= 0 ) stop 'error allocating array kappa_kl'
! anisotropic kernels
- if ( ANISOTROPIC_KL ) then
- allocate(cijkl_kl(21,NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
- if( ier /= 0 ) stop 'error allocating array cijkl_kl'
- else
- allocate(cijkl_kl(1,1,1,1,1),stat=ier)
- if( ier /= 0 ) stop 'error allocating array cijkl_kl'
- endif
+!! DK DK commented this out for now; must be made optional
+! allocate(cijkl_kl(21,NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+! if( ier /= 0 ) stop 'error allocating array cijkl_kl'
+!! DK DK added this for now
+ allocate(cijkl_kl(1,1,1,1,1),stat=ier)
! derived kernels
! density prime kernel
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2013-05-07 19:00:15 UTC (rev 21997)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2013-05-07 19:56:46 UTC (rev 21998)
@@ -90,7 +90,7 @@
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: Veloc_dsm_boundary,Tract_dsm_boundary
! attenuation
- integer :: NSPEC_ATTENUATION_AB,NSPEC_ATTENUATION_AB_kappa !ZN
+ integer :: NSPEC_ATTENUATION_AB,NSPEC_ATTENUATION_AB_kappa
character(len=256) prname_Q
! additional mass matrix for ocean load
@@ -296,17 +296,17 @@
implicit none
! memory variables and standard linear solids for attenuation
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: one_minus_sum_beta,one_minus_sum_beta_kappa !ZN
- real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: factor_common,factor_common_kappa !ZN
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: one_minus_sum_beta,one_minus_sum_beta_kappa
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: factor_common,factor_common_kappa
real(kind=CUSTOM_REAL), dimension(N_SLS) :: tau_sigma
real(kind=CUSTOM_REAL) :: min_resolved_period
real(kind=CUSTOM_REAL), dimension(N_SLS) :: &
alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: &
- R_trace,R_xx,R_yy,R_xy,R_xz,R_yz !ZN
+ R_trace,R_xx,R_yy,R_xy,R_xz,R_yz
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
- epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz !ZN
+ epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: epsilon_trace_over_3
! displacement, velocity, acceleration
@@ -353,9 +353,9 @@
real(kind=CUSTOM_REAL), dimension(N_SLS) :: &
b_alphaval, b_betaval, b_gammaval
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: &
- b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz !ZN
+ b_R_trace,b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
- b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz !ZN
+ b_epsilondev_trace,b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: b_epsilon_trace_over_3
! adjoint kernels
@@ -396,12 +396,13 @@
implicit none
! potential
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_acoustic, &
- potential_dot_acoustic,potential_dot_dot_acoustic
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_acoustic,potential_dot_acoustic, &
+ potential_dot_dot_acoustic,potential_dot_dot_acoustic_interface
real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_acoustic_adj_coupling
! mass matrix
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_acoustic
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_acoustic_interface
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassz_acoustic
! acoustic-elastic coupling surface
More information about the CIG-COMMITS
mailing list