[cig-commits] r21845 - in seismo/3D/SPECFEM3D/trunk/src: decompose_mesh generate_databases meshfem3D shared specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Fri Apr 12 15:24:17 PDT 2013
Author: dkomati1
Date: 2013-04-12 15:24:17 -0700 (Fri, 12 Apr 2013)
New Revision: 21845
Modified:
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/program_decompose_mesh.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/meshfem3D/check_mesh_quality.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90
seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90
seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90
seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90
seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_viscoelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
Log:
committed Vadim's modifications to the trunk (I also changed the calls to read_parameter_file() in all the other routines)
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/decompose_mesh.F90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -150,7 +150,7 @@
APPROXIMATE_OCEAN_LOAD,TOPOGRAPHY,USE_FORCE_POINT_SOURCE,FULL_ATTENUATION_SOLID
logical :: STACEY_ABSORBING_CONDITIONS,SAVE_FORWARD,STACEY_INSTEAD_OF_FREE_SURFACE
logical :: ANISOTROPY,SAVE_MESH_FILES,USE_RICKER_TIME_FUNCTION,PRINT_SOURCE_TIME_FUNCTION
- character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH
+ character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH,TRAC_PATH
integer :: IMODEL
contains
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/program_decompose_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/program_decompose_mesh.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh/program_decompose_mesh.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -43,7 +43,7 @@
STACEY_ABSORBING_CONDITIONS,SAVE_FORWARD,STACEY_INSTEAD_OF_FREE_SURFACE, &
ANISOTROPY,SAVE_MESH_FILES,USE_RICKER_TIME_FUNCTION,PRINT_SOURCE_TIME_FUNCTION, &
LOCAL_PATH,TOMOGRAPHY_PATH,PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE, &
- f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID
+ f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRAC_PATH
implicit none
@@ -81,7 +81,7 @@
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE, &
USE_RICKER_TIME_FUNCTION,OLSEN_ATTENUATION_RATIO,PML_CONDITIONS, &
- PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID)
+ PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRAC_PATH)
! reads in (CUBIT) mesh files: mesh_file,nodes_coord_file, ...
call read_mesh_files()
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -265,7 +265,7 @@
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE, &
USE_RICKER_TIME_FUNCTION,OLSEN_ATTENUATION_RATIO,PML_CONDITIONS, &
- PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID)
+ PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRAC_PATH)
! check that the code is running with the requested nb of processes
if(sizeprocs /= NPROC) then
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases_par.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -67,7 +67,7 @@
logical :: USE_RICKER_TIME_FUNCTION,PRINT_SOURCE_TIME_FUNCTION
logical :: MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT,USE_HIGHRES_FOR_MOVIES
- character(len=256) OUTPUT_FILES,LOCAL_PATH,TOMOGRAPHY_PATH
+ character(len=256) OUTPUT_FILES,LOCAL_PATH,TOMOGRAPHY_PATH,TRAC_PATH
! parameters deduced from parameters read from file
integer :: NPROC
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/check_mesh_quality.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/check_mesh_quality.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/check_mesh_quality.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -107,7 +107,7 @@
logical STACEY_ABSORBING_CONDITIONS,SAVE_FORWARD,STACEY_INSTEAD_OF_FREE_SURFACE
logical ANISOTROPY,SAVE_MESH_FILES,USE_RICKER_TIME_FUNCTION,PRINT_SOURCE_TIME_FUNCTION
logical PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE,FULL_ATTENUATION_SOLID
- character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH
+ character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH,TRAC_PATH
integer NPROC
integer MOVIE_TYPE,IMODEL
@@ -132,7 +132,7 @@
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE, &
USE_RICKER_TIME_FUNCTION,OLSEN_ATTENUATION_RATIO,PML_CONDITIONS, &
- PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID)
+ PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRAC_PATH)
if(NGNOD /= 8) stop 'error: check_mesh_quality only supports NGNOD == 8 for now'
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -341,7 +341,7 @@
logical ANISOTROPY,SAVE_MESH_FILES,USE_RICKER_TIME_FUNCTION,PRINT_SOURCE_TIME_FUNCTION
logical PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE,FULL_ATTENUATION_SOLID
integer MOVIE_TYPE,IMODEL
- character(len=256) OUTPUT_FILES,LOCAL_PATH,TOMOGRAPHY_PATH
+ character(len=256) OUTPUT_FILES,LOCAL_PATH,TOMOGRAPHY_PATH,TRAC_PATH
! ************** PROGRAM STARTS HERE **************
@@ -382,7 +382,7 @@
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE, &
USE_RICKER_TIME_FUNCTION,OLSEN_ATTENUATION_RATIO,PML_CONDITIONS, &
- PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID)
+ PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRAC_PATH)
! read the mesh parameter file
! nullify(subregions,material_properties)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -100,7 +100,7 @@
logical :: STACEY_ABSORBING_CONDITIONS,SAVE_FORWARD,STACEY_INSTEAD_OF_FREE_SURFACE
logical :: ANISOTROPY,SAVE_MESH_FILES,USE_RICKER_TIME_FUNCTION,PRINT_SOURCE_TIME_FUNCTION
logical :: PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE,FULL_ATTENUATION_SOLID
- character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH
+ character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH,TRAC_PATH
integer :: IMODEL
! checks given arguments
@@ -184,7 +184,7 @@
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE, &
USE_RICKER_TIME_FUNCTION,OLSEN_ATTENUATION_RATIO,PML_CONDITIONS, &
- PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID)
+ PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRAC_PATH)
print *, 'Slice list: '
print *, node_list(1:num_node)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2013-04-12 22:24:17 UTC (rev 21845)
@@ -442,3 +442,7 @@
integer, parameter :: IMODEL_IPATI = 10
integer, parameter :: IMODEL_IPATI_WATER = 11
+! Vadim tests
+ logical, parameter :: USE_VADIM = .false.
+ integer, parameter :: IIN_veloc_dsm = 51, IIN_tract_dsm = 52, Ntime_step_dsm = 100
+
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -139,7 +139,7 @@
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE, &
USE_RICKER_TIME_FUNCTION,OLSEN_ATTENUATION_RATIO,PML_CONDITIONS, &
- PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID)
+ PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRAC_PATH)
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', OUTPUT_FILES_PATH(1:len_trim(OUTPUT_FILES_PATH)))
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/read_parameter_file.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -34,7 +34,7 @@
SIMULATION_TYPE,SAVE_FORWARD,NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE, &
USE_RICKER_TIME_FUNCTION,OLSEN_ATTENUATION_RATIO,PML_CONDITIONS, &
- PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID)
+ PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRAC_PATH)
implicit none
@@ -53,7 +53,7 @@
logical USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE,USE_RICKER_TIME_FUNCTION
logical PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE,FULL_ATTENUATION_SOLID
- character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH,CMTSOLUTION,FORCESOLUTION
+ character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH,CMTSOLUTION,FORCESOLUTION,TRAC_PATH ! VM VM adds TRAC_PATH
! local variables
integer ::ios,icounter,isource,idummy,nproc_eta_old,nproc_xi_old
@@ -172,6 +172,11 @@
if(err_occurred() /= 0) return
call read_value_logical(PRINT_SOURCE_TIME_FUNCTION, 'solver.PRINT_SOURCE_TIME_FUNCTION')
if(err_occurred() /= 0) return
+ !! VM VM read the traction path directory
+ if (USE_VADIM) then
+ call read_value_string(TRAC_PATH, 'TRAC_PATH')
+ if(err_occurred() /= 0) return
+ endif
! close parameter file
call close_parameter_file()
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -108,7 +108,7 @@
logical :: STACEY_ABSORBING_CONDITIONS,SAVE_FORWARD,STACEY_INSTEAD_OF_FREE_SURFACE
logical :: ANISOTROPY,SAVE_MESH_FILES,USE_RICKER_TIME_FUNCTION,PRINT_SOURCE_TIME_FUNCTION
logical :: PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE,FULL_ATTENUATION_SOLID
- character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH
+ character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH,TRAC_PATH
integer :: MOVIE_TYPE,IMODEL
! smoothing parameters
@@ -222,7 +222,7 @@
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE, &
USE_RICKER_TIME_FUNCTION,OLSEN_ATTENUATION_RATIO,PML_CONDITIONS, &
- PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID)
+ PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRAC_PATH)
! checks if number of MPI process as specified
if (sizeprocs /= NPROC) then
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/sum_kernels.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -115,7 +115,7 @@
logical :: STACEY_ABSORBING_CONDITIONS,SAVE_FORWARD,STACEY_INSTEAD_OF_FREE_SURFACE
logical :: ANISOTROPY,SAVE_MESH_FILES,USE_RICKER_TIME_FUNCTION,PRINT_SOURCE_TIME_FUNCTION
logical :: PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE,FULL_ATTENUATION_SOLID
- character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH
+ character(len=256) LOCAL_PATH,TOMOGRAPHY_PATH,TRAC_PATH
! ============ program starts here =====================
! initialize the MPI communicator and start the NPROCTOT MPI processes
@@ -161,7 +161,7 @@
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE, &
USE_RICKER_TIME_FUNCTION,OLSEN_ATTENUATION_RATIO,PML_CONDITIONS, &
- PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID)
+ PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRAC_PATH)
! checks if number of MPI process as specified
if (sizeprocs /= NPROC) then
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-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -242,7 +242,7 @@
endif
! adds source term (single-force/moment-tensor solution)
- call compute_add_sources_viscoelastic( NSPEC_AB,NGLOB_AB,accel, &
+ if (.not. USE_VADIM) call compute_add_sources_viscoelastic( NSPEC_AB,NGLOB_AB,accel, &
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
hdur,hdur_gaussian,tshift_src,dt,t0,sourcearrays, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_viscoelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_viscoelastic.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_viscoelastic.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -82,8 +82,26 @@
! local parameters
real(kind=CUSTOM_REAL) vx,vy,vz,nx,ny,nz,tx,ty,tz,vn,jacobianw
integer :: ispec,iglob,i,j,k,iface,igll
- !integer:: reclen1,reclen2
+! VM VM for new method
+!! DK DK for Vadim: this MUST be declared in the main program (i.e. in the calling program) and sent
+!! DK DK to this subroutine as an argument, otherwise it is allocated and deallocated every time the code
+!! DK DK enters this subroutine, thus this will be extremely slow, and also what the array contains
+!! DK DK will be lost between two calls
+ real(kind=CUSTOM_REAL) :: Veloc_dsm_boundary(3,Ntime_step_dsm,NGLLSQUARE,num_abs_boundary_faces)
+ real(kind=CUSTOM_REAL) :: Tract_dsm_boundary(3,Ntime_step_dsm,NGLLSQUARE,num_abs_boundary_faces)
+
+!! DK DK for Vadim: I had to add this missing declaration; but then of course now it is declared but undefined / unassigned
+ integer :: it_dsm
+
+ if (USE_VADIM) then
+ if ( phase_is_inner .eqv. .false. ) then
+ if (mod(it_dsm,Ntime_step_dsm+1) == 0 .or. it == 1) then
+ call read_dsm_file(Veloc_dsm_boundary,Tract_dsm_boundary,num_abs_boundary_faces,it_dsm)
+ end if
+ end if
+ end if
+
! checks if anything to do
if( num_abs_boundary_faces == 0 ) return
@@ -126,7 +144,11 @@
vx=veloc(1,iglob)
vy=veloc(2,iglob)
vz=veloc(3,iglob)
-
+ if (USE_VADIM) then
+ vx = vx - Veloc_dsm_boundary(1,it_dsm,igll,iface)
+ vy = vy - Veloc_dsm_boundary(2,it_dsm,igll,iface)
+ vz = vz - Veloc_dsm_boundary(3,it_dsm,igll,iface)
+ end if
! gets associated normal
nx = abs_boundary_normal(1,igll,iface)
ny = abs_boundary_normal(2,igll,iface)
@@ -140,6 +162,12 @@
ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(vy-vn*ny)
tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(vz-vn*nz)
+ if (USE_VADIM) then
+ tx = -Tract_dsm_boundary(1,it_dsm,igll,iface) + tx
+ ty = -Tract_dsm_boundary(2,it_dsm,igll,iface) + ty
+ tz = -Tract_dsm_boundary(3,it_dsm,igll,iface) + tz
+ end if
+
! gets associated, weighted jacobian
jacobianw = abs_boundary_jacobian2Dw(igll,iface)
@@ -182,3 +210,37 @@
end subroutine compute_stacey_viscoelastic
+!---------------------------------------------------------------------------------------
+
+ subroutine read_dsm_file(Veloc_dsm_boundary,Tract_dsm_boundary,num_abs_boundary_faces)
+
+ implicit none
+
+ include "constants.h"
+
+ integer igll,it_dsm
+ integer iface,num_abs_boundary_faces,i,j
+ real(kind=CUSTOM_REAL) :: Veloc_dsm_boundary(3,Ntime_step_dsm,NGLLSQUARE,num_abs_boundary_faces)
+ real(kind=CUSTOM_REAL) :: Tract_dsm_boundary(3,Ntime_step_dsm,NGLLSQUARE,num_abs_boundary_faces)
+
+!! DK DK why use 5 and not NGLLX here? (I assume 5 means 5 GLL points here?)
+ real(kind=CUSTOM_REAL) :: dsm_boundary_tmp(3,100,5,5) !!! warning: hardwired
+
+ it_dsm = 1
+ write(*,*) 'read dsm files',it_dsm
+ do iface=1,num_abs_boundary_faces
+
+ igll = 0
+ do j=1,5 !! DK DK why use 5 and not NGLLY here? (I assume 5 means 5 GLL points here?)
+ do i=1,5 !! DK DK why use 5 and not NGLLX here? (I assume 5 means 5 GLL points here?)
+ igll = igll + 1
+ read(IIN_veloc_dsm) dsm_boundary_tmp(:,:,i,j)
+ Veloc_dsm_boundary(:,:,igll,iface) = dsm_boundary_tmp(:,:,i,j)
+ read(IIN_tract_dsm) dsm_boundary_tmp(:,:,i,j)
+ Tract_dsm_boundary(:,:,igll,iface) = dsm_boundary_tmp(:,:,i,j)
+ enddo
+ enddo
+ enddo
+
+ end subroutine read_dsm_file
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -50,7 +50,7 @@
NTSTEP_BETWEEN_READ_ADJSRC,NOISE_TOMOGRAPHY, &
USE_FORCE_POINT_SOURCE,STACEY_INSTEAD_OF_FREE_SURFACE, &
USE_RICKER_TIME_FUNCTION,OLSEN_ATTENUATION_RATIO,PML_CONDITIONS, &
- PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID)
+ PML_INSTEAD_OF_FREE_SURFACE,f0_FOR_PML,IMODEL,FULL_ATTENUATION_SOLID,TRAC_PATH)
!! DK DK added this for now (March 2013) because CPML is not yet implemented for USE_DEVILLE_PRODUCTS;
!! DK DK we will soon add it (in a month or so)
@@ -133,8 +133,9 @@
write(IMAIN,*)
endif
- ! reads in numbers of spectral elements and points for this process' domain
+ ! reads in numbers of spectral elements and points for the part of the mesh handled by this process
call create_name_database(prname,myrank,LOCAL_PATH)
+ if (USE_VADIM) call create_name_database(dsmname,myrank,TRAC_PATH) !! VM VM
open(unit=27,file=prname(1:len_trim(prname))//'external_mesh.bin',status='old',&
action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -378,6 +378,18 @@
abs_boundary_jacobian2Dw(NGLLSQUARE,num_abs_boundary_faces), &
abs_boundary_normal(NDIM,NGLLSQUARE,num_abs_boundary_faces),stat=ier)
if( ier /= 0 ) stop 'error allocating array abs_boundary_ispec etc.'
+
+ if (USE_VADIM) then
+ ! VM for new method
+!! DK DK for Vadim: these two arrays are undeclared, thus I comment them out for now otherwise the code does not compile
+! allocate(Veloc_dsm_boundary(3,Ntime_step_dsm,NGLLSQUARE,num_abs_boundary_faces))
+! allocate(Tract_dsm_boundary(3,Ntime_step_dsm,NGLLSQUARE,num_abs_boundary_faces))
+ open(unit=IIN_veloc_dsm,file=dsmname(1:len_trim(dsmname))//'vel.bin',status='old', &
+ action='read',form='unformatted',iostat=ier)
+ open(unit=IIN_tract_dsm,file=dsmname(1:len_trim(dsmname))//'tract.bin',status='old', &
+ action='read',form='unformatted',iostat=ier)
+ endif
+
if( num_abs_boundary_faces > 0 ) then
read(27) abs_boundary_ispec
read(27) abs_boundary_ijk
@@ -716,8 +728,11 @@
if( ier /= 0 ) stop 'error allocating array kappa_kl'
! anisotropic kernels
- allocate(cijkl_kl(21,NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
- if( ier /= 0 ) stop 'error allocating array cijkl_kl'
+!! DK DK commented this out for now; must be made optional
+! allocate(cijkl_kl(21,NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+! if( ier /= 0 ) stop 'error allocating array cijkl_kl'
+!! DK DK added this for now
+ allocate(cijkl_kl(1,1,1,1,1),stat=ier)
! derived kernels
! density prime kernel
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2013-04-12 20:39:24 UTC (rev 21844)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2013-04-12 22:24:17 UTC (rev 21845)
@@ -194,7 +194,7 @@
integer :: NPROC_XI,NPROC_ETA
double precision :: LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX
- character(len=256) OUTPUT_FILES,LOCAL_PATH,TOMOGRAPHY_PATH,prname
+ character(len=256) OUTPUT_FILES,LOCAL_PATH,TOMOGRAPHY_PATH,prname,dsmname,TRAC_PATH
! names of the data files for all the processors in MPI
character(len=256) outputname
More information about the CIG-COMMITS
mailing list