[cig-commits] r22691 - in seismo/3D/SPECFEM3D/trunk/src: generate_databases shared specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Jul 31 09:48:32 PDT 2013
Author: dkomati1
Date: 2013-07-31 09:48:31 -0700 (Wed, 31 Jul 2013)
New Revision: 22691
Modified:
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.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/shared/check_mesh_resolution.f90
seismo/3D/SPECFEM3D/trunk/src/shared/gll_library.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.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_Dev.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/finalize_simulation.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.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:
removed useless white spaces
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -115,14 +115,14 @@
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
+ 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
+ do ispec=1,nspec
if( ispec_is_acoustic(ispec) ) then
do k=1,NGLLZ
do j=1,NGLLY
@@ -480,7 +480,7 @@
use generate_databases_par, only: NGLLX,NGLLY,NGLLZ,SIZE_REAL,CUSTOM_REAL,DT,&
is_CPML,CPML_regions,d_store_x,d_store_y,d_store_z, &
- K_store_x,K_store_y,K_store_z,nspec_cpml,CPML_to_spec,&
+ K_store_x,K_store_y,K_store_z,nspec_cpml,CPML_to_spec,&
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 : rmass,rhostore,jacobianstore,wxgll,wygll,wzgll,ispec_is_elastic
@@ -737,7 +737,7 @@
use generate_databases_par, only: NGLLX,NGLLY,NGLLZ,SIZE_REAL,CUSTOM_REAL,DT,&
is_CPML,CPML_regions,d_store_x,d_store_y,d_store_z, &
- K_store_x,K_store_y,K_store_z,nspec_cpml,CPML_to_spec,&
+ K_store_x,K_store_y,K_store_z,nspec_cpml,CPML_to_spec,&
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 : rmass_acoustic,kappastore,jacobianstore,wxgll,wygll,wzgll,ispec_is_acoustic
@@ -980,7 +980,7 @@
enddo
else
stop 'error in PML mesh file'
- endif
+ endif
endif
enddo ! do ispec_CPML=1,nspec_cpml
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -391,14 +391,14 @@
endif
write(IMAIN,*)
- if(STACEY_ABSORBING_CONDITIONS) then
+ if(STACEY_ABSORBING_CONDITIONS) then
write(IMAIN,*) 'incorporating Stacey absorbing conditions'
else
- if(PML_CONDITIONS) then
- write(IMAIN,*) 'incorporating absorbing conditions of perfectly matched layer'
- else
- write(IMAIN,*) 'no absorbing condition'
- endif
+ if(PML_CONDITIONS) then
+ write(IMAIN,*) 'incorporating absorbing conditions of perfectly matched layer'
+ else
+ write(IMAIN,*) 'no absorbing condition'
+ endif
endif
write(IMAIN,*)
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -208,9 +208,9 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass,rmass_acoustic,&
rmass_solid_poroelastic,rmass_fluid_poroelastic
-! mass matrix for interface
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_acoustic_interface
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_elastic_interface
+! mass matrix for interface
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_acoustic_interface
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_elastic_interface
! mass matrix contributions
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx,rmassy,rmassz
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-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/pml_set_local_dampingcoeff.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -70,8 +70,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,&
- vp_max,vp_max_all
+ CPML_width_z_top_max_all,CPML_width_z_bottom_max_all,&
+ vp_max,vp_max_all
! stores damping profiles
allocate(d_store_x(NGLLX,NGLLY,NGLLZ,nspec_cpml),stat=ier)
@@ -235,7 +235,7 @@
ispec = CPML_to_spec(ispec_CPML)
do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
vp = rho_vp(i,j,k,ispec)/rhostore(i,j,k,ispec)
- if(vp .ge. vp_max)then
+ if(vp >= vp_max)then
vp_max = vp
endif
enddo; enddo; enddo
@@ -276,12 +276,12 @@
! calculates P-velocity
if( ispec_is_acoustic(ispec) ) then
! vp = rho_vp(i,j,k,ispec)/rhostore(i,j,k,ispec)
-! For convenience only, when computing the damping profile inside PML,
+! 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 = rho_vp(i,j,k,ispec)/rhostore(i,j,k,ispec)
-! For convenience only, when computing the damping profile inside PML,
+! 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
@@ -1342,7 +1342,7 @@
do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
mask_ibool_interior_domain(ibool(i,j,k,ispec))=.true.
enddo; enddo; enddo
- endif
+ endif
enddo
!------------------------------------------------------
@@ -1350,13 +1350,13 @@
!------------------------------------------------------
nglob_interface_PML_acoustic = 0
- if(ACOUSTIC_SIMULATION) then
+ if(ACOUSTIC_SIMULATION) then
do ispec=1,nspec
if( ispec_is_acoustic(ispec) .and. is_CPML(ispec)) then
do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
if(mask_ibool_interior_domain(ibool(i,j,k,ispec)))then
- nglob_interface_PML_acoustic = nglob_interface_PML_acoustic + 1
+ nglob_interface_PML_acoustic = nglob_interface_PML_acoustic + 1
endif
enddo; enddo; enddo
endif
@@ -1366,17 +1366,17 @@
allocate(points_interface_PML_acoustic(nglob_interface_PML_acoustic),stat=ier)
if(ier /= 0) stop 'error allocating array points_interface_PML_acoustic'
points_interface_PML_acoustic = 0
- nglob_interface_PML_acoustic = 0
+ nglob_interface_PML_acoustic = 0
do ispec=1,nspec
if( ispec_is_acoustic(ispec) .and. is_CPML(ispec)) then
do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
if(mask_ibool_interior_domain(ibool(i,j,k,ispec)))then
- nglob_interface_PML_acoustic = nglob_interface_PML_acoustic + 1
+ nglob_interface_PML_acoustic = nglob_interface_PML_acoustic + 1
points_interface_PML_acoustic(nglob_interface_PML_acoustic) = ibool(i,j,k,ispec)
endif
enddo; enddo; enddo
endif
- enddo
+ enddo
endif
endif
@@ -1389,13 +1389,13 @@
!------------------------------------------------------
nglob_interface_PML_elastic = 0
- if(ELASTIC_SIMULATION) then
+ if(ELASTIC_SIMULATION) then
do ispec=1,nspec
if( ispec_is_elastic(ispec) .and. is_CPML(ispec)) then
do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
if(mask_ibool_interior_domain(ibool(i,j,k,ispec)))then
- nglob_interface_PML_elastic = nglob_interface_PML_elastic + 1
+ nglob_interface_PML_elastic = nglob_interface_PML_elastic + 1
endif
enddo; enddo; enddo
endif
@@ -1410,12 +1410,12 @@
if( ispec_is_elastic(ispec) .and. is_CPML(ispec)) then
do k=1,NGLLZ; do j=1,NGLLY; do i=1,NGLLX
if(mask_ibool_interior_domain(ibool(i,j,k,ispec)))then
- nglob_interface_PML_elastic = nglob_interface_PML_elastic + 1
+ nglob_interface_PML_elastic = nglob_interface_PML_elastic + 1
points_interface_PML_elastic(nglob_interface_PML_elastic) = ibool(i,j,k,ispec)
endif
enddo; enddo; enddo
endif
- enddo
+ enddo
endif
endif
@@ -1450,7 +1450,7 @@
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._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
+ ! 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._CUSTOM_REAL * delta)) &
* dist**(1.2_CUSTOM_REAL*NPOWER)
else
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -1004,7 +1004,7 @@
iglob_b = ibool(i,j+1,k,ispec)
dist = ( xstore(iglob_a) - xstore(iglob_b) )**2 &
+ ( ystore(iglob_a) - ystore(iglob_b) )**2 &
- + ( zstore(iglob_a) - zstore(iglob_b) )**2
+ + ( zstore(iglob_a) - zstore(iglob_b) )**2
if( dist < distance_min) distance_min = dist
if( dist > distance_max) distance_max = dist
enddo
@@ -1019,7 +1019,7 @@
iglob_b = ibool(i,j,k+1,ispec)
dist = ( xstore(iglob_a) - xstore(iglob_b) )**2 &
+ ( ystore(iglob_a) - ystore(iglob_b) )**2 &
- + ( zstore(iglob_a) - zstore(iglob_b) )**2
+ + ( zstore(iglob_a) - zstore(iglob_b) )**2
if( dist < distance_min) distance_min = dist
if( dist > distance_max) distance_max = dist
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/gll_library.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/gll_library.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/gll_library.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -243,7 +243,7 @@
! looks for index with minimum value
do j=i,np
- ! note: some compilers (cray) might be too aggressive in optimizing this loop,
+ ! note: some compilers (cray) might be too aggressive in optimizing this loop,
! thus we need this temporary array value x to store and compare values
x = xjac(j)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -141,7 +141,7 @@
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
! the sign is negative because pressure p = - Chi_dot_dot therefore we need
! to add minus the source to Chi_dot_dot to get plus the source in pressure:
-
+
! acoustic source for pressure gets divided by kappa
! source contribution
do k=1,NGLLZ
@@ -486,7 +486,7 @@
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
! the sign is negative because pressure p = - Chi_dot_dot therefore we need
! to add minus the source to Chi_dot_dot to get plus the source in pressure:
-
+
! acoustic source for pressure gets divided by kappa
! source contribution
do k=1,NGLLZ
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -32,10 +32,10 @@
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner,&
+ ispec_is_inner,phase_is_inner,&
PML_CONDITIONS,spec_to_CPML,is_CPML,&
potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ,&
- SIMULATION_TYPE,backward_simulation,accel_interface)
+ SIMULATION_TYPE,backward_simulation,accel_interface)
! returns the updated pressure array: potential_dot_dot_acoustic
@@ -67,10 +67,10 @@
logical, dimension(NSPEC_AB) :: ispec_is_inner
logical :: phase_is_inner
-! CPML
- logical :: PML_CONDITIONS
- integer :: spec_to_CPML(NSPEC_AB)
- logical :: is_CPML(NSPEC_AB)
+! CPML
+ logical :: PML_CONDITIONS
+ integer :: spec_to_CPML(NSPEC_AB)
+ logical :: is_CPML(NSPEC_AB)
! local parameters
real(kind=CUSTOM_REAL) :: displ_x,displ_y,displ_z,displ_n
@@ -98,7 +98,7 @@
iglob = ibool(i,j,k,ispec)
! elastic displacement on global point
- if(PML_CONDITIONS)then
+ if(PML_CONDITIONS)then
if(.not. backward_simulation)then
if(is_CPML(ispec))then
if(SIMULATION_TYPE == 1)then
@@ -120,8 +120,8 @@
displ_z = displ(3,iglob)
endif
else
- if(is_CPML(ispec))then
-! left blank, since no operation needed
+ if(is_CPML(ispec))then
+! left blank, since no operation needed
else
displ_x = displ(1,iglob)
displ_y = displ(2,iglob)
@@ -158,8 +158,8 @@
! calculating the coupling on the elastic side for the acceleration...
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + jacobianw*displ_n
- if(PML_CONDITIONS)then
- if(is_CPML(ispec))then
+ if(PML_CONDITIONS)then
+ if(is_CPML(ispec))then
potential_dot_dot_acoustic_interface(iglob) = potential_dot_dot_acoustic_interface(iglob) &
+ jacobianw*displ_n
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_viscoelastic_ac.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -36,14 +36,14 @@
PML_CONDITIONS,is_CPML,potential_dot_dot_acoustic_interface,&
SIMULATION_TYPE,backward_simulation,accel_interface,&
rmemory_coupling_el_ac_potential,spec_to_CPML, &
- potential_acoustic,potential_dot_acoustic)
+ potential_acoustic,potential_dot_acoustic)
! returns the updated acceleration array: accel
implicit none
include 'constants.h'
- integer :: NSPEC_AB,NGLOB_AB,SIMULATION_TYPE
+ integer :: NSPEC_AB,NGLOB_AB,SIMULATION_TYPE
logical :: backward_simulation
! displacement and pressure
@@ -72,16 +72,16 @@
integer :: iface,igll,ispec,iglob
integer :: i,j,k
-! CPML
+! CPML
integer :: ispec_CPML
integer :: spec_to_CPML(NSPEC_AB)
- logical :: PML_CONDITIONS
- logical :: is_CPML(NSPEC_AB)
+ logical :: PML_CONDITIONS
+ logical :: is_CPML(NSPEC_AB)
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: accel_interface
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2) :: rmemory_coupling_el_ac_potential
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic_interface
-
+
! loops on all coupling faces
do iface = 1,num_coupling_ac_el_faces
@@ -103,7 +103,7 @@
iglob = ibool(i,j,k,ispec)
! acoustic pressure on global point
- if(PML_CONDITIONS)then
+ if(PML_CONDITIONS)then
if(.not. backward_simulation)then
if(is_CPML(ispec))then
if(SIMULATION_TYPE == 1)then
@@ -111,17 +111,17 @@
endif
if(SIMULATION_TYPE == 3)then
- ispec_CPML = spec_to_CPML(ispec)
- call pml_compute_memory_variables_elastic_acoustic(ispec_CPML,iface,iglob,i,j,k,&
- pressure,potential_acoustic,potential_dot_acoustic,&
- num_coupling_ac_el_faces,rmemory_coupling_el_ac_potential)
+ ispec_CPML = spec_to_CPML(ispec)
+ call pml_compute_memory_variables_elastic_acoustic(ispec_CPML,iface,iglob,i,j,k,&
+ pressure,potential_acoustic,potential_dot_acoustic,&
+ num_coupling_ac_el_faces,rmemory_coupling_el_ac_potential)
endif
else
pressure = - potential_dot_dot_acoustic(iglob)
endif
else
if(is_CPML(ispec))then
-! left blank, since no operation needed
+! left blank, since no operation needed
else
pressure = - potential_dot_dot_acoustic(iglob)
endif
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-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -62,7 +62,7 @@
implicit none
! local parameters
- integer:: iphase,iface,ispec,iglob,igll,i,j,k
+ integer:: iphase,iface,ispec,iglob,igll,i,j,k
logical:: phase_is_inner
! enforces free surface (zeroes potentials at free surface)
@@ -108,7 +108,7 @@
rhostore,jacobian,ibool,deltat, &
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
phase_ispec_inner_acoustic,ELASTIC_SIMULATION,&
- .false.,potential_dot_dot_acoustic_interface)
+ .false.,potential_dot_dot_acoustic_interface)
endif
! ! Stacey absorbing boundary conditions
@@ -133,10 +133,10 @@
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner,&
+ ispec_is_inner,phase_is_inner,&
PML_CONDITIONS,spec_to_CPML,is_CPML,&
potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ,&
- SIMULATION_TYPE,.false.,accel_interface)
+ SIMULATION_TYPE,.false.,accel_interface)
else
! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
@@ -147,10 +147,10 @@
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner,&
+ ispec_is_inner,phase_is_inner,&
PML_CONDITIONS,spec_to_CPML,is_CPML,&
potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ,&
- SIMULATION_TYPE,.false.,accel_interface)
+ SIMULATION_TYPE,.false.,accel_interface)
endif
endif
endif
@@ -206,11 +206,11 @@
! 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(:)
+ potential_dot_dot_acoustic_interface(:) = potential_dot_dot_acoustic_interface(:) * &
+ rmass_acoustic_interface(:)
endif
endif
@@ -221,7 +221,7 @@
! 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
+! 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
@@ -243,7 +243,7 @@
potential_dot_dot_acoustic(iglob) = 0.0
potential_dot_acoustic(iglob) = 0.0
potential_acoustic(iglob) = 0.0
- if(ELASTIC_SIMULATION ) then
+ if(ELASTIC_SIMULATION ) then
potential_dot_dot_acoustic_interface(iglob) = 0.0
endif
enddo
@@ -375,7 +375,7 @@
rhostore,jacobian,ibool,deltat, &
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
phase_ispec_inner_acoustic,ELASTIC_SIMULATION,&
- .true.,potential_dot_dot_acoustic_interface)
+ .true.,potential_dot_dot_acoustic_interface)
endif
endif
@@ -401,10 +401,10 @@
coupling_ac_el_ispec,coupling_ac_el_ijk, &
coupling_ac_el_normal, &
coupling_ac_el_jacobian2Dw, &
- ispec_is_inner,phase_is_inner,&
+ ispec_is_inner,phase_is_inner,&
PML_CONDITIONS,spec_to_CPML,is_CPML,&
potential_dot_dot_acoustic_interface,veloc,rmemory_coupling_ac_el_displ,&
- SIMULATION_TYPE,.true.,accel_interface)
+ SIMULATION_TYPE,.true.,accel_interface)
endif
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -56,7 +56,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
@@ -83,7 +83,7 @@
! CPML fluid-solid interface
logical :: ELASTIC_SIMULATION
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic_interface
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic_interface
! CPML adjoint
logical :: backward_simulation
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_Dev.F90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_Dev.F90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -35,11 +35,11 @@
kappastore,mustore,jacobian,ibool, &
ATTENUATION,deltat, &
one_minus_sum_beta,factor_common,&
- one_minus_sum_beta_kappa,factor_common_kappa,&
- alphaval,betaval,gammaval,&
+ one_minus_sum_beta_kappa,factor_common_kappa,&
+ alphaval,betaval,gammaval,&
NSPEC_ATTENUATION_AB,NSPEC_ATTENUATION_AB_Kappa, &
- R_trace,R_xx,R_yy,R_xy,R_xz,R_yz, &
- epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy, &
+ R_trace,R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilondev_trace,epsilondev_xx,epsilondev_yy,epsilondev_xy, &
epsilondev_xz,epsilondev_yz,epsilon_trace_over_3, &
ANISOTROPY,NSPEC_ANISO, &
c11store,c12store,c13store,c14store,c15store,c16store,&
@@ -94,17 +94,17 @@
integer :: NSPEC_ATTENUATION_AB,NSPEC_ATTENUATION_AB_Kappa
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB) :: one_minus_sum_beta
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB) :: factor_common
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_Kappa) :: one_minus_sum_beta_kappa
- real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_Kappa) :: factor_common_kappa
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_Kappa) :: one_minus_sum_beta_kappa
+ real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_Kappa) :: factor_common_kappa
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS) :: &
R_xx,R_yy,R_xy,R_xz,R_yz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_Kappa,N_SLS) :: R_trace
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_Kappa,N_SLS) :: R_trace
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY) :: &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_Kappa) :: epsilondev_trace
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB_Kappa) :: epsilondev_trace
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: epsilon_trace_over_3
! anisotropy
@@ -198,10 +198,10 @@
equivalence(tempz3_att,C3_mxm_m2_m1_5points_att)
! local attenuation parameters
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_trace_loc,epsilondev_xx_loc, &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_trace_loc,epsilondev_xx_loc, &
epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
- real(kind=CUSTOM_REAL) R_xx_val1,R_yy_val1,R_xx_val2,R_yy_val2,R_xx_val3,R_yy_val3, &
- R_trace_val1,R_trace_val2,R_trace_val3
+ real(kind=CUSTOM_REAL) R_xx_val1,R_yy_val1,R_xx_val2,R_yy_val2,R_xx_val3,R_yy_val3, &
+ R_trace_val1,R_trace_val2,R_trace_val3
real(kind=CUSTOM_REAL) factor_loc,alphaval_loc,betaval_loc,gammaval_loc
real(kind=CUSTOM_REAL) Sn,Snp1
real(kind=CUSTOM_REAL) templ
@@ -563,7 +563,7 @@
if(ATTENUATION) then
! use unrelaxed parameters if attenuation
mul = mul * one_minus_sum_beta(i,j,k,ispec)
- if(FULL_ATTENUATION_SOLID) kappal = kappal * one_minus_sum_beta_kappa(i,j,k,ispec)
+ if(FULL_ATTENUATION_SOLID) kappal = kappal * one_minus_sum_beta_kappa(i,j,k,ispec)
endif
! full anisotropic case, stress calculations
@@ -823,13 +823,13 @@
gammaval_loc = gammaval(i_sls)
if(FULL_ATTENUATION_SOLID) then
- ! term in trace
- factor_loc = kappastore(i,j,k,ispec) * factor_common_kappa(i_sls,i,j,k,ispec)
+ ! term in trace
+ factor_loc = kappastore(i,j,k,ispec) * factor_common_kappa(i_sls,i,j,k,ispec)
- Sn = factor_loc * epsilondev_trace(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_trace_loc(i,j,k)
- R_trace(i,j,k,ispec,i_sls) = alphaval_loc * R_trace(i,j,k,ispec,i_sls) + &
- betaval_loc * Sn + gammaval_loc * Snp1
+ Sn = factor_loc * epsilondev_trace(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_trace_loc(i,j,k)
+ R_trace(i,j,k,ispec,i_sls) = alphaval_loc * R_trace(i,j,k,ispec,i_sls) + &
+ betaval_loc * Sn + gammaval_loc * Snp1
endif
! term in xx yy zz xy xz yz
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-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -115,13 +115,13 @@
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,&
SIMULATION_TYPE,.false.,accel_interface,&
rmemory_coupling_el_ac_potential,spec_to_CPML,&
- potential_acoustic,potential_dot_acoustic)
-
+ potential_acoustic,potential_dot_acoustic)
+
else
! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
! adoint definition: pressure^\dagger=potential^\dagger
@@ -131,12 +131,12 @@
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,&
SIMULATION_TYPE,.false.,accel_interface,&
rmemory_coupling_el_ac_potential,spec_to_CPML,&
- potential_acoustic,potential_dot_acoustic)
-
+ potential_acoustic,potential_dot_acoustic)
+
endif
endif ! num_coupling_ac_el_faces
@@ -292,7 +292,7 @@
!
!=====================================================================
-! elastic solver for back
+! elastic solver for back
subroutine compute_forces_viscoelastic_bpwf()
@@ -384,12 +384,12 @@
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,&
SIMULATION_TYPE,.true.,accel_interface,&
rmemory_coupling_el_ac_potential,spec_to_CPML,&
- potential_acoustic,potential_dot_acoustic)
-
+ potential_acoustic,potential_dot_acoustic)
+
endif ! num_coupling_ac_el_faces
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -736,7 +736,7 @@
! because array is_CPML() is unallocated when PML_CONDITIONS is false
if(is_CPML(ispec)) then
! sets C-PML elastic memory variables to compute stress sigma and form dot product with test vector
- call pml_compute_memory_variables_elastic(ispec,ispec_CPML,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
+ call pml_compute_memory_variables_elastic(ispec,ispec_CPML,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
tempx3,tempy3,tempz3, &
rmemory_dux_dxl_x, rmemory_duy_dyl_x, rmemory_duz_dzl_x, &
rmemory_dux_dyl_x, rmemory_dux_dzl_x, rmemory_duz_dxl_x, rmemory_duy_dxl_x, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -187,7 +187,7 @@
deallocate(alpha_store)
deallocate(spec_to_CPML)
deallocate(CPML_type)
-
+
if( ELASTIC_SIMULATION ) then
deallocate(PML_dux_dxl)
deallocate(PML_dux_dyl)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -152,15 +152,15 @@
if( ATTENUATION ) then
!pll
NSPEC_ATTENUATION_AB = NSPEC_AB
- if(FULL_ATTENUATION_SOLID) then
- NSPEC_ATTENUATION_AB_kappa = NSPEC_AB
- else
- NSPEC_ATTENUATION_AB_kappa = 1
- endif
+ if(FULL_ATTENUATION_SOLID) then
+ NSPEC_ATTENUATION_AB_kappa = NSPEC_AB
+ else
+ NSPEC_ATTENUATION_AB_kappa = 1
+ endif
else
! if attenuation is off, set dummy size of arrays to one
NSPEC_ATTENUATION_AB = 1
- NSPEC_ATTENUATION_AB_kappa = 1
+ NSPEC_ATTENUATION_AB_kappa = 1
endif
! needed for attenuation and/or kernel computations
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.F90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -817,7 +817,7 @@
if (SIMULATION_TYPE == 3 .and. .NOT. GPU_MODE) then
! acoustic backward fields
if( ACOUSTIC_SIMULATION ) then
- if(PML_CONDITIONS)then
+ if(PML_CONDITIONS)then
if(nglob_interface_PML_acoustic > 0)then
call read_potential_on_pml_interface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic,b_potential_acoustic,&
nglob_interface_PML_acoustic,b_PML_potential,b_reclen_PML_potential)
@@ -833,7 +833,7 @@
! elastic backward fields
if( ELASTIC_SIMULATION ) then
- if(PML_CONDITIONS)then
+ if(PML_CONDITIONS)then
if(nglob_interface_PML_elastic > 0)then
call read_field_on_pml_interface(b_accel,b_veloc,b_displ,nglob_interface_PML_elastic,&
b_PML_field,b_reclen_PML_field)
@@ -1057,13 +1057,13 @@
size(epsilondev_xx))
endif
- if(FULL_ATTENUATION_SOLID) write(27) R_trace
+ if(FULL_ATTENUATION_SOLID) write(27) R_trace
write(27) R_xx
write(27) R_yy
write(27) R_xy
write(27) R_xz
write(27) R_yz
- if(FULL_ATTENUATION_SOLID) write(27) epsilondev_trace
+ if(FULL_ATTENUATION_SOLID) write(27) epsilondev_trace
write(27) epsilondev_xx
write(27) epsilondev_yy
write(27) epsilondev_xy
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -32,7 +32,7 @@
use specfem_par, only: NSPEC_AB,PML_CONDITIONS,SIMULATION_TYPE,SAVE_FORWARD,NSTEP,myrank,prname !
use constants, only: NDIM,NGLLX,NGLLY,NGLLZ
use specfem_par_acoustic, only: ACOUSTIC_SIMULATION,num_coupling_ac_el_faces
- use specfem_par_elastic, only: ELASTIC_SIMULATION
+ use specfem_par_elastic, only: ELASTIC_SIMULATION
implicit none
@@ -177,16 +177,16 @@
endif
! stores C-PML contribution on elastic/acoustic interface
- if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
- allocate(rmemory_coupling_ac_el_displ(3,NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_coupling_ac_el_displ array'
- endif
+ if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
+ allocate(rmemory_coupling_ac_el_displ(3,NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_coupling_ac_el_displ array'
+ endif
if(SIMULATION_TYPE == 3)then
- if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
- allocate(rmemory_coupling_el_ac_potential(NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_coupling_el_ac_potential array'
- endif
+ if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
+ allocate(rmemory_coupling_el_ac_potential(NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2),stat=ier)
+ if(ier /= 0) stop 'error allocating rmemory_coupling_el_ac_potential array'
+ endif
endif
spec_to_CPML(:) = 0
@@ -259,14 +259,14 @@
potential_dot_dot_acoustic_CPML(:,:,:) = 0._CUSTOM_REAL
endif
- if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
- rmemory_coupling_ac_el_displ(:,:,:,:,:,:) = 0._CUSTOM_REAL
- endif
+ if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
+ rmemory_coupling_ac_el_displ(:,:,:,:,:,:) = 0._CUSTOM_REAL
+ endif
if(SIMULATION_TYPE == 3)then
- if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
- rmemory_coupling_el_ac_potential(:,:,:,:,:) = 0._CUSTOM_REAL
- endif
+ if(ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION)then
+ rmemory_coupling_el_ac_potential(:,:,:,:,:) = 0._CUSTOM_REAL
+ endif
endif
! fields on PML interface will be reconstructed for adjoint simulations
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -37,7 +37,7 @@
use specfem_par, only: NGLOB_AB,it,deltat,wgll_cube,jacobian,ibool,rhostore
use pml_par, only: CPML_regions,d_store_x,d_store_y,d_store_z,K_store_x,K_store_y,K_store_z,&
- alpha_store,NSPEC_CPML,accel_elastic_CPML
+ alpha_store,NSPEC_CPML,accel_elastic_CPML
use constants, only: CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ,CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY, &
CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ
@@ -426,8 +426,8 @@
! A3 = temp_A3 + (it+0.0)*deltat*temp_A4 + ((it+0.0)*deltat)**2*temp_A5
! A4 = -temp_A4 -2.0*(it+0.0)*deltat*temp_A5
! A5 = temp_A5
-!!! the full experssion of A3,A4,A5 are given by above equation, here we use reduced
-!!! exprssion of A3,A4,A5 in order to stabilized the code.
+!!! the full experssion of A3,A4,A5 are given by above equation, here we use reduced
+!!! exprssion of A3,A4,A5 in order to stabilized the code.
A3 = temp_A3
A4 = 0.0
@@ -464,7 +464,7 @@
!
!=====================================================================
!
-!
+!
subroutine pml_compute_accel_contribution_acoustic(ispec,ispec_CPML,potential_acoustic,&
potential_dot_acoustic,rmemory_potential_acoustic)
@@ -485,7 +485,7 @@
implicit none
integer, intent(in) :: ispec,ispec_CPML
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB), intent(in) :: potential_acoustic,potential_dot_acoustic
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB), intent(in) :: potential_acoustic,potential_dot_acoustic
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,3) :: rmemory_potential_acoustic
@@ -786,8 +786,8 @@
! A4 = -temp_A4 -2.0*it*deltat*temp_A5
! A5 = temp_A5
-!!! the full experssion of A3,A4,A5 are given by above equation, here we use reduced
-!!! exprssion of A3,A4,A5 in order to stabilized the code.
+!!! the full experssion of A3,A4,A5 are given by above equation, here we use reduced
+!!! exprssion of A3,A4,A5 in order to stabilized the code.
A3 = temp_A3
A4 = 0.0
@@ -823,17 +823,17 @@
integer :: iglob
do iglob = 1, nglob_interface_PML_elastic
- b_PML_field(1,iglob) = displ(1,iglob)
- b_PML_field(2,iglob) = displ(2,iglob)
+ b_PML_field(1,iglob) = displ(1,iglob)
+ b_PML_field(2,iglob) = displ(2,iglob)
b_PML_field(3,iglob) = displ(3,iglob)
- b_PML_field(4,iglob) = veloc(1,iglob)
- b_PML_field(5,iglob) = veloc(2,iglob)
- b_PML_field(6,iglob) = veloc(3,iglob)
+ b_PML_field(4,iglob) = veloc(1,iglob)
+ b_PML_field(5,iglob) = veloc(2,iglob)
+ b_PML_field(6,iglob) = veloc(3,iglob)
b_PML_field(7,iglob) = accel(1,iglob)
- b_PML_field(8,iglob) = accel(2,iglob)
- b_PML_field(9,iglob) = accel(3,iglob)
+ b_PML_field(8,iglob) = accel(2,iglob)
+ b_PML_field(9,iglob) = accel(3,iglob)
enddo
call write_abs(0,b_PML_field,b_reclen_PML_field,it)
@@ -861,25 +861,25 @@
do i = 1, NGLLX; do j = 1, NGLLY; do k = 1, NGLLZ
iglob = ibool(i,j,k,ispec)
b_displ(:,iglob) = 0._CUSTOM_REAL
- b_veloc(:,iglob) = 0._CUSTOM_REAL
+ b_veloc(:,iglob) = 0._CUSTOM_REAL
b_accel(:,iglob) = 0._CUSTOM_REAL
enddo; enddo; enddo
- enddo
+ enddo
call read_abs(0,b_PML_field,b_reclen_PML_field,NSTEP-it+1)
do iglob = 1, nglob_interface_PML_elastic
- b_displ(1,iglob) = b_PML_field(1,iglob)
- b_displ(2,iglob) = b_PML_field(2,iglob)
+ b_displ(1,iglob) = b_PML_field(1,iglob)
+ b_displ(2,iglob) = b_PML_field(2,iglob)
b_displ(3,iglob) = b_PML_field(3,iglob)
- b_veloc(1,iglob) = b_PML_field(4,iglob)
- b_veloc(2,iglob) = b_PML_field(5,iglob)
- b_veloc(3,iglob) = b_PML_field(6,iglob)
+ b_veloc(1,iglob) = b_PML_field(4,iglob)
+ b_veloc(2,iglob) = b_PML_field(5,iglob)
+ b_veloc(3,iglob) = b_PML_field(6,iglob)
b_accel(1,iglob) = b_PML_field(7,iglob)
- b_accel(2,iglob) = b_PML_field(8,iglob)
- b_accel(3,iglob) = b_PML_field(9,iglob)
+ b_accel(2,iglob) = b_PML_field(8,iglob)
+ b_accel(3,iglob) = b_PML_field(9,iglob)
enddo
end subroutine read_field_on_pml_interface
@@ -898,11 +898,11 @@
real(kind=CUSTOM_REAL), dimension(3,nglob_interface_PML_acoustic) :: b_PML_potential
integer :: iglob
-
+
do iglob = 1, nglob_interface_PML_acoustic
b_PML_potential(1,iglob) = potential_acoustic(iglob)
- b_PML_potential(2,iglob) = potential_dot_acoustic(iglob)
- b_PML_potential(3,iglob) = potential_dot_dot_acoustic(iglob)
+ b_PML_potential(2,iglob) = potential_dot_acoustic(iglob)
+ b_PML_potential(3,iglob) = potential_dot_dot_acoustic(iglob)
enddo
call write_abs(1,b_PML_potential,b_reclen_PML_potential,it)
@@ -930,17 +930,17 @@
do i = 1, NGLLX; do j = 1, NGLLY; do k = 1, NGLLZ
iglob = ibool(i,j,k,ispec)
b_potential_acoustic(iglob) = 0._CUSTOM_REAL
- b_potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
- b_potential_dot_dot_acoustic(iglob) = 0._CUSTOM_REAL
+ b_potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
+ b_potential_dot_dot_acoustic(iglob) = 0._CUSTOM_REAL
enddo; enddo; enddo
- enddo
+ enddo
call read_abs(1,b_PML_potential,b_reclen_PML_potential,NSTEP-it+1)
do iglob = 1, nglob_interface_PML_acoustic
- b_potential_acoustic(iglob) = b_PML_potential(1,iglob)
- b_potential_dot_acoustic(iglob) = b_PML_potential(2,iglob)
- b_potential_dot_dot_acoustic(iglob) = b_PML_potential(3,iglob)
+ b_potential_acoustic(iglob) = b_PML_potential(1,iglob)
+ b_potential_dot_acoustic(iglob) = b_PML_potential(2,iglob)
+ b_potential_dot_dot_acoustic(iglob) = b_PML_potential(3,iglob)
enddo
end subroutine read_potential_on_pml_interface
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -63,7 +63,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,2) :: &
rmemory_dux_dxl_x, rmemory_dux_dyl_x, rmemory_dux_dzl_x, &
rmemory_duy_dxl_y, rmemory_duy_dyl_y, rmemory_duy_dzl_y, &
- rmemory_duz_dxl_z, rmemory_duz_dyl_z, rmemory_duz_dzl_z
+ rmemory_duz_dxl_z, rmemory_duz_dyl_z, rmemory_duz_dzl_z
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CPML) :: &
rmemory_duy_dyl_x, rmemory_duz_dzl_x, rmemory_duz_dxl_x, rmemory_duy_dxl_x, &
@@ -2028,7 +2028,7 @@
integer, intent(in) :: ispec_CPML,iface,iglob,num_coupling_ac_el_faces
real(kind=CUSTOM_REAL) :: displ_x,displ_y,displ_z
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB), intent(in) :: displ,veloc
- real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2) :: &
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2) :: &
rmemory_coupling_ac_el_displ
! local parameters
@@ -2052,7 +2052,7 @@
rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = 0.d0
displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
- + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
! displ_y
A9 = k_store_x(i,j,k,ispec_CPML)
@@ -2075,7 +2075,7 @@
rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = 0.d0
displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
- + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
! displ_z
A12 = k_store_x(i,j,k,ispec_CPML)
@@ -2098,7 +2098,7 @@
rmemory_coupling_ac_el_displ(3,i,j,k,iface,2) = 0.d0
displ_z = A12 * displ(3,iglob) + A13 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) &
- + A14 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2)
+ + A14 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2)
else if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
@@ -2127,7 +2127,7 @@
rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = 0.d0
displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
- + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
! displ_y
A9 = 1.d0
@@ -2138,7 +2138,7 @@
rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = 0.d0
displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
- + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
! displ_z
A12 = k_store_y(i,j,k,ispec_CPML)
@@ -2162,7 +2162,7 @@
rmemory_coupling_ac_el_displ(3,i,j,k,iface,2) = 0.d0
displ_z = A12 * displ(3,iglob) + A13 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,1) &
- + A14 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2)
+ + A14 * rmemory_coupling_ac_el_displ(3,i,j,k,iface,2)
else if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
@@ -2192,7 +2192,7 @@
rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = 0.d0
displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
- + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
! displ_y
A9 = k_store_z(i,j,k,ispec_CPML)
@@ -2216,7 +2216,7 @@
rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = 0.d0
displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
- + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
! displ_z
A12 = 1.d0
@@ -2256,7 +2256,7 @@
rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = 0.d0
displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
- + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
! displ_y
A9 = k_store_x(i,j,k,ispec_CPML)
@@ -2279,7 +2279,7 @@
rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = 0.d0
displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
- + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
! displ_z
A12 = k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML)
@@ -2341,7 +2341,7 @@
rmemory_coupling_ac_el_displ(1,i,j,k,iface,2) = 0.d0
displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
- + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
! displ_y
A9 = k_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML)
@@ -2373,7 +2373,7 @@
+ (displ(2,iglob)) * it*deltat * coef2_2
displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
- + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
! displ_z
A12 = k_store_x(i,j,k,ispec_CPML)
@@ -2434,7 +2434,7 @@
+ (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_2 + (displ(1,iglob)) * coef2_2
displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
- + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
! displ_y
A9 = k_store_z(i,j,k,ispec_CPML)
@@ -2457,7 +2457,7 @@
rmemory_coupling_ac_el_displ(2,i,j,k,iface,2) = 0.d0
displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
- + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
! displ_z
A12 = k_store_y(i,j,k,ispec_CPML)
@@ -2516,7 +2516,7 @@
+ (displ(1,iglob) + deltat * veloc(1,iglob)) * coef1_2 + (displ(1,iglob)) * coef2_2
displ_x = A6 * displ(1,iglob) + A7 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,1) &
- + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
+ + A8 * rmemory_coupling_ac_el_displ(1,i,j,k,iface,2)
! displ_y
A9 = k_store_x(i,j,k,ispec_CPML) * k_store_z(i,j,k,ispec_CPML)
@@ -2548,7 +2548,7 @@
+ (displ(2,iglob)) * it*deltat * coef2_2
displ_y = A9 * displ(2,iglob) + A10 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,1) &
- + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
+ + A11 * rmemory_coupling_ac_el_displ(2,i,j,k,iface,2)
! displ_z
A12 = k_store_x(i,j,k,ispec_CPML) * k_store_y(i,j,k,ispec_CPML)
@@ -2609,7 +2609,7 @@
integer, intent(in) :: ispec_CPML,iface,iglob,num_coupling_ac_el_faces
real(kind=CUSTOM_REAL) :: pressure
real(kind=CUSTOM_REAL), dimension(NGLOB_AB), intent(in) :: potential_acoustic,potential_dot_acoustic
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,num_coupling_ac_el_faces,2) :: &
rmemory_coupling_el_ac_potential
! local parameters
@@ -2646,7 +2646,7 @@
rmemory_coupling_el_ac_potential(i,j,k,iface,2) = 0.d0
pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
- + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
else if( CPML_regions(ispec_CPML) == CPML_Y_ONLY ) then
@@ -2676,7 +2676,7 @@
rmemory_coupling_el_ac_potential(i,j,k,iface,2) = 0.d0
pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
- + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
else if( CPML_regions(ispec_CPML) == CPML_Z_ONLY ) then
@@ -2693,7 +2693,7 @@
rmemory_coupling_el_ac_potential(i,j,k,iface,2) = 0.d0
pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
- + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
else if( CPML_regions(ispec_CPML) == CPML_XY_ONLY ) then
@@ -2731,7 +2731,7 @@
+ (potential_acoustic(iglob)) * it*deltat * coef2_2
pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
- + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
else if( CPML_regions(ispec_CPML) == CPML_XZ_ONLY ) then
@@ -2762,7 +2762,7 @@
rmemory_coupling_el_ac_potential(i,j,k,iface,2) = 0.d0
pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
- + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
else if( CPML_regions(ispec_CPML) == CPML_YZ_ONLY ) then
@@ -2793,7 +2793,7 @@
rmemory_coupling_el_ac_potential(i,j,k,iface,2) = 0.d0
pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
- + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
else if( CPML_regions(ispec_CPML) == CPML_XYZ ) then
@@ -2831,7 +2831,7 @@
+ (potential_acoustic(iglob)) * it*deltat * coef2_2
pressure = A12 * potential_acoustic(iglob) + A13 * rmemory_coupling_el_ac_potential(i,j,k,iface,1) &
- + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
+ + A14 * rmemory_coupling_el_ac_potential(i,j,k,iface,2)
else
stop 'wrong PML flag in PML memory variable calculation routine'
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_par.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -115,10 +115,10 @@
! for adjoint tomography
! need the array stored the points on interface between PML and interior computational domain
! --------------------------------------------------------------------------------------------
- integer :: nglob_interface_PML_acoustic,nglob_interface_PML_elastic
- integer, dimension(:), allocatable :: points_interface_PML_acoustic, points_interface_PML_elastic
+ integer :: nglob_interface_PML_acoustic,nglob_interface_PML_elastic
+ integer, dimension(:), allocatable :: points_interface_PML_acoustic, points_interface_PML_elastic
- integer :: b_reclen_PML_field,b_reclen_PML_potential
+ integer :: b_reclen_PML_field,b_reclen_PML_potential
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: b_PML_field,b_PML_potential
end module pml_par
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -261,8 +261,8 @@
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(:)
+ 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
@@ -303,14 +303,14 @@
rmassz(:) = 1._CUSTOM_REAL / rmassz(:)
if(PML_CONDITIONS)then
- if(ACOUSTIC_SIMULATION)then
- call assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,rmass_elastic_interface, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ if(ACOUSTIC_SIMULATION)then
+ call assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,rmass_elastic_interface, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
my_neighbours_ext_mesh)
where(rmass_elastic_interface <= 0._CUSTOM_REAL) rmass_elastic_interface = 1._CUSTOM_REAL
rmass_elastic_interface(:) = 1._CUSTOM_REAL / rmass_elastic_interface(:)
- endif
+ endif
endif
! ocean load
@@ -485,7 +485,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
+ 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
@@ -499,11 +499,11 @@
if( ier /= 0 ) call exit_mpi(myrank,'error allocation scale_factor')
scale_factor(:,:,:,:) = 1._CUSTOM_REAL
- 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
+ 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', &
@@ -522,11 +522,11 @@
read(27) factor_common
read(27) scale_factor
- if(FULL_ATTENUATION_SOLID)then
- read(27) one_minus_sum_beta_kappa
- read(27) factor_common_kappa
- read(27) scale_factor_kappa
- endif
+ 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)
@@ -560,11 +560,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
+ 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
+ endif
enddo
enddo
@@ -572,7 +572,7 @@
enddo
deallocate(scale_factor)
- deallocate(scale_factor_kappa)
+ deallocate(scale_factor_kappa)
! statistics
! user output
@@ -590,14 +590,14 @@
! clear memory variables if attenuation
! initialize memory variables for attenuation
- epsilondev_trace(:,:,:,:) = 0._CUSTOM_REAL
+ 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
+ R_trace(:,:,:,:,:) = 0._CUSTOM_REAL
R_xx(:,:,:,:,:) = 0._CUSTOM_REAL
R_yy(:,:,:,:,:) = 0._CUSTOM_REAL
R_xy(:,:,:,:,:) = 0._CUSTOM_REAL
@@ -605,7 +605,7 @@
R_yz(:,:,:,:,:) = 0._CUSTOM_REAL
if(FIX_UNDERFLOW_PROBLEM) then
- R_trace(:,:,:,:,:) = VERYSMALLVAL
+ R_trace(:,:,:,:,:) = VERYSMALLVAL
R_xx(:,:,:,:,:) = VERYSMALLVAL
R_yy(:,:,:,:,:) = VERYSMALLVAL
R_xy(:,:,:,:,:) = VERYSMALLVAL
@@ -863,13 +863,13 @@
! memory variables if attenuation
if( ATTENUATION ) then
- b_R_trace = 0._CUSTOM_REAL
+ 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
+ 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-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -92,8 +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'
+ 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'
@@ -101,10 +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'
+ 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
+ read(27) rmass_acoustic_interface
! initializes mass matrix contribution
allocate(rmassz_acoustic(NGLOB_AB),stat=ier)
@@ -140,24 +140,24 @@
! allocates mass matrix
allocate(rmass(NGLOB_AB),stat=ier)
if(PML_CONDITIONS)then ! need to be optimized
- if(ACOUSTIC_SIMULATION)then
- allocate(rmass_elastic_interface(NGLOB_AB),stat=ier)
- if( ier /= 0 ) stop 'error allocating array rmass_elastic_interface'
+ if(ACOUSTIC_SIMULATION)then
+ allocate(rmass_elastic_interface(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rmass_elastic_interface'
rmass_elastic_interface(:) = 0._CUSTOM_REAL
if(SIMULATION_TYPE == 3)then
- allocate(accel_interface(NDIM,NGLOB_AB),stat=ier)
+ allocate(accel_interface(NDIM,NGLOB_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array accel_interface'
accel_interface(:,:) = 0._CUSTOM_REAL
else
- allocate(accel_interface(NDIM,1),stat=ier)
+ allocate(accel_interface(NDIM,1),stat=ier)
if( ier /= 0 ) stop 'error allocating array accel_interface'
accel_interface(:,:) = 0._CUSTOM_REAL
endif
else
allocate(rmass_elastic_interface(1),stat=ier)
allocate(accel_interface(NDIM,1),stat=ier)
- endif
- endif
+ endif
+ endif
if( ier /= 0 ) stop 'error allocating array rmass'
! initializes mass matrix contributions
allocate(rmassx(NGLOB_AB), &
@@ -235,11 +235,11 @@
if( ier /= 0 ) stop 'error reading in array rmass'
if(PML_CONDITIONS)then !need to be optimized
- if(ACOUSTIC_SIMULATION)then
- read(27,iostat=ier) rmass_elastic_interface
- if( ier /= 0 ) stop 'error reading in array rmass_elastic_interface'
- endif
- endif
+ if(ACOUSTIC_SIMULATION)then
+ read(27,iostat=ier) rmass_elastic_interface
+ if( ier /= 0 ) stop 'error reading in array rmass_elastic_interface'
+ endif
+ endif
if( APPROXIMATE_OCEAN_LOAD ) then
! ocean mass matrix
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2013-07-31 15:02:27 UTC (rev 22690)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2013-07-31 16:48:31 UTC (rev 22691)
@@ -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
+ 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
- real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: factor_common,factor_common_kappa
+ 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
+ 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
+ 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
@@ -317,8 +317,8 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass
! PML on fluid-solid interface
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_interface
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_elastic_interface
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_interface
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_elastic_interface
! Stacey
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassx,rmassy,rmassz
@@ -357,9 +357,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
+ 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
+ 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
@@ -401,12 +401,12 @@
! potential
real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_acoustic,potential_dot_acoustic, &
- potential_dot_dot_acoustic,potential_dot_dot_acoustic_interface
+ 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 :: rmass_acoustic_interface
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmassz_acoustic
! acoustic-elastic coupling surface
More information about the CIG-COMMITS
mailing list