[cig-commits] r22779 - in seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER: . setup src/create_header_file src/cuda src/meshfem3D src/shared src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Sep 9 04:33:20 PDT 2013
Author: danielpeter
Date: 2013-09-09 04:33:19 -0700 (Mon, 09 Sep 2013)
New Revision: 22779
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/create_header_file/rules.mk
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/finalize_mesher.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_par.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/adios_method_stubs.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_seismograms.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time_undoatt.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/multiply_arrays_source.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_topography_bathymetry.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/rules.mk
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_forward_arrays.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/update_displacement_Newmark.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_output_ASCII.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_output_SAC.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90
Log:
updates undoatt routines, seismogram outputs and force vectorization routines
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in 2013-09-09 11:33:19 UTC (rev 22779)
@@ -75,7 +75,7 @@
CUDA_LIB_LOCATION = @CUDA_LIB@
CUDA_LINK = $(CUDA_LIB_LOCATION) $(CUDA_LIBS)
-CUDA_INC = @CUDA_INC@ -I../../setup -I../../
+CUDA_INC = @CUDA_INC@ -I${SETUP}
MPI_INC = @MPI_INC@
@COND_CUDA_TRUE at NVCC = nvcc
@@ -119,11 +119,11 @@
@COND_VTK_TRUE at VTK = yes
@COND_VTK_FALSE at VTK = no
- at COND_VTK_TRUE@CPPFLAGS =${CPPFLAGS} -DWITH_VTK
- at COND_VTK_TRUE@FCCOMPILE_CHECK =${FCCOMPILE_CHECK} -DWITH_VTK
-#@COND_VTK_TRUE at FCCOMPILE_NO_CHECK =${FCCOMPILE_NO_CHECK} -DWITH_VTK
- at COND_VTK_TRUE@MPIFCCOMPILE_CHECK =${MPIFCCOMPILE_CHECK} -DWITH_VTK
-#@COND_VTK_TRUE at MPIFCCOMPILE_NO_CHECK =${MPIFCCOMPILE_NO_CHECK} -DWITH_VTK
+ at COND_VTK_TRUE@CPPFLAGS += -DWITH_VTK
+ at COND_VTK_TRUE@FCCOMPILE_CHECK += -DWITH_VTK
+#@COND_VTK_TRUE at FCCOMPILE_NO_CHECK += -DWITH_VTK
+ at COND_VTK_TRUE@MPIFCCOMPILE_CHECK += -DWITH_VTK
+#@COND_VTK_TRUE at MPIFCCOMPILE_NO_CHECK += -DWITH_VTK
#######################################
####
@@ -145,11 +145,11 @@
@COND_VECTORIZATION_TRUE at FORCE_VECTORIZATION = yes
@COND_VECTORIZATION_FALSE at FORCE_VECTORIZATION = no
- at COND_VTK_TRUE@CPPFLAGS =${CPPFLAGS} -DFORCE_VECTORIZATION
- at COND_VTK_TRUE@FCCOMPILE_CHECK =${FCCOMPILE_CHECK} -DFORCE_VECTORIZATION
-#@COND_VTK_TRUE at FCCOMPILE_NO_CHECK =${FCCOMPILE_NO_CHECK} -DFORCE_VECTORIZATION
- at COND_VTK_TRUE@MPIFCCOMPILE_CHECK =${MPIFCCOMPILE_CHECK} -DFORCE_VECTORIZATION
-#@COND_VTK_TRUE at MPIFCCOMPILE_NO_CHECK =${MPIFCCOMPILE_NO_CHECK} -DFORCE_VECTORIZATION
+ at COND_VECTORIZATION_TRUE@CPPFLAGS += -DFORCE_VECTORIZATION
+ at COND_VECTORIZATION_TRUE@FCCOMPILE_CHECK += -DFORCE_VECTORIZATION
+#@COND_VECTORIZATION_TRUE at FCCOMPILE_NO_CHECK += -DFORCE_VECTORIZATION
+ at COND_VECTORIZATION_TRUE@MPIFCCOMPILE_CHECK += -DFORCE_VECTORIZATION
+#@COND_VECTORIZATION_TRUE at MPIFCCOMPILE_NO_CHECK += -DFORCE_VECTORIZATION
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in 2013-09-09 11:33:19 UTC (rev 22779)
@@ -236,7 +236,7 @@
! the mesher) and use them for the computation of boundary kernel (in the solver)
logical, parameter :: SAVE_BOUNDARY_MESH = .false.
-!*********************************************************************************************************
+!************************************************************
! added these parameters for the GPU version of the solver with mesh coloring
! add mesh coloring for the GPU + MPI implementation
@@ -276,12 +276,10 @@
!! (uncomment desired resolution)
! old version
! old version 5.1.5 uses full 3d attenuation arrays (set to .true.), custom_real for attenuation arrays (Qmu_store, tau_e_store)
- logical, parameter :: USE_VERSION_5_1_5 = .true.
-! integer, parameter :: CUSTOM_REAL_ATT = CUSTOM_REAL
+! logical, parameter :: USE_VERSION_5_1_5 = .true.
! new version
-! new version uses full 3d attenuation, double precision for attenuation arrays (Qmu_store, tau_e_store)
-!! logical, parameter :: USE_VERSION_5_1_5 = .false.
-!! integer,parameter :: CUSTOM_REAL_ATT = SIZE_DOUBLE
+! new version uses optional full 3d attenuation
+ logical, parameter :: USE_VERSION_5_1_5 = .false.
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/create_header_file/rules.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/create_header_file/rules.mk 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/create_header_file/rules.mk 2013-09-09 11:33:19 UTC (rev 22779)
@@ -67,7 +67,11 @@
${OUTPUT}/values_from_mesher.h: $E/xcreate_header_file $B/DATA/Par_file
@-rm -f $@
+ @echo ""
+ @echo "running xcreate_header_file..."
+ @echo ""
$E/xcreate_header_file
+ @echo ""
@test -f $@
#######################################
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu 2013-09-09 11:33:19 UTC (rev 22779)
@@ -608,7 +608,6 @@
COMPUTE_COUPLING_OCEAN_CUDA)(long* Mesh_pointer_f,
int* NCHUNKS_VAL,
int* exact_mass_matrix_for_rotation,
- int* use_lddrk,
int* FORWARD_OR_ADJOINT) {
TRACE("compute_coupling_ocean_cuda");
@@ -628,8 +627,7 @@
dim3 grid(num_blocks_x,num_blocks_y);
dim3 threads(blocksize,1,1);
- if( ( *NCHUNKS_VAL != 6 && mp->absorbing_conditions || (mp->rotation && *exact_mass_matrix_for_rotation)) &&
- ! *use_lddrk ){
+ if( ( *NCHUNKS_VAL != 6 && mp->absorbing_conditions) || (mp->rotation && *exact_mass_matrix_for_rotation) ){
// uses corrected mass matrices
if( *FORWARD_OR_ADJOINT == 1 ){
compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_accel_crust_mantle,
@@ -643,9 +641,9 @@
}else if( *FORWARD_OR_ADJOINT == 3){
// for backward/reconstructed potentials
compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_b_accel_crust_mantle,
- mp->d_rmassx_crust_mantle,
- mp->d_rmassy_crust_mantle,
- mp->d_rmassz_crust_mantle,
+ mp->d_b_rmassx_crust_mantle,
+ mp->d_b_rmassy_crust_mantle,
+ mp->d_b_rmassz_crust_mantle,
mp->d_rmass_ocean_load,
mp->npoin_oceans,
mp->d_ibool_ocean_load,
@@ -665,9 +663,9 @@
}else if( *FORWARD_OR_ADJOINT == 3){
// for backward/reconstructed potentials
compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_b_accel_crust_mantle,
- mp->d_rmassz_crust_mantle,
- mp->d_rmassz_crust_mantle,
- mp->d_rmassz_crust_mantle,
+ mp->d_b_rmassz_crust_mantle,
+ mp->d_b_rmassz_crust_mantle,
+ mp->d_b_rmassz_crust_mantle,
mp->d_rmass_ocean_load,
mp->npoin_oceans,
mp->d_ibool_ocean_load,
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h 2013-09-09 11:33:19 UTC (rev 22779)
@@ -237,10 +237,16 @@
realw* d_kappavstore_crust_mantle; realw* d_muvstore_crust_mantle;
realw* d_kappahstore_crust_mantle; realw* d_muhstore_crust_mantle;
realw* d_eta_anisostore_crust_mantle;
+
+ // mass matrices
realw* d_rmassx_crust_mantle;
realw* d_rmassy_crust_mantle;
realw* d_rmassz_crust_mantle;
+ realw* d_b_rmassx_crust_mantle;
+ realw* d_b_rmassy_crust_mantle;
+ realw* d_b_rmassz_crust_mantle;
+
// global indexing
int* d_ibool_crust_mantle;
int* d_ispec_is_tiso_crust_mantle;
@@ -348,7 +354,9 @@
// model parameters
realw* d_rhostore_outer_core; realw* d_kappavstore_outer_core;
+
realw* d_rmass_outer_core;
+ realw* d_b_rmass_outer_core;
// global indexing
int* d_ibool_outer_core;
@@ -411,9 +419,13 @@
// model parameters
realw* d_rhostore_inner_core;
realw* d_kappavstore_inner_core; realw* d_muvstore_inner_core;
+
realw* d_rmassx_inner_core;
realw* d_rmassy_inner_core;
realw* d_rmassz_inner_core;
+ realw* d_b_rmassx_inner_core;
+ realw* d_b_rmassy_inner_core;
+ realw* d_b_rmassz_inner_core;
// global indexing
int* d_ibool_inner_core;
@@ -553,6 +565,7 @@
int anisotropic_3D_mantle;
int gravity;
int rotation;
+ int exact_mass_matrix_for_rotation;
int oceans;
int anisotropic_inner_core;
int save_boundary_mesh;
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2013-09-09 11:33:19 UTC (rev 22779)
@@ -134,7 +134,7 @@
int* ABSORBING_CONDITIONS_f,
int* OCEANS_f,
int* GRAVITY_f,
- int* ROTATION_f,
+ int* ROTATION_f,int* EXACT_MASS_MATRIX_FOR_ROTATION_f,
int* ATTENUATION_f,int* UNDO_ATTENUATION_f,
int* USE_ATTENUATION_MIMIC_f,
int* USE_3D_ATTENUATION_ARRAYS_f,
@@ -236,6 +236,7 @@
mp->oceans = *OCEANS_f;
mp->gravity = *GRAVITY_f;
mp->rotation = *ROTATION_f;
+ mp->exact_mass_matrix_for_rotation = *EXACT_MASS_MATRIX_FOR_ROTATION_f;
mp->attenuation = *ATTENUATION_f;
mp->undo_attenuation = *UNDO_ATTENUATION_f;
@@ -1085,9 +1086,8 @@
realw* h_kappav, realw* h_muv,
realw* h_kappah, realw* h_muh,
realw* h_eta_aniso,
- realw* h_rmassx,
- realw* h_rmassy,
- realw* h_rmassz,
+ realw* h_rmassx,realw* h_rmassy,realw* h_rmassz,
+ realw* h_b_rmassx,realw* h_b_rmassy,
int* h_ibool,
realw* h_xstore, realw* h_ystore, realw* h_zstore,
int* h_ispec_is_tiso,
@@ -1361,7 +1361,7 @@
// mass matrices
copy_todevice_realw((void**)&mp->d_rmassz_crust_mantle,h_rmassz,size_glob);
- if( *NCHUNKS_VAL != 6 && mp->absorbing_conditions){
+ if( (*NCHUNKS_VAL != 6 && mp->absorbing_conditions) || ( mp->rotation && mp->exact_mass_matrix_for_rotation)){
copy_todevice_realw((void**)&mp->d_rmassx_crust_mantle,h_rmassx,size_glob);
copy_todevice_realw((void**)&mp->d_rmassy_crust_mantle,h_rmassy,size_glob);
}else{
@@ -1369,8 +1369,18 @@
mp->d_rmassy_crust_mantle = mp->d_rmassz_crust_mantle;
}
- // kernels
+ // kernel simulations
if( mp->simulation_type == 3 ){
+ mp->d_b_rmassz_crust_mantle = mp->d_rmassz_crust_mantle;
+ if( mp->rotation && mp->exact_mass_matrix_for_rotation ){
+ copy_todevice_realw((void**)&mp->d_b_rmassx_crust_mantle,h_b_rmassx,size_glob);
+ copy_todevice_realw((void**)&mp->d_b_rmassy_crust_mantle,h_b_rmassy,size_glob);
+ }else{
+ mp->d_b_rmassx_crust_mantle = mp->d_rmassx_crust_mantle;
+ mp->d_b_rmassy_crust_mantle = mp->d_rmassy_crust_mantle;
+ }
+
+ // kernels
size = NGLL3*(mp->NSPEC_CRUST_MANTLE);
// density kernel
@@ -1575,8 +1585,12 @@
// mass matrix
copy_todevice_realw((void**)&mp->d_rmass_outer_core,h_rmass,size_glob);
- // kernels
+ // kernel simulations
if( mp->simulation_type == 3 ){
+ // mass matrix
+ mp->d_b_rmass_outer_core = mp->d_rmass_outer_core;
+
+ //kernels
size = NGLL3*(mp->NSPEC_OUTER_CORE);
// density kernel
@@ -1615,6 +1629,7 @@
realw* h_gammax, realw* h_gammay, realw* h_gammaz,
realw* h_rho, realw* h_kappav, realw* h_muv,
realw* h_rmassx,realw* h_rmassy,realw* h_rmassz,
+ realw* h_b_rmassx,realw* h_b_rmassy,
int* h_ibool,
realw* h_xstore, realw* h_ystore, realw* h_zstore,
realw *c11store,realw *c12store,realw *c13store,
@@ -1791,12 +1806,28 @@
#endif
// mass matrix
- copy_todevice_realw((void**)&mp->d_rmassx_inner_core,h_rmassx,size_glob);
- copy_todevice_realw((void**)&mp->d_rmassy_inner_core,h_rmassy,size_glob);
copy_todevice_realw((void**)&mp->d_rmassz_inner_core,h_rmassz,size_glob);
+ if( mp->rotation && mp->exact_mass_matrix_for_rotation ){
+ copy_todevice_realw((void**)&mp->d_rmassx_inner_core,h_rmassx,size_glob);
+ copy_todevice_realw((void**)&mp->d_rmassy_inner_core,h_rmassy,size_glob);
+ }else{
+ mp->d_rmassx_inner_core = mp->d_rmassz_inner_core;
+ mp->d_rmassy_inner_core = mp->d_rmassz_inner_core;
+ }
- // kernels
+ // kernel simulations
if( mp->simulation_type == 3 ){
+ // mass matrices
+ mp->d_b_rmassz_inner_core = mp->d_rmassz_inner_core;
+ if( mp->rotation && mp->exact_mass_matrix_for_rotation ){
+ copy_todevice_realw((void**)&mp->d_b_rmassx_inner_core,h_b_rmassx,size_glob);
+ copy_todevice_realw((void**)&mp->d_b_rmassy_inner_core,h_b_rmassy,size_glob);
+ }else{
+ mp->d_b_rmassx_inner_core = mp->d_rmassx_inner_core;
+ mp->d_b_rmassy_inner_core = mp->d_rmassy_inner_core;
+ }
+
+ // kernels
size = NGLL3*(mp->NSPEC_INNER_CORE);
// density kernel
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/save_and_compare_cpu_vs_gpu.c 2013-09-09 11:33:19 UTC (rev 22779)
@@ -30,6 +30,7 @@
#include <stdlib.h>
#include <math.h>
#include <errno.h>
+#include <string.h>
#ifdef WITH_MPI
#include <mpi.h>
@@ -109,6 +110,7 @@
float* vector = (float*)malloc(nodes_per_iteration*sizeof(float));
float max_val;
int i;
+ max_val = 0.0;
for(it=0;it<*NSTEP;it++) {
int pos = (sizeof(float)*nodes_per_iteration)*(it);
fseek(fp,pos,SEEK_SET);
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c 2013-09-09 11:33:19 UTC (rev 22779)
@@ -174,7 +174,6 @@
COMPUTE_COUPLING_OCEAN_CUDA)(long* Mesh_pointer_f,
int* NCHUNKS_VAL,
int* exact_mass_matrix_for_rotation,
- int* use_lddrk,
int* FORWARD_OR_ADJOINT) {}
@@ -370,7 +369,7 @@
int* ABSORBING_CONDITIONS_f,
int* OCEANS_f,
int* GRAVITY_f,
- int* ROTATION_f,
+ int* ROTATION_f,int* EXACT_MASS_MATRIX_FOR_ROTATION_f,
int* ATTENUATION_f,int* UNDO_ATTENUATION_f,
int* USE_ATTENUATION_MIMIC_f,
int* USE_3D_ATTENUATION_ARRAYS_f,
@@ -536,9 +535,8 @@
realw* h_kappav, realw* h_muv,
realw* h_kappah, realw* h_muh,
realw* h_eta_aniso,
- realw* h_rmassx,
- realw* h_rmassy,
- realw* h_rmassz,
+ realw* h_rmassx,realw* h_rmassy,realw* h_rmassz,
+ realw* h_b_rmassx,realw* h_b_rmassy,
int* h_ibool,
realw* h_xstore, realw* h_ystore, realw* h_zstore,
int* h_ispec_is_tiso,
@@ -592,6 +590,7 @@
realw* h_gammax, realw* h_gammay, realw* h_gammaz,
realw* h_rho, realw* h_kappav, realw* h_muv,
realw* h_rmassx,realw* h_rmassy,realw* h_rmassz,
+ realw* h_b_rmassx,realw* h_b_rmassy,
int* h_ibool,
realw* h_xstore, realw* h_ystore, realw* h_zstore,
realw *c11store,realw *c12store,realw *c13store,
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -48,7 +48,7 @@
use meshfem3D_par,only: &
NCHUNKS,ABSORBING_CONDITIONS, &
- ROTATION,EXACT_MASS_MATRIX_FOR_ROTATION
+ ROTATION,EXACT_MASS_MATRIX_FOR_ROTATION,INCLUDE_CENTRAL_CUBE
use create_regions_mesh_par,only: &
wxgll,wygll,wzgll
@@ -186,7 +186,7 @@
endif
! check that mass matrix is positive
- if( iregion_code == IREGION_INNER_CORE ) then
+ if( iregion_code == IREGION_INNER_CORE .and. INCLUDE_CENTRAL_CUBE ) then
! note: in fictitious elements mass matrix is still zero
if(minval(rmassz(:)) < 0._CUSTOM_REAL) call exit_MPI(myrank,'negative rmassz matrix term')
! check that the additional mass matrices are positive, if they exist
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/finalize_mesher.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/finalize_mesher.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/finalize_mesher.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -35,6 +35,7 @@
! local parameters
! timing
+ double precision :: tCPU
double precision, external :: wtime
!--- print number of points and elements in the mesh for each region
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_par.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_par.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -169,7 +169,7 @@
integer :: numelem_total
! timer MPI
- double precision :: time_start,tCPU
+ double precision :: time_start
! addressing for all the slices
integer, dimension(:), allocatable :: ichunk_slice,iproc_xi_slice,iproc_eta_slice
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/adios_method_stubs.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/adios_method_stubs.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/adios_method_stubs.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -36,6 +36,8 @@
subroutine adios_setup()
end subroutine
+! for xmeshfem3D compilation
+
subroutine crm_save_mesh_files_adios()
end subroutine
@@ -45,6 +47,9 @@
subroutine save_arrays_solver_adios()
end subroutine
+ subroutine save_arrays_solver_boundary_adios()
+ end subroutine
+
subroutine save_mpi_arrays_adios()
end subroutine
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -254,6 +254,9 @@
! produces simulations compatible with old globe version 5.1.5
if( USE_VERSION_5_1_5 ) then
+ print*
+ print*,'**************'
+ print*,'using globe version 5.1.5 compatible simulation parameters'
if( .not. ATTENUATION_1D_WITH_3D_STORAGE ) then
print*,'setting ATTENUATION_1D_WITH_3D_STORAGE to .true. for compatibility with globe version 5.1.5 '
ATTENUATION_1D_WITH_3D_STORAGE = .true.
@@ -270,6 +273,8 @@
print*,'setting EXACT_MASS_MATRIX_FOR_ROTATION to .false. for compatibility with globe version 5.1.5 '
EXACT_MASS_MATRIX_FOR_ROTATION = .false.
endif
+ print*,'**************'
+ print*
endif
!----------------------------------------------
@@ -296,8 +301,8 @@
if( EXACT_MASS_MATRIX_FOR_ROTATION .and. GPU_MODE ) &
stop 'EXACT_MASS_MATRIX_FOR_ROTATION support not implemented yet for GPU simulations'
- if( UNDO_ATTENUATION ) &
- stop 'UNDO_ATTENUATION support not implemented yet'
+ !if( UNDO_ATTENUATION ) &
+ ! stop 'UNDO_ATTENUATION support not implemented yet'
if( UNDO_ATTENUATION .and. NOISE_TOMOGRAPHY > 0 ) &
stop 'UNDO_ATTENUATION support not implemented yet for noise simulations'
if( UNDO_ATTENUATION .and. MOVIE_VOLUME .and. MOVIE_VOLUME_TYPE == 4 ) &
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/rules.mk 2013-09-09 11:33:19 UTC (rev 22779)
@@ -67,6 +67,7 @@
$O/save_header_file.shared.o \
$O/spline_routines.shared.o \
$O/write_c_binary.cc.o \
+ $O/write_VTK_file.shared.o \
$(EMPTY_MACRO)
ADIOS_OBJECTS = \
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -25,6 +25,7 @@
!
!=====================================================================
+
subroutine check_stability()
! computes the maximum of the norm of the displacement
@@ -59,7 +60,10 @@
real(kind=CUSTOM_REAL) Strain_norm,Strain_norm_all,Strain2_norm,Strain2_norm_all
real(kind=CUSTOM_REAL) b_Usolidnorm,b_Usolidnorm_all,b_Ufluidnorm,b_Ufluidnorm_all
! timer MPI
- double precision :: tCPU,t_remain,t_total
+ double precision :: tCPU
+ double precision, external :: wtime
+ double precision :: time
+ double precision :: t_remain,t_total
integer :: ihours,iminutes,iseconds,int_tCPU, &
ihours_remain,iminutes_remain,iseconds_remain,int_t_remain, &
ihours_total,iminutes_total,iseconds_total,int_t_total
@@ -78,8 +82,6 @@
integer :: year,mon,day,hr,minutes,timestamp,julian_day_number,day_of_week, &
timestamp_remote,year_remote,mon_remote,day_remote,hr_remote,minutes_remote,day_of_week_remote
integer, external :: idaywk
- ! timing
- double precision, external :: wtime
double precision,parameter :: scale_displ = R_EARTH
@@ -193,9 +195,12 @@
iminutes_total = (int_t_total - 3600*ihours_total) / 60
iseconds_total = int_t_total - 3600*ihours_total - 60*iminutes_total
+ ! current time (in seconds)
+ time = dble(it-1)*DT - t0
+
! user output
write(IMAIN,*) 'Time step # ',it
- write(IMAIN,*) 'Time: ',sngl(((it-1)*DT-t0)/60.d0),' minutes'
+ write(IMAIN,*) 'Time: ',sngl((time)/60.d0),' minutes'
! rescale maximum displacement to correct dimensions
Usolidnorm_all = Usolidnorm_all * sngl(scale_displ)
@@ -363,6 +368,133 @@
!------------------------------------------------------------------------------------------------------------------
!
+ subroutine check_stability_backward()
+
+! only for backward/reconstructed wavefield
+
+ use constants,only: CUSTOM_REAL,IMAIN,R_EARTH, &
+ ADD_TIME_ESTIMATE_ELSEWHERE,HOURS_TIME_DIFFERENCE,MINUTES_TIME_DIFFERENCE, &
+ STABILITY_THRESHOLD
+
+ use specfem_par,only: &
+ GPU_MODE,Mesh_pointer, &
+ SIMULATION_TYPE,time_start,DT,t0, &
+ NSTEP,it,it_begin,it_end,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN, &
+ myrank
+
+ use specfem_par_crustmantle,only: b_displ_crust_mantle
+ use specfem_par_innercore,only: b_displ_inner_core
+ use specfem_par_outercore,only: b_displ_outer_core
+
+ implicit none
+
+ ! local parameters
+ ! maximum of the norm of the displacement and of the potential in the fluid
+ real(kind=CUSTOM_REAL) b_Usolidnorm,b_Usolidnorm_all,b_Ufluidnorm,b_Ufluidnorm_all
+ ! timer MPI
+ double precision :: tCPU
+ double precision, external :: wtime
+ double precision :: time
+ integer :: ihours,iminutes,iseconds,int_tCPU
+
+ integer :: it_run,nstep_run
+ logical :: SHOW_SEPARATE_RUN_INFORMATION
+
+ double precision,parameter :: scale_displ = R_EARTH
+
+ ! checks if anything to do
+ if( SIMULATION_TYPE /= 3 ) return
+
+ ! compute maximum of norm of displacement in each slice
+ if( .not. GPU_MODE) then
+ ! on CPU
+ b_Usolidnorm = max( &
+ maxval(sqrt(b_displ_crust_mantle(1,:)**2 + &
+ b_displ_crust_mantle(2,:)**2 + b_displ_crust_mantle(3,:)**2)), &
+ maxval(sqrt(b_displ_inner_core(1,:)**2 &
+ + b_displ_inner_core(2,:)**2 &
+ + b_displ_inner_core(3,:)**2)))
+
+ b_Ufluidnorm = maxval(abs(b_displ_outer_core))
+ else
+ ! on GPU
+ call check_norm_elastic_from_device(b_Usolidnorm,Mesh_pointer,3)
+ call check_norm_acoustic_from_device(b_Ufluidnorm,Mesh_pointer,3)
+ endif
+
+ if(b_Usolidnorm > STABILITY_THRESHOLD .or. b_Usolidnorm < 0) &
+ call exit_MPI(myrank,'backward simulation became unstable and blew up in the solid')
+ if(b_Ufluidnorm > STABILITY_THRESHOLD .or. b_Ufluidnorm < 0) &
+ call exit_MPI(myrank,'backward simulation became unstable and blew up in the fluid')
+
+ ! compute the maximum of the maxima for all the slices using an MPI reduction
+ call max_all_cr(b_Usolidnorm,b_Usolidnorm_all)
+ call max_all_cr(b_Ufluidnorm,b_Ufluidnorm_all)
+
+ if(myrank == 0) then
+
+ ! this is in the case of restart files, when a given run consists of several partial runs
+ ! information about the current run only
+ SHOW_SEPARATE_RUN_INFORMATION = ( NUMBER_OF_RUNS > 1 .and. NUMBER_OF_THIS_RUN < NUMBER_OF_RUNS )
+ it_run = it - it_begin + 1
+ nstep_run = it_end - it_begin + 1
+
+ ! elapsed time since beginning of the simulation
+ tCPU = wtime() - time_start
+
+ int_tCPU = int(tCPU)
+ ihours = int_tCPU / 3600
+ iminutes = (int_tCPU - 3600*ihours) / 60
+ iseconds = int_tCPU - 3600*ihours - 60*iminutes
+
+ ! no further time estimation since only partially computed solution yet...
+
+ ! current time (in seconds)
+ time = dble(it-1)*DT - t0
+
+ ! user output
+ write(IMAIN,*) 'Time step for back propagation # ',it
+ write(IMAIN,*) 'Time: ',sngl((time)/60.d0),' minutes'
+
+ ! rescale maximum displacement to correct dimensions
+ b_Usolidnorm_all = b_Usolidnorm_all * sngl(scale_displ)
+ write(IMAIN,*) 'Max norm displacement vector U in solid in all slices for back prop.(m) = ',b_Usolidnorm_all
+ write(IMAIN,*) 'Max non-dimensional potential Ufluid in fluid in all slices for back prop.= ',b_Ufluidnorm_all
+
+ write(IMAIN,*) 'Elapsed time in seconds = ',tCPU
+ write(IMAIN,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
+ write(IMAIN,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
+
+ if (SHOW_SEPARATE_RUN_INFORMATION) then
+ write(IMAIN,*) 'Time steps done for this run = ',it_run,' out of ',nstep_run
+ write(IMAIN,*) 'Time steps done in total = ',it,' out of ',NSTEP
+ write(IMAIN,*) 'Time steps remaining for this run = ',it_end - it
+ write(IMAIN,*) 'Time steps remaining for all runs = ',NSTEP - it
+ else
+ write(IMAIN,*) 'Time steps done = ',it,' out of ',NSTEP
+ write(IMAIN,*) 'Time steps remaining = ',NSTEP - it
+ endif
+ write(IMAIN,*)
+
+ ! flushes file buffer for main output file (IMAIN)
+ call flush_IMAIN()
+
+ ! check stability of the code, exit if unstable
+ ! negative values can occur with some compilers when the unstable value is greater
+ ! than the greatest possible floating-point number of the machine
+ if(b_Usolidnorm_all > STABILITY_THRESHOLD .or. b_Usolidnorm_all < 0) &
+ call exit_MPI(myrank,'backward simulation became unstable and blew up in the solid')
+ if(b_Ufluidnorm_all > STABILITY_THRESHOLD .or. b_Ufluidnorm_all < 0) &
+ call exit_MPI(myrank,'backward simulation became unstable and blew up in the fluid')
+
+ endif
+
+ end subroutine check_stability_backward
+
+!
+!------------------------------------------------------------------------------------------------------------------
+!
+
subroutine write_timestamp_file(Usolidnorm_all,Ufluidnorm_all,b_Usolidnorm_all,b_Ufluidnorm_all, &
tCPU,ihours,iminutes,iseconds, &
t_remain,ihours_remain,iminutes_remain,iseconds_remain, &
@@ -491,3 +623,37 @@
close(IOUT)
end subroutine write_timestamp_file
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine print_elapsed_time()
+
+ use specfem_par,only: time_start,IMAIN,myrank
+ implicit none
+
+ ! local parameters
+ integer :: ihours,iminutes,iseconds,int_tCPU
+ ! timing
+ double precision :: tCPU
+ double precision, external :: wtime
+
+ if(myrank == 0) then
+ ! elapsed time since beginning of the simulation
+ tCPU = wtime() - time_start
+
+ int_tCPU = int(tCPU)
+ ihours = int_tCPU / 3600
+ iminutes = (int_tCPU - 3600*ihours) / 60
+ iseconds = int_tCPU - 3600*ihours - 60*iminutes
+ write(IMAIN,*) 'Time-Loop Complete. Timing info:'
+ write(IMAIN,*) 'Total elapsed time in seconds = ',tCPU
+ write(IMAIN,"(' Total elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
+ call flush_IMAIN()
+ endif
+
+ end subroutine print_elapsed_time
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -359,11 +359,10 @@
rmassx_crust_mantle, rmassy_crust_mantle, rmassz_crust_mantle, &
rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load,NGLOB_XY, &
+ updated_dof_ocean_load, &
nspec_top)
use constants_solver
-! use specfem_par,only: ABSORBING_CONDITIONS,EXACT_MASS_MATRIX_FOR_ROTATION
implicit none
@@ -376,10 +375,9 @@
! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
!
! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
- ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
- integer :: NGLOB_XY
- real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassx_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLOB_XY) :: rmassy_crust_mantle
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be pointers to it
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmassx_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmassy_crust_mantle
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmassz_crust_mantle
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
@@ -395,110 +393,57 @@
! local parameters
real(kind=CUSTOM_REAL) :: force_normal_comp
real(kind=CUSTOM_REAL) :: additional_term_x,additional_term_y,additional_term_z
-! real(kind=CUSTOM_REAL) :: additional_term
real(kind=CUSTOM_REAL) :: nx,ny,nz
integer :: i,j,k,ispec,ispec2D,iglob
! initialize the updates
updated_dof_ocean_load(:) = .false.
-!daniel debug: check not needed anymore
-! if( (NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) .or. &
-! (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION) ) then
+ ! for surface elements exactly at the top of the crust (ocean bottom)
+ do ispec2D = 1,nspec_top !NSPEC2D_TOP(IREGION_CRUST_MANTLE)
- ! for surface elements exactly at the top of the crust (ocean bottom)
- do ispec2D = 1,nspec_top !NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+ ispec = ibelm_top_crust_mantle(ispec2D)
- ispec = ibelm_top_crust_mantle(ispec2D)
+ ! only for DOFs exactly at the top of the crust (ocean bottom)
+ k = NGLLZ
- ! only for DOFs exactly at the top of the crust (ocean bottom)
- k = NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
- do j = 1,NGLLY
- do i = 1,NGLLX
+ ! get global point number
+ iglob = ibool_crust_mantle(i,j,k,ispec)
- ! get global point number
- iglob = ibool_crust_mantle(i,j,k,ispec)
+ ! only update once
+ if(.not. updated_dof_ocean_load(iglob)) then
- ! only update once
- if(.not. updated_dof_ocean_load(iglob)) then
+ ! get normal
+ nx = normal_top_crust_mantle(1,i,j,ispec2D)
+ ny = normal_top_crust_mantle(2,i,j,ispec2D)
+ nz = normal_top_crust_mantle(3,i,j,ispec2D)
- ! get normal
- nx = normal_top_crust_mantle(1,i,j,ispec2D)
- ny = normal_top_crust_mantle(2,i,j,ispec2D)
- nz = normal_top_crust_mantle(3,i,j,ispec2D)
+ ! make updated component of right-hand side
+ ! we divide by rmass_crust_mantle() which is 1 / M
+ ! we use the total force which includes the Coriolis term above
+ force_normal_comp = accel_crust_mantle(1,iglob)*nx / rmassx_crust_mantle(iglob) &
+ + accel_crust_mantle(2,iglob)*ny / rmassy_crust_mantle(iglob) &
+ + accel_crust_mantle(3,iglob)*nz / rmassz_crust_mantle(iglob)
- ! make updated component of right-hand side
- ! we divide by rmass_crust_mantle() which is 1 / M
- ! we use the total force which includes the Coriolis term above
- force_normal_comp = accel_crust_mantle(1,iglob)*nx / rmassx_crust_mantle(iglob) &
- + accel_crust_mantle(2,iglob)*ny / rmassy_crust_mantle(iglob) &
- + accel_crust_mantle(3,iglob)*nz / rmassz_crust_mantle(iglob)
+ additional_term_x = (rmass_ocean_load(iglob) - rmassx_crust_mantle(iglob)) * force_normal_comp
+ additional_term_y = (rmass_ocean_load(iglob) - rmassy_crust_mantle(iglob)) * force_normal_comp
+ additional_term_z = (rmass_ocean_load(iglob) - rmassz_crust_mantle(iglob)) * force_normal_comp
- additional_term_x = (rmass_ocean_load(iglob) - rmassx_crust_mantle(iglob)) * force_normal_comp
- additional_term_y = (rmass_ocean_load(iglob) - rmassy_crust_mantle(iglob)) * force_normal_comp
- additional_term_z = (rmass_ocean_load(iglob) - rmassz_crust_mantle(iglob)) * force_normal_comp
+ accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + additional_term_x * nx
+ accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + additional_term_y * ny
+ accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + additional_term_z * nz
- accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + additional_term_x * nx
- accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + additional_term_y * ny
- accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + additional_term_z * nz
+ ! done with this point
+ updated_dof_ocean_load(iglob) = .true.
- ! done with this point
- updated_dof_ocean_load(iglob) = .true.
+ endif
- endif
+ enddo
+ enddo
+ enddo
- enddo
- enddo
- enddo
-
-! else
-!
-! ! for surface elements exactly at the top of the crust (ocean bottom)
-! do ispec2D = 1,nspec_top !NSPEC2D_TOP(IREGION_CRUST_MANTLE)
-!
-! ispec = ibelm_top_crust_mantle(ispec2D)
-!
-! ! only for DOFs exactly at the top of the crust (ocean bottom)
-! k = NGLLZ
-!
-! do j = 1,NGLLY
-! do i = 1,NGLLX
-!
-! ! get global point number
-! iglob = ibool_crust_mantle(i,j,k,ispec)
-!
-! ! only update once
-! if(.not. updated_dof_ocean_load(iglob)) then
-!
-! ! get normal
-! nx = normal_top_crust_mantle(1,i,j,ispec2D)
-! ny = normal_top_crust_mantle(2,i,j,ispec2D)
-! nz = normal_top_crust_mantle(3,i,j,ispec2D)
-!
-! ! make updated component of right-hand side
-! ! we divide by rmass_crust_mantle() which is 1 / M
-! ! we use the total force which includes the Coriolis term above
-! force_normal_comp = (accel_crust_mantle(1,iglob)*nx + &
-! accel_crust_mantle(2,iglob)*ny + &
-! accel_crust_mantle(3,iglob)*nz) / rmassz_crust_mantle(iglob)
-!
-! additional_term = (rmass_ocean_load(iglob) - rmassz_crust_mantle(iglob)) * force_normal_comp
-!
-! accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + additional_term * nx
-! accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + additional_term * ny
-! accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + additional_term * nz
-!
-! ! done with this point
-! updated_dof_ocean_load(iglob) = .true.
-!
-! endif
-!
-! enddo
-! enddo
-! enddo
-!
-! endif
-!
end subroutine compute_coupling_ocean
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -204,7 +204,7 @@
enddo ! iphase
! multiply by the inverse of the mass matrix
- call it_multiply_accel_acoustic(NGLOB_OUTER_CORE,accel_outer_core,rmass_outer_core)
+ call multiply_accel_acoustic(NGLOB_OUTER_CORE,accel_outer_core,rmass_outer_core)
! time schemes
if( USE_LDDRK ) then
@@ -425,7 +425,7 @@
enddo ! iphase
! multiply by the inverse of the mass matrix
- call it_multiply_accel_acoustic(NGLOB_OUTER_CORE_ADJOINT,b_accel_outer_core,b_rmass_outer_core)
+ call multiply_accel_acoustic(NGLOB_OUTER_CORE_ADJOINT,b_accel_outer_core,b_rmass_outer_core)
! time schemes
if( USE_LDDRK ) then
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.F90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.F90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -95,7 +95,6 @@
integer :: i,j,k
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
real(kind=CUSTOM_REAL) :: dpotentialdxl,dpotentialdyl,dpotentialdzl
- real(kind=CUSTOM_REAL) :: sum_terms
! manually inline the calls to the Deville et al. (2002) routines
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc
@@ -124,6 +123,13 @@
integer :: num_elements,ispec_p
integer :: iphase
+#ifdef FORCE_VECTORIZATION
+ integer :: ijk
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: sum_terms
+#else
+ real(kind=CUSTOM_REAL) :: sum_terms
+#endif
+
! ****************************************************
! big loop over all spectral elements in the fluid
! ****************************************************
@@ -144,7 +150,44 @@
ispec = phase_ispec_inner(ispec_p,iphase)
! only compute element which belong to current phase (inner or outer elements)
+#ifdef FORCE_VECTORIZATION
+ do ijk=1,NGLLCUBE
+ iglob = ibool(ijk,1,1,ispec)
+ ! get a local copy of the potential field
+ dummyx_loc(ijk,1,1) = displfluid(iglob)
+
+ ! pre-computes factors
+ ! use mesh coordinates to get theta and phi
+ ! x y z contain r theta phi
+ radius = dble(xstore(iglob))
+ theta = dble(ystore(iglob))
+ phi = dble(zstore(iglob))
+
+ cos_theta = dcos(theta)
+ sin_theta = dsin(theta)
+ cos_phi = dcos(phi)
+ sin_phi = dsin(phi)
+
+ int_radius = nint(radius * R_EARTH_KM * 10.d0)
+
+ if( .not. GRAVITY_VAL ) then
+ ! grad(rho)/rho in Cartesian components
+ displ_times_grad_x_ln_rho(ijk,1,1) = dummyx_loc(ijk,1,1) &
+ * sin_theta * cos_phi * d_ln_density_dr_table(int_radius)
+ displ_times_grad_y_ln_rho(ijk,1,1) = dummyx_loc(ijk,1,1) &
+ * sin_theta * sin_phi * d_ln_density_dr_table(int_radius)
+ displ_times_grad_z_ln_rho(ijk,1,1) = dummyx_loc(ijk,1,1) &
+ * cos_theta * d_ln_density_dr_table(int_radius)
+ else
+ ! Cartesian components of the gravitational acceleration
+ ! integrate and multiply by rho / Kappa
+ temp_gxl(ijk,1,1) = sin_theta*cos_phi
+ temp_gyl(ijk,1,1) = sin_theta*sin_phi
+ temp_gzl(ijk,1,1) = cos_theta
+ endif
+ enddo
+#else
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -186,6 +229,7 @@
enddo
enddo
enddo
+#endif
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
@@ -401,8 +445,8 @@
enddo
enddo
+ ! sum contributions from each element to the global mesh and add gravity term
#ifdef FORCE_VECTORIZATION
- ! sum contributions from each element to the global mesh and add gravity term
do ijk=1,NGLLCUBE
sum_terms(ijk,1,1) = - (wgllwgll_yz_3D(ijk,1,1)*newtempx1(ijk,1,1) &
+ wgllwgll_xz_3D(ijk,1,1)*newtempx2(ijk,1,1) &
@@ -426,10 +470,25 @@
enddo
! update rotation term with Euler scheme
if(ROTATION_VAL) then
- do ijk = 1,NGLLCUBE
- A_array_rotation(ijk,1,1,ispec) = A_array_rotation(ijk,1,1,ispec) + source_euler_A(ijk,1,1)
- B_array_rotation(ijk,1,1,ispec) = B_array_rotation(ijk,1,1,ispec) + source_euler_B(ijk,1,1)
- enddo
+ if(USE_LDDRK) then
+ ! use the source saved above
+ do ijk = 1,NGLLCUBE
+ A_array_rotation_lddrk(ijk,1,1,ispec) = ALPHA_LDDRK(istage) * A_array_rotation_lddrk(ijk,1,1,ispec) &
+ + source_euler_A(ijk,1,1)
+ A_array_rotation(ijk,1,1,ispec) = A_array_rotation(ijk,1,1,ispec) &
+ + BETA_LDDRK(istage) * A_array_rotation_lddrk(ijk,1,1,ispec)
+
+ B_array_rotation_lddrk(ijk,1,1,ispec) = ALPHA_LDDRK(istage) * B_array_rotation_lddrk(ijk,1,1,ispec) &
+ + source_euler_B(ijk,1,1)
+ B_array_rotation(ijk,1,1,ispec) = B_array_rotation(ijk,1,1,ispec) &
+ + BETA_LDDRK(istage) * B_array_rotation_lddrk(ijk,1,1,ispec)
+ enddo
+ else
+ do ijk = 1,NGLLCUBE
+ A_array_rotation(ijk,1,1,ispec) = A_array_rotation(ijk,1,1,ispec) + source_euler_A(ijk,1,1)
+ B_array_rotation(ijk,1,1,ispec) = B_array_rotation(ijk,1,1,ispec) + source_euler_B(ijk,1,1)
+ enddo
+ endif
endif
#else
! sum contributions from each element to the global mesh and add gravity term
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -327,13 +327,13 @@
if(.NOT. GPU_MODE) then
! on CPU
! crust/mantle region
- call it_multiply_accel_elastic(NGLOB_CRUST_MANTLE,veloc_crust_mantle,accel_crust_mantle, &
- two_omega_earth, &
- rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle)
+ call multiply_accel_elastic(NGLOB_CRUST_MANTLE,veloc_crust_mantle,accel_crust_mantle, &
+ two_omega_earth, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle)
! inner core region
- call it_multiply_accel_elastic(NGLOB_INNER_CORE,veloc_inner_core,accel_inner_core, &
- two_omega_earth, &
- rmassx_inner_core,rmassy_inner_core,rmassz_inner_core)
+ call multiply_accel_elastic(NGLOB_INNER_CORE,veloc_inner_core,accel_inner_core, &
+ two_omega_earth, &
+ rmassx_inner_core,rmassy_inner_core,rmassz_inner_core)
else
! on GPU
! includes FORWARD_OR_ADJOINT == 1
@@ -349,11 +349,11 @@
rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load,NGLOB_XY_CM, &
+ updated_dof_ocean_load, &
NSPEC2D_TOP(IREGION_CRUST_MANTLE) )
else
! on GPU
- call compute_coupling_ocean_cuda(Mesh_pointer,NCHUNKS_VAL,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK, &
+ call compute_coupling_ocean_cuda(Mesh_pointer,NCHUNKS_VAL,EXACT_MASS_MATRIX_FOR_ROTATION, &
1) ! <- 1 == forward arrays
endif
endif
@@ -706,13 +706,13 @@
! on CPU
! adjoint / kernel runs
! crust/mantle region
- call it_multiply_accel_elastic(NGLOB_CRUST_MANTLE_ADJOINT,b_veloc_crust_mantle,b_accel_crust_mantle, &
- b_two_omega_earth, &
- b_rmassx_crust_mantle,b_rmassy_crust_mantle,b_rmassz_crust_mantle)
+ call multiply_accel_elastic(NGLOB_CRUST_MANTLE_ADJOINT,b_veloc_crust_mantle,b_accel_crust_mantle, &
+ b_two_omega_earth, &
+ b_rmassx_crust_mantle,b_rmassy_crust_mantle,b_rmassz_crust_mantle)
! inner core region
- call it_multiply_accel_elastic(NGLOB_INNER_CORE,b_veloc_inner_core,b_accel_inner_core, &
- b_two_omega_earth, &
- b_rmassx_inner_core,b_rmassy_inner_core,b_rmassz_inner_core)
+ call multiply_accel_elastic(NGLOB_INNER_CORE,b_veloc_inner_core,b_accel_inner_core, &
+ b_two_omega_earth, &
+ b_rmassx_inner_core,b_rmassy_inner_core,b_rmassz_inner_core)
else
! on GPU
! includes FORWARD_OR_ADJOINT == 3
@@ -728,11 +728,11 @@
b_rmassx_crust_mantle,b_rmassy_crust_mantle,b_rmassz_crust_mantle, &
rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
- updated_dof_ocean_load,NGLOB_XY_CM, &
+ updated_dof_ocean_load, &
NSPEC2D_TOP(IREGION_CRUST_MANTLE) )
else
! on GPU
- call compute_coupling_ocean_cuda(Mesh_pointer,NCHUNKS_VAL,EXACT_MASS_MATRIX_FOR_ROTATION,USE_LDDRK, &
+ call compute_coupling_ocean_cuda(Mesh_pointer,NCHUNKS_VAL,EXACT_MASS_MATRIX_FOR_ROTATION, &
3) ! <- 3 == backward/reconstructed arrays
endif
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_seismograms.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_seismograms.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -25,131 +25,43 @@
!
!=====================================================================
- subroutine compute_seismograms(nrec_local,nrec,displ_crust_mantle, &
- nu,hxir_store,hetar_store,hgammar_store, &
- scale_displ,ibool_crust_mantle, &
- ispec_selected_rec,number_receiver_global, &
- seismo_current,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
- seismograms)
- use constants_solver
+ subroutine compute_seismograms(nglob,displ,seismo_current,seismograms)
- implicit none
+ use constants_solver,only: &
+ CUSTOM_REAL,SIZE_REAL,ZERO,NGLLX,NGLLY,NGLLZ, &
+ NDIM
- integer nrec_local,nrec
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
- displ_crust_mantle
+ use specfem_par,only: &
+ NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+ nrec,nrec_local, &
+ nu,hxir_store,hetar_store,hgammar_store, &
+ ispec_selected_rec,number_receiver_global, &
+ scale_displ
- double precision, dimension(NDIM,NDIM,nrec) :: nu
+ use specfem_par_crustmantle,only: ibool_crust_mantle
- double precision, dimension(nrec_local,NGLLX) :: hxir_store
- double precision, dimension(nrec_local,NGLLY) :: hetar_store
- double precision, dimension(nrec_local,NGLLZ) :: hgammar_store
-
- double precision scale_displ
-
- integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
-
- integer, dimension(nrec) :: ispec_selected_rec
- integer, dimension(nrec_local) :: number_receiver_global
-
- integer :: seismo_current
- integer :: NTSTEP_BETWEEN_OUTPUT_SEISMOS
-
- real(kind=CUSTOM_REAL), dimension(NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS) :: &
- seismograms
-
- ! local parameters
- double precision :: uxd,uyd,uzd,hlagrange
- integer :: i,j,k,iglob,irec_local,irec
-
- do irec_local = 1,nrec_local
-
- ! get global number of that receiver
- irec = number_receiver_global(irec_local)
-
- ! perform the general interpolation using Lagrange polynomials
- uxd = ZERO
- uyd = ZERO
- uzd = ZERO
-
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
-
- iglob = ibool_crust_mantle(i,j,k,ispec_selected_rec(irec))
-
- hlagrange = hxir_store(irec_local,i)*hetar_store(irec_local,j)*hgammar_store(irec_local,k)
-
- uxd = uxd + dble(displ_crust_mantle(1,iglob))*hlagrange
- uyd = uyd + dble(displ_crust_mantle(2,iglob))*hlagrange
- uzd = uzd + dble(displ_crust_mantle(3,iglob))*hlagrange
-
- enddo
- enddo
- enddo
- ! store North, East and Vertical components
-
- ! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- seismograms(:,irec_local,seismo_current) = sngl(scale_displ*(nu(:,1,irec)*uxd + &
- nu(:,2,irec)*uyd + nu(:,3,irec)*uzd))
- else
- seismograms(:,irec_local,seismo_current) = scale_displ*(nu(:,1,irec)*uxd + &
- nu(:,2,irec)*uyd + nu(:,3,irec)*uzd)
- endif
-
- enddo
-
- end subroutine compute_seismograms
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
- subroutine compute_seismograms_backward(nrec_local,nrec,b_displ_crust_mantle, &
- nu,hxir_store,hetar_store,hgammar_store, &
- scale_displ,ibool_crust_mantle, &
- ispec_selected_rec,number_receiver_global, &
- seismo_current,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
- seismograms)
-
- use constants_solver
-
implicit none
- integer nrec_local,nrec
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE_ADJOINT) :: &
- b_displ_crust_mantle
+ integer,intent(in) :: nglob
+ real(kind=CUSTOM_REAL), dimension(NDIM,nglob),intent(in) :: displ
- double precision, dimension(NDIM,NDIM,nrec) :: nu
+ integer,intent(in) :: seismo_current
- double precision, dimension(nrec_local,NGLLX) :: hxir_store
- double precision, dimension(nrec_local,NGLLY) :: hetar_store
- double precision, dimension(nrec_local,NGLLZ) :: hgammar_store
-
- double precision scale_displ
-
- integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
-
- integer, dimension(nrec) :: ispec_selected_rec
- integer, dimension(nrec_local) :: number_receiver_global
-
- integer :: seismo_current
- integer :: NTSTEP_BETWEEN_OUTPUT_SEISMOS
-
- real(kind=CUSTOM_REAL), dimension(NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS) :: &
+ real(kind=CUSTOM_REAL), dimension(NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),intent(out) :: &
seismograms
! local parameters
double precision :: uxd,uyd,uzd,hlagrange
- integer :: i,j,k,iglob,irec_local,irec
+ integer :: i,j,k,ispec,iglob,irec_local,irec
do irec_local = 1,nrec_local
! get global number of that receiver
irec = number_receiver_global(irec_local)
+ ispec = ispec_selected_rec(irec)
+
! perform the general interpolation using Lagrange polynomials
uxd = ZERO
uyd = ZERO
@@ -159,13 +71,13 @@
do j = 1,NGLLY
do i = 1,NGLLX
- iglob = ibool_crust_mantle(i,j,k,ispec_selected_rec(irec))
+ iglob = ibool_crust_mantle(i,j,k,ispec)
hlagrange = hxir_store(irec_local,i)*hetar_store(irec_local,j)*hgammar_store(irec_local,k)
- uxd = uxd + dble(b_displ_crust_mantle(1,iglob))*hlagrange
- uyd = uyd + dble(b_displ_crust_mantle(2,iglob))*hlagrange
- uzd = uzd + dble(b_displ_crust_mantle(3,iglob))*hlagrange
+ uxd = uxd + dble(displ(1,iglob))*hlagrange
+ uyd = uyd + dble(displ(2,iglob))*hlagrange
+ uzd = uzd + dble(displ(3,iglob))*hlagrange
enddo
enddo
@@ -184,86 +96,68 @@
enddo
- end subroutine compute_seismograms_backward
+ end subroutine compute_seismograms
!
!-------------------------------------------------------------------------------------------------
!
- subroutine compute_seismograms_adjoint(NSOURCES,nrec_local,displ_crust_mantle, &
- eps_trace_over_3_crust_mantle, &
- epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
- epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
- nu_source,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
- hxir_store,hetar_store,hgammar_store, &
- hpxir_store,hpetar_store,hpgammar_store, &
- tshift_cmt,hdur_gaussian,DT,t0,scale_displ, &
- hprime_xx,hprime_yy,hprime_zz, &
- xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
- etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
- gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
- moment_der,sloc_der,stshift_der,shdur_der,&
- NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismograms,deltat, &
- ibool_crust_mantle,ispec_selected_source,number_receiver_global, &
- NSTEP,it,nit_written)
+ subroutine compute_seismograms_adjoint(displ_crust_mantle, &
+ eps_trace_over_3_crust_mantle, &
+ epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
+ epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
+ nit_written, &
+ moment_der,sloc_der,stshift_der,shdur_der, &
+ seismograms)
- use constants_solver
- use specfem_par,only: UNDO_ATTENUATION
+ use constants_solver,only: &
+ CUSTOM_REAL,SIZE_REAL,ZERO,ONE,PI,GRAV,RHOAV,NGLLX,NGLLY,NGLLZ, &
+ NDIM,NGLOB_CRUST_MANTLE,NSPEC_CRUST_MANTLE, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_CRUST_MANTLE_STR_OR_ATT
- implicit none
+ use specfem_par,only: &
+ NSTEP,NTSTEP_BETWEEN_OUTPUT_SEISMOS,UNDO_ATTENUATION, &
+ NSOURCES,nrec_local, &
+ nu_source,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
+ hxir_store,hpxir_store,hetar_store,hpetar_store,hgammar_store,hpgammar_store, &
+ tshift_cmt,hdur_gaussian, &
+ DT,t0,deltat,it, &
+ scale_displ, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ ispec_selected_source,number_receiver_global
- integer NSOURCES,nrec_local
+ use specfem_par_crustmantle,only: ibool_crust_mantle, &
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
+ implicit none
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE),intent(in) :: &
displ_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY),intent(in) :: &
eps_trace_over_3_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT),intent(in) :: &
epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle
- double precision, dimension(NDIM,NDIM,NSOURCES) :: nu_source
- double precision, dimension(NSOURCES) :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
+ integer,intent(in) :: nit_written
- double precision, dimension(nrec_local,NGLLX) :: hxir_store,hpxir_store
- double precision, dimension(nrec_local,NGLLY) :: hetar_store,hpetar_store
- double precision, dimension(nrec_local,NGLLZ) :: hgammar_store,hpgammar_store
+ real(kind=CUSTOM_REAL), dimension(NDIM,NDIM,nrec_local),intent(inout) :: moment_der
+ real(kind=CUSTOM_REAL), dimension(NDIM,nrec_local),intent(inout) :: sloc_der
+ real(kind=CUSTOM_REAL), dimension(nrec_local),intent(inout) :: stshift_der, shdur_der
- double precision, dimension(NSOURCES) :: tshift_cmt,hdur_gaussian
- double precision :: DT,t0
- double precision :: scale_displ, scale_t
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
- real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
- real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: &
- xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
- etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
- gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle
-
- real(kind=CUSTOM_REAL), dimension(NDIM,NDIM,nrec_local) :: moment_der
- real(kind=CUSTOM_REAL), dimension(NDIM,nrec_local) :: sloc_der
- real(kind=CUSTOM_REAL), dimension(nrec_local) :: stshift_der, shdur_der
-
- integer NTSTEP_BETWEEN_OUTPUT_SEISMOS
-
- real(kind=CUSTOM_REAL), dimension(NDIM*NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS) :: &
+ real(kind=CUSTOM_REAL), dimension(NDIM*NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),intent(out) :: &
seismograms
- real(kind=CUSTOM_REAL) :: deltat
- integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
-
- integer,dimension(NSOURCES) :: ispec_selected_source
- integer, dimension(nrec_local) :: number_receiver_global
- integer :: NSTEP,it,nit_written
-
! local parameters
double precision :: uxd,uyd,uzd,hlagrange
double precision :: eps_trace,dxx,dyy,dxy,dxz,dyz
double precision :: eps_loc(NDIM,NDIM), eps_loc_new(NDIM,NDIM)
double precision :: stf
+ double precision :: scale_t
+
real(kind=CUSTOM_REAL) :: displ_s(NDIM,NGLLX,NGLLY,NGLLZ)
real(kind=CUSTOM_REAL) :: eps_s(NDIM,NDIM), eps_m_s, &
eps_m_l_s(NDIM), stf_deltat, Kp_deltat, Hp_deltat
@@ -393,6 +287,7 @@
sloc_der(:,irec_local) = sloc_der(:,irec_local) + eps_m_l_s(:) * stf_deltat
scale_t = ONE/dsqrt(PI*GRAV*RHOAV)
+
Kp_deltat= -1.0d0/sqrt(PI)/hdur_gaussian(irec)*exp(-((dble(NSTEP-it)*DT-t0-tshift_cmt(irec))/hdur_gaussian(irec))**2) &
* deltat * scale_t
Hp_deltat= (dble(NSTEP-it)*DT-t0-tshift_cmt(irec))/hdur_gaussian(irec)*Kp_deltat
@@ -410,58 +305,59 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine compute_seismograms_undoatt(seismo_current,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismograms)
+! unused...
+!
+! subroutine compute_seismograms_undoatt()
+!
+!! re-orders seismogram entries
+!
+! use specfem_par,only: CUSTOM_REAL,NDIM, &
+! NT_DUMP_ATTENUATION,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+! nrec_local,myrank, &
+! seismo_current,seismograms
+!
+! implicit none
+!
+! ! local parameters
+! integer :: i,j,k,irec_local
+! real(kind=CUSTOM_REAL), dimension(3) :: seismograms_temp
+!
+! ! checks if anything to do
+! if( nrec_local == 0 ) return
+!
+! if(mod(NT_DUMP_ATTENUATION,2) == 0)then
+!
+! do irec_local = 1,nrec_local
+! do i = 1,seismo_current/NT_DUMP_ATTENUATION
+! do j = 1,NT_DUMP_ATTENUATION/2
+! do k = 1,NDIM
+! seismograms_temp(k) = seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + j)
+!
+! seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + j) = &
+! seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + (NT_DUMP_ATTENUATION-j+1))
+!
+! seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + (NT_DUMP_ATTENUATION-j+1)) = seismograms_temp(k)
+! enddo
+! enddo
+! enddo
+! enddo
+!
+! else
+!
+! do irec_local = 1,nrec_local
+! do i = 1,seismo_current/NT_DUMP_ATTENUATION
+! do j = 1,(NT_DUMP_ATTENUATION-1)/2
+! do k = 1,NDIM
+! seismograms_temp(k) = seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + j)
+! seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + j) = &
+! seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + (NT_DUMP_ATTENUATION-j+1))
+! seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + (NT_DUMP_ATTENUATION-j+1)) = seismograms_temp(k)
+! enddo
+! enddo
+! enddo
+! enddo
+!
+! endif
+!
+! end subroutine compute_seismograms_undoatt
-! re-orders seismogram entries
-
- use specfem_par,only: CUSTOM_REAL,NDIM,NT_DUMP_ATTENUATION
-
- implicit none
-
- integer :: seismo_current
- integer :: nrec_local
- integer :: NTSTEP_BETWEEN_OUTPUT_SEISMOS
-
- real(kind=CUSTOM_REAL), dimension(NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS) :: &
- seismograms
-
- ! local parameters
- integer :: i,j,k,irec_local
- real(kind=CUSTOM_REAL), dimension(3) :: seismograms_temp
-
- if(mod(NT_DUMP_ATTENUATION,2) == 0)then
-
- do irec_local = 1,nrec_local
- do i = 1,seismo_current/NT_DUMP_ATTENUATION
- do j = 1,NT_DUMP_ATTENUATION/2
- do k = 1,NDIM
- seismograms_temp(k) = seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + j)
-
- seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + j) = &
- seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + (NT_DUMP_ATTENUATION-j+1))
-
- seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + (NT_DUMP_ATTENUATION-j+1)) = seismograms_temp(k)
- enddo
- enddo
- enddo
- enddo
-
- else
-
- do irec_local = 1,nrec_local
- do i = 1,seismo_current/NT_DUMP_ATTENUATION
- do j = 1,(NT_DUMP_ATTENUATION-1)/2
- do k = 1,NDIM
- seismograms_temp(k) = seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + j)
- seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + j) = &
- seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + (NT_DUMP_ATTENUATION-j+1))
- seismograms(k,irec_local,(i-1)*NT_DUMP_ATTENUATION + (NT_DUMP_ATTENUATION-j+1)) = seismograms_temp(k)
- enddo
- enddo
- enddo
- enddo
-
- endif
-
- end subroutine compute_seismograms_undoatt
-
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -183,7 +183,6 @@
! mass matrices
! crust/mantle
- deallocate(rmassz_crust_mantle)
if( (NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) .or. &
(ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION) ) then
deallocate(rmassx_crust_mantle,rmassy_crust_mantle)
@@ -200,11 +199,9 @@
endif
! outer core
- deallocate(rmass_outer_core)
if(SIMULATION_TYPE == 3 ) nullify(b_rmass_outer_core)
! inner core
- deallocate(rmassz_inner_core)
if( ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION ) then
deallocate(rmassx_inner_core,rmassy_inner_core)
else
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -71,13 +71,13 @@
do it = it_begin,it_end
+ ! simulation status output and stability check
+ if( mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin + 4 .or. it == it_end ) then
+ call check_stability()
+ endif
+
do istage = 1, NSTAGE_TIME_SCHEME ! is equal to 1 if Newmark because only one stage then
- ! simulation status output and stability check
- if((mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin + 4 .or. it == it_end) .and. istage == 1) then
- call check_stability()
- endif
-
if(USE_LDDRK)then
! update displacement using runge-kutta time scheme
call update_displacement_lddrk()
@@ -93,10 +93,16 @@
! elastic solver for crust/mantle and inner core
call compute_forces_viscoelastic()
- ! kernel simulations (forward and adjoint wavefields)
- if( SIMULATION_TYPE == 3 ) then
- ! reconstructs forward wavefields based on last stored wavefield data
+ enddo ! end of very big external loop on istage for all the stages of the LDDRK time scheme (only one stage if Newmark)
+ ! kernel simulations (forward and adjoint wavefields)
+ if( SIMULATION_TYPE == 3 ) then
+
+ ! note: we step back in time (using time steps - DT ), i.e. wavefields b_displ_..() are time-reversed here
+
+ ! reconstructs forward wavefields based on last stored wavefield data
+ do istage = 1, NSTAGE_TIME_SCHEME ! is equal to 1 if Newmark because only one stage then
+
if(USE_LDDRK)then
! update displacement using runge-kutta time scheme
call update_displacement_lddrk_backward()
@@ -111,12 +117,9 @@
! elastic solver for crust/mantle and inner core
call compute_forces_viscoelastic_backward()
- endif
- enddo ! end of very big external loop on istage for all the stages of the LDDRK time scheme (only one stage if Newmark)
+ enddo
- ! kernel simulations (forward and adjoint wavefields)
- if( SIMULATION_TYPE == 3 ) then
! restores last time snapshot saved for backward/reconstruction of wavefields
! note: this is done here after the Newmark time scheme, otherwise the indexing for sources
! and adjoint sources will become more complicated
@@ -127,8 +130,9 @@
! adjoint simulations: kernels
call compute_kernels()
- endif
+ endif ! kernel simulations
+
! write the seismograms with time shift
if( nrec_local > 0 .or. ( WRITE_SEISMOGRAMS_BY_MASTER .and. myrank == 0 ) ) then
call write_seismograms()
@@ -154,130 +158,18 @@
!---- end of time iteration loop
!
- call it_print_elapsed_time()
+ call print_elapsed_time()
! Transfer fields from GPU card to host for further analysis
if(GPU_MODE) call it_transfer_from_GPU()
end subroutine iterate_time
-!
-!-------------------------------------------------------------------------------------------------
-!
- subroutine it_multiply_accel_elastic(NGLOB,veloc,accel, &
- two_omega_earth, &
- rmassx,rmassy,rmassz)
-
-! multiplies acceleration with inverse of mass matrices in crust/mantle,solid inner core region
-
- use constants_solver,only: CUSTOM_REAL,NDIM
-
- implicit none
-
- integer :: NGLOB
-
- ! velocity & acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: veloc,accel
-
- real(kind=CUSTOM_REAL) :: two_omega_earth
-
- ! mass matrices
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: rmassx,rmassy,rmassz
-
- ! local parameters
- integer :: i
-
- ! note: mass matrices
- !
- ! in the case of Stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
- ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
- ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
- !
- ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
- ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
-
- ! updates acceleration w/ rotation in elastic region
-
- ! see input call, differs for corrected mass matrices for rmassx,rmassy,rmassz
- do i=1,NGLOB
- accel(1,i) = accel(1,i)*rmassx(i) + two_omega_earth*veloc(2,i)
- accel(2,i) = accel(2,i)*rmassy(i) - two_omega_earth*veloc(1,i)
- accel(3,i) = accel(3,i)*rmassz(i)
- enddo
-
- end subroutine it_multiply_accel_elastic
-
!
!-------------------------------------------------------------------------------------------------
!
- subroutine it_multiply_accel_acoustic(NGLOB,accel,rmass)
-
-! multiplies acceleration with inverse of mass matrix in outer core region
-
- use constants_solver,only: CUSTOM_REAL
-
- implicit none
-
- integer :: NGLOB
-
- ! velocity & acceleration
- ! crust/mantle region
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: accel
-
- ! mass matrices
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: rmass
-
- ! local parameters
- integer :: i
-
- ! note: mass matrices for fluid region has no stacey or rotation correction
- ! it is also the same for forward and backward/reconstructed wavefields
-
- do i=1,NGLOB
- accel(i) = accel(i)*rmass(i)
- enddo
-
- end subroutine it_multiply_accel_acoustic
-
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-
- subroutine it_print_elapsed_time()
-
- use specfem_par
- implicit none
-
- ! local parameters
- integer :: ihours,iminutes,iseconds,int_tCPU
- ! timing
- double precision, external :: wtime
-
- if(myrank == 0) then
- ! elapsed time since beginning of the simulation
- tCPU = wtime() - time_start
-
- int_tCPU = int(tCPU)
- ihours = int_tCPU / 3600
- iminutes = (int_tCPU - 3600*ihours) / 60
- iseconds = int_tCPU - 3600*ihours - 60*iminutes
- write(IMAIN,*) 'Time-Loop Complete. Timing info:'
- write(IMAIN,*) 'Total elapsed time in seconds = ',tCPU
- write(IMAIN,"(' Total elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
- call flush_IMAIN()
- endif
-
- end subroutine it_print_elapsed_time
-
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
subroutine it_transfer_from_GPU()
! transfers fields on GPU back onto CPU
@@ -354,9 +246,12 @@
end subroutine it_transfer_from_GPU
-!=====================================================================
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine it_update_vtkwindow()
use specfem_par
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time_undoatt.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time_undoatt.F90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time_undoatt.F90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -52,28 +52,8 @@
! checks
if( .not. UNDO_ATTENUATION ) return
- ! synchronize all processes to make sure everybody is ready to start time loop
- call sync_all()
-
- if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) 'Starting time iteration loop in undoing attenuation...'
- write(IMAIN,*)
- call flush_IMAIN()
- endif
-
- ! create an empty file to monitor the start of the simulation
- if(myrank == 0) then
- open(unit=IOUT,file=trim(OUTPUT_FILES)//'/starttimeloop_undoatt.txt',status='unknown',action='write')
- write(IOUT,*) 'hello, starting time loop'
- close(IOUT)
- endif
-
+ ! allocates buffers
if( SIMULATION_TYPE == 3 ) then
- ! to switch between simulation type 1 mode and simulation type 3 mode
- ! in exact undoing of attenuation
- undo_att_sim_type_3 = .true.
-
!! DK DK to Daniel, July 2013: in the case of GPU_MODE it will be *crucial* to leave these arrays on the host
!! i.e. on the CPU, in order to be able to use all the (unused) memory of the host for them, since they are
!! (purposely) huge and designed to use almost all the memory available (by carefully optimizing the
@@ -93,6 +73,23 @@
if( ier /= 0 ) call exit_MPI(myrank,'error allocating b_displ_inner_core_store_buffer')
endif
+ ! synchronize all processes to make sure everybody is ready to start time loop
+ call sync_all()
+
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'Starting time iteration loop in undoing attenuation...'
+ write(IMAIN,*)
+ call flush_IMAIN()
+ endif
+
+ ! create an empty file to monitor the start of the simulation
+ if(myrank == 0) then
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/starttimeloop_undoatt.txt',status='unknown',action='write')
+ write(IOUT,*) 'hello, starting time loop'
+ close(IOUT)
+ endif
+
! initialize variables for writing seismograms
seismo_offset = it_begin-1
seismo_current = 0
@@ -104,7 +101,6 @@
! ************* MAIN LOOP OVER THE TIME STEPS *************
! *********************************************************
-
it = 0
do iteration_on_subset = 1, NSTEP / NT_DUMP_ATTENUATION
@@ -127,21 +123,23 @@
seismo_current_temp = seismo_current
endif
- ! forward and adjoint simulations
- if(SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 2) then
+ ! time loop within this iteration subset
+ select case( SIMULATION_TYPE )
+ case( 1, 2 )
+ ! forward and adjoint simulations
! subset loop
do it_of_this_subset = 1, NT_DUMP_ATTENUATION
it = it + 1
+ ! simulation status output and stability check
+ if( mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin + 4 .or. it == it_end ) then
+ call check_stability()
+ endif
+
do istage = 1, NSTAGE_TIME_SCHEME ! is equal to 1 if Newmark because only one stage then
- ! simulation status output and stability check
- if((mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin + 4 .or. it == it_end) .and. istage == 1) then
- call check_stability()
- endif
-
if(USE_LDDRK)then
! update displacement using runge-kutta time scheme
call update_displacement_lddrk()
@@ -180,22 +178,26 @@
enddo ! subset loop
- else if( SIMULATION_TYPE == 3 ) then
+ case( 3 )
! kernel simulations
! reconstructs forward wavefield based on last stored wavefield data
+ !
+ ! note: we step forward in time here, starting from last snapshot.
+ ! the newly computed, reconstructed forward wavefields (b_displ_..) get stored in buffers.
+
! subset loop
do it_of_this_subset = 1, NT_DUMP_ATTENUATION
it = it + 1
+ ! simulation status output and stability check
+ if( mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin + 4 .or. it == it_end ) then
+ call check_stability_backward()
+ endif
+
do istage = 1, NSTAGE_TIME_SCHEME ! is equal to 1 if Newmark because only one stage then
- ! simulation status output and stability check
- if((mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin + 4 .or. it == it_end) .and. istage == 1) then
- call check_stability()
- endif
-
if(USE_LDDRK)then
! update displacement using runge-kutta time scheme
call update_displacement_lddrk_backward()
@@ -213,11 +215,6 @@
enddo ! istage
- ! write the seismograms with time shift
- if( nrec_local > 0 .or. ( WRITE_SEISMOGRAMS_BY_MASTER .and. myrank == 0 ) ) then
- call write_seismograms()
- endif
-
! stores wavefield in buffers
b_displ_crust_mantle_store_buffer(:,:,it_of_this_subset) = b_displ_crust_mantle(:,:)
b_displ_outer_core_store_buffer(:,it_of_this_subset) = b_displ_outer_core(:)
@@ -236,6 +233,7 @@
do it_of_this_subset = 1, NT_DUMP_ATTENUATION
! reads backward/reconstructed wavefield from buffers
+ ! note: uses wavefield at corresponding time (NSTEP - it + 1 ), i.e. we have now time-reversed wavefields
! crust/mantle
do i = 1, NDIM
do j =1,NGLOB_CRUST_MANTLE_ADJOINT
@@ -256,13 +254,13 @@
it = it + 1
+ ! simulation status output and stability check
+ if( mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin + 4 .or. it == it_end ) then
+ call check_stability()
+ endif
+
do istage = 1, NSTAGE_TIME_SCHEME ! is equal to 1 if Newmark because only one stage then
- ! simulation status output and stability check
- if((mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin + 4 .or. it == it_end) .and. istage == 1) then
- call check_stability()
- endif
-
if(USE_LDDRK)then
! update displacement using runge-kutta time scheme
call update_displacement_lddrk()
@@ -291,7 +289,7 @@
enddo ! subset loop
- endif ! SIMULATION_TYPE == 3
+ end select ! SIMULATION_TYPE
enddo ! end of main time loop
@@ -307,7 +305,7 @@
b_displ_inner_core_store_buffer)
endif
- call it_print_elapsed_time()
+ call print_elapsed_time()
! Transfer fields from GPU card to host for further analysis
if(GPU_MODE) call it_transfer_from_GPU()
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -104,7 +104,7 @@
double precision gammax,gammay,gammaz
! timer MPI
- double precision time_start,tCPU
+ double precision :: time_start,tCPU
! use dynamic allocation
double precision, dimension(:), allocatable :: final_distance
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -132,7 +132,7 @@
double precision :: sec
! timer MPI
- double precision time_start,tCPU
+ double precision :: time_start,tCPU
double precision, external :: wtime
character(len=150) :: OUTPUT_FILES
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/multiply_arrays_source.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/multiply_arrays_source.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/multiply_arrays_source.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -25,12 +25,110 @@
!
!=====================================================================
+
! we put these multiplications in a separate routine AND IN A SEPARATE FILE because otherwise
! some compilers try to unroll the six loops above and take forever to compile
! we leave this in a separate file otherwise many compilers perform subroutine inlining when
! two subroutines are in the same file and one calls the other
+
+!-------------------------------------------------------------------------------------------------
+!
+! elastic domains
+!
+!-------------------------------------------------------------------------------------------------
+
+
+ subroutine multiply_accel_elastic(NGLOB,veloc,accel, &
+ two_omega_earth, &
+ rmassx,rmassy,rmassz)
+
+! multiplies acceleration with inverse of mass matrices in crust/mantle,solid inner core region
+
+ use constants_solver,only: CUSTOM_REAL,NDIM
+
+ implicit none
+
+ integer :: NGLOB
+
+ ! velocity & acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: veloc,accel
+
+ real(kind=CUSTOM_REAL) :: two_omega_earth
+
+ ! mass matrices
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: rmassx,rmassy,rmassz
+
+ ! local parameters
+ integer :: i
+
+ ! note: mass matrices
+ !
+ ! in the case of Stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+
+ ! updates acceleration w/ rotation in elastic region
+
+ ! see input call, differs for corrected mass matrices for rmassx,rmassy,rmassz
+ do i=1,NGLOB
+ accel(1,i) = accel(1,i)*rmassx(i) + two_omega_earth*veloc(2,i)
+ accel(2,i) = accel(2,i)*rmassy(i) - two_omega_earth*veloc(1,i)
+ accel(3,i) = accel(3,i)*rmassz(i)
+ enddo
+
+ end subroutine multiply_accel_elastic
+
+
+
+!-------------------------------------------------------------------------------------------------
+!
+! acoustic/fluid domains
+!
+!-------------------------------------------------------------------------------------------------
+
+
+ subroutine multiply_accel_acoustic(NGLOB,accel,rmass)
+
+! multiplies acceleration with inverse of mass matrix in outer core region
+
+ use constants_solver,only: CUSTOM_REAL
+
+ implicit none
+
+ integer :: NGLOB
+
+ ! velocity & acceleration
+ ! crust/mantle region
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: accel
+
+ ! mass matrices
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: rmass
+
+ ! local parameters
+ integer :: i
+
+ ! note: mass matrices for fluid region has no stacey or rotation correction
+ ! it is also the same for forward and backward/reconstructed wavefields
+
+ do i=1,NGLOB
+ accel(i) = accel(i)*rmass(i)
+ enddo
+
+ end subroutine multiply_accel_acoustic
+
+
+
+!-------------------------------------------------------------------------------------------------
+!
+! interpolated source arrays
+!
+!-------------------------------------------------------------------------------------------------
+
subroutine multiply_arrays_source(sourcearrayd,G11,G12,G13,G21,G22,G23, &
G31,G32,G33,hxis,hpxis,hetas,hpetas,hgammas,hpgammas,k,l,m)
@@ -78,7 +176,9 @@
end subroutine multiply_arrays_source
-!================================================================
+!
+!-------------------------------------------------------------------------------------------------
+!
subroutine multiply_arrays_adjoint(sourcearrayd,hxir,hetar,hgammar,adj_src_ud)
@@ -104,3 +204,4 @@
enddo
end subroutine multiply_arrays_adjoint
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -34,6 +34,7 @@
! local parameters
! timing
+ double precision :: tCPU
double precision, external :: wtime
! get MPI starting time
@@ -353,8 +354,8 @@
nibool_interfaces_crust_mantle,ibool_interfaces_crust_mantle,&
my_neighbours_crust_mantle)
- if( ((NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) .or. &
- (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION)) .and. NGLOB_CRUST_MANTLE > 0 ) then
+ if( (NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) .or. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION) ) then
call assemble_MPI_scalar(NPROCTOT_VAL,NGLOB_CRUST_MANTLE, &
rmassx_crust_mantle, &
num_interfaces_crust_mantle,max_nibool_interfaces_cm, &
@@ -369,7 +370,7 @@
endif
if( SIMULATION_TYPE == 3 ) then
- if( (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION ) .and. NGLOB_XY_CM > 0)then
+ if( ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION )then
call assemble_MPI_scalar(NPROCTOT_VAL,NGLOB_XY_CM, &
b_rmassx_crust_mantle, &
num_interfaces_crust_mantle,max_nibool_interfaces_cm, &
@@ -1192,9 +1193,6 @@
if( ier /= 0 ) call exit_MPI(myrank,'error allocating b_epsilondev*** arrays for inner core')
endif
endif
- ! to switch between simulation type 1 mode and simulation type 3 mode
- ! in exact undoing of attenuation
- undo_att_sim_type_3 = .false.
! inner core
eps_trace_over_3_inner_core(:,:,:,:) = init_value
@@ -1861,7 +1859,7 @@
SAVE_FORWARD,ABSORBING_CONDITIONS, &
OCEANS_VAL, &
GRAVITY_VAL, &
- ROTATION_VAL, &
+ ROTATION_VAL,EXACT_MASS_MATRIX_FOR_ROTATION, &
ATTENUATION_VAL,UNDO_ATTENUATION, &
PARTIAL_PHYS_DISPERSION_ONLY,USE_3D_ATTENUATION_ARRAYS, &
COMPUTE_AND_STORE_STRAIN, &
@@ -2117,9 +2115,8 @@
kappavstore_crust_mantle,muvstore_crust_mantle, &
kappahstore_crust_mantle,muhstore_crust_mantle, &
eta_anisostore_crust_mantle, &
- rmassx_crust_mantle, &
- rmassy_crust_mantle, &
- rmassz_crust_mantle, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ b_rmassx_crust_mantle,b_rmassy_crust_mantle, &
ibool_crust_mantle, &
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
ispec_is_tiso_crust_mantle, &
@@ -2171,6 +2168,7 @@
gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
rhostore_inner_core,kappavstore_inner_core,muvstore_inner_core, &
rmassx_inner_core,rmassy_inner_core,rmassz_inner_core, &
+ b_rmassx_inner_core,b_rmassy_inner_core, &
ibool_inner_core, &
xstore_inner_core,ystore_inner_core,zstore_inner_core, &
c11store_inner_core,c12store_inner_core,c13store_inner_core, &
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -40,9 +40,6 @@
integer :: ier
character(len=150) outputname
- ! checks if anything to do
- if( UNDO_ATTENUATION ) return
-
! checks run/checkpoint number
if(NUMBER_OF_RUNS < 1 .or. NUMBER_OF_RUNS > NSTEP) &
stop 'number of restart runs can not be less than 1 or greater than NSTEP'
@@ -61,6 +58,10 @@
it_end = NSTEP
endif
+ ! checks if anything to do
+ ! undoing attenuation doesn't support the following checkpointing
+ if( UNDO_ATTENUATION ) return
+
! read files back from local disk or MT tape system if restart file
if(NUMBER_OF_THIS_RUN > 1) then
if( ADIOS_ENABLED .and. ADIOS_FOR_FORWARD_ARRAYS ) then
@@ -264,7 +265,13 @@
! reads in saved wavefield
write(outputname,'(a,i6.6,a,i6.6,a)') 'proc',myrank,'_save_frame_at',iteration_on_subset_tmp,'.bin'
- open(unit=IIN,file=trim(LOCAL_PATH)//'/'//outputname,status='old',action='read',form='unformatted',iostat=ier)
+
+ ! debug
+ !if(myrank == 0 ) print*,'reading in: ',trim(LOCAL_PATH)//'/'//outputname, NSTEP/NT_DUMP_ATTENUATION,iteration_on_subset
+
+ ! opens corresponding snapshot file for reading
+ open(unit=IIN,file=trim(LOCAL_PATH)//'/'//outputname, &
+ status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error opening file proc***_save_frame_at** for reading')
read(IIN) b_displ_crust_mantle
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.F90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.F90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -36,6 +36,7 @@
! local parameters
! timing
+ double precision :: tCPU
double precision, external :: wtime
! get MPI starting time
@@ -148,6 +149,10 @@
! sets number of top elements for surface movies & noise tomography
NSPEC_TOP = NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+ ! allocates dummy array
+ allocate(dummy_idoubling(NSPEC_CRUST_MANTLE),stat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error allocating dummy idoubling in crust_mantle')
+
! allocates mass matrices in this slice (will be fully assembled in the solver)
!
! in the case of Stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
@@ -157,18 +162,11 @@
! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
- ! allocates dummy array
- allocate(dummy_idoubling(NSPEC_CRUST_MANTLE),stat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error allocating dummy idoubling in crust_mantle')
-
! allocates mass matrices
allocate(rmassx_crust_mantle(NGLOB_XY_CM), &
rmassy_crust_mantle(NGLOB_XY_CM),stat=ier)
if(ier /= 0) stop 'error allocating rmassx, rmassy in crust_mantle'
- allocate(rmassz_crust_mantle(NGLOB_CRUST_MANTLE),stat=ier)
- if(ier /= 0) stop 'error allocating rmassz in crust_mantle'
-
! b_rmassx and b_rmassy will be different to rmassx and rmassy
! needs new arrays
allocate(b_rmassx_crust_mantle(NGLOB_XY_CM), &
@@ -231,8 +229,11 @@
deallocate(dummy_idoubling)
! mass matrix corrections
- if( .not. ((NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) .or. &
- (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION) ) ) then
+ if( (NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) .or. &
+ (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION) ) then
+ ! mass matrices differ for rmassx,rmassy
+ ! continue
+ else
! uses single mass matrix without correction
! frees pointer memory
deallocate(rmassx_crust_mantle,rmassy_crust_mantle)
@@ -246,7 +247,10 @@
! associates mass matrix used for backward/reconstructed wavefields
b_rmassz_crust_mantle => rmassz_crust_mantle
! checks if we can take rmassx and rmassy (only differs for rotation correction)
- if( .not. (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION) ) then
+ if( ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION ) then
+ ! mass matrices differ for b_rmassx,b_rmassy
+ ! continue
+ else
! frees pointer memory
deallocate(b_rmassx_crust_mantle,b_rmassy_crust_mantle)
! re-associates with corresponding rmassx,rmassy
@@ -256,6 +260,7 @@
else
! b_rmassx,b_rmassy not used anymore
deallocate(b_rmassx_crust_mantle,b_rmassy_crust_mantle)
+ nullify(b_rmassz_crust_mantle)
endif
@@ -303,16 +308,7 @@
stat=ier)
if(ier /= 0) stop 'error allocating dummy rmass and dummy ispec/idoubling in outer core'
- ! allocates mass matrices in this slice (will be fully assembled in the solver)
- !
- ! in the case of Stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
- ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
- ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
- !
- ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
- ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
- allocate(rmass_outer_core(NGLOB_OUTER_CORE),stat=ier)
- if(ier /= 0) stop 'error allocating rmass in outer core'
+ ! reads in mesh arrays
if( ADIOS_ENABLED .and. ADIOS_FOR_ARRAYS_SOLVER ) then
call read_arrays_solver_adios(IREGION_OUTER_CORE,myrank, &
@@ -333,7 +329,7 @@
dummy_array,dummy_array,dummy_array, &
dummy_array,dummy_array,dummy_array, &
ibool_outer_core,dummy_idoubling_outer_core,dummy_ispec_is_tiso, &
- dummy_rmass,dummy_rmass,rmass_outer_core,rmass_ocean_load, &
+ dummy_rmass,dummy_rmass,rmass_outer_core,dummy_array, &
READ_KAPPA_MU,READ_TISO, &
dummy_rmass,dummy_rmass)
else
@@ -355,7 +351,7 @@
dummy_array,dummy_array,dummy_array, &
dummy_array,dummy_array,dummy_array, &
ibool_outer_core,dummy_idoubling_outer_core,dummy_ispec_is_tiso, &
- dummy_rmass,dummy_rmass,rmass_outer_core,rmass_ocean_load, &
+ dummy_rmass,dummy_rmass,rmass_outer_core,dummy_array, &
READ_KAPPA_MU,READ_TISO, &
dummy_rmass,dummy_rmass)
endif
@@ -375,6 +371,8 @@
if( SIMULATION_TYPE == 3 ) then
! associates mass matrix used for backward/reconstructed wavefields
b_rmass_outer_core => rmass_outer_core
+ else
+ nullify(b_rmass_outer_core)
endif
end subroutine read_mesh_databases_OC
@@ -428,9 +426,6 @@
rmassy_inner_core(NGLOB_XY_IC),stat=ier)
if(ier /= 0) stop 'error allocating rmassx, rmassy in inner_core'
- allocate(rmassz_inner_core(NGLOB_INNER_CORE),stat=ier)
- if(ier /= 0) stop 'error allocating rmass in inner core'
-
! b_rmassx and b_rmassy maybe different to rmassx,rmassy
allocate(b_rmassx_inner_core(NGLOB_XY_IC), &
b_rmassy_inner_core(NGLOB_XY_IC),stat=ier)
@@ -456,7 +451,7 @@
c44store_inner_core,dummy_array,dummy_array, &
dummy_array,dummy_array,dummy_array, &
ibool_inner_core,idoubling_inner_core,dummy_ispec_is_tiso, &
- rmassx_inner_core,rmassy_inner_core,rmassz_inner_core,rmass_ocean_load, &
+ rmassx_inner_core,rmassy_inner_core,rmassz_inner_core,dummy_array, &
READ_KAPPA_MU,READ_TISO, &
b_rmassx_inner_core,b_rmassy_inner_core)
else
@@ -478,7 +473,7 @@
c44store_inner_core,dummy_array,dummy_array, &
dummy_array,dummy_array,dummy_array, &
ibool_inner_core,idoubling_inner_core,dummy_ispec_is_tiso, &
- rmassx_inner_core,rmassy_inner_core,rmassz_inner_core,rmass_ocean_load, &
+ rmassx_inner_core,rmassy_inner_core,rmassz_inner_core,dummy_array, &
READ_KAPPA_MU,READ_TISO, &
b_rmassx_inner_core,b_rmassy_inner_core)
endif
@@ -490,7 +485,10 @@
call exit_MPI(myrank,'incorrect global numbering: iboolmax does not equal nglob in inner core')
! mass matrix corrections
- if( .not. ( ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION ) ) then
+ if( ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION ) then
+ ! uses corrected mass matrices
+ ! continue
+ else
! uses single mass matrix without correction
! frees pointer memory
deallocate(rmassx_inner_core,rmassy_inner_core)
@@ -504,7 +502,10 @@
! associates mass matrix used for backward/reconstructed wavefields
b_rmassz_inner_core => rmassz_inner_core
! checks if we can take rmassx and rmassy (only differs for rotation correction)
- if( .not. (ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION) ) then
+ if( ROTATION_VAL .and. EXACT_MASS_MATRIX_FOR_ROTATION ) then
+ ! uses corrected mass matrices
+ ! continue
+ else
! frees pointer memory
deallocate(b_rmassx_inner_core,b_rmassy_inner_core)
! re-associates with corresponding rmassx,rmassy
@@ -514,6 +515,12 @@
else
! b_rmassx,b_rmassy not used anymore
deallocate(b_rmassx_inner_core,b_rmassy_inner_core)
+ ! use dummy pointers used for passing as function arguments
+ ! associates mass matrix used for backward/reconstructed wavefields
+ !b_rmassz_inner_core => rmassz_inner_core
+ !b_rmassx_inner_core => rmassz_inner_core
+ !b_rmassy_inner_core => rmassz_inner_core
+ nullify(b_rmassz_inner_core)
endif
end subroutine read_mesh_databases_IC
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_topography_bathymetry.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_topography_bathymetry.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_topography_bathymetry.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -37,6 +37,7 @@
! local parameters
integer :: ier
! timing
+ double precision :: tCPU
double precision, external :: wtime
! get MPI starting time
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/rules.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/rules.mk 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/rules.mk 2013-09-09 11:33:19 UTC (rev 22779)
@@ -43,11 +43,8 @@
$O/get_backazimuth.solver.o \
$O/get_cmt.solver.o \
$O/get_event_info.solver.o \
- $O/multiply_arrays_source.solver.o \
$O/netlib_specfun_erf.solver.o \
$O/recompute_jacobian.solver.o \
- $O/write_output_ASCII.solver.o \
- $O/write_output_SAC.solver.o \
$(EMPTY_MACRO)
# solver objects with statically allocated arrays; dependent upon
@@ -81,6 +78,7 @@
$O/locate_receivers.solverstatic.o \
$O/locate_regular_points.solverstatic.o \
$O/locate_sources.solverstatic.o \
+ $O/multiply_arrays_source.solverstatic.o \
$O/noise_tomography.solverstatic.o \
$O/prepare_timerun.solverstatic.o \
$O/read_arrays_solver.solverstatic.o \
@@ -98,6 +96,8 @@
$O/write_movie_output.solverstatic.o \
$O/write_movie_volume.solverstatic.o \
$O/write_movie_surface.solverstatic.o \
+ $O/write_output_ASCII.solverstatic.o \
+ $O/write_output_SAC.solverstatic.o \
$O/write_seismograms.solverstatic.o \
$(EMPTY_MACRO)
@@ -135,6 +135,7 @@
$O/rthetaphi_xyz.shared.o \
$O/spline_routines.shared.o \
$O/write_c_binary.cc.o \
+ $O/write_VTK_file.shared.o \
$(EMPTY_MACRO)
###
@@ -237,14 +238,18 @@
## cuda 5 version
${E}/xspecfem3D: $(XSPECFEM_OBJECTS)
+ @echo ""
@echo "building xspecfem3D with CUDA 5 support"
+ @echo ""
${NVCCLINK} -o $(cuda_DEVICE_OBJ) $(cuda_OBJECTS)
${FCLINK} -o ${E}/xspecfem3D $(XSPECFEM_OBJECTS) $(cuda_DEVICE_OBJ) $(MPILIBS) $(CUDA_LINK)
else
## cuda 4 version
${E}/xspecfem3D: $(XSPECFEM_OBJECTS)
+ @echo ""
@echo "building xspecfem3D with CUDA 4 support"
+ @echo ""
${FCLINK} -o ${E}/xspecfem3D $(XSPECFEM_OBJECTS) $(MPILIBS) $(CUDA_LINK)
endif
@@ -252,7 +257,9 @@
## non-cuda version
${E}/xspecfem3D: $(XSPECFEM_OBJECTS)
+ @echo ""
@echo "building xspecfem3D without CUDA support"
+ @echo ""
## use MPI here
## DK DK add OpenMP compiler flag here if needed
# ${MPIFCCOMPILE_CHECK} -qsmp=omp -o ${E}/xspecfem3D $(XSPECFEM_OBJECTS) $(MPILIBS)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_forward_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_forward_arrays.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_forward_arrays.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -186,7 +186,12 @@
! saves frame of the forward simulation
write(outputname,'(a,i6.6,a,i6.6,a)') 'proc',myrank,'_save_frame_at',iteration_on_subset_tmp,'.bin'
- open(unit=IOUT,file=trim(LOCAL_PATH)//'/'//outputname,status='unknown',form='unformatted',action='write',iostat=ier)
+
+ ! debug
+ !if(myrank == 0 ) print*,'saving in: ',trim(LOCAL_PATH)//'/'//outputname, NSTEP/NT_DUMP_ATTENUATION
+
+ open(unit=IOUT,file=trim(LOCAL_PATH)//'/'//outputname, &
+ status='unknown',form='unformatted',action='write',iostat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error opening file proc***_save_frame_at** for writing')
write(IOUT) displ_crust_mantle
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -716,19 +716,19 @@
integer :: ier
! define local to global receiver numbering mapping
- ! needs to be allocate for subroutine calls (even if nrec_local == 0)
+ ! needs to be allocated for subroutine calls (even if nrec_local == 0)
allocate(number_receiver_global(nrec_local),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating global receiver numbering')
! allocates receiver interpolators
if (nrec_local > 0) then
- ! allocate Lagrange interpolators for receivers
+ ! allocates Lagrange interpolators for receivers
allocate(hxir_store(nrec_local,NGLLX), &
hetar_store(nrec_local,NGLLY), &
hgammar_store(nrec_local,NGLLZ),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating receiver interpolators')
- ! define and store Lagrange interpolators at all the receivers
+ ! defines and stores Lagrange interpolators at all the receivers
if (SIMULATION_TYPE == 2) then
nadj_hprec_local = nrec_local
else
@@ -750,7 +750,7 @@
hxir_store,hetar_store,hgammar_store, &
nadj_hprec_local,hpxir_store,hpetar_store,hpgammar_store)
- ! allocate seismogram array
+ ! allocates seismogram array
if (SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then
allocate(seismograms(NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
if(ier /= 0) stop 'error while allocating seismograms'
@@ -758,7 +758,7 @@
! adjoint seismograms
allocate(seismograms(NDIM*NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
if(ier /= 0) stop 'error while allocating adjoint seismograms'
- ! allocate Frechet derivatives array
+ ! allocates Frechet derivatives array
allocate(moment_der(NDIM,NDIM,nrec_local),sloc_der(NDIM,nrec_local), &
stshift_der(nrec_local),shdur_der(nrec_local),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating frechet derivatives arrays')
@@ -769,11 +769,11 @@
shdur_der(:) = 0._CUSTOM_REAL
endif
- ! initialize seismograms
+ ! initializes seismograms
seismograms(:,:,:) = 0._CUSTOM_REAL
nit_written = 0
else
- ! allocate dummy array since we need it to pass as argument e.g. in write_seismograms() routine
+ ! allocates dummy array since we need it to pass as argument e.g. in write_seismograms() routine
! note: nrec_local is zero, fortran 90/95 should allow zero-sized array allocation...
allocate(seismograms(NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
if( ier /= 0) stop 'error while allocating zero seismograms'
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -211,7 +211,7 @@
integer :: ichunk ! needed for Stacey boundaries
! time loop timing
- double precision :: time_start,tCPU
+ double precision :: time_start
!-----------------------------------------------------------------
! assembly
@@ -280,9 +280,6 @@
real(kind=CUSTOM_REAL),dimension(N_SLS) :: tau_sigma_CUSTOM_REAL
! undo_attenuation
- ! to switch between simulation type 1 mode and simulation type 3 mode
- ! in exact undoing of attenuation
- logical :: undo_att_sim_type_3
integer :: iteration_on_subset,it_of_this_subset
! serial i/o mesh reading
@@ -343,12 +340,11 @@
!
! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE),target :: rmassz_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: b_rmassz_crust_mantle
real(kind=CUSTOM_REAL), dimension(:), pointer :: rmassx_crust_mantle,rmassy_crust_mantle
real(kind=CUSTOM_REAL), dimension(:), pointer :: b_rmassx_crust_mantle,b_rmassy_crust_mantle
- real(kind=CUSTOM_REAL), dimension(:), allocatable,target :: rmassz_crust_mantle
- real(kind=CUSTOM_REAL), dimension(:), pointer :: b_rmassz_crust_mantle
-
! displacement, velocity, acceleration
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
displ_crust_mantle,veloc_crust_mantle,accel_crust_mantle
@@ -517,9 +513,8 @@
integer, dimension(NSPEC_INNER_CORE) :: idoubling_inner_core
! mass matrix
- real(kind=CUSTOM_REAL), dimension(:), allocatable,target :: rmassz_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE), target :: rmassz_inner_core
real(kind=CUSTOM_REAL), dimension(:), pointer :: b_rmassz_inner_core
-
real(kind=CUSTOM_REAL), dimension(:), pointer :: rmassx_inner_core,rmassy_inner_core
real(kind=CUSTOM_REAL), dimension(:), pointer :: b_rmassx_inner_core,b_rmassy_inner_core
@@ -619,7 +614,7 @@
rhostore_outer_core,kappavstore_outer_core
! mass matrix
- real(kind=CUSTOM_REAL), dimension(:), allocatable,target :: rmass_outer_core
+ real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE), target :: rmass_outer_core
real(kind=CUSTOM_REAL), dimension(:), pointer :: b_rmass_outer_core
! velocity potential
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/update_displacement_Newmark.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/update_displacement_Newmark.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/update_displacement_Newmark.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -160,17 +160,21 @@
! Newmark time scheme update
if(FORCE_VECTORIZATION_VAL) then
+
do i=1,NGLOB * NDIM
displ(i,1) = displ(i,1) + deltat * veloc(i,1) + deltatsqover2 * accel(i,1)
veloc(i,1) = veloc(i,1) + deltatover2 * accel(i,1)
accel(i,1) = 0._CUSTOM_REAL
enddo
+
else
+
do i=1,NGLOB
displ(:,i) = displ(:,i) + deltat * veloc(:,i) + deltatsqover2 * accel(:,i)
veloc(:,i) = veloc(:,i) + deltatover2 * accel(:,i)
accel(:,i) = 0._CUSTOM_REAL
enddo
+
endif
end subroutine update_displ_elastic
@@ -393,7 +397,7 @@
if(FORCE_VECTORIZATION_VAL) then
do i=1,NGLOB_CM * NDIM
- veloc_crust_mantle(i,i) = veloc_crust_mantle(i,i) + deltatover2*accel_crust_mantle(i,i)
+ veloc_crust_mantle(i,1) = veloc_crust_mantle(i,1) + deltatover2*accel_crust_mantle(i,1)
enddo
else
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -46,11 +46,12 @@
iNIT = 1
endif
- ipoints_3dmovie=0
+ ! initializes movie points
+ ipoints_3dmovie = 0
num_ibool_3dmovie(:) = -99
ispecel_3dmovie = 0
- mask_ibool(:)=.false.
- mask_3dmovie(:,:,:,:)=.false.
+ mask_ibool(:) = .false.
+ mask_3dmovie(:,:,:,:) = .false.
! create name of database
write(prname,'(a,i6.6,a)') trim(LOCAL_TMP_PATH)//'/'//'proc',myrank,'_'
@@ -72,23 +73,28 @@
( (phival < MOVIE_EAST .and. phival > MOVIE_WEST) .or. &
( (MOVIE_EAST < MOVIE_WEST) .and. (phival >MOVIE_EAST .or. phival < MOVIE_WEST) ) ) ) then
ispecel_3dmovie=ispecel_3dmovie+1
+
do k=1,NGLLZ,iNIT
do j=1,NGLLY,iNIT
do i=1,NGLLX,iNIT
iglob = ibool_crust_mantle(i,j,k,ispec)
+
if(.not. mask_ibool(iglob)) then
ipoints_3dmovie = ipoints_3dmovie + 1
+
mask_ibool(iglob)=.true.
mask_3dmovie(i,j,k,ispec)=.true.
num_ibool_3dmovie(iglob)= ipoints_3dmovie
endif
+
enddo !i
enddo !j
enddo !k
+
endif !in region
enddo !ispec
- npoints_3dmovie=ipoints_3dmovie
- nspecel_3dmovie=ispecel_3dmovie
+ npoints_3dmovie = ipoints_3dmovie
+ nspecel_3dmovie = ispecel_3dmovie
write(IOUT,*) npoints_3dmovie, nspecel_3dmovie
close(IOUT)
@@ -163,42 +169,55 @@
if(NDIM /= 3) stop 'movie volume output requires NDIM = 3'
+ ! stepping
if(MOVIE_COARSE) then
iNIT = NGLLX-1
else
iNIT = 1
endif
- ipoints_3dmovie=0
- do ispec=1,NSPEC_CRUST_MANTLE
- do k=1,NGLLZ,iNIT
- do j=1,NGLLY,iNIT
- do i=1,NGLLX,iNIT
- if(mask_3dmovie(i,j,k,ispec)) then
- ipoints_3dmovie=ipoints_3dmovie+1
- iglob= ibool_crust_mantle(i,j,k,ispec)
- rval = xstore_crust_mantle(iglob)
+ ! loops over all elements
+ ipoints_3dmovie = 0
+ do ispec = 1,NSPEC_CRUST_MANTLE
+
+ do k = 1,NGLLZ,iNIT
+ do j = 1,NGLLY,iNIT
+ do i = 1,NGLLX,iNIT
+
+ ! only store points once
+ if( mask_3dmovie(i,j,k,ispec) ) then
+ ! point increment
+ ipoints_3dmovie = ipoints_3dmovie + 1
+
+ ! gets point position
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+
+ rval = xstore_crust_mantle(iglob)
thetaval = ystore_crust_mantle(iglob)
- phival = zstore_crust_mantle(iglob)
+ phival = zstore_crust_mantle(iglob)
+
!x,y,z store have been converted to r theta phi already, need to revert back for xyz output
call rthetaphi_2_xyz(xval,yval,zval,rval,thetaval,phival)
- store_val3D_x(ipoints_3dmovie)=xval
- store_val3D_y(ipoints_3dmovie)=yval
- store_val3D_z(ipoints_3dmovie)=zval
- store_val3D_mu(ipoints_3dmovie)=muvstore_crust_mantle_3dmovie(i,j,k,ispec)
+
+ store_val3D_x(ipoints_3dmovie) = xval
+ store_val3D_y(ipoints_3dmovie) = yval
+ store_val3D_z(ipoints_3dmovie) = zval
+ store_val3D_mu(ipoints_3dmovie) = muvstore_crust_mantle_3dmovie(i,j,k,ispec)
+
st = sin(thetaval)
ct = cos(thetaval)
sp = sin(phival)
cp = cos(phival)
- nu_3dmovie(1,1,ipoints_3dmovie)=-ct*cp
- nu_3dmovie(1,2,ipoints_3dmovie)=-ct*sp
- nu_3dmovie(1,3,ipoints_3dmovie)=st
- nu_3dmovie(2,1,ipoints_3dmovie)=-sp
- nu_3dmovie(2,2,ipoints_3dmovie)=cp
- nu_3dmovie(2,3,ipoints_3dmovie)=0.d0
- nu_3dmovie(3,1,ipoints_3dmovie)=st*cp
- nu_3dmovie(3,2,ipoints_3dmovie)=st*sp
- nu_3dmovie(3,3,ipoints_3dmovie)=ct
+
+ nu_3dmovie(1,1,ipoints_3dmovie) = -ct*cp
+ nu_3dmovie(1,2,ipoints_3dmovie) = -ct*sp
+ nu_3dmovie(1,3,ipoints_3dmovie) = st
+ nu_3dmovie(2,1,ipoints_3dmovie) = -sp
+ nu_3dmovie(2,2,ipoints_3dmovie) = cp
+ nu_3dmovie(2,3,ipoints_3dmovie) = 0.d0
+ nu_3dmovie(3,1,ipoints_3dmovie) = st*cp
+ nu_3dmovie(3,2,ipoints_3dmovie) = st*sp
+ nu_3dmovie(3,3,ipoints_3dmovie) = ct
endif !mask_3dmovie
enddo !i
enddo !j
@@ -208,29 +227,30 @@
write(prname,'(a,i6.6,a)') trim(LOCAL_TMP_PATH)//'/'//'proc',myrank,'_'
open(unit=IOUT,file=trim(prname)//'movie3D_x.bin',status='unknown',form='unformatted')
- if(npoints_3dmovie>0) then
+ if( npoints_3dmovie > 0 ) then
write(IOUT) store_val3D_x(1:npoints_3dmovie)
endif
close(IOUT)
open(unit=IOUT,file=trim(prname)//'movie3D_y.bin',status='unknown',form='unformatted')
- if(npoints_3dmovie>0) then
+ if( npoints_3dmovie > 0 ) then
write(IOUT) store_val3D_y(1:npoints_3dmovie)
endif
close(IOUT)
open(unit=IOUT,file=trim(prname)//'movie3D_z.bin',status='unknown',form='unformatted')
- if(npoints_3dmovie>0) then
+ if( npoints_3dmovie > 0 ) then
write(IOUT) store_val3D_z(1:npoints_3dmovie)
endif
close(IOUT)
open(unit=IOUT,file=trim(prname)//'ascii_output.txt',status='unknown')
- if(npoints_3dmovie>0) then
+ if( npoints_3dmovie > 0 ) then
do i=1,npoints_3dmovie
write(IOUT,*) store_val3D_x(i),store_val3D_y(i),store_val3D_z(i),store_val3D_mu(i)
enddo
endif
close(IOUT)
+
open(unit=IOUT,file=trim(prname)//'movie3D_elements.bin',status='unknown',form='unformatted')
ispecele=0
! open(unit=IOUT,file=trim(prname)//'movie3D_elements.txt',status='unknown')
@@ -263,6 +283,7 @@
n7 = num_ibool_3dmovie(iglob7)-1
n8 = num_ibool_3dmovie(iglob8)-1
write(IOUT) n1,n2,n3,n4,n5,n6,n7,n8
+ ! text output
! write(57,*) n1,n2,n3,n4,n5,n6,n7,n8
! endif !mask3dmovie
enddo !i
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_output_ASCII.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_output_ASCII.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_output_ASCII.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -25,11 +25,7 @@
!
!=====================================================================
- subroutine write_output_ASCII(seismogram_tmp, &
- DT,hdur,OUTPUT_FILES, &
- NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current, &
- SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,myrank, &
- iorientation,sisname,sisname_big_file)
+ subroutine write_output_ASCII(seismogram_tmp,iorientation,sisname,sisname_big_file)
! save seismograms in text format with no subsampling.
! Because we do not subsample the output, this can result in large files
@@ -37,36 +33,35 @@
! here would result in a loss of accuracy when one later convolves
! the results with the source time function
- use constants
+ use constants,only: CUSTOM_REAL,SIZE_REAL,IOUT,C_LDDRK
+ use specfem_par,only: &
+ DT,t0,NSTEP, &
+ seismo_offset,seismo_current, &
+ NTSTEP_BETWEEN_OUTPUT_SEISMOS,OUTPUT_FILES,SIMULATION_TYPE, &
+ SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE, &
+ myrank
+
implicit none
- integer :: seismo_offset, seismo_current, NTSTEP_BETWEEN_OUTPUT_SEISMOS
-
real(kind=CUSTOM_REAL), dimension(5,NTSTEP_BETWEEN_OUTPUT_SEISMOS) :: seismogram_tmp
- integer myrank
- double precision hdur,DT
+ integer :: iorientation
+ character(len=256) :: sisname,sisname_big_file
- integer iorientation
+ ! local parameters
+ integer :: it
+ integer :: ier,isample
+ double precision :: value
+ double precision :: time
+ character(len=256) :: sisname_2
- character(len=256) sisname,sisname_big_file
- character(len=150) OUTPUT_FILES
+ ! add .ascii extension to seismogram file name for ASCII seismograms
+ write(sisname_2,"('/',a,'.ascii')") trim(sisname)
! save all seismograms in one large combined file instead of one file per seismogram
! to avoid overloading shared non-local file systems such as GPFS for instance
- logical SAVE_ALL_SEISMOS_IN_ONE_FILE
- logical USE_BINARY_FOR_LARGE_FILE
- ! local parameters
- integer ier,isample
- character(len=256) sisname_2
- double precision value
-
-
- ! add .ascii extension to seismogram file name for ASCII seismograms
- write(sisname_2,"('/',a,'.ascii')") trim(sisname)
-
! create one large file instead of one small file per station to avoid file system overload
if(SAVE_ALL_SEISMOS_IN_ONE_FILE) then
if(USE_BINARY_FOR_LARGE_FILE) then
@@ -87,21 +82,34 @@
! subtract half duration of the source to make sure travel time is correct
do isample = 1,seismo_current
+
+ ! seismogram value
value = dble(seismogram_tmp(iorientation,isample))
+ ! current time increment
+ it = seismo_offset + isample
+
+ ! current time
+ if( SIMULATION_TYPE == 3 ) then
+ time = dble(NSTEP-it)*DT - t0
+ else
+ time = dble(it-1)*DT - t0
+ endif
+
+ ! writes out to file
if(SAVE_ALL_SEISMOS_IN_ONE_FILE .and. USE_BINARY_FOR_LARGE_FILE) then
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- write(IOUT) sngl(dble(seismo_offset+isample-1)*DT - hdur),sngl(value)
+ write(IOUT) sngl(time),sngl(value)
else
- write(IOUT) dble(seismo_offset+isample-1)*DT - hdur,value
+ write(IOUT) time,value
endif
else
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
- write(IOUT,*) sngl(dble(seismo_offset+isample-1)*DT - hdur),' ',sngl(value)
+ write(IOUT,*) sngl(time),' ',sngl(value)
else
- write(IOUT,*) dble(seismo_offset+isample-1)*DT - hdur,' ',value
+ write(IOUT,*) time,' ',value
endif
endif
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_output_SAC.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_output_SAC.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_output_SAC.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -25,60 +25,49 @@
!
!=====================================================================
- subroutine write_output_SAC(seismogram_tmp,irec, &
- station_name,network_name,stlat,stlon,stele,stbur,nrec, &
- ANGULAR_WIDTH_XI_IN_DEGREES,NEX_XI,DT,hdur,it_end, &
- yr,jda,ho,mi,sec,tshift_cmt,t_shift,&
- elat,elon,depth,event_name,cmt_lat,cmt_lon,cmt_depth,cmt_hdur, &
- OUTPUT_FILES, &
- OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, MODEL, &
- NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current, &
- iorientation,phi,chn,sisname)
+ subroutine write_output_SAC(seismogram_tmp,irec,iorientation,sisname,chn,phi)
! SAC headers have new format
! by Ebru
use constants
- implicit none
+ use specfem_par,only: &
+ ANGULAR_WIDTH_XI_IN_DEGREES,NEX_XI, &
+ myrank,nrec, &
+ station_name,network_name,stlat,stlon,stele,stbur, &
+ DT,t0, &
+ seismo_offset,seismo_current,it_end, &
+ OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
+ NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+ MODEL,OUTPUT_FILES
- integer nrec,it_end
+ use specfem_par,only: &
+ yr=>yr_SAC,jda=>jda_SAC,ho=>ho_SAC,mi=>mi_SAC,sec=>sec_SAC, &
+ tshift_cmt=>t_cmt_SAC,t_shift=>t_shift_SAC, &
+ elat=>elat_SAC,elon=>elon_SAC,depth=>depth_SAC, &
+ event_name=>event_name_SAC,cmt_lat=>cmt_lat_SAC,cmt_lon=>cmt_lon_SAC,&
+ cmt_depth=>cmt_depth_SAC,cmt_hdur=>cmt_hdur_SAC
- integer :: seismo_offset, seismo_current, NTSTEP_BETWEEN_OUTPUT_SEISMOS
+ implicit none
+
real(kind=CUSTOM_REAL), dimension(5,NTSTEP_BETWEEN_OUTPUT_SEISMOS) :: seismogram_tmp
- integer NEX_XI
- double precision ANGULAR_WIDTH_XI_IN_DEGREES
+ integer :: irec
+ integer :: iorientation
- double precision hdur,DT
+ character(len=256) :: sisname
+ character(len=4) :: chn
- character(len=MAX_LENGTH_STATION_NAME), dimension(nrec) :: station_name
- character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec) :: network_name
+ double precision :: phi
- integer irec
- integer iorientation
-
- character(len=4) chn
- character(len=256) sisname
- character(len=150) OUTPUT_FILES,MODEL
-
- double precision tshift_cmt,t_shift,elat,elon,depth
- double precision cmt_lat,cmt_lon,cmt_depth,cmt_hdur
- double precision, dimension(nrec) :: stlat,stlon,stele,stbur
- integer yr,jda,ho,mi
- double precision sec
- character(len=20) event_name
-
- ! flags to determine seismogram type
- logical OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY
-
- real(kind=CUSTOM_REAL) phi
- double precision :: phi_dble
-
! local parameters
- integer time_sec,isample
- character(len=256) sisname_2
+ double precision :: btime
+ integer :: time_sec,isample
+ character(len=256) :: sisname_2
+
+ ! SAC header variables
real DELTA
real DEPMIN
real DEPMAX
@@ -122,8 +111,8 @@
real BYSAC
! end SAC header variables
- double precision shortest_period
- double precision value1,value2, value3,value4,value5
+ double precision :: shortest_period
+ double precision :: value1,value2, value3,value4,value5
logical, external :: is_leap_year
!----------------------------------------------------------------
@@ -173,10 +162,14 @@
SCALE_F= 1000000000 ! factor for y-value, set to 10e9, so that values are in nm
ODELTA = undef ! increment from delta
- B = sngl((seismo_offset)*DT-hdur + tshift_cmt) ! [REQUIRED]
+ ! begin time
+ btime = (seismo_offset)*DT - t0 + tshift_cmt
+
+ B = sngl(btime) ! [REQUIRED]
E = BYSAC ! [REQUIRED]
O = 0 !
A = undef !###
+
!station values:
STLA = stlat(irec)
STLO = stlon(irec)
@@ -194,7 +187,7 @@
! by Ebru
! SAC headers will have new format
- USER0 = cmt_hdur !half duration from CMT file if not changed to hdur=0.d0 (point source)
+ USER0 = cmt_hdur !half duration from CMT file if not changed to t0=0.d0 (point source)
! USER1 and USER2 slots are used for the shortest and longest periods at which
! simulations are accurate, respectively.
@@ -216,7 +209,7 @@
!USER0 = elat
!USER1 = elon
!USER2 = depth
- !USER3 = cmt_hdur !half duration from CMT if not changed to hdur=0.d0 (point source)
+ !USER3 = cmt_hdur !half duration from CMT if not changed to t0=0.d0 (point source)
! just to avoid compiler warning
value1 = elat
@@ -253,26 +246,24 @@
CMPAZ = 0.00
CMPINC = 0.00
else if(iorientation == 4) then !R
- phi_dble = phi
- CMPAZ = sngl(modulo(phi_dble,360.d0)) ! phi is calculated above (see call distaz())
+ CMPAZ = sngl(modulo(phi,360.d0)) ! phi is calculated above (see call distaz())
CMPINC =90.00
else if(iorientation == 5) then !T
- phi_dble = phi
- CMPAZ = sngl(modulo(phi_dble+90.d0,360.d0)) ! phi is calculated above (see call distaz())
+ CMPAZ = sngl(modulo(phi+90.d0,360.d0)) ! phi is calculated above (see call distaz())
CMPINC =90.00
endif
!----------------end format G15.7--------
! date and time:
- NZYEAR =yr
- NZJDAY =jda
- NZHOUR =ho
- NZMIN =mi
+ NZYEAR = yr
+ NZJDAY = jda
+ NZHOUR = ho
+ NZMIN = mi
! adds time-shift to get the CMT time in the headers as origin time of events
! by Ebru
- NZSEC =int(sec+t_shift)
- NZMSEC =int((sec+t_shift-int(sec+t_shift))*1000)
+ NZSEC = int(sec+t_shift)
+ NZMSEC = int((sec+t_shift-int(sec+t_shift))*1000)
!NZSEC =int(sec)
!NZMSEC =int((sec-int(sec))*1000)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90 2013-09-07 01:03:26 UTC (rev 22778)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90 2013-09-09 11:33:19 UTC (rev 22779)
@@ -54,60 +54,31 @@
! computes traces at interpolated receiver locations
select case( SIMULATION_TYPE )
case( 1 )
- call compute_seismograms(nrec_local,nrec,displ_crust_mantle, &
- nu,hxir_store,hetar_store,hgammar_store, &
- scale_displ,ibool_crust_mantle, &
- ispec_selected_rec,number_receiver_global, &
- seismo_current,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
- seismograms)
-
+ call compute_seismograms(NGLOB_CRUST_MANTLE,displ_crust_mantle, &
+ seismo_current,seismograms)
case( 2 )
- call compute_seismograms_adjoint(NSOURCES,nrec_local,displ_crust_mantle, &
- eps_trace_over_3_crust_mantle, &
- epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
- epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
- nu_source,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
- hxir_store,hetar_store,hgammar_store, &
- hpxir_store,hpetar_store,hpgammar_store, &
- tshift_cmt,hdur_gaussian,DT,t0,scale_displ, &
- hprime_xx,hprime_yy,hprime_zz, &
- xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
- etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
- gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
- moment_der,sloc_der,stshift_der,shdur_der, &
- NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismograms,deltat, &
- ibool_crust_mantle,ispec_selected_source,number_receiver_global, &
- NSTEP,it,nit_written)
-
+ call compute_seismograms_adjoint(displ_crust_mantle, &
+ eps_trace_over_3_crust_mantle, &
+ epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
+ epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
+ nit_written, &
+ moment_der,sloc_der,stshift_der,shdur_der, &
+ seismograms)
case( 3 )
- call compute_seismograms_backward(nrec_local,nrec,b_displ_crust_mantle, &
- nu,hxir_store,hetar_store,hgammar_store, &
- scale_displ,ibool_crust_mantle, &
- ispec_selected_rec,number_receiver_global, &
- seismo_current,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
- seismograms)
-
+ call compute_seismograms(NGLOB_CRUST_MANTLE_ADJOINT,b_displ_crust_mantle, &
+ seismo_current,seismograms)
end select
endif ! nrec_local
-
! write the current or final seismograms
if(seismo_current == NTSTEP_BETWEEN_OUTPUT_SEISMOS .or. it == it_end) then
! writes out seismogram files
if (SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then
- ! stores seismograms in right order
- if( UNDO_ATTENUATION .and. SIMULATION_TYPE == 3) then
- call compute_seismograms_undoatt(seismo_current,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismograms)
- endif
+ call write_seismograms_to_file()
- ! writes out seismogram files
- if( .not. undo_att_sim_type_3 ) then
- call write_seismograms_to_file()
- endif
-
! user output
if(myrank==0) then
write(IMAIN,*)
@@ -118,14 +89,11 @@
else if( SIMULATION_TYPE == 2 ) then
if( nrec_local > 0 ) &
call write_adj_seismograms(nit_written)
- nit_written = it
+ nit_written = it
endif
! resets current seismogram position
- if( .not. undo_att_sim_type_3 ) then
- seismo_offset = seismo_offset + seismo_current
- endif
-
+ seismo_offset = seismo_offset + seismo_current
seismo_current = 0
endif
@@ -159,7 +127,6 @@
integer :: iproc,sender,irec_local,irec,ier,receiver
integer :: nrec_local_received
- integer :: nrec_tot_found
integer :: total_seismos,total_seismos_local
integer,dimension(:),allocatable:: islice_num_rec_local
character(len=256) :: sisname
@@ -170,17 +137,14 @@
allocate(one_seismogram(NDIM,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
if(ier /= 0) call exit_mpi(myrank,'error while allocating one temporary seismogram')
- ! check that the sum of the number of receivers in each slice is nrec
- call sum_all_i(nrec_local,nrec_tot_found)
- if(myrank == 0 .and. nrec_tot_found /= nrec) &
- call exit_MPI(myrank,'total number of receivers is incorrect')
-
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
- ! all the processes write their local seismograms themselves
+ ! writes out seismograms
if(.not. WRITE_SEISMOGRAMS_BY_MASTER) then
+ ! all the processes write their local seismograms themselves
+
write_time_begin = wtime()
if(OUTPUT_SEISMOS_ASCII_TEXT .and. SAVE_ALL_SEISMOS_IN_ONE_FILE) then
@@ -234,10 +198,11 @@
call flush_IMAIN()
endif
- ! now only the master process does the writing of seismograms and
- ! collects the data from all other processes
else ! WRITE_SEISMOGRAMS_BY_MASTER
+ ! now only the master process does the writing of seismograms and
+ ! collects the data from all other processes
+
write_time_begin = wtime()
if(myrank == 0) then ! on the master, gather all the seismograms
@@ -285,16 +250,18 @@
! loop on all the slices
do iproc = 0,NPROCTOT_VAL-1
- ! communicates only with processes which contain local receivers
+ ! communicates only with processes which contain local receivers (to minimize mpi chatter)
if( islice_num_rec_local(iproc) == 0 ) cycle
! receive except from proc 0, which is me and therefore I already have this value
sender = iproc
- if(iproc /= 0) then
- call recv_singlei(nrec_local_received,sender,itag)
- if(nrec_local_received < 0) call exit_MPI(myrank,'error while receiving local number of receivers')
- else
+ if(iproc == 0) then
+ ! master is current slice
nrec_local_received = nrec_local
+ else
+ ! receives info from slave processes
+ call recv_singlei(nrec_local_received,sender,itag)
+ if(nrec_local_received <= 0) call exit_MPI(myrank,'error while receiving local number of receivers')
endif
if (nrec_local_received > 0) then
do irec_local = 1,nrec_local_received
@@ -304,6 +271,7 @@
irec = number_receiver_global(irec_local)
one_seismogram(:,:) = seismograms(:,irec_local,:)
else
+ ! receives info from slave processes
call recv_singlei(irec,sender,itag)
if(irec < 1 .or. irec > nrec) call exit_MPI(myrank,'error while receiving global receiver number')
call recv_cr(one_seismogram,NDIM*seismo_current,sender,itag)
@@ -329,8 +297,9 @@
else ! on the nodes, send the seismograms to the master
receiver = 0
- call send_singlei(nrec_local,receiver,itag)
+ ! only sends if this slice contains receiver stations
if (nrec_local > 0) then
+ call send_singlei(nrec_local,receiver,itag)
do irec_local = 1,nrec_local
! get global number of that receiver
irec = number_receiver_global(irec_local)
@@ -368,18 +337,15 @@
ANGULAR_WIDTH_XI_IN_DEGREES,NEX_XI, &
myrank,nrec, &
station_name,network_name,stlat,stlon,stele,stbur, &
- DT,seismo_offset,seismo_current,it_end, &
+ DT,t0, &
+ seismo_offset,seismo_current,it_end, &
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM, &
OUTPUT_SEISMOS_SAC_BINARY,ROTATE_SEISMOGRAMS_RT,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE, &
MODEL,OUTPUT_FILES
use specfem_par,only: &
- hdur=>t0,yr=>yr_SAC,jda=>jda_SAC,ho=>ho_SAC,mi=>mi_SAC,sec=>sec_SAC, &
- tshift_cmt=>t_cmt_SAC,t_shift=>t_shift_SAC, &
- elat=>elat_SAC,elon=>elon_SAC,depth=>depth_SAC, &
- event_name=>event_name_SAC,cmt_lat=>cmt_lat_SAC,cmt_lon=>cmt_lon_SAC,&
- cmt_depth=>cmt_depth_SAC,cmt_hdur=>cmt_hdur_SAC
+ cmt_lat=>cmt_lat_SAC,cmt_lon=>cmt_lon_SAC
implicit none
@@ -389,16 +355,22 @@
! local parameters
real(kind=CUSTOM_REAL), dimension(5,NTSTEP_BETWEEN_OUTPUT_SEISMOS) :: seismogram_tmp
integer :: iorientation,length_station_name,length_network_name
+
character(len=4) :: chn
character(len=256) :: sisname,sisname_big_file
character(len=2) :: bic
+
! variables used for calculation of backazimuth and
! rotation of components if ROTATE_SEISMOGRAMS=.true.
integer :: ior_start,ior_end
double precision :: backaz
- real(kind=CUSTOM_REAL) :: phi,cphi,sphi
+ double precision :: phi
+ real(kind=CUSTOM_REAL) :: cphi,sphi
integer :: isample
+ ! initializes
+ seismogram_tmp(:,:) = 0.0_CUSTOM_REAL
+
! get band code
call band_instrument_code(DT,bic)
@@ -412,24 +384,25 @@
do iorientation = ior_start,ior_end ! BS BS changed according to ROTATE_SEISMOGRAMS_RT
- if(iorientation == 1) then
+ select case( iorientation )
+ case( 1 )
!chn = 'LHN'
chn = bic(1:2)//'N'
- else if(iorientation == 2) then
+ case( 2 )
!chn = 'LHE'
chn = bic(1:2)//'E'
- else if(iorientation == 3) then
+ case( 3 )
!chn = 'LHZ'
chn = bic(1:2)//'Z'
- else if(iorientation == 4) then
+ case( 4 )
!chn = 'LHR'
chn = bic(1:2)//'R'
- else if(iorientation == 5) then
+ case( 5 )
!chn = 'LHT'
chn = bic(1:2)//'T'
- else
+ case default
call exit_MPI(myrank,'incorrect channel value')
- endif
+ end select
if (iorientation == 4 .or. iorientation == 5) then ! LMU BS BS
@@ -439,11 +412,11 @@
call get_backazimuth(cmt_lat,cmt_lon,stlat(irec),stlon(irec),backaz)
phi = backaz
- if (phi>180.) then
- phi = phi-180.
- else if (phi<180.) then
- phi = phi+180.
- else if (phi==180.) then
+ if (phi>180.d0) then
+ phi = phi-180.d0
+ else if (phi<180.d0) then
+ phi = phi+180.d0
+ else if (phi==180.d0) then
phi = backaz
endif
@@ -492,26 +465,12 @@
network_name(irec)(1:length_network_name),chn
! SAC output format
- if (OUTPUT_SEISMOS_SAC_ALPHANUM .or. OUTPUT_SEISMOS_SAC_BINARY) then
- call write_output_SAC(seismogram_tmp,irec, &
- station_name,network_name,stlat,stlon,stele,stbur,nrec, &
- ANGULAR_WIDTH_XI_IN_DEGREES,NEX_XI,DT,hdur,it_end, &
- yr,jda,ho,mi,sec,tshift_cmt,t_shift,&
- elat,elon,depth,event_name,cmt_lat,cmt_lon,cmt_depth,cmt_hdur, &
- OUTPUT_FILES, &
- OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY,MODEL, &
- NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current, &
- iorientation,phi,chn,sisname)
- endif ! OUTPUT_SEISMOS_SAC_ALPHANUM .or. OUTPUT_SEISMOS_SAC_BINARY
+ if( OUTPUT_SEISMOS_SAC_ALPHANUM .or. OUTPUT_SEISMOS_SAC_BINARY ) &
+ call write_output_SAC(seismogram_tmp,irec,iorientation,sisname,chn,phi)
! ASCII output format
- if(OUTPUT_SEISMOS_ASCII_TEXT) then
- call write_output_ASCII(seismogram_tmp, &
- DT,hdur,OUTPUT_FILES, &
- NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current, &
- SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,myrank, &
- iorientation,sisname,sisname_big_file)
- endif ! OUTPUT_SEISMOS_ASCII_TEXT
+ if(OUTPUT_SEISMOS_ASCII_TEXT) &
+ call write_output_ASCII(seismogram_tmp,iorientation,sisname,sisname_big_file)
enddo ! do iorientation
@@ -620,16 +579,18 @@
!
subroutine band_instrument_code(DT,bic)
- ! This subroutine is to choose the appropriate band and instrument codes for channel names of seismograms
- ! based on the IRIS convention (first two letters of channel codes which were LH(Z/E/N) previously).
- ! For consistency with observed data, we now use the IRIS convention for band codes (first letter in channel codes)of
- ! SEM seismograms governed by their sampling rate.
- ! Instrument code (second letter in channel codes) is fixed to "X" which is assigned by IRIS for synthetic seismograms.
- ! See the manual for further explanations!
- ! Ebru, November 2010
+
+! This subroutine is to choose the appropriate band and instrument codes for channel names of seismograms
+! based on the IRIS convention (first two letters of channel codes which were LH(Z/E/N) previously).
+! For consistency with observed data, we now use the IRIS convention for band codes (first letter in channel codes)of
+! SEM seismograms governed by their sampling rate.
+! Instrument code (second letter in channel codes) is fixed to "X" which is assigned by IRIS for synthetic seismograms.
+! See the manual for further explanations!
+! Ebru, November 2010
+
implicit none
- double precision DT
- character(len=2) bic
+ double precision :: DT
+ character(len=2) :: bic
if (DT >= 1.0d0) bic = 'LX'
if (DT < 1.0d0 .and. DT > 0.1d0) bic = 'MX'
More information about the CIG-COMMITS
mailing list