[cig-commits] r19586 - in seismo/3D/SPECFEM3D/trunk/src: decompose_mesh_SCOTCH generate_databases shared specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Fri Feb 3 06:39:55 PST 2012
Author: danielpeter
Date: 2012-02-03 06:39:54 -0800 (Fri, 03 Feb 2012)
New Revision: 19586
Modified:
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_arrays_source.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_po.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_ac.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_po.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_ac.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_el.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_pot.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_gradient.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
Log:
updates smooth_vol_data.f90; cleans trailing spaces in source code files
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -120,7 +120,7 @@
implicit none
character(len=256) :: line
logical :: use_poroelastic_file
-
+
! sets number of nodes per element
ngnod = esize
@@ -283,7 +283,7 @@
! mufr : frame shear modulus
open(unit=97, file=localpath_name(1:len_trim(localpath_name))//'/nummaterial_poroelastic_file', &
status='old', form='formatted', iostat=ier)
- ! checks if we can use file
+ ! checks if we can use file
if( ier /= 0 ) then
use_poroelastic_file = .false.
!stop 'error opening nummaterial_poroelastic_file'
@@ -292,7 +292,7 @@
print*, ' poroelastic material file found'
endif
ier = 0
-
+
! note: entries in nummaterial_velocity_file can be an unsorted list of all
! defined materials (material_id > 0) and undefined materials (material_id < 0 )
do imat=1,count_def_mat
@@ -316,7 +316,7 @@
! checks material_id bounds
if(num_mat < 1 .or. num_mat > count_def_mat) stop "ERROR : Invalid nummaterial_velocity_file file."
- if(idomain_id == 1 .or. idomain_id == 2) then
+ if(idomain_id == 1 .or. idomain_id == 2) then
! material is elastic or acoustic
!read(98,*) num_mat, mat_prop(1,num_mat),mat_prop(2,num_mat),&
@@ -328,10 +328,10 @@
mat_prop(5,num_mat) = aniso_flag
mat_prop(6,num_mat) = idomain_id
- else
+ else
! material is poroelastic
if( use_poroelastic_file .eqv. .false. ) stop 'error poroelastic material requires nummaterial_poroelastic_file'
-
+
read(97,*) rhos,rhof,phi,tort,kxx,kxy,kxz,kyy,kyz,kzz,kappas,kappaf,kappafr,eta,mufr
mat_prop(1,num_mat) = rhos
mat_prop(2,num_mat) = rhof
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -416,7 +416,7 @@
integer,dimension(:,:,:),allocatable :: tmp_ijk
integer,dimension(:),allocatable :: tmp_ispec
- integer,dimension(NGNOD2D) :: iglob_corners_ref
+ integer,dimension(NGNOD2D) :: iglob_corners_ref
integer :: ispec,i,j,k,igll,ier,iglob
integer :: inum,iface_ref,icorner
integer :: count_poroelastic,count_acoustic
@@ -852,7 +852,7 @@
iface_ref_el = iface_el ![CM]: for some reason this shows a wrong orientation
! but the calcul is ok.
ispec_ref_el = ispec_el
-
+
! gets face GLL points i,j,k indices from poroelastic element face
call get_element_face_gll_indices(iface_ref,ijk_face_po,NGLLX,NGLLY)
! gets face GLL points i,j,k indices from elastic element face
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -307,7 +307,7 @@
write(IOUT) coupling_el_po_normal
endif
-!MPI interfaces
+ !MPI interfaces
max_nibool_interfaces_ext_mesh = maxval(nibool_interfaces_ext_mesh(:))
allocate(ibool_interfaces_ext_mesh_dummy(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -151,7 +151,7 @@
cmax_glob = max(cmax_glob,cmax)
! debug: for vtk output
- if( SAVE_MESH_FILES ) tmp1(ispec) = cmax
+ if( SAVE_MESH_FILES ) tmp1(ispec) = cmax
endif
! suggested timestep
@@ -331,7 +331,7 @@
! debug: for vtk output
if( SAVE_MESH_FILES ) then
- call create_name_database(prname,myrank,LOCAL_PATH)
+ call create_name_database(prname,myrank,LOCAL_PATH)
! courant number
if( DT_PRESENT ) then
filename = trim(prname)//'res_courant_number'
@@ -345,7 +345,7 @@
xstore,ystore,zstore,ibool, &
tmp2,filename)
deallocate(tmp1,tmp2)
- endif
+ endif
end subroutine check_mesh_resolution
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -72,6 +72,8 @@
real(kind=CUSTOM_REAL), dimension(:),allocatable :: dummy_1
real(kind=CUSTOM_REAL), dimension(:,:),allocatable :: dummy_2
real(kind=CUSTOM_REAL), dimension(:,:,:),allocatable :: dummy_3
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:),allocatable :: dummy_5
+
integer, dimension(:),allocatable :: idummy
logical, dimension(:),allocatable :: ldummy
integer, dimension(:,:,:),allocatable :: idummy_3
@@ -88,7 +90,6 @@
integer :: i,j,k,ios,it,iglob,ier,ispec2,ispec
integer :: iproc, node_list(300)
- integer :: np
character(len=256) :: arg(7), filename, indir, outdir
character(len=256) :: prname, prname_lp
@@ -145,7 +146,7 @@
call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
- if (myrank == 0) print*,"smooth:"
+ if (myrank == 0) print*,"smooth_vol_data:"
call mpi_barrier(MPI_COMM_WORLD,ier)
! reads arguments
@@ -220,7 +221,18 @@
NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY)
- if (sizeprocs /= NPROC) call exit_mpi(myrank,'Error total number of slices')
+ ! checks if number of MPI process as specified
+ if (sizeprocs /= NPROC) then
+ if( myrank == 0 ) then
+ print*,''
+ print*,'error: run xsmooth_vol_data with the same number of MPI processes '
+ print*,' as specified in Par_file by NPROC when slices were created'
+ print*,''
+ print*,'for example: mpirun -np ',NPROC,' ./xsmooth_vol_data ...'
+ print*,''
+ endif
+ call exit_mpi(myrank,'Error total number of slices')
+ endif
call mpi_barrier(MPI_COMM_WORLD,ier)
! GLL points weights
@@ -235,29 +247,33 @@
enddo
enddo
-! ---------------------
+ ! user output formatting
+ if( myrank == 0 ) then
+ print *, ' '
+ endif
-! reads mesh file
-
+ ! initializes flags
ACOUSTIC_SIMULATION = .false.
ELASTIC_SIMULATION = .false.
POROELASTIC_SIMULATION = .false.
- np = 0
- iproc = myrank
- if( myrank == 0 ) then
- print *, ' '
- !print *, 'Reading partition ', iproc
- endif
+! ---------------------
- ! gets number of elements and global points for this partition
- write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'/proc',iproc,'_'//'external_mesh.bin'
- !if( myrank == 0 ) print*,trim(prname_lp)
-
+ ! reads mesh file
+ !
+ ! needs to get point locations, jacobians and MPI neighbours
+
+ ! opens external mesh file
+ write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'/proc',myrank,'_'//'external_mesh.bin'
open(unit=27,file=trim(prname_lp),&
status='old',action='read',form='unformatted',iostat=ios)
- if( ios /= 0 ) call exit_mpi(myrank, 'error reading external mesh file')
+ if( ier /= 0 ) then
+ print*,'error: could not open database '
+ print*,'path: ',trim(prname_lp)
+ call exit_mpi(myrank, 'error reading external mesh file')
+ endif
+ ! gets number of elements and global points for this partition
read(27) NSPEC_AB
read(27) NGLOB_AB
@@ -273,10 +289,11 @@
read(27) ystore
read(27) zstore
- ! reads in jacobian
allocate(dummy(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
jacobian(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array dummy and jacobian'
+
+ ! needs jacobian
read(27) dummy ! xix
read(27) dummy ! xiy
read(27) dummy ! xiz
@@ -288,6 +305,8 @@
read(27) dummy ! gammaz
read(27) jacobian
+ ! now skips all until MPI section can be read in
+
! reads in partiton neighbors
read(27) dummy ! kappastore
read(27) dummy ! mustore
@@ -304,13 +323,16 @@
call any_all_l( ANY(ldummy), POROELASTIC_SIMULATION )
deallocate(ldummy)
-
allocate(dummy_1(NGLOB_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array dummy_1'
+
+ ! acoustic
if( ACOUSTIC_SIMULATION ) then
read(27) dummy_1 ! rmass_acoustic
read(27) dummy ! rhostore
endif
+
+ ! elastic
if( ELASTIC_SIMULATION ) then
read(27) dummy_1 ! rmass
if( OCEANS ) then
@@ -319,9 +341,40 @@
read(27) dummy ! rho_vp
read(27) dummy ! rho_vs
endif
+
+ ! poroelastic
+ if( POROELASTIC_SIMULATION ) then
+ read(27) dummy_1 ! rmass_solid_poroelastic
+ read(27) dummy_1 ! rmass_fluid_poroelastic
+ allocate(dummy_5(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ read(27) dummy_5 ! rhoarraystore
+ deallocate(dummy_5)
+ allocate(dummy_5(3,NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ read(27) dummy_5 ! kappaarraystore
+ deallocate(dummy_5)
+ read(27) dummy ! etastore
+ read(27) dummy ! tortstore
+ allocate(dummy_5(6,NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ read(27) dummy_5 ! permstore
+ deallocate(dummy_5)
+ read(27) dummy ! phistore
+ read(27) dummy ! rho_vpI
+ read(27) dummy ! rho_vpII
+ read(27) dummy ! rho_vsI
+ endif
+
deallocate(dummy_1)
deallocate(dummy)
+ ! checks simulation types are valid
+ if( (.not. ACOUSTIC_SIMULATION ) .and. &
+ (.not. ELASTIC_SIMULATION ) .and. &
+ (.not. POROELASTIC_SIMULATION ) ) then
+ close(27)
+ call exit_mpi(myrank,'error no simulation type defined')
+ endif
+
+ ! absorbing boundary surface
read(27) idummy_a ! num_abs_boundary_faces
allocate(idummy(idummy_a), &
idummy_3(3,NGLLSQUARE,idummy_a), &
@@ -334,39 +387,85 @@
read(27) dummy_3 ! abs_boundary_normal
deallocate( idummy,idummy_3,dummy_2,dummy_3)
+ ! free surface
read(27) idummy_a ! num_free_surface_faces
allocate(idummy(idummy_a), &
idummy_3(3,NGLLSQUARE,idummy_a), &
dummy_2(NGLLSQUARE,idummy_a), &
dummy_3(NDIM,NGLLSQUARE,idummy_a),stat=ier)
if( ier /= 0 ) stop 'error allocating array idummy etc.'
- read(27) idummy ! free_surface_ispec
+ read(27) idummy ! free_surface_ispec
read(27) idummy_3 ! free_surface_ijk
- read(27) dummy_2 ! free_surface_jacobian2Dw
- read(27) dummy_3 ! free_surface_normal
+ read(27) dummy_2 ! free_surface_jacobian2Dw
+ read(27) dummy_3 ! free_surface_normal
deallocate( idummy,idummy_3,dummy_2,dummy_3)
+ ! acoustic-elastic coupling surface
read(27) idummy_a ! num_coupling_ac_el_faces
- allocate(idummy(idummy_a), &
- idummy_3(3,NGLLSQUARE,idummy_a), &
- dummy_2(NGLLSQUARE,idummy_a), &
- dummy_3(NDIM,NGLLSQUARE,idummy_a),stat=ier)
- if( ier /= 0 ) stop 'error allocating array idummy etc.'
- read(27) idummy ! coupling_ac_el_ispec
- read(27) idummy_3 ! coupling_ac_el_ijk
- read(27) dummy_2 ! coupling_ac_el_jacobian2Dw
- read(27) dummy_3 ! coupling_ac_el_normal
- deallocate( idummy,idummy_3,dummy_2,dummy_3)
+ if( idummy_a > 0 ) then
+ allocate(idummy(idummy_a), &
+ idummy_3(3,NGLLSQUARE,idummy_a), &
+ dummy_2(NGLLSQUARE,idummy_a), &
+ dummy_3(NDIM,NGLLSQUARE,idummy_a),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array idummy etc.'
+ read(27) idummy ! coupling_ac_el_ispec
+ read(27) idummy_3 ! coupling_ac_el_ijk
+ read(27) dummy_2 ! coupling_ac_el_jacobian2Dw
+ read(27) dummy_3 ! coupling_ac_el_normal
+ deallocate( idummy,idummy_3,dummy_2,dummy_3)
+ endif
+ ! acoustic-poroelastic coupling surface
+ read(27) idummy_a ! num_coupling_ac_po_faces
+ if( idummy_a > 0 ) then
+ allocate(idummy(idummy_a), &
+ idummy_3(3,NGLLSQUARE,idummy_a), &
+ dummy_2(NGLLSQUARE,idummy_a), &
+ dummy_3(NDIM,NGLLSQUARE,idummy_a),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array idummy etc.'
+ read(27) idummy ! coupling_ac_po_ispec
+ read(27) idummy_3 ! coupling_ac_po_ijk
+ read(27) dummy_2 ! coupling_ac_po_jacobian2Dw
+ read(27) dummy_3 ! coupling_ac_po_normal
+ deallocate( idummy,idummy_3,dummy_2,dummy_3)
+ endif
+
+ ! elastic-poroelastic coupling surface
+ read(27) idummy_a ! num_coupling_el_po_faces
+ if( idummy_a > 0 ) then
+ allocate(idummy(idummy_a), &
+ idummy_3(3,NGLLSQUARE,idummy_a), &
+ dummy_2(NGLLSQUARE,idummy_a), &
+ dummy_3(NDIM,NGLLSQUARE,idummy_a),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array idummy etc.'
+ read(27) idummy ! coupling_el_po_ispec
+ read(27) idummy ! coupling_po_el_ispec
+ read(27) idummy_3 ! coupling_el_po_ijk
+ read(27) idummy_3 ! coupling_po_el_ijk
+ read(27) dummy_2 ! coupling_el_po_jacobian2Dw
+ read(27) dummy_3 ! coupling_el_po_normal
+ deallocate( idummy,idummy_3,dummy_2,dummy_3)
+ endif
+
+ ! needs MPI neighbours
+ ! MPI interfaces
read(27) num_interfaces_ext_mesh ! num_interfaces_ext_mesh
- read(27) idummy_a ! max_nibool_interfaces_ext_mesh
- allocate(my_neighbours_ext_mesh(num_interfaces_ext_mesh),stat=ier)
- if( ier /= 0 ) stop 'error allocating array my_neighbours_ext_mesh'
- read(27) my_neighbours_ext_mesh
-
+ if( num_interfaces_ext_mesh > 0 ) then
+ read(27) idummy_a ! max_nibool_interfaces_ext_mesh
+ allocate(my_neighbours_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array my_neighbours_ext_mesh'
+ read(27) my_neighbours_ext_mesh
+ ! no more information is needed from external mesh files
+ endif
+
+ ! we're done reading in mesh arrays
close(27)
- ! get the location of the center of the elements and local points
+! ---------------------
+
+ ! for smoothing, we use cell centers to find and locate nearby elements
+ !
+ ! sets the location of the center of the elements and local points
allocate(xl(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
yl(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
zl(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
@@ -375,6 +474,7 @@
cz0(NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array xl etc.'
+ ! sets element center location
do ispec = 1, nspec_AB
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -395,6 +495,7 @@
deallocate(jacobian)
! sets up slices to process
+ node_list(:) = -1
do i=1,num_interfaces_ext_mesh
! adds neighbors
node_list(i) = my_neighbours_ext_mesh(i)
@@ -486,6 +587,7 @@
cz(NSPEC_N),stat=ier)
if( ier /= 0 ) stop 'error allocating array xx etc.'
+ ! sets element center location
do ispec = 1, nspec_N
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -527,7 +629,6 @@
! finds closest elements for smoothing
!if(myrank==0) print*, ' start looping over elements and points for smoothing ...'
-
! loop over elements to be smoothed in the current slice
do ispec = 1, nspec_AB
! --- only double loop over the elements in the search radius ---
@@ -611,7 +712,7 @@
!------------------
-! file output
+ ! file output
! smoothed kernel file name
write(ks_file,'(a,i6.6,a)') trim(outdir)//'/proc',myrank,'_'//trim(filename)//'_smooth.bin'
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -70,9 +70,6 @@
! mask ibool needed for time marching
logical,dimension(:),allocatable :: PML_mask_ibool
- ! PML damping flag
- logical:: PML = .false.
-
end module PML_par
!--------
@@ -144,9 +141,6 @@
real(kind=CUSTOM_REAL) :: dominant_wavelength,hdur_max
integer :: count,ilayer,sign,ier
- ! sets flag
- PML = .true.
-
! user output
if( myrank == 0 ) then
write(IMAIN,*)
@@ -157,6 +151,7 @@
! PML element type array: 1 = face, 2 = edge, 3 = corner
allocate(ispec_is_PML_inum(NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array ispec_is_PML_inum'
+ ispec_is_PML_inum(:) = 0
num_PML_ispec = 0
! PML interface points between PML and "regular" region
@@ -439,6 +434,8 @@
allocate(PML_ispec(num_PML_ispec), &
PML_normal(NDIM,num_PML_ispec),stat=ier)
if( ier /= 0 ) stop 'error allocating array PML_ispec'
+ PML_normal(:,:) = 0.0_CUSTOM_REAL
+ PML_ispec(:) = 0
! stores PML elements flags and normals
ispecPML = 0
@@ -888,7 +885,6 @@
iglob_pml_normal(NDIM,NGLOB_AB), &
ispec_pml_normal(NDIM,NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array is_pml_elem etc.'
-
iglob_pml_normal(:,:) = 0._CUSTOM_REAL
ispec_pml_normal(:,:) = 0._CUSTOM_REAL
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -147,7 +147,7 @@
stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
! source encoding
- stf_used = stf_used * pm1_source_encoding(isource)
+ stf_used = stf_used * pm1_source_encoding(isource)
! 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
@@ -212,12 +212,12 @@
!
! backward/reconstructed wavefields:
! time for b_potential( it ) would correspond to (NSTEP - it - 1 )*DT - t0
-! if we read in saved wavefields b_potential() before Newark time scheme
+! if we read in saved wavefields b_potential() before Newmark time scheme
! (see sources for simulation_type 1 and seismograms)
-! since at the beginning of the time loop, the numerical Newark time scheme updates
+! since at the beginning of the time loop, the numerical Newmark time scheme updates
! the wavefields, that is b_potential( it=1) would correspond to time (NSTEP -1 - 1)*DT - t0
!
-! b_potential is now read in after Newark time scheme:
+! b_potential is now read in after Newmark time scheme:
! we read the backward/reconstructed wavefield at the end of the first time loop,
! such that b_potential(it=1) corresponds to -t0 + (NSTEP-1)*DT.
! assuming that until that end the backward/reconstructed wavefield and adjoint fields
@@ -348,7 +348,7 @@
endif ! it
endif
-! note: b_potential() is read in after Newark time scheme, thus
+! note: b_potential() is read in after Newmark time scheme, thus
! b_potential(it=1) corresponds to -t0 + (NSTEP-1)*DT.
! thus indexing is NSTEP - it , instead of NSTEP - it - 1
@@ -390,7 +390,7 @@
stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
! source encoding
- stf_used = stf_used * pm1_source_encoding(isource)
+ stf_used = stf_used * pm1_source_encoding(isource)
! 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
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -190,12 +190,12 @@
!
! backward/reconstructed wavefields:
! time for b_displ( it ) would correspond to (NSTEP - it - 1 )*DT - t0
-! if we read in saved wavefields b_displ() before Newark time scheme
+! if we read in saved wavefields b_displ() before Newmark time scheme
! (see sources for simulation_type 1 and seismograms)
-! since at the beginning of the time loop, the numerical Newark time scheme updates
+! since at the beginning of the time loop, the numerical Newmark time scheme updates
! the wavefields, that is b_displ( it=1) would correspond to time (NSTEP -1 - 1)*DT - t0
!
-! b_displ is now read in after Newark time scheme:
+! b_displ is now read in after Newmark time scheme:
! we read the backward/reconstructed wavefield at the end of the first time loop,
! such that b_displ(it=1) corresponds to -t0 + (NSTEP-1)*DT.
! assuming that until that end the backward/reconstructed wavefield and adjoint fields
@@ -339,7 +339,7 @@
endif !adjoint
-! note: b_displ() is read in after Newark time scheme, thus
+! note: b_displ() is read in after Newmark time scheme, thus
! b_displ(it=1) corresponds to -t0 + (NSTEP-1)*DT.
! thus indexing is NSTEP - it , instead of NSTEP - it - 1
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -212,7 +212,7 @@
! since we start with saved wavefields b_displ( 0 ) = displ( NSTEP ) which correspond
! to a time (NSTEP - 1)*DT - t0
! (see sources for simulation_type 1 and seismograms)
-! now, at the beginning of the time loop, the numerical Newark time scheme updates
+! now, at the beginning of the time loop, the numerical Newmark time scheme updates
! the wavefields, that is b_displ( it=1) corresponds now to time (NSTEP -1 - 1)*DT - t0
!
! let's define the start time t to (1-1)*DT - t0 = -t0, and the end time T to (NSTEP-1)*DT - t0
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_arrays_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_arrays_source.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_arrays_source.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -319,7 +319,7 @@
!!! ! reads in adjoint source trace
!!! do itime = 1, NSTEP
!!!
-!!! ! things become a bit tricky because of the Newark time scheme at
+!!! ! things become a bit tricky because of the Newmark time scheme at
!!! ! the very beginning of the time loop. however, when we read in the backward/reconstructed
!!! ! wavefields at the end of the first time loop, we can use the adjoint source index from 1 to NSTEP
!!! ! (and then access it in reverse NSTEP-it+1 down to 1, for it=1,..NSTEP; see compute_add_sources*.f90).
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_el.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -106,7 +106,7 @@
! continuity of pressure and normal displacement on global point
!
- ! note: newark time scheme together with definition of scalar potential:
+ ! note: Newmark time scheme together with definition of scalar potential:
! pressure = - chi_dot_dot
! requires that this coupling term uses the updated displacement at time step [t+delta_t],
! which is done at the very beginning of the time loop
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_po.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_po.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_po.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -104,7 +104,7 @@
! continuity of pressure and normal displacement on global point
!
- ! note: newark time scheme together with definition of scalar potential:
+ ! note: Newmark time scheme together with definition of scalar potential:
! pressure = - chi_dot_dot
! requires that this coupling term uses the updated displacement at time step [t+delta_t],
! which is done at the very beginning of the time loop
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_ac.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_ac.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_ac.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -101,7 +101,7 @@
! continuity of displacement and pressure on global point
!
- ! note: newark time scheme together with definition of scalar potential:
+ ! note: Newmark time scheme together with definition of scalar potential:
! pressure = - chi_dot_dot
! requires that this coupling term uses the *UPDATED* pressure (chi_dot_dot), i.e.
! pressure at time step [t + delta_t]
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_po.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_po.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_elastic_po.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -52,11 +52,11 @@
integer :: NSPEC_AB,NGLOB_AB
-! displacements, etc
+! displacements, etc
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displs_poroelastic,&
displw_poroelastic
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
-
+
! global indexing
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -67,12 +67,12 @@
! elastic-poroelastic coupling surface
integer :: num_coupling_el_po_faces
- real(kind=CUSTOM_REAL) :: coupling_el_po_normal(NDIM,NGLLSQUARE,num_coupling_el_po_faces)
+ real(kind=CUSTOM_REAL) :: coupling_el_po_normal(NDIM,NGLLSQUARE,num_coupling_el_po_faces)
real(kind=CUSTOM_REAL) :: coupling_el_po_jacobian2Dw(NGLLSQUARE,num_coupling_el_po_faces)
integer :: coupling_el_po_ijk(3,NGLLSQUARE,num_coupling_el_po_faces)
integer :: coupling_po_el_ijk(3,NGLLSQUARE,num_coupling_el_po_faces)
- integer :: coupling_el_po_ispec(num_coupling_el_po_faces)
- integer :: coupling_po_el_ispec(num_coupling_el_po_faces)
+ integer :: coupling_el_po_ispec(num_coupling_el_po_faces)
+ integer :: coupling_po_el_ispec(num_coupling_el_po_faces)
! properties
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -101,7 +101,7 @@
logical :: phase_is_inner
! local parameters
- real(kind=CUSTOM_REAL) :: sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+ real(kind=CUSTOM_REAL) :: sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
real(kind=CUSTOM_REAL) :: rhol_s,rhol_f,phil,tortl,rhol_bar
real(kind=CUSTOM_REAL) :: nx,ny,nz,jacobianw
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
@@ -115,7 +115,7 @@
! local anisotropy parameters
real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
-
+
integer :: iface,igll,ispec_po,ispec_el,iglob,iglob_el,iglob_po
integer :: i,j,k,l
@@ -137,7 +137,7 @@
! Jacobian matrix and determinant
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
-
+
! loops on all coupling faces
do iface = 1,num_coupling_el_po_faces
@@ -151,7 +151,7 @@
do igll = 1, NGLLSQUARE
!-----------------------
- ! from the poroelastic side
+ ! from the poroelastic side
!-----------------------
i = coupling_el_po_ijk(1,igll,iface)
j = coupling_el_po_ijk(2,igll,iface)
@@ -329,7 +329,7 @@
!-----------------------
- ! from the elastic side
+ ! from the elastic side
!-----------------------
i = coupling_po_el_ijk(1,igll,iface)
j = coupling_po_el_ijk(2,igll,iface)
@@ -469,7 +469,7 @@
ny = coupling_el_po_normal(2,igll,iface)
nz = coupling_el_po_normal(3,igll,iface)
- ! gets associated, weighted 2D jacobian
+ ! gets associated, weighted 2D jacobian
! (note: should be the same for poroelastic and elastic element)
jacobianw = coupling_el_po_jacobian2Dw(igll,iface)
@@ -484,11 +484,11 @@
accel(3,iglob_el) = accel(3,iglob_el) - jacobianw* &
( sigma_xz*nx + sigma_yz*ny + sigma_zz*nz )/2.d0
-
+
enddo ! igll
endif
-
+
enddo ! iface
end subroutine compute_coupling_elastic_po
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_ac.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_ac.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_ac.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -114,7 +114,7 @@
! continuity of displacement and pressure on global point
!
- ! note: newark time scheme together with definition of scalar potential:
+ ! note: Newmark time scheme together with definition of scalar potential:
! pressure = - chi_dot_dot
! requires that this coupling term uses the *UPDATED* pressure (chi_dot_dot), i.e.
! pressure at time step [t + delta_t]
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_el.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_el.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_el.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -52,11 +52,11 @@
integer :: NSPEC_AB,NGLOB_AB
-! displacements, etc
+! displacements, etc
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displs_poroelastic,accels_poroelastic,&
displw_poroelastic
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ
-
+
! global indexing
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -67,12 +67,12 @@
! elastic-poroelastic coupling surface
integer :: num_coupling_el_po_faces
- real(kind=CUSTOM_REAL) :: coupling_el_po_normal(NDIM,NGLLSQUARE,num_coupling_el_po_faces)
+ real(kind=CUSTOM_REAL) :: coupling_el_po_normal(NDIM,NGLLSQUARE,num_coupling_el_po_faces)
real(kind=CUSTOM_REAL) :: coupling_el_po_jacobian2Dw(NGLLSQUARE,num_coupling_el_po_faces)
integer :: coupling_el_po_ijk(3,NGLLSQUARE,num_coupling_el_po_faces)
integer :: coupling_po_el_ijk(3,NGLLSQUARE,num_coupling_el_po_faces)
- integer :: coupling_el_po_ispec(num_coupling_el_po_faces)
- integer :: coupling_po_el_ispec(num_coupling_el_po_faces)
+ integer :: coupling_el_po_ispec(num_coupling_el_po_faces)
+ integer :: coupling_po_el_ispec(num_coupling_el_po_faces)
! properties
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -101,7 +101,7 @@
logical :: phase_is_inner
! local parameters
- real(kind=CUSTOM_REAL) :: sigmap,sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+ real(kind=CUSTOM_REAL) :: sigmap,sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
real(kind=CUSTOM_REAL) :: rhol_s,rhol_f,phil,tortl,rhol_bar
real(kind=CUSTOM_REAL) :: nx,ny,nz,jacobianw
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
@@ -115,7 +115,7 @@
! local anisotropy parameters
real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
-
+
integer :: iface,igll,ispec_po,ispec_el,iglob,iglob_po,iglob_el
integer :: i,j,k,l
@@ -137,7 +137,7 @@
! Jacobian matrix and determinant
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
-
+
! loops on all coupling faces
do iface = 1,num_coupling_el_po_faces
@@ -151,12 +151,12 @@
do igll = 1, NGLLSQUARE
!-----------------------
- ! from the elastic side
+ ! from the elastic side
!-----------------------
i = coupling_po_el_ijk(1,igll,iface)
j = coupling_po_el_ijk(2,igll,iface)
k = coupling_po_el_ijk(3,igll,iface)
-
+
! gets global index of this common GLL point
! (note: should be the same as for corresponding i',j',k',ispec_poroelastic or ispec_elastic )
iglob_el = ibool(i,j,k,ispec_el)
@@ -286,7 +286,7 @@
!-----------------------
- ! from the poroelastic side
+ ! from the poroelastic side
!-----------------------
i = coupling_el_po_ijk(1,igll,iface)
j = coupling_el_po_ijk(2,igll,iface)
@@ -473,12 +473,12 @@
! (note convention: pointing outwards of poroelastic element)
nx = coupling_el_po_normal(1,igll,iface)
ny = coupling_el_po_normal(2,igll,iface)
- nz = coupling_el_po_normal(3,igll,iface)
-
- ! gets associated, weighted 2D jacobian
+ nz = coupling_el_po_normal(3,igll,iface)
+
+ ! gets associated, weighted 2D jacobian
! (note: should be the same for poroelastic and elastic element)
jacobianw = coupling_el_po_jacobian2Dw(igll,iface)
-
+
! continuity of displacement and traction on global point
!
! note: continuity of displacement is enforced after the velocity update
@@ -493,11 +493,11 @@
( (sigma_xz*nx + sigma_yz*ny + sigma_zz*nz)/2.d0 - phil/tortl*sigmap*nz )
! contribution to the fluid phase
! w = 0
-
+
enddo ! igll
endif
-
+
enddo ! iface
end subroutine compute_coupling_poroelastic_el
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -79,7 +79,8 @@
num_free_surface_faces,ispec_is_acoustic)
- if(PML) call PML_acoustic_enforce_free_srfc(NSPEC_AB,NGLOB_AB, &
+ if(ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
+ call PML_acoustic_enforce_free_srfc(NSPEC_AB,NGLOB_AB, &
potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, &
ibool,free_surface_ijk,free_surface_ispec, &
num_free_surface_faces,ispec_is_acoustic, &
@@ -88,6 +89,7 @@
chi1_dot,chi2_t_dot,chi3_dot,chi4_dot,&
chi1_dot_dot,chi2_t_dot_dot,&
chi3_dot_dot,chi4_dot_dot)
+ endif
! distinguishes two runs: for points on MPI interfaces, and points within the partitions
do iphase=1,2
@@ -122,7 +124,7 @@
num_phase_ispec_acoustic,nspec_inner_acoustic,nspec_outer_acoustic,&
phase_ispec_inner_acoustic )
- if(PML) then
+ if(ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
call compute_forces_acoustic_PML(NSPEC_AB,NGLOB_AB, &
ibool,ispec_is_inner,phase_is_inner, &
rhostore,ispec_is_acoustic,potential_acoustic, &
@@ -147,9 +149,10 @@
! absorbing boundaries
if(ABSORBING_CONDITIONS) then
- if( PML .and. PML_USE_SOMMERFELD ) then
- ! adds a Sommerfeld condition on the domain's absorbing boundaries
- call PML_acoustic_abs_boundaries(phase_is_inner,NSPEC_AB,NGLOB_AB,&
+ if(ABSORB_USE_PML) then
+ if( PML_USE_SOMMERFELD ) then
+ ! adds a Sommerfeld condition on the domain's absorbing boundaries
+ call PML_acoustic_abs_boundaries(phase_is_inner,NSPEC_AB,NGLOB_AB,&
abs_boundary_jacobian2Dw,abs_boundary_ijk,abs_boundary_ispec, &
num_abs_boundary_faces, &
kappastore,ibool,ispec_is_inner, &
@@ -158,6 +161,7 @@
num_PML_ispec,PML_ispec,ispec_is_PML_inum,&
chi1_dot,chi2_t,chi2_t_dot,chi3_dot,chi4_dot,&
chi1_dot_dot,chi3_dot_dot,chi4_dot_dot)
+ endif
else
call compute_stacey_acoustic(NSPEC_AB,NGLOB_AB, &
potential_dot_dot_acoustic,potential_dot_acoustic, &
@@ -218,10 +222,10 @@
coupling_ac_po_jacobian2Dw, &
ispec_is_inner,phase_is_inner)
else
- stop 'not implemented yet'
+ stop 'not implemented yet'
endif
if( SIMULATION_TYPE == 3 ) &
- stop 'not implemented yet'
+ stop 'not implemented yet'
endif
endif
@@ -282,14 +286,14 @@
b_potential_dot_dot_acoustic(:) = b_potential_dot_dot_acoustic(:) * rmass_acoustic(:)
- if(PML) then
+ if(ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
! divides local contributions with mass term
call PML_acoustic_mass_update(NSPEC_AB,NGLOB_AB,&
ispec_is_acoustic,rmass_acoustic,ibool,&
num_PML_ispec,PML_ispec,&
chi1_dot_dot,chi2_t_dot_dot,chi3_dot_dot,chi4_dot_dot)
- ! Newark time scheme corrector terms
+ ! Newmark time scheme corrector terms
call PML_acoustic_time_corrector(NSPEC_AB,ispec_is_acoustic,deltatover2,&
num_PML_ispec,PML_ispec,PML_damping_d,&
chi1_dot,chi2_t_dot,chi3_dot,chi4_dot,&
@@ -298,7 +302,7 @@
! update velocity
-! note: Newark finite-difference time scheme with acoustic domains:
+! note: Newmark finite-difference time scheme with acoustic domains:
! (see e.g. Hughes, 1987; Chaljub et al., 2003)
!
! chi(t+delta_t) = chi(t) + delta_t chi_dot(t) + 1/2 delta_t**2 chi_dot_dot(t)
@@ -320,7 +324,8 @@
b_potential_dot_acoustic(:) = b_potential_dot_acoustic(:) + b_deltatover2*b_potential_dot_dot_acoustic(:)
! updates potential_dot_acoustic and potential_dot_dot_acoustic inside PML region for plotting seismograms/movies
- if(PML) call PML_acoustic_update_potentials(NGLOB_AB,NSPEC_AB, &
+ if(ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
+ call PML_acoustic_update_potentials(NGLOB_AB,NSPEC_AB, &
ibool,ispec_is_acoustic, &
potential_dot_acoustic,potential_dot_dot_acoustic,&
num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
@@ -331,6 +336,7 @@
chi1,chi2,chi2_t,chi3,&
chi1_dot,chi2_t_dot,chi3_dot,chi4_dot,&
chi1_dot_dot,chi3_dot_dot,chi4_dot_dot)
+ endif
! enforces free surface (zeroes potentials at free surface)
call acoustic_enforce_free_surface(NSPEC_AB,NGLOB_AB, &
@@ -343,7 +349,7 @@
+ deltat * potential_dot_acoustic(:) &
+ deltatsqover2 * potential_dot_dot_acoustic(:)
endif
-
+
! adjoint simulations
if (SIMULATION_TYPE == 3) &
call acoustic_enforce_free_surface(NSPEC_AB,NGLOB_ADJOINT, &
@@ -352,7 +358,8 @@
num_free_surface_faces,ispec_is_acoustic)
- if(PML) call PML_acoustic_enforce_free_srfc(NSPEC_AB,NGLOB_AB, &
+ if(ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
+ call PML_acoustic_enforce_free_srfc(NSPEC_AB,NGLOB_AB, &
potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, &
ibool,free_surface_ijk,free_surface_ispec, &
num_free_surface_faces, &
@@ -362,6 +369,7 @@
chi1_dot,chi2_t_dot,chi3_dot,chi4_dot,&
chi1_dot_dot,chi2_t_dot_dot,&
chi3_dot_dot,chi4_dot_dot)
+ endif
end subroutine compute_forces_acoustic
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_pot.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_pot.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_pot.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -41,24 +41,25 @@
! note that pressure is defined as:
! p = - Chi_dot_dot
!
- use constants,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,TINYVAL_SNGL
- use PML_par,only:PML,ispec_is_PML_inum
+ use specfem_par,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,TINYVAL_SNGL,ABSORB_USE_PML,ABSORBING_CONDITIONS
+ use PML_par,only:ispec_is_PML_inum
+
implicit none
!include "constants.h"
integer :: NSPEC_AB,NGLOB_AB
-! acoustic potentials
+ ! acoustic potentials
real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: &
potential_acoustic,potential_dot_dot_acoustic
-! arrays with mesh parameters per slice
+ ! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
rhostore,jacobian
-! array with derivatives of Lagrange polynomials and precalculated products
+ ! array with derivatives of Lagrange polynomials and precalculated products
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
@@ -94,106 +95,107 @@
num_elements = nspec_inner_acoustic
endif
-! loop over spectral elements
+ ! loop over spectral elements
do ispec_p = 1,num_elements
!if ( (ispec_is_inner(ispec) .eqv. phase_is_inner) ) then
- ispec = phase_ispec_inner_acoustic(ispec_p,iphase)
+ ispec = phase_ispec_inner_acoustic(ispec_p,iphase)
- ! only elements outside PML, inside "regular" domain
- if( PML ) then
- if( ispec_is_PML_inum(ispec) > 0 ) then
- cycle
- endif
- endif
+ ! only elements outside PML, inside "regular" domain
+ if( ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
+ if( ispec_is_PML_inum(ispec) > 0 ) cycle
+ endif
! if( ispec_is_acoustic(ispec) ) then
- ! gets values for element
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- chi_elem(i,j,k) = potential_acoustic(ibool(i,j,k,ispec))
- enddo
- enddo
+ ! gets values for element
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ chi_elem(i,j,k) = potential_acoustic(ibool(i,j,k,ispec))
enddo
- ! would check if anything to do, but might lower accuracy of computation
- !if( maxval( abs( chi_elem ) ) < TINYVAL_SNGL ) cycle
+ enddo
+ enddo
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ ! would check if anything to do, but might lower accuracy of computation
+ !if( maxval( abs( chi_elem ) ) < TINYVAL_SNGL ) cycle
- ! density (reciproc)
- rho_invl = 1.0_CUSTOM_REAL / rhostore(i,j,k,ispec)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
- ! derivative along x, y, z
- ! first double loop over GLL points to compute and store gradients
- ! we can merge the loops because NGLLX == NGLLY == NGLLZ
- temp1l = 0._CUSTOM_REAL
- temp2l = 0._CUSTOM_REAL
- temp3l = 0._CUSTOM_REAL
- do l = 1,NGLLX
- temp1l = temp1l + chi_elem(l,j,k)*hprime_xx(i,l)
- temp2l = temp2l + chi_elem(i,l,k)*hprime_yy(j,l)
- temp3l = temp3l + chi_elem(i,j,l)*hprime_zz(k,l)
- enddo
+ ! derivative along x, y, z
+ ! first double loop over GLL points to compute and store gradients
+ ! we can merge the loops because NGLLX == NGLLY == NGLLZ
+ temp1l = 0._CUSTOM_REAL
+ temp2l = 0._CUSTOM_REAL
+ temp3l = 0._CUSTOM_REAL
+ do l = 1,NGLLX
+ temp1l = temp1l + chi_elem(l,j,k)*hprime_xx(i,l)
+ temp2l = temp2l + chi_elem(i,l,k)*hprime_yy(j,l)
+ temp3l = temp3l + chi_elem(i,j,l)*hprime_zz(k,l)
+ enddo
- ! get derivatives of potential with respect to x, y and z
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
- jacobianl = jacobian(i,j,k,ispec)
+ ! get derivatives of potential with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
- ! derivatives of potential
- dpotentialdxl = xixl*temp1l + etaxl*temp2l + gammaxl*temp3l
- dpotentialdyl = xiyl*temp1l + etayl*temp2l + gammayl*temp3l
- dpotentialdzl = xizl*temp1l + etazl*temp2l + gammazl*temp3l
+ ! derivatives of potential
+ dpotentialdxl = xixl*temp1l + etaxl*temp2l + gammaxl*temp3l
+ dpotentialdyl = xiyl*temp1l + etayl*temp2l + gammayl*temp3l
+ dpotentialdzl = xizl*temp1l + etazl*temp2l + gammazl*temp3l
- ! for acoustic medium
- ! also add GLL integration weights
- temp1(i,j,k) = rho_invl * wgllwgll_yz(j,k) * jacobianl* &
- (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- temp2(i,j,k) = rho_invl * wgllwgll_xz(i,k) * jacobianl* &
- (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- temp3(i,j,k) = rho_invl * wgllwgll_xy(i,j) * jacobianl* &
- (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
- enddo
- enddo
+ jacobianl = jacobian(i,j,k,ispec)
+
+ ! density (reciproc)
+ rho_invl = 1.0_CUSTOM_REAL / rhostore(i,j,k,ispec)
+
+ ! for acoustic medium
+ ! also add GLL integration weights
+ temp1(i,j,k) = rho_invl * wgllwgll_yz(j,k) * jacobianl* &
+ (xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
+ temp2(i,j,k) = rho_invl * wgllwgll_xz(i,k) * jacobianl* &
+ (etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
+ temp3(i,j,k) = rho_invl * wgllwgll_xy(i,j) * jacobianl* &
+ (gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
+
enddo
+ enddo
+ enddo
- ! second double-loop over GLL to compute all the terms
- do k = 1,NGLLZ
- do j = 1,NGLLZ
- do i = 1,NGLLX
+ ! second double-loop over GLL to compute all the terms
+ do k = 1,NGLLZ
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
- ! along x,y,z direction
- ! and assemble the contributions
- !!! can merge these loops because NGLLX = NGLLY = NGLLZ
- temp1l = 0._CUSTOM_REAL
- temp2l = 0._CUSTOM_REAL
- temp3l = 0._CUSTOM_REAL
- do l=1,NGLLX
- temp1l = temp1l + temp1(l,j,k) * hprimewgll_xx(l,i)
- temp2l = temp2l + temp2(i,l,k) * hprimewgll_yy(l,j)
- temp3l = temp3l + temp3(i,j,l) * hprimewgll_zz(l,k)
- enddo
+ ! along x,y,z direction
+ ! and assemble the contributions
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ temp1l = 0._CUSTOM_REAL
+ temp2l = 0._CUSTOM_REAL
+ temp3l = 0._CUSTOM_REAL
+ do l=1,NGLLX
+ temp1l = temp1l + temp1(l,j,k) * hprimewgll_xx(l,i)
+ temp2l = temp2l + temp2(i,l,k) * hprimewgll_yy(l,j)
+ temp3l = temp3l + temp3(i,j,l) * hprimewgll_zz(l,k)
+ enddo
- ! sum contributions from each element to the global values
- iglob = ibool(i,j,k,ispec)
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - ( temp1l + temp2l + temp3l )
+ ! sum contributions from each element to the global values
+ iglob = ibool(i,j,k,ispec)
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - ( temp1l + temp2l + temp3l )
- enddo
- enddo
enddo
+ enddo
+ enddo
! endif ! end of test if acoustic element
! endif ! ispec_is_inner
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -142,7 +142,7 @@
coupling_ac_el_jacobian2Dw, &
ispec_is_inner,phase_is_inner)
else
- ! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
+ ! handles adjoint runs coupling between adjoint potential and adjoint elastic wavefield
! adoint definition: pressure^\dagger=potential^\dagger
call compute_coupling_elastic_ac(NSPEC_AB,NGLOB_AB, &
ibool,accel,-potential_acoustic_adj_coupling, &
@@ -167,7 +167,7 @@
! poroelastic coupling
- !print *,'entering poro counpling'
+ !print *,'entering poro counpling'
if( POROELASTIC_SIMULATION ) &
call compute_coupling_elastic_po(NSPEC_AB,NGLOB_AB,ibool,&
displs_poroelastic,displw_poroelastic,&
@@ -188,7 +188,7 @@
coupling_el_po_normal, &
coupling_el_po_jacobian2Dw, &
ispec_is_inner,phase_is_inner)
- !print *,'ok poro counpling'
+ !print *,'ok poro counpling'
! adds source term (single-force/moment-tensor solution)
call compute_add_sources_elastic( NSPEC_AB,NGLOB_AB,accel, &
@@ -271,7 +271,7 @@
endif
! updates velocities
-! Newark finite-difference time scheme with elastic domains:
+! Newmark finite-difference time scheme with elastic domains:
! (see e.g. Hughes, 1987; Chaljub et al., 2003)
!
! u(t+delta_t) = u(t) + delta_t v(t) + 1/2 delta_t**2 a(t)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -133,7 +133,7 @@
coupling_el_po_normal, &
coupling_el_po_jacobian2Dw, &
ispec_is_inner,phase_is_inner)
- !print *,'ok elastic counpling'
+ !print *,'ok elastic counpling'
! adjoint simulations
!chris: TO DO
@@ -244,7 +244,7 @@
endif !adjoint
! updates velocities
-! Newark finite-difference time scheme with elastic domains:
+! Newmark finite-difference time scheme with elastic domains:
! (see e.g. Hughes, 1987; Chaljub et al., 2003)
!
! u(t+delta_t) = u(t) + delta_t v(t) + 1/2 delta_t**2 a(t)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_gradient.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_gradient.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_gradient.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -25,7 +25,7 @@
!=====================================================================
-subroutine compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
+ subroutine compute_gradient(ispec,NSPEC_AB,NGLOB_AB, &
scalar_field, vector_field_element,&
hprime_xx,hprime_yy,hprime_zz, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
@@ -111,6 +111,6 @@
enddo
enddo
-end subroutine compute_gradient
+ end subroutine compute_gradient
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_acoustic.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -72,7 +72,7 @@
integer :: ispec,iglob,i,j,k,iface,igll
!integer:: reclen1,reclen2
-! adjoint simulations:
+ ! adjoint simulations:
if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
! reads in absorbing boundary array when first phase is running
if( phase_is_inner .eqv. .false. ) then
@@ -86,7 +86,7 @@
endif
endif !adjoint
-! absorbs absorbing-boundary surface using Sommerfeld condition (vanishing field in the outer-space)
+ ! absorbs absorbing-boundary surface using Sommerfeld condition (vanishing field in the outer-space)
do iface=1,num_abs_boundary_faces
ispec = abs_boundary_ispec(iface)
@@ -118,7 +118,6 @@
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- absorbl
-
! adjoint simulations
if (SIMULATION_TYPE == 3) then
b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_elastic.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -84,7 +84,7 @@
if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
! reads in absorbing boundary array when first phase is running
if( phase_is_inner .eqv. .false. ) then
- ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newark scheme
+ ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newmark scheme
! uses fortran routine
!read(IOABS,rec=NSTEP-it+1) reclen1,b_absorb_field,reclen2
!if (reclen1 /= b_reclen_field .or. reclen1 /= reclen2) &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -37,28 +37,28 @@
! type = 2 : velocity V_y component
! type = 3 : velocity V_z component
! type = 4 : velocity V norm
- integer,parameter:: IMAGE_TYPE = 2
+ integer,parameter:: IMAGE_TYPE = 3
! cross-section surface
! cross-section origin point
real(kind=CUSTOM_REAL),parameter:: section_xorg = 67000.0
real(kind=CUSTOM_REAL),parameter:: section_yorg = 0.0
- real(kind=CUSTOM_REAL),parameter:: section_zorg = 0.0
+ real(kind=CUSTOM_REAL),parameter:: section_zorg = -1000.0
! cross-section surface normal
- real(kind=CUSTOM_REAL),parameter:: section_nx = 1.0
+ real(kind=CUSTOM_REAL),parameter:: section_nx = 0.0
real(kind=CUSTOM_REAL),parameter:: section_ny = 0.0
- real(kind=CUSTOM_REAL),parameter:: section_nz = 0.0
+ real(kind=CUSTOM_REAL),parameter:: section_nz = 1.0
! cross-section (in-plane) horizontal-direction
- real(kind=CUSTOM_REAL),parameter:: section_hdirx = 0.0
- real(kind=CUSTOM_REAL),parameter:: section_hdiry = 1.0
+ real(kind=CUSTOM_REAL),parameter:: section_hdirx = 1.0
+ real(kind=CUSTOM_REAL),parameter:: section_hdiry = 0.0
real(kind=CUSTOM_REAL),parameter:: section_hdirz = 0.0
! cross-section (in-plane) vertical-direction
real(kind=CUSTOM_REAL),parameter:: section_vdirx = 0.0
- real(kind=CUSTOM_REAL),parameter:: section_vdiry = 0.0
- real(kind=CUSTOM_REAL),parameter:: section_vdirz = 1.0
+ real(kind=CUSTOM_REAL),parameter:: section_vdiry = 1.0
+ real(kind=CUSTOM_REAL),parameter:: section_vdirz = 0.0
! non linear display to enhance small amplitudes in color images
real(kind=CUSTOM_REAL), parameter :: POWER_DISPLAY_COLOR = 0.30_CUSTOM_REAL
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -66,11 +66,11 @@
do it = 1,NSTEP
! simulation status output and stability check
- if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5) then
+ if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5 .or. it == NSTEP) then
call it_check_stability()
endif
- ! update displacement using Newark time scheme
+ ! update displacement using Newmark time scheme
call it_update_displacement_scheme()
! acoustic solver
@@ -85,7 +85,7 @@
if( POROELASTIC_SIMULATION ) call compute_forces_poroelastic()
! restores last time snapshot saved for backward/reconstruction of wavefields
- ! note: this must be read in after the Newark time scheme
+ ! note: this must be read in after the Newmark time scheme
if( SIMULATION_TYPE == 3 .and. it == 1 ) then
call it_read_foward_arrays()
endif
@@ -282,7 +282,7 @@
subroutine it_update_displacement_scheme()
-! explicit Newark time scheme with acoustic & elastic domains:
+! explicit Newmark time scheme with acoustic & elastic domains:
! (see e.g. Hughes, 1987; Chaljub et al., 2003)
!
! chi(t+delta_t) = chi(t) + delta_t chi_dot(t) + 1/2 delta_t**2 chi_dot_dot(t)
@@ -329,7 +329,8 @@
potential_dot_dot_acoustic(:) = 0._CUSTOM_REAL
! time marching potentials
- if(PML) call PML_acoustic_time_march(NSPEC_AB,NGLOB_AB,ibool,&
+ if(ABSORB_USE_PML .and. ABSORBING_CONDITIONS) then
+ call PML_acoustic_time_march(NSPEC_AB,NGLOB_AB,ibool,&
potential_acoustic,potential_dot_acoustic,&
deltat,deltatsqover2,deltatover2,&
num_PML_ispec,PML_ispec,PML_damping_d,&
@@ -341,6 +342,7 @@
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
my_neighbours_ext_mesh,NPROC,&
ispec_is_acoustic)
+ endif
endif
! updates elastic displacement and velocity
@@ -470,7 +472,7 @@
! note backward/reconstructed wavefields:
! storing wavefield displ() at time step it, corresponds to time (it-1)*DT - t0 (see routine write_seismograms_to_file )
! reconstucted wavefield b_displ() at it corresponds to time (NSTEP-it-1)*DT - t0
-! we read in the reconstructed wavefield at the end of the time iteration loop, i.e. after the Newark scheme,
+! we read in the reconstructed wavefield at the end of the time iteration loop, i.e. after the Newmark scheme,
! thus, indexing is NSTEP-it (rather than something like NSTEP-(it-1) )
if (SIMULATION_TYPE == 3 .and. mod(NSTEP-it,NSTEP_Q_SAVE) == 0) then
! reads files content
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -76,11 +76,11 @@
integer, dimension(3,NGLLSQUARE,num_free_surface_faces) :: free_surface_ijk
! local parameters
-
+
! for surface locating and normal computing with external mesh
real(kind=CUSTOM_REAL), dimension(3) :: u_vector,v_vector,w_vector
integer :: pt0_ix,pt0_iy,pt0_iz,pt1_ix,pt1_iy,pt1_iz,pt2_ix,pt2_iy,pt2_iz
-
+
integer, allocatable, dimension(:) :: ix_initial_guess,iy_initial_guess,iz_initial_guess
integer iprocloop
integer ios
@@ -145,7 +145,7 @@
! get MPI starting time
time_start = wtime()
-
+
! user output
if(myrank == 0) then
write(IMAIN,*)
@@ -161,7 +161,7 @@
xmin=minval(xstore(:)); xmax=maxval(xstore(:))
ymin=minval(ystore(:)); ymax=maxval(ystore(:))
zmin=minval(zstore(:)); zmax=maxval(zstore(:))
-
+
! loop limits
if(FASTER_RECEIVERS_POINTS_ONLY) then
imin_temp = 1; imax_temp = NGLLX
@@ -192,13 +192,13 @@
do irec=1,nrec
read(1,*,iostat=ios) station_name(irec),network_name(irec),llat,llon,lele,lbur
if (ios /= 0) call exit_mpi(myrank, 'Error reading station file '//trim(rec_filename))
- enddo
- close(1)
+ enddo
+ close(1)
! master reads in available station information
if( myrank == 0 ) then
open(unit=IOUT_SU,file=trim(OUTPUT_FILES)//'/SU_stations_info.bin', &
- status='old',action='read',form='unformatted',iostat=ios)
- if (ios /= 0) call exit_mpi(myrank,'error opening file '//trim(rec_filename))
+ status='old',action='read',form='unformatted',iostat=ios)
+ if (ios /= 0) call exit_mpi(myrank,'error opening file '//trim(rec_filename))
write(IMAIN,*) 'station details from SU_stations_info.bin'
allocate(x_found(nrec),y_found(nrec),z_found(nrec))
! reads in station infos
@@ -206,7 +206,7 @@
read(IOUT_SU) xi_receiver,eta_receiver,gamma_receiver
read(IOUT_SU) x_found,y_found,z_found
read(IOUT_SU) nu
- close(IOUT_SU)
+ close(IOUT_SU)
! write the locations of stations, so that we can load them and write them to SU headers later
open(unit=IOUT_SU,file=trim(OUTPUT_FILES)//'/output_list_stations.txt', &
status='unknown',action='write',iostat=ios)
@@ -214,7 +214,7 @@
do irec=1,nrec
write(IOUT_SU,*) x_found(irec),y_found(irec),z_found(irec)
enddo
- close(IOUT_SU)
+ close(IOUT_SU)
deallocate(x_found,y_found,z_found)
endif
! main process broadcasts the results to all the slices
@@ -223,7 +223,7 @@
call bcast_all_dp(xi_receiver,nrec)
call bcast_all_dp(eta_receiver,nrec)
call bcast_all_dp(gamma_receiver,nrec)
- call bcast_all_dp(nu,NDIM*NDIM*nrec)
+ call bcast_all_dp(nu,NDIM*NDIM*nrec)
call sync_all()
! user output
if( myrank == 0 ) then
@@ -233,9 +233,9 @@
write(IMAIN,*) 'Elapsed time for receiver detection in seconds = ',tCPU
write(IMAIN,*)
write(IMAIN,*) 'End of receiver detection - done'
- write(IMAIN,*)
+ write(IMAIN,*)
endif
- ! everything done
+ ! everything done
return
endif
endif
@@ -292,7 +292,7 @@
' horizontal distance: ',sngl(horiz_dist(irec)),' km'
endif
endif
-
+
! get approximate topography elevation at source long/lat coordinates
! set distance to huge initial value
distmin = HUGEVAL
@@ -872,10 +872,10 @@
write(IMAIN,*) 'error locating station # ',irec,' ',station_name(irec),network_name(irec)
call exit_MPI(myrank,'error locating receiver')
endif
-
+
! limits user output if too many receivers
if( nrec < 1000 .and. (.not. SU_FORMAT ) ) then
-
+
write(IMAIN,*)
write(IMAIN,*) 'station # ',irec,' ',station_name(irec),network_name(irec)
@@ -974,7 +974,7 @@
! stores station infos for later runs
if( SU_FORMAT ) then
open(unit=IOUT_SU,file=trim(OUTPUT_FILES)//'/SU_stations_info.bin', &
- status='unknown',action='write',form='unformatted',iostat=ios)
+ status='unknown',action='write',form='unformatted',iostat=ios)
if( ios == 0 ) then
write(IOUT_SU) islice_selected_rec,ispec_selected_rec
write(IOUT_SU) xi_receiver,eta_receiver,gamma_receiver
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -238,7 +238,7 @@
read(27) rho_vsI
endif
-! checks simulation types are valid
+ ! checks simulation types are valid
if( (.not. ACOUSTIC_SIMULATION ) .and. &
(.not. ELASTIC_SIMULATION ) .and. &
(.not. POROELASTIC_SIMULATION ) ) then
@@ -246,7 +246,7 @@
call exit_mpi(myrank,'error no simulation type defined')
endif
-! absorbing boundary surface
+ ! absorbing boundary surface
read(27) num_abs_boundary_faces
allocate(abs_boundary_ispec(num_abs_boundary_faces), &
abs_boundary_ijk(3,NGLLSQUARE,num_abs_boundary_faces), &
@@ -260,7 +260,7 @@
read(27) abs_boundary_normal
endif
-! free surface
+ ! free surface
read(27) num_free_surface_faces
allocate(free_surface_ispec(num_free_surface_faces), &
free_surface_ijk(3,NGLLSQUARE,num_free_surface_faces), &
@@ -274,7 +274,7 @@
read(27) free_surface_normal
endif
-! acoustic-elastic coupling surface
+ ! acoustic-elastic coupling surface
read(27) num_coupling_ac_el_faces
allocate(coupling_ac_el_normal(NDIM,NGLLSQUARE,num_coupling_ac_el_faces), &
coupling_ac_el_jacobian2Dw(NGLLSQUARE,num_coupling_ac_el_faces), &
@@ -288,7 +288,7 @@
read(27) coupling_ac_el_normal
endif
-! acoustic-poroelastic coupling surface
+ ! acoustic-poroelastic coupling surface
read(27) num_coupling_ac_po_faces
allocate(coupling_ac_po_normal(NDIM,NGLLSQUARE,num_coupling_ac_po_faces), &
coupling_ac_po_jacobian2Dw(NGLLSQUARE,num_coupling_ac_po_faces), &
@@ -302,7 +302,7 @@
read(27) coupling_ac_po_normal
endif
-! elastic-poroelastic coupling surface
+ ! elastic-poroelastic coupling surface
read(27) num_coupling_el_po_faces
allocate(coupling_el_po_normal(NDIM,NGLLSQUARE,num_coupling_el_po_faces), &
coupling_el_po_jacobian2Dw(NGLLSQUARE,num_coupling_el_po_faces), &
@@ -320,7 +320,7 @@
read(27) coupling_el_po_normal
endif
-! MPI interfaces
+ ! MPI interfaces
read(27) num_interfaces_ext_mesh
allocate(my_neighbours_ext_mesh(num_interfaces_ext_mesh), &
nibool_interfaces_ext_mesh(num_interfaces_ext_mesh),stat=ier)
@@ -361,7 +361,7 @@
read(27) c66store
endif
-! inner / outer elements
+ ! inner / outer elements
allocate(ispec_is_inner(NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array ispec_is_inner'
read(27) ispec_is_inner
@@ -408,7 +408,7 @@
if( myrank == 0 ) then
write(IMAIN,*) 'total poroelastic elements :',inum
endif
-
+
! debug
!call sum_all_i(num_interfaces_ext_mesh,inum)
!if(myrank == 0) then
@@ -416,7 +416,7 @@
! write(IMAIN,*)
!endif
-! MPI communications
+ ! MPI communications
allocate(buffer_send_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
buffer_recv_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
buffer_send_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
@@ -435,7 +435,7 @@
request_recv_vector_ext_mesh_w(num_interfaces_ext_mesh),stat=ier)
if( ier /= 0 ) stop 'error allocating array buffer_send_vector_ext_mesh etc.'
-! gets model dimensions
+ ! gets model dimensions
minl = minval( xstore )
maxl = maxval( xstore )
call min_all_all_cr(minl,min_all)
@@ -483,7 +483,7 @@
deallocate(rho_vp,rho_vs)
endif
-! reads adjoint parameters
+ ! reads adjoint parameters
call read_mesh_databases_adjoint()
end subroutine read_mesh_databases
@@ -505,7 +505,7 @@
integer :: ier
-! allocates adjoint arrays for elastic simulations
+ ! allocates adjoint arrays for elastic simulations
if( ELASTIC_SIMULATION .and. SIMULATION_TYPE == 3 ) then
! backward displacement,velocity,acceleration fields
allocate(b_displ(NDIM,NGLOB_ADJOINT),stat=ier)
@@ -593,7 +593,7 @@
endif
-! allocates adjoint arrays for acoustic simulations
+ ! allocates adjoint arrays for acoustic simulations
if( ACOUSTIC_SIMULATION .and. SIMULATION_TYPE == 3 ) then
! backward potentials
@@ -650,8 +650,8 @@
endif
-! ADJOINT moho
-! moho boundary
+ ! ADJOINT moho
+ ! moho boundary
if( ELASTIC_SIMULATION ) then
allocate( is_moho_top(NSPEC_BOUN),is_moho_bot(NSPEC_BOUN),stat=ier)
if( ier /= 0 ) stop 'error allocating array is_moho_top etc.'
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -56,13 +56,13 @@
if(NSOURCES > 1) write(IMAIN,*) 'Using ',NSOURCES,' point sources'
endif
-end subroutine setup_sources_receivers
+ end subroutine setup_sources_receivers
!
!-------------------------------------------------------------------------------------------------
!
-subroutine setup_sources()
+ subroutine setup_sources()
use specfem_par
use specfem_par_acoustic
@@ -96,7 +96,7 @@
if( ier /= 0 ) stop 'error allocating arrays for sources'
! for source encoding (acoustic sources so far only)
- allocate(pm1_source_encoding(NSOURCES),stat=ier)
+ allocate(pm1_source_encoding(NSOURCES),stat=ier)
if( ier /= 0 ) stop 'error allocating arrays for sources'
! locate sources in the mesh
@@ -214,14 +214,14 @@
! checks if source is in an acoustic element and exactly on the free surface because pressure is zero there
call setup_sources_check_acoustic()
-end subroutine setup_sources
+ end subroutine setup_sources
!
!-------------------------------------------------------------------------------------------------
!
-subroutine setup_sources_check_acoustic()
+ subroutine setup_sources_check_acoustic()
! checks if source is in an acoustic element and exactly on the free surface because pressure is zero there
@@ -321,13 +321,13 @@
enddo ! num_free_surface_faces
-end subroutine setup_sources_check_acoustic
+ end subroutine setup_sources_check_acoustic
!
!-------------------------------------------------------------------------------------------------
!
-subroutine setup_receivers()
+ subroutine setup_receivers()
use specfem_par
use specfem_par_acoustic
@@ -417,14 +417,14 @@
! checks if acoustic receiver is exactly on the free surface because pressure is zero there
call setup_receivers_check_acoustic()
-end subroutine setup_receivers
+ end subroutine setup_receivers
!
!-------------------------------------------------------------------------------------------------
!
-subroutine setup_receivers_check_acoustic()
+ subroutine setup_receivers_check_acoustic()
! checks if acoustic receiver is exactly on the free surface because pressure is zero there
@@ -492,7 +492,7 @@
call any_all_l( is_on, is_on_all )
if( myrank == 0 .and. is_on_all ) then
! limits user output if too many receivers
- if( nrec < 1000 .and. (.not. SU_FORMAT ) ) then
+ if( nrec < 1000 .and. (.not. SU_FORMAT ) ) then
write(IMAIN,*) '**********************************************************************'
write(IMAIN,*) '*** station:',irec,' ***'
write(IMAIN,*) '*** Warning: acoustic receiver located exactly on the free surface ***'
@@ -504,14 +504,14 @@
enddo ! num_free_surface_faces
-end subroutine setup_receivers_check_acoustic
+ end subroutine setup_receivers_check_acoustic
!
!-------------------------------------------------------------------------------------------------
!
-subroutine setup_sources_precompute_arrays()
+ subroutine setup_sources_precompute_arrays()
use specfem_par
use specfem_par_elastic
@@ -567,10 +567,10 @@
factor_source = factor_source * sqrt(2.0) / sqrt(3.0)
! source encoding
- ! determines factor +/-1 depending on sign of moment tensor
+ ! determines factor +/-1 depending on sign of moment tensor
! (see e.g. Krebs et al., 2009. Fast full-wavefield seismic inversion using encoded sources,
! Geophysics, 74 (6), WCC177-WCC188.)
- pm1_source_encoding(isource) = sign(1.0,Mxx(isource))
+ pm1_source_encoding(isource) = sign(1.0,Mxx(isource))
! source array interpolated on all element gll points (only used for non point sources)
call compute_arrays_source_acoustic(xi_source(isource),eta_source(isource),gamma_source(isource),&
@@ -665,14 +665,14 @@
endif
-end subroutine setup_sources_precompute_arrays
+ end subroutine setup_sources_precompute_arrays
!
!-------------------------------------------------------------------------------------------------
!
-subroutine setup_receivers_precompute_intp()
+ subroutine setup_receivers_precompute_intp()
use specfem_par
implicit none
@@ -747,12 +747,12 @@
endif
-end subroutine setup_receivers_precompute_intp
+ end subroutine setup_receivers_precompute_intp
!
!-------------------------------------------------------------------------------------------------
!
-subroutine setup_sources_receivers_VTKfile()
+ subroutine setup_sources_receivers_VTKfile()
use specfem_par
implicit none
@@ -905,6 +905,4 @@
endif
-
-
-end subroutine setup_sources_receivers_VTKfile
+ end subroutine setup_sources_receivers_VTKfile
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -99,7 +99,7 @@
real(kind=CUSTOM_REAL) :: stf_used_total
integer :: NSOURCES
! source encoding
- ! for acoustic sources: takes +/- 1 sign, depending on sign(Mxx)[ = sign(Myy) = sign(Mzz)
+ ! for acoustic sources: takes +/- 1 sign, depending on sign(Mxx)[ = sign(Myy) = sign(Mzz)
! since they have to equal in the acoustic setting]
real(kind=CUSTOM_REAL), dimension(:), allocatable :: pm1_source_encoding
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2012-02-03 00:53:03 UTC (rev 19585)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2012-02-03 14:39:54 UTC (rev 19586)
@@ -1065,7 +1065,7 @@
if( ier /= 0 ) stop 'error allocating arrays for movie div and curl'
! calculates divergence and curl of velocity field
call wmo_movie_div_curl(NSPEC_AB,NGLOB_AB,velocs_poroelastic, &
- div_glob,curl_glob,valency, &
+ div_glob,curl_glob,valency, &
div,curl_x,curl_y,curl_z, &
velocity_x,velocity_y,velocity_z, &
ibool,ispec_is_poroelastic, &
@@ -1106,7 +1106,7 @@
!write(27) veloc
!close(27)
- endif
+ endif
if( ACOUSTIC_SIMULATION .or. ELASTIC_SIMULATION .or. POROELASTIC_SIMULATION) then
More information about the CIG-COMMITS
mailing list