[cig-commits] r22968 - in seismo/3D/SPECFEM3D_GLOBE/trunk: . setup src/cuda src/meshfem3D src/shared src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Tue Oct 22 07:57:01 PDT 2013
Author: danielpeter
Date: 2013-10-22 07:57:01 -0700 (Tue, 22 Oct 2013)
New Revision: 22968
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_method_stubs.c
Removed:
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_method_stubs.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/configure.ac
seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/compute_forces_crust_mantle_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/mesh_constants_cuda.h
seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/prepare_mesh_constants_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh_adios.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_absorb_adios.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_aniso_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_atten3D_QRFSI12.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_eucrust.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_gapp2.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp3d.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_ppm.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_sea99_s.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/rules.mk
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver_adios.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/memory_eval.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/parallel.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/param_reader.c
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/rules.mk
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_coupling.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element_att_memory.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element_strain.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_noDev.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/iterate_time_undoatt.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver_adios.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D_par.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_output.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90
Log:
bug fix for receivers outside of chunk; add vectorization to element routines; rewrites some routines to compile with ifort and -warn all without warnings/errors
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/configure.ac
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/configure.ac 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/configure.ac 2013-10-22 14:57:01 UTC (rev 22968)
@@ -83,9 +83,9 @@
AC_ARG_ENABLE([vectorization],
[AC_HELP_STRING([--enable-vectorization],
- [build FORCE_VECTORIZATION enabled version @<:@default=yes@:>@])],
+ [build FORCE_VECTORIZATION enabled version @<:@default=no@:>@])],
[want_vec="$enableval"],
- [want_vec=yes])
+ [want_vec=no])
AM_CONDITIONAL([COND_VECTORIZATION], [test x"$want_vec" != xno])
############################################################
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2013-10-22 14:57:01 UTC (rev 22968)
@@ -253,6 +253,7 @@
! postprocess the colors to balance them if Droux (1993) algorithm is not used
logical, parameter :: BALANCE_COLORS_SIMPLE_ALGO = .true.
+
!!-----------------------------------------------------------
!!
!! ADIOS Related values
@@ -278,6 +279,7 @@
real, parameter :: KL_REG_MIN_LAT = -90.0
real, parameter :: KL_REG_MAX_LAT = +90.0
+
!!-----------------------------------------------------------
!!
!! new version compatibility
@@ -292,7 +294,6 @@
logical, parameter :: USE_VERSION_5_1_5 = .false.
-
!!-----------------------------------------------------------
!!
!! time stamp information
@@ -311,6 +312,18 @@
!!-----------------------------------------------------------
!!
+!! movie outputs
+!!
+!!-----------------------------------------------------------
+
+! runs external bash script after storing movie files
+! (useful for postprocessing during simulation time)
+ logical, parameter :: RUN_EXTERNAL_MOVIE_SCRIPT = .false.
+ character(len=256) :: MOVIE_SCRIPT_NAME = "./tar_movie_files.sh"
+
+
+!!-----------------------------------------------------------
+!!
!! debugging flags
!!
!!-----------------------------------------------------------
@@ -660,7 +673,11 @@
integer, parameter :: NSPEC_DOUBLING_BASICBRICK = 8
integer, parameter :: NGLOB_DOUBLING_BASICBRICK = 27
-! for LDDRK high-order time scheme
+!!-----------------------------------------------------------
+!!
+!! for LDDRK high-order time scheme
+!!
+!!-----------------------------------------------------------
integer, parameter :: NSTAGE = 6
real(kind=CUSTOM_REAL), parameter, dimension(NSTAGE) :: ALPHA_LDDRK = &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/compute_forces_crust_mantle_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/compute_forces_crust_mantle_cuda.cu 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/compute_forces_crust_mantle_cuda.cu 2013-10-22 14:57:01 UTC (rev 22968)
@@ -39,8 +39,32 @@
#ifdef USE_TEXTURES_FIELDS
-texture<realw, cudaTextureType1D, cudaReadModeElementType> d_displ_cm_tex;
-texture<realw, cudaTextureType1D, cudaReadModeElementType> d_accel_cm_tex;
+//forward
+realw_texture d_displ_cm_tex;
+realw_texture d_accel_cm_tex;
+//backward/reconstructed
+realw_texture d_b_displ_cm_tex;
+realw_texture d_b_accel_cm_tex;
+
+//note: texture variables are implicitly static, and cannot be passed as arguments to cuda kernels;
+// thus, 1) we thus use if-statements (FORWARD_OR_ADJOINT) to determine from which texture to fetch from
+// 2) we use templates
+// since if-statements are a bit slower as the variable is only known at runtime, we use option 2)
+
+// templates definitions
+template<int FORWARD_OR_ADJOINT> __device__ float texfetch_displ_cm(int x);
+template<int FORWARD_OR_ADJOINT> __device__ float texfetch_veloc_cm(int x);
+template<int FORWARD_OR_ADJOINT> __device__ float texfetch_accel_cm(int x);
+
+// templates for texture fetching
+// FORWARD_OR_ADJOINT == 1 <- forward arrays
+template<> __device__ float texfetch_displ_cm<1>(int x) { return tex1Dfetch(d_displ_cm_tex, x); }
+template<> __device__ float texfetch_veloc_cm<1>(int x) { return tex1Dfetch(d_veloc_cm_tex, x); }
+template<> __device__ float texfetch_accel_cm<1>(int x) { return tex1Dfetch(d_accel_cm_tex, x); }
+// FORWARD_OR_ADJOINT == 3 <- backward/reconstructed arrays
+template<> __device__ float texfetch_displ_cm<3>(int x) { return tex1Dfetch(d_b_displ_cm_tex, x); }
+template<> __device__ float texfetch_veloc_cm<3>(int x) { return tex1Dfetch(d_b_veloc_cm_tex, x); }
+template<> __device__ float texfetch_accel_cm<3>(int x) { return tex1Dfetch(d_b_accel_cm_tex, x); }
#endif
#ifdef USE_TEXTURES_CONSTANTS
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/mesh_constants_cuda.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/mesh_constants_cuda.h 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/mesh_constants_cuda.h 2013-10-22 14:57:01 UTC (rev 22968)
@@ -99,31 +99,6 @@
/* ----------------------------------------------------------------------------------------------- */
-// type of "working" variables: see also CUSTOM_REAL
-// double precision temporary variables leads to 10% performance decrease
-// in Kernel_2_impl (not very much..)
-typedef float realw;
-
-
-/* ----------------------------------------------------------------------------------------------- */
-
-// utility functions: defined in check_fields_cuda.cu
-
-/* ----------------------------------------------------------------------------------------------- */
-
-double get_time();
-void get_free_memory(double* free_db, double* used_db, double* total_db);
-void print_CUDA_error_if_any(cudaError_t err, int num);
-void pause_for_debugger(int pause);
-void exit_on_cuda_error(char* kernel_name);
-void exit_on_error(char* info);
-void synchronize_cuda();
-void synchronize_mpi();
-void get_blocks_xy(int num_blocks,int* num_blocks_x,int* num_blocks_y);
-realw get_device_array_maximum_value(realw* array,int size);
-
-/* ----------------------------------------------------------------------------------------------- */
-
// cuda constant arrays
/* ----------------------------------------------------------------------------------------------- */
@@ -157,7 +132,6 @@
// uncomment line below for PREM with oceans
//#define R_EARTH_KM 6368.0f
-
/* ----------------------------------------------------------------------------------------------- */
// (optional) pre-processing directive used in kernels: if defined check that it is also set in src/shared/constants.h:
@@ -243,6 +217,36 @@
/* ----------------------------------------------------------------------------------------------- */
+// type of "working" variables: see also CUSTOM_REAL
+// double precision temporary variables leads to 10% performance decrease
+// in Kernel_2_impl (not very much..)
+typedef float realw;
+
+// textures
+typedef texture<float, cudaTextureType1D, cudaReadModeElementType> realw_texture;
+
+
+/* ----------------------------------------------------------------------------------------------- */
+
+// utility functions: defined in check_fields_cuda.cu
+
+/* ----------------------------------------------------------------------------------------------- */
+
+double get_time();
+void get_free_memory(double* free_db, double* used_db, double* total_db);
+void print_CUDA_error_if_any(cudaError_t err, int num);
+void pause_for_debugger(int pause);
+void exit_on_cuda_error(char* kernel_name);
+void exit_on_error(char* info);
+void synchronize_cuda();
+void synchronize_mpi();
+void get_blocks_xy(int num_blocks,int* num_blocks_x,int* num_blocks_y);
+realw get_device_array_maximum_value(realw* array,int size);
+
+
+
+/* ----------------------------------------------------------------------------------------------- */
+
// mesh pointer wrapper structure
/* ----------------------------------------------------------------------------------------------- */
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/prepare_mesh_constants_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/prepare_mesh_constants_cuda.cu 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/cuda/prepare_mesh_constants_cuda.cu 2013-10-22 14:57:01 UTC (rev 22968)
@@ -40,20 +40,31 @@
#ifdef USE_OLDER_CUDA4_GPU
#else
#ifdef USE_TEXTURES_FIELDS
-extern texture<realw, cudaTextureType1D, cudaReadModeElementType> d_displ_cm_tex;
-extern texture<realw, cudaTextureType1D, cudaReadModeElementType> d_accel_cm_tex;
+ // forward
+ extern realw_texture d_displ_cm_tex;
+ extern realw_texture d_accel_cm_tex;
-extern texture<realw, cudaTextureType1D, cudaReadModeElementType> d_displ_oc_tex;
-extern texture<realw, cudaTextureType1D, cudaReadModeElementType> d_accel_oc_tex;
+ extern realw_texture d_displ_oc_tex;
+ extern realw_texture d_accel_oc_tex;
-extern texture<realw, cudaTextureType1D, cudaReadModeElementType> d_displ_ic_tex;
-extern texture<realw, cudaTextureType1D, cudaReadModeElementType> d_accel_ic_tex;
+ extern realw_texture d_displ_ic_tex;
+ extern realw_texture d_accel_ic_tex;
+
+ // backward/reconstructed
+ extern realw_texture d_b_displ_cm_tex;
+ extern realw_texture d_b_accel_cm_tex;
+
+ extern realw_texture d_b_displ_oc_tex;
+ extern realw_texture d_b_accel_oc_tex;
+
+ extern realw_texture d_b_displ_ic_tex;
+ extern realw_texture d_b_accel_ic_tex;
#endif
#ifdef USE_TEXTURES_CONSTANTS
-extern texture<realw, cudaTextureType1D, cudaReadModeElementType> d_hprime_xx_cm_tex;
-extern texture<realw, cudaTextureType1D, cudaReadModeElementType> d_hprime_xx_oc_tex;
-extern texture<realw, cudaTextureType1D, cudaReadModeElementType> d_hprime_xx_ic_tex;
+ extern realw_texture d_hprime_xx_cm_tex;
+ extern realw_texture d_hprime_xx_oc_tex;
+ extern realw_texture d_hprime_xx_ic_tex;
#endif
#endif
@@ -168,35 +179,29 @@
#ifdef USE_TEXTURES_CONSTANTS
{
#ifdef USE_OLDER_CUDA4_GPU
+ cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<realw>();
const textureReference* d_hprime_xx_cm_tex_ptr;
print_CUDA_error_if_any(cudaGetTextureReference(&d_hprime_xx_cm_tex_ptr, "d_hprime_xx_cm_tex"), 1101);
- cudaChannelFormatDesc channelDesc1 = cudaCreateChannelDesc<realw>();
print_CUDA_error_if_any(cudaBindTexture(0, d_hprime_xx_cm_tex_ptr, mp->d_hprime_xx,
- &channelDesc1, sizeof(realw)*(NGLL2)), 1102);
+ &channelDesc, sizeof(realw)*(NGLL2)), 1102);
const textureReference* d_hprime_xx_oc_tex_ptr;
print_CUDA_error_if_any(cudaGetTextureReference(&d_hprime_xx_oc_tex_ptr, "d_hprime_xx_oc_tex"), 1103);
- cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<realw>();
print_CUDA_error_if_any(cudaBindTexture(0, d_hprime_xx_oc_tex_ptr, mp->d_hprime_xx,
- &channelDesc2, sizeof(realw)*(NGLL2)), 1104);
+ &channelDesc, sizeof(realw)*(NGLL2)), 1104);
const textureReference* d_hprime_xx_ic_tex_ptr;
print_CUDA_error_if_any(cudaGetTextureReference(&d_hprime_xx_ic_tex_ptr, "d_hprime_xx_ic_tex"), 1105);
- cudaChannelFormatDesc channelDesc3 = cudaCreateChannelDesc<realw>();
print_CUDA_error_if_any(cudaBindTexture(0, d_hprime_xx_ic_tex_ptr, mp->d_hprime_xx,
- &channelDesc3, sizeof(realw)*(NGLL2)), 1106);
+ &channelDesc, sizeof(realw)*(NGLL2)), 1106);
#else
- cudaChannelFormatDesc channelDesc1 = cudaCreateChannelDesc<float>();
+ cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
print_CUDA_error_if_any(cudaBindTexture(0, &d_hprime_xx_cm_tex, mp->d_hprime_xx,
- &channelDesc1, sizeof(realw)*(NGLL2)), 1102);
-
- cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<float>();
+ &channelDesc, sizeof(realw)*(NGLL2)), 1102);
print_CUDA_error_if_any(cudaBindTexture(0, &d_hprime_xx_oc_tex, mp->d_hprime_xx,
- &channelDesc2, sizeof(realw)*(NGLL2)), 1104);
-
- cudaChannelFormatDesc channelDesc3 = cudaCreateChannelDesc<float>();
+ &channelDesc, sizeof(realw)*(NGLL2)), 1104);
print_CUDA_error_if_any(cudaBindTexture(0, &d_hprime_xx_ic_tex, mp->d_hprime_xx,
- &channelDesc3, sizeof(realw)*(NGLL2)), 1106);
+ &channelDesc, sizeof(realw)*(NGLL2)), 1106);
#endif
}
#endif
@@ -1344,25 +1349,42 @@
#ifdef USE_TEXTURES_FIELDS
{
#ifdef USE_OLDER_CUDA4_GPU
+ cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<realw>();
const textureReference* d_displ_cm_tex_ref_ptr;
- print_CUDA_error_if_any(cudaGetTextureReference(&mp->d_displ_cm_tex_ref_ptr, "d_displ_cm_tex"), 4021);
- cudaChannelFormatDesc channelDesc1 = cudaCreateChannelDesc<realw>();
- print_CUDA_error_if_any(cudaBindTexture(0, mp->d_displ_cm_tex_ref_ptr, mp->d_displ_crust_mantle,
- &channelDesc1, sizeof(realw)*size), 4021);
+ print_CUDA_error_if_any(cudaGetTextureReference(&d_displ_cm_tex_ref_ptr, "d_displ_cm_tex"), 4021);
+ print_CUDA_error_if_any(cudaBindTexture(0, d_displ_cm_tex_ref_ptr, mp->d_displ_crust_mantle,
+ &channelDesc, sizeof(realw)*size), 4021);
const textureReference* d_accel_cm_tex_ref_ptr;
print_CUDA_error_if_any(cudaGetTextureReference(&d_accel_cm_tex_ref_ptr, "d_accel_cm_tex"), 4023);
- cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<realw>();
print_CUDA_error_if_any(cudaBindTexture(0, d_accel_cm_tex_ref_ptr, mp->d_accel_crust_mantle,
- &channelDesc2, sizeof(realw)*size), 4023);
+ &channelDesc, sizeof(realw)*size), 4023);
+
+ // backward/reconstructed wavefields
+ if( mp->simulation_type == 3 ){
+ const textureReference* d_b_displ_cm_tex_ref_ptr;
+ print_CUDA_error_if_any(cudaGetTextureReference(&d_b_displ_cm_tex_ref_ptr, "d_b_displ_cm_tex"), 4021);
+ print_CUDA_error_if_any(cudaBindTexture(0, d_b_displ_cm_tex_ref_ptr, mp->d_b_displ_crust_mantle,
+ &channelDesc, sizeof(realw)*size), 4021);
+
+ const textureReference* d_b_accel_cm_tex_ref_ptr;
+ print_CUDA_error_if_any(cudaGetTextureReference(&d_b_accel_cm_tex_ref_ptr, "d_b_accel_cm_tex"), 4023);
+ print_CUDA_error_if_any(cudaBindTexture(0, d_b_accel_cm_tex_ref_ptr, mp->d_b_accel_crust_mantle,
+ &channelDesc, sizeof(realw)*size), 4023);
+ }
#else
- cudaChannelFormatDesc channelDesc1 = cudaCreateChannelDesc<float>();
+ cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
print_CUDA_error_if_any(cudaBindTexture(0, &d_displ_cm_tex, mp->d_displ_crust_mantle,
- &channelDesc1, sizeof(realw)*size), 4021);
-
- cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<float>();
+ &channelDesc, sizeof(realw)*size), 4021);
print_CUDA_error_if_any(cudaBindTexture(0, &d_accel_cm_tex, mp->d_accel_crust_mantle,
- &channelDesc2, sizeof(realw)*size), 4023);
+ &channelDesc, sizeof(realw)*size), 4023);
+ // backward/reconstructed wavefields
+ if( mp->simulation_type == 3 ){
+ print_CUDA_error_if_any(cudaBindTexture(0, &d_b_displ_cm_tex, mp->d_b_displ_crust_mantle,
+ &channelDesc, sizeof(realw)*size), 4021);
+ print_CUDA_error_if_any(cudaBindTexture(0, &d_b_accel_cm_tex, mp->d_b_accel_crust_mantle,
+ &channelDesc, sizeof(realw)*size), 4023);
+ }
#endif
}
#endif
@@ -1574,25 +1596,41 @@
#ifdef USE_TEXTURES_FIELDS
{
#ifdef USE_OLDER_CUDA4_GPU
+ cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<realw>();
const textureReference* d_displ_oc_tex_ref_ptr;
- print_CUDA_error_if_any(cudaGetTextureReference(&mp->d_displ_oc_tex_ref_ptr, "d_displ_oc_tex"), 5021);
- cudaChannelFormatDesc channelDesc1 = cudaCreateChannelDesc<realw>();
- print_CUDA_error_if_any(cudaBindTexture(0, mp->d_displ_oc_tex_ref_ptr, mp->d_displ_outer_core,
- &channelDesc1, sizeof(realw)*size_glob), 5021);
+ print_CUDA_error_if_any(cudaGetTextureReference(&d_displ_oc_tex_ref_ptr, "d_displ_oc_tex"), 5021);
+ print_CUDA_error_if_any(cudaBindTexture(0, d_displ_oc_tex_ref_ptr, mp->d_displ_outer_core,
+ &channelDesc, sizeof(realw)*size_glob), 5021);
const textureReference* d_accel_oc_tex_ref_ptr;
print_CUDA_error_if_any(cudaGetTextureReference(&d_accel_oc_tex_ref_ptr, "d_accel_oc_tex"), 5023);
- cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<realw>();
print_CUDA_error_if_any(cudaBindTexture(0, d_accel_oc_tex_ref_ptr, mp->d_accel_outer_core,
- &channelDesc2, sizeof(realw)*size_glob), 5023);
+ &channelDesc, sizeof(realw)*size_glob), 5023);
+ // backward/reconstructed wavefields
+ if( mp->simulation_type == 3 ){
+ const textureReference* d_b_displ_oc_tex_ref_ptr;
+ print_CUDA_error_if_any(cudaGetTextureReference(&d_b_displ_oc_tex_ref_ptr, "d_b_displ_oc_tex"), 5021);
+ print_CUDA_error_if_any(cudaBindTexture(0, d_b_displ_oc_tex_ref_ptr, mp->d_b_displ_outer_core,
+ &channelDesc, sizeof(realw)*size_glob), 5021);
+
+ const textureReference* d_b_accel_oc_tex_ref_ptr;
+ print_CUDA_error_if_any(cudaGetTextureReference(&d_b_accel_oc_tex_ref_ptr, "d_b_accel_oc_tex"), 5023);
+ print_CUDA_error_if_any(cudaBindTexture(0, d_b_accel_oc_tex_ref_ptr, mp->d_b_accel_outer_core,
+ &channelDesc, sizeof(realw)*size_glob), 5023);
+ }
#else
- cudaChannelFormatDesc channelDesc1 = cudaCreateChannelDesc<float>();
+ cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
print_CUDA_error_if_any(cudaBindTexture(0, &d_displ_oc_tex, mp->d_displ_outer_core,
- &channelDesc1, sizeof(realw)*size_glob), 5021);
-
- cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<float>();
+ &channelDesc, sizeof(realw)*size_glob), 5021);
print_CUDA_error_if_any(cudaBindTexture(0, &d_accel_oc_tex, mp->d_accel_outer_core,
- &channelDesc2, sizeof(realw)*size_glob), 5023);
+ &channelDesc, sizeof(realw)*size_glob), 5023);
+ // backward/reconstructed wavefields
+ if( mp->simulation_type == 3 ){
+ print_CUDA_error_if_any(cudaBindTexture(0, &d_b_displ_oc_tex, mp->d_b_displ_outer_core,
+ &channelDesc, sizeof(realw)*size_glob), 5021);
+ print_CUDA_error_if_any(cudaBindTexture(0, &d_b_accel_oc_tex, mp->d_b_accel_outer_core,
+ &channelDesc, sizeof(realw)*size_glob), 5023);
+ }
#endif
}
#endif
@@ -1804,28 +1842,43 @@
#ifdef USE_TEXTURES_FIELDS
{
#ifdef USE_OLDER_CUDA4_GPU
+ cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<realw>();
const textureReference* d_displ_ic_tex_ref_ptr;
- print_CUDA_error_if_any(cudaGetTextureReference(&mp->d_displ_ic_tex_ref_ptr, "d_displ_ic_tex"), 6021);
- cudaChannelFormatDesc channelDesc1 = cudaCreateChannelDesc<realw>();
- print_CUDA_error_if_any(cudaBindTexture(0, mp->d_displ_ic_tex_ref_ptr, mp->d_displ_inner_core,
- &channelDesc1, sizeof(realw)*size), 6021);
+ print_CUDA_error_if_any(cudaGetTextureReference(&d_displ_ic_tex_ref_ptr, "d_displ_ic_tex"), 6021);
+ print_CUDA_error_if_any(cudaBindTexture(0, d_displ_ic_tex_ref_ptr, mp->d_displ_inner_core,
+ &channelDesc, sizeof(realw)*size), 6021);
const textureReference* d_accel_ic_tex_ref_ptr;
print_CUDA_error_if_any(cudaGetTextureReference(&d_accel_ic_tex_ref_ptr, "d_accel_ic_tex"), 6023);
- cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<realw>();
print_CUDA_error_if_any(cudaBindTexture(0, d_accel_ic_tex_ref_ptr, mp->d_accel_inner_core,
- &channelDesc2, sizeof(realw)*size), 6023);
+ &channelDesc, sizeof(realw)*size), 6023);
+ // backward/reconstructed wavefields
+ if( mp->simulation_type == 3 ){
+ const textureReference* d_b_displ_ic_tex_ref_ptr;
+ print_CUDA_error_if_any(cudaGetTextureReference(&d_b_displ_ic_tex_ref_ptr, "d_b_displ_ic_tex"), 6021);
+ print_CUDA_error_if_any(cudaBindTexture(0, d_b_displ_ic_tex_ref_ptr, mp->d_b_displ_inner_core,
+ &channelDesc, sizeof(realw)*size), 6021);
+
+ const textureReference* d_b_accel_ic_tex_ref_ptr;
+ print_CUDA_error_if_any(cudaGetTextureReference(&d_b_accel_ic_tex_ref_ptr, "d_b_accel_ic_tex"), 6023);
+ print_CUDA_error_if_any(cudaBindTexture(0, d_b_accel_ic_tex_ref_ptr, mp->d_b_accel_inner_core,
+ &channelDesc, sizeof(realw)*size), 6023);
+
+ }
#else
- cudaChannelFormatDesc channelDesc1 = cudaCreateChannelDesc<float>();
+ cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
print_CUDA_error_if_any(cudaBindTexture(0, &d_displ_ic_tex, mp->d_displ_inner_core,
- &channelDesc1, sizeof(realw)*size), 6021);
-
- cudaChannelFormatDesc channelDesc2 = cudaCreateChannelDesc<float>();
+ &channelDesc, sizeof(realw)*size), 6021);
print_CUDA_error_if_any(cudaBindTexture(0, &d_accel_ic_tex, mp->d_accel_inner_core,
- &channelDesc2, sizeof(realw)*size), 6023);
+ &channelDesc, sizeof(realw)*size), 6023);
+ // backward/reconstructed wavefields
+ if( mp->simulation_type == 3 ){
+ print_CUDA_error_if_any(cudaBindTexture(0, &d_b_displ_ic_tex, mp->d_b_displ_inner_core,
+ &channelDesc, sizeof(realw)*size), 6021);
+ print_CUDA_error_if_any(cudaBindTexture(0, &d_b_accel_ic_tex, mp->d_b_accel_inner_core,
+ &channelDesc, sizeof(realw)*size), 6023);
+ }
#endif
-
-
}
#endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh_adios.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh_adios.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh_adios.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -56,10 +56,11 @@
! number of spectral elements in each block
integer,intent(in) :: nspec,npointot,iregion_code
- ! local parameters
! arrays used for AVS or DX files
integer, dimension(npointot), intent(inout) :: num_ibool_AVS_DX
logical, dimension(npointot), intent(inout) :: mask_ibool
+
+ ! local parameters
! structures used for ADIOS AVS/DX files
type(avs_dx_global_t) :: avs_dx_global_vars
type(avs_dx_global_faces_t) :: avs_dx_global_faces_vars
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_absorb_adios.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_absorb_adios.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_absorb_adios.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -66,6 +66,7 @@
implicit none
integer :: myrank
+ integer :: iregion
integer :: NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX
integer,dimension(2,NSPEC2DMAX_YMIN_YMAX) :: nimin,nimax
@@ -73,8 +74,9 @@
integer,dimension(2,NSPEC2DMAX_XMIN_XMAX) :: nkmin_xi
integer,dimension(2,NSPEC2DMAX_YMIN_YMAX) :: nkmin_eta
+ ! local parameters
character(len=150) :: outputname, group_name
- integer :: comm, local_dim, iregion
+ integer :: comm, local_dim
integer(kind=8) :: group_size_inc
! ADIOS variables
integer :: adios_err
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_aniso_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_aniso_mantle.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_aniso_mantle.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -77,7 +77,7 @@
if(myrank == 0) call read_aniso_mantle_model()
! broadcast the information read on the master to the nodes
- call bcast_all_i(AMM_V_npar1,1)
+ call bcast_all_singlei(AMM_V_npar1)
call bcast_all_dp(AMM_V_beta,14*34*37*73)
call bcast_all_dp(AMM_V_pro,47)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_atten3D_QRFSI12.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_atten3D_QRFSI12.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_atten3D_QRFSI12.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -80,10 +80,10 @@
! allocates model arrays
allocate(QRFSI12_Q_dqmu(NKQ,NSQ), &
- QRFSI12_Q_spknt(NKQ), &
- QRFSI12_Q_refdepth(NDEPTHS_REFQ), &
- QRFSI12_Q_refqmu(NDEPTHS_REFQ), &
- stat=ier)
+ QRFSI12_Q_spknt(NKQ), &
+ QRFSI12_Q_refdepth(NDEPTHS_REFQ), &
+ QRFSI12_Q_refqmu(NDEPTHS_REFQ), &
+ stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating QRFSI12_Q model arrays')
! master process reads in file values
@@ -203,7 +203,8 @@
real(kind=4) :: splpts(NKQ),splcon(NKQ),splcond(NKQ)
real(kind=4) :: depth,ylat,xlon
real(kind=4) :: shdep(NSQ)
- real(kind=4) :: xlmvec(NSQ),wk1(NSQ),wk2(NSQ),wk3(NSQ)
+ real(kind=4) :: wk1(NSQ),wk2(NSQ),wk3(NSQ)
+ real(kind=4) :: xlmvec(NSQ**2)
double precision, parameter :: rmoho_prem = 6371.0-24.4
double precision, parameter :: rcmb = 3480.0
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -102,9 +102,9 @@
if(myrank == 0) call read_attenuation_model(MIN_ATTENUATION_PERIOD, MAX_ATTENUATION_PERIOD, AM_V)
! broadcasts to all others
- call bcast_all_dp(AM_V%min_period, 1)
- call bcast_all_dp(AM_V%max_period, 1)
- call bcast_all_dp(AM_V%QT_c_source, 1)
+ call bcast_all_singledp(AM_V%min_period)
+ call bcast_all_singledp(AM_V%max_period)
+ call bcast_all_singledp(AM_V%QT_c_source)
call bcast_all_dp(AM_V%Qtau_s, N_SLS)
end subroutine model_attenuation_broadcast
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_eucrust.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_eucrust.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_eucrust.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -63,7 +63,7 @@
if( myrank == 0 ) call read_EuCrust()
! broadcast the information read on the master to the nodes
- call bcast_all_i(num_eucrust,1)
+ call bcast_all_singlei(num_eucrust)
! allocates on all other processes
if( myrank /= 0 ) then
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_gapp2.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_gapp2.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_gapp2.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -37,8 +37,8 @@
module gapp2_mantle_model_constants
! data file resolution
integer, parameter :: ma=288,mo=576,mr=32,mr1=64
- integer no,na,nnr,nr1
- real dela,delo
+ integer :: no,na,nnr,nr1
+ real :: dela,delo
! allocatable model arrays
real,dimension(:),allocatable :: dep,dep1,vp1
real,dimension(:,:,:),allocatable :: vp3
@@ -71,18 +71,18 @@
if(myrank == 0) call read_mantle_gapmodel()
! master process broadcasts data to all processes
- call bcast_all_r( dep,mr+1)
+ call bcast_all_r(dep,mr+1)
call bcast_all_r(dep1,mr1+1)
- call bcast_all_r( vp1,mr1+1)
- call bcast_all_r( vp3,ma*mo*mr)
+ call bcast_all_r(vp1,mr1+1)
+ call bcast_all_r(vp3,ma*mo*mr)
- call bcast_all_i( nnr,1)
- call bcast_all_i( nr1,1)
- call bcast_all_i( no,1)
- call bcast_all_i( na,1)
+ call bcast_all_singlei(nnr)
+ call bcast_all_singlei(nr1)
+ call bcast_all_singlei(no)
+ call bcast_all_singlei(na)
- call bcast_all_r( dela,1)
- call bcast_all_r( delo,1)
+ call bcast_all_singler(dela)
+ call bcast_all_singler(delo)
end subroutine model_gapp2_broadcast
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp3d.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp3d.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp3d.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -153,12 +153,12 @@
if(myrank == 0) call read_jp3d_iso_zhao_model()
! JP3DM_V
- call bcast_all_i(JP3DM_NPA,1)
- call bcast_all_i(JP3DM_NRA,1)
- call bcast_all_i(JP3DM_NHA,1)
- call bcast_all_i(JP3DM_NPB,1)
- call bcast_all_i(JP3DM_NRB,1)
- call bcast_all_i(JP3DM_NHB,1)
+ call bcast_all_singlei(JP3DM_NPA)
+ call bcast_all_singlei(JP3DM_NRA)
+ call bcast_all_singlei(JP3DM_NHA)
+ call bcast_all_singlei(JP3DM_NPB)
+ call bcast_all_singlei(JP3DM_NRB)
+ call bcast_all_singlei(JP3DM_NHB)
call bcast_all_dp(JP3DM_PNA,MPA)
call bcast_all_dp(JP3DM_RNA,MRA)
@@ -181,33 +181,33 @@
call bcast_all_i(JP3DM_IRLOCB,MKB)
call bcast_all_i(JP3DM_IHLOCB,MKB)
- call bcast_all_dp(JP3DM_PLA,1)
- call bcast_all_dp(JP3DM_RLA,1)
- call bcast_all_dp(JP3DM_HLA,1)
- call bcast_all_dp(JP3DM_PLB,1)
- call bcast_all_dp(JP3DM_RLB,1)
- call bcast_all_dp(JP3DM_HLB,1)
+ call bcast_all_singledp(JP3DM_PLA)
+ call bcast_all_singledp(JP3DM_RLA)
+ call bcast_all_singledp(JP3DM_HLA)
+ call bcast_all_singledp(JP3DM_PLB)
+ call bcast_all_singledp(JP3DM_RLB)
+ call bcast_all_singledp(JP3DM_HLB)
- call bcast_all_i(JP3DM_IP,1)
- call bcast_all_i(JP3DM_JP,1)
- call bcast_all_i(JP3DM_KP,1)
- call bcast_all_i(JP3DM_IP1,1)
- call bcast_all_i(JP3DM_JP1,1)
- call bcast_all_i(JP3DM_KP1,1)
+ call bcast_all_singlei(JP3DM_IP)
+ call bcast_all_singlei(JP3DM_JP)
+ call bcast_all_singlei(JP3DM_KP)
+ call bcast_all_singlei(JP3DM_IP1)
+ call bcast_all_singlei(JP3DM_JP1)
+ call bcast_all_singlei(JP3DM_KP1)
call bcast_all_dp(JP3DM_WV,8)
- call bcast_all_dp(JP3DM_P,1)
- call bcast_all_dp(JP3DM_R,1)
- call bcast_all_dp(JP3DM_H,1)
- call bcast_all_dp(JP3DM_PF,1)
- call bcast_all_dp(JP3DM_RF,1)
- call bcast_all_dp(JP3DM_HF,1)
- call bcast_all_dp(JP3DM_PF1,1)
- call bcast_all_dp(JP3DM_RF1,1)
- call bcast_all_dp(JP3DM_HF1,1)
- call bcast_all_dp(JP3DM_PD,1)
- call bcast_all_dp(JP3DM_RD,1)
- call bcast_all_dp(JP3DM_HD,1)
+ call bcast_all_singledp(JP3DM_P)
+ call bcast_all_singledp(JP3DM_R)
+ call bcast_all_singledp(JP3DM_H)
+ call bcast_all_singledp(JP3DM_PF)
+ call bcast_all_singledp(JP3DM_RF)
+ call bcast_all_singledp(JP3DM_HF)
+ call bcast_all_singledp(JP3DM_PF1)
+ call bcast_all_singledp(JP3DM_RF1)
+ call bcast_all_singledp(JP3DM_HF1)
+ call bcast_all_singledp(JP3DM_PD)
+ call bcast_all_singledp(JP3DM_RD)
+ call bcast_all_singledp(JP3DM_HD)
call bcast_all_dp(JP3DM_VP,29)
call bcast_all_dp(JP3DM_VS,29)
call bcast_all_dp(JP3DM_RA,29)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_ppm.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_ppm.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_ppm.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -107,9 +107,10 @@
if(myrank == 0) call read_model_ppm()
! broadcast the information read on the master to the nodes
- call bcast_all_i(PPM_num_v,1)
- call bcast_all_i(PPM_num_latperlon,1)
- call bcast_all_i(PPM_num_lonperdepth,1)
+ call bcast_all_singlei(PPM_num_v)
+ call bcast_all_singlei(PPM_num_latperlon)
+ call bcast_all_singlei(PPM_num_lonperdepth)
+
if( myrank /= 0 ) then
allocate(PPM_lat(PPM_num_v),PPM_lon(PPM_num_v),PPM_depth(PPM_num_v),PPM_dvs(PPM_num_v))
endif
@@ -118,15 +119,16 @@
call bcast_all_dp(PPM_lat(1:PPM_num_v),PPM_num_v)
call bcast_all_dp(PPM_lon(1:PPM_num_v),PPM_num_v)
call bcast_all_dp(PPM_depth(1:PPM_num_v),PPM_num_v)
- call bcast_all_dp(PPM_maxlat,1)
- call bcast_all_dp(PPM_minlat,1)
- call bcast_all_dp(PPM_maxlon,1)
- call bcast_all_dp(PPM_minlon,1)
- call bcast_all_dp(PPM_maxdepth,1)
- call bcast_all_dp(PPM_mindepth,1)
- call bcast_all_dp(PPM_dlat,1)
- call bcast_all_dp(PPM_dlon,1)
- call bcast_all_dp(PPM_ddepth,1)
+
+ call bcast_all_singledp(PPM_maxlat)
+ call bcast_all_singledp(PPM_minlat)
+ call bcast_all_singledp(PPM_maxlon)
+ call bcast_all_singledp(PPM_minlon)
+ call bcast_all_singledp(PPM_maxdepth)
+ call bcast_all_singledp(PPM_mindepth)
+ call bcast_all_singledp(PPM_dlat)
+ call bcast_all_singledp(PPM_dlon)
+ call bcast_all_singledp(PPM_ddepth)
end subroutine model_ppm_broadcast
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -115,9 +115,10 @@
if(myrank == 0) call read_model_s362ani(THREE_D_MODEL,THREE_D_MODEL_S362ANI,THREE_D_MODEL_S362WMANI, &
THREE_D_MODEL_S362ANI_PREM,THREE_D_MODEL_S29EA)
- call bcast_all_i(numker,1)
- call bcast_all_i(numhpa,1)
- call bcast_all_i(ihpa,1)
+ call bcast_all_singlei(numker)
+ call bcast_all_singlei(numhpa)
+ call bcast_all_singlei(ihpa)
+
call bcast_all_i(lmxhpa,maxhpa)
call bcast_all_i(itypehpa,maxhpa)
call bcast_all_i(ihpakern,maxker)
@@ -1773,8 +1774,6 @@
subroutine ylm(XLAT,XLON,LMAX,Y,WK1,WK2,WK3)
- use model_s362ani_par,only : maxl
-
implicit none
complex TEMP,FAC,DFAC
@@ -1786,9 +1785,9 @@
!
! WK1,WK2,WK3 SHOULD BE DIMENSIONED AT LEAST (LMAX+1)*4
!
- real(kind=4) WK1(LMAX+1),WK2(LMAX+1),WK3(LMAX+1)
- real(kind=4) XLAT,XLON
- real(kind=4),dimension((maxl+1)**2) :: Y !! Y should go at least from 1 to fac(LMAX)
+ real(kind=4),dimension(LMAX+1) :: WK1,WK2,WK3
+ real(kind=4) :: XLAT,XLON
+ real(kind=4),dimension((LMAX+1)**2) :: Y !! Y should go at least from 1 to fac(LMAX)
real(kind=4), parameter :: RADIAN = 57.2957795
@@ -1806,8 +1805,10 @@
! index L goes from 0 to LMAX
L=IL1-1
+
+ !! see legndr(): WK1,WK2,WK3 should go from 1 to L+1
!CALL legndr(THETA,L,L,WK1,WK2,WK3)
- CALL legndr(THETA,L,L,WK1(1:L+1),WK2(1:L+1),WK3(1:L+1)) !! see legndr(): WK1,WK2,WK3 should go from 1 to L+1
+ CALL legndr(THETA,L,L,WK1(1:L+1),WK2(1:L+1),WK3(1:L+1))
FAC=(1.,0.)
DFAC=CEXP(CMPLX(0.,PHI))
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_sea99_s.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_sea99_s.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_sea99_s.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -74,25 +74,27 @@
! allocates model arrays
allocate(sea99_vs(100,100,100), &
- sea99_depth(100), &
- stat=ier)
+ sea99_depth(100), &
+ stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating sea99 arrays')
! master proc reads in values
if(myrank == 0) call read_sea99_s_model()
! broadcast the information read on the master to the nodes
- call bcast_all_i(sea99_ndep,1)
- call bcast_all_i(sea99_nlat,1)
- call bcast_all_i(sea99_nlon,1)
- call bcast_all_dp(sea99_ddeg,1)
- call bcast_all_dp(alatmin,1)
- call bcast_all_dp(alatmax,1)
- call bcast_all_dp(alonmin,1)
- call bcast_all_i(alonmax,1)
- call bcast_all_i(sea99_vs,100*100*100)
- call bcast_all_i(sea99_depth,100)
+ call bcast_all_singlei(sea99_ndep)
+ call bcast_all_singlei(sea99_nlat)
+ call bcast_all_singlei(sea99_nlon)
+ call bcast_all_singledp(sea99_ddeg)
+ call bcast_all_singledp(alatmin)
+ call bcast_all_singledp(alatmax)
+ call bcast_all_singledp(alonmin)
+ call bcast_all_singledp(alonmax)
+
+ call bcast_all_dp(sea99_vs,100*100*100)
+ call bcast_all_dp(sea99_depth,100)
+
end subroutine model_sea99_s_broadcast
!
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/rules.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/rules.mk 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/rules.mk 2013-10-22 14:57:01 UTC (rev 22968)
@@ -214,7 +214,7 @@
$(EMPTY_MACRO)
adios_meshfem3D_SHARED_STUBS = \
- $O/adios_method_stubs.shared.o \
+ $O/adios_method_stubs.cc.o \
$(EMPTY_MACRO)
# conditional adios linking
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver_adios.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver_adios.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver_adios.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -1243,7 +1243,7 @@
!! \param phase_ispec_inner
!! \param num_colors_inner Number of colors for GPU computing in the inner core.
!! \param num_colors_outer Number of colors for GPU computing in the outer core.
-subroutine save_MPI_arrays_adios(myrank,iregion_code,LOCAL_PATH, &
+subroutine save_mpi_arrays_adios(myrank,iregion_code,LOCAL_PATH, &
num_interfaces,max_nibool_interfaces, my_neighbours,nibool_interfaces, &
ibool_interfaces, nspec_inner,nspec_outer, num_phase_ispec, &
phase_ispec_inner, num_colors_outer,num_colors_inner, num_elem_colors)
@@ -1434,7 +1434,7 @@
num_regions_written = num_regions_written + 1
-end subroutine save_MPI_arrays_adios
+end subroutine save_mpi_arrays_adios
!===============================================================================
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_method_stubs.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_method_stubs.c (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_method_stubs.c 2013-10-22 14:57:01 UTC (rev 22968)
@@ -0,0 +1,113 @@
+/*
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 6 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! August 2013
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "config.h"
+
+typedef float realw;
+
+// placeholders for non-adios compilation
+
+void warn_no_adios() {
+ fprintf(stderr,"ERROR: ADIOS enabled without ADIOS Support. To enable ADIOS support, reconfigure with --with-adios flag.\n");
+ exit(1);
+}
+
+void FC_FUNC_(adios_setup,ADIOS_SETUP)(){ warn_no_adios(); }
+
+void FC_FUNC_(adios_cleanup,ADIOS_CLEANUP)(){}
+
+
+
+// for xmeshfem3D compilation
+
+void FC_FUNC_(crm_save_mesh_files_adios,CRM_SAVE_MESH_FILES_ADIOS)(int* nspec, int* npointot, int* iregion_code,
+ int* num_ibool_AVS_DX, int* mask_ibool){}
+
+void FC_FUNC_(get_absorb_adios,GET_ABSORB_ADIOS)(int* myrank, int* iregion,
+ int* nimin, int* nimax, int* njmin, int* njmax, int* nkmin_xi, int* nkmin_eta,
+ int* NSPEC2DMAX_XMIN_XMAX, int* NSPEC2DMAX_YMIN_YMAX){}
+
+void FC_FUNC_(save_arrays_solver_adios,SAVE_ARRAYS_SOLVER_ADIOS)(int* myrank, int* nspec, int* nglob,
+ int* idoubling, int* ibool,
+ int* iregion_code,
+ realw* xstore, realw* ystore, realw* zstore,
+ int* NSPEC2DMAX_XMIN_XMAX, int* NSPEC2DMAX_YMIN_YMAX,
+ int* NSPEC2D_TOP, int* NSPEC2D_BOTTOM){}
+
+void FC_FUNC_(save_arrays_solver_meshfiles_adios,SAVE_ARRAYS_SOLVER_MESHFILES_ADIOS)(){}
+
+void FC_FUNC_(save_arrays_solver_boundary_adios,SAVE_ARRAYS_SOLVER_BOUNDARY_ADIOS)(){}
+
+void FC_FUNC_(save_mpi_arrays_adios,SAVE_MPI_ARRAYS_ADIOS)(){}
+
+void FC_FUNC_(read_gll_model_adios,READ_GLL_MODEL_ADIOS)(){}
+
+// for xspecfem3D compilation
+
+void FC_FUNC_(read_arrays_solver_adios,READ_ARRAYS_SOLVER_ADIOS)(){}
+
+void FC_FUNC_(read_attenuation_adios,READ_ATTENUATION_ADIOS)(){}
+
+void FC_FUNC_(read_forward_arrays_adios,READ_FORWARD_ARRAYS_ADIOS)(){}
+
+void FC_FUNC_(read_intermediate_forward_arrays_adios,READ_INTERMEDIATE_FORWARD_ARRAYS_ADIOS)(){}
+
+void FC_FUNC_(read_mesh_databases_coupling_adios,READ_MESH_DATABASES_COUPLING_ADIOS)(){}
+
+void FC_FUNC_(read_mesh_databases_mpi_cm_adios,READ_MESH_DATABASES_MPI_CM_ADIOS)(){}
+
+void FC_FUNC_(read_mesh_databases_mpi_ic_adios,READ_MESH_DATABASES_MPI_IC_ADIOS)(){}
+
+void FC_FUNC_(read_mesh_databases_mpi_oc_adios,READ_MESH_DATABASES_MPI_OC_ADIOS)(){}
+
+void FC_FUNC_(read_mesh_databases_stacey_adios,READ_MESH_DATABASES_STACEY_ADIOS)(){}
+
+void FC_FUNC_(save_forward_arrays_adios,SAVE_FORWARD_ARRAYS_ADIOS)(){}
+
+void FC_FUNC_(save_intermediate_forward_arrays_adios,SAVE_INTERMEDIATE_FORWARD_ARRAYS_ADIOS)(){}
+
+void FC_FUNC_(perform_write_adios_kernels,PERFORM_WRITE_ADIOS_KERNELS)(){}
+
+void FC_FUNC_(define_kernel_adios_variables,DEFINE_KERNEL_ADIOS_VARIABLES)(){}
+
+void FC_FUNC_(write_kernels_crust_mantle_adios,WRITE_KERNELS_CRUST_MANTLE_ADIOS)(){}
+
+void FC_FUNC_(write_kernels_outer_core_adios,WRITE_KERNELS_OUTER_CORE_ADIOS)(){}
+
+void FC_FUNC_(write_kernels_inner_core_adios,WRITE_KERNELS_INNER_CORE_ADIOS)(){}
+
+void FC_FUNC_(write_kernels_boundary_kl_adios,WRITE_KERNELS_BOUNDARY_KL_ADIOS)(){}
+
+void FC_FUNC_(write_kernels_source_derivatives_adios,WRITE_KERNELS_SOURCE_DERIVATIVES_ADIOS)(){}
+
+void FC_FUNC_(write_kernels_hessian_adios,WRITE_KERNELS_HESSIAN_ADIOS)(){}
+
Deleted: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_method_stubs.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_method_stubs.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/adios_method_stubs.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -1,219 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 6 . 0
-! --------------------------------------------------
-!
-! Main authors: Dimitri Komatitsch and Jeroen Tromp
-! Princeton University, USA
-! and CNRS / INRIA / University of Pau, France
-! (c) Princeton University and CNRS / INRIA / University of Pau
-! August 2013
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
-
-! placeholders for non-adios compilation
-
- subroutine warn_no_adios()
- stop 'Trying to use ADIOS while not configured for it.'
- end subroutine warn_no_adios
-
-! modules
-
-! module AVS_DX_global_mod
-! type avs_dx_global_t
-! end type avs_dx_global_t
-! end module AVS_DX_global_mod
-
-! module AVS_DX_global_faces_mod
-! type avs_dx_global_faces_t
-! end type avs_dx_global_faces_t
-! end module AVS_DX_global_faces_mod
-
-! module AVS_DX_global_chunks_mod
-! type avs_dx_global_chunks_t
-! end type avs_dx_global_chunks_t
-! end module AVS_DX_global_chunks_mod
-
-! module AVS_DX_surface_mod
-! type avs_dx_surface_t
-! end type avs_dx_surface_t
-! end module AVS_DX_surface_mod
-
-! for both xmeshfem3D/xspecfem3D compilation
-
- subroutine adios_setup()
- call warn_no_adios()
- end subroutine
-
- subroutine adios_cleanup()
- end subroutine
-
-! for xmeshfem3D compilation
-
- subroutine crm_save_mesh_files_adios()
- end subroutine
-
- subroutine get_absorb_adios()
- end subroutine
-
- subroutine save_arrays_solver_adios()
- end subroutine
-
- subroutine save_arrays_solver_meshfiles_adios()
- end subroutine
-
- subroutine save_arrays_solver_boundary_adios()
- end subroutine
-
- subroutine save_MPI_arrays_adios()
- end subroutine
-
- subroutine read_gll_model_adios
- end subroutine
-
-! subroutine prepare_AVS_DX_global_chunks_data_adios()
-! end subroutine
-!
-! subroutine write_AVS_DX_global_chunks_data_adios()
-! end subroutine
-!
-! subroutine free_AVS_DX_global_chunks_data_adios()
-! end subroutine
-!
-! subroutine define_AVS_DX_global_data_adios()
-! end subroutine
-!
-! subroutine prepare_AVS_DX_global_data_adios()
-! end subroutine
-!
-! subroutine write_AVS_DX_global_data_adios()
-! end subroutine
-!
-! subroutine free_AVS_DX_global_data_adios()
-! end subroutine
-!
-! subroutine define_AVS_DX_global_faces_data_adios()
-! end subroutine
-!
-! subroutine prepare_AVS_DX_global_faces_data_adios()
-! end subroutine
-!
-! subroutine write_AVS_DX_global_faces_data_adios()
-! end subroutine
-!
-! subroutine free_AVS_DX_global_faces_data_adios()
-! end subroutine
-!
-! subroutine define_AVS_DX_surfaces_data_adios()
-! end subroutine
-!
-! subroutine prepare_AVS_DX_surfaces_data_adios()
-! end subroutine
-!
-! subroutine write_AVS_DX_surfaces_data_adios()
-! end subroutine
-!
-! subroutine free_AVS_DX_surfaces_data_adios()
-! end subroutine
-!
-
-! for xspecfem3D compilation
-
- subroutine read_arrays_solver_adios()
- end subroutine
-
- subroutine read_attenuation_adios()
- end subroutine
-
- subroutine read_forward_arrays_adios()
- end subroutine
-
- subroutine read_intermediate_forward_arrays_adios()
- end subroutine
-
- subroutine read_mesh_databases_coupling_adios()
- end subroutine
-
- subroutine read_mesh_databases_mpi_cm_adios()
- end subroutine
-
- subroutine read_mesh_databases_mpi_ic_adios()
- end subroutine
-
- subroutine read_mesh_databases_mpi_oc_adios()
- end subroutine
-
- subroutine read_mesh_databases_stacey_adios()
- end subroutine
-
- subroutine save_forward_arrays_adios()
- end subroutine
-
- subroutine save_intermediate_forward_arrays_adios()
- end subroutine
-
-! subroutine define_common_forward_arrays_adios()
-! end subroutine
-!
-! subroutine define_rotation_forward_arrays_adios()
-! end subroutine
-!
-! subroutine define_attenuation_forward_arrays_adios()
-! end subroutine
-!
-! subroutine write_common_forward_arrays_adios()
-! end subroutine
-!
-! subroutine write_rotation_forward_arrays_adios()
-! end subroutine
-!
-! subroutine write_attenuation_forward_arrays_adios()
-! end subroutine
-!
-! subroutine write_1D_global_array_adios_dims()
-! end subroutine
-!
-
- subroutine perform_write_adios_kernels()
- end subroutine
-
- subroutine define_kernel_adios_variables()
- end subroutine
-
- subroutine write_kernels_crust_mantle_adios()
- end subroutine
-
- subroutine write_kernels_outer_core_adios()
- end subroutine
-
- subroutine write_kernels_inner_core_adios()
- end subroutine
-
- subroutine write_kernels_boundary_kl_adios()
- end subroutine
-
- subroutine write_kernels_source_derivatives_adios()
- end subroutine
-
- subroutine write_kernels_hessian_adios()
- end subroutine
-
-!
-! subroutine write_specfem_header_adios()
-! end subroutine
-
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/memory_eval.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/memory_eval.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/memory_eval.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -47,13 +47,15 @@
static_memory_size)
use constants
- use shared_parameters,only: ATT1,ATT2,ATT3, &
+ use shared_parameters,only: &
+ ATT1,ATT2,ATT3, &
APPROXIMATE_HESS_KL,ANISOTROPIC_KL,NOISE_TOMOGRAPHY, &
EXACT_MASS_MATRIX_FOR_ROTATION, &
OCEANS,ABSORBING_CONDITIONS,ATTENUATION,ANISOTROPIC_3D_MANTLE, &
TRANSVERSE_ISOTROPY,ANISOTROPIC_INNER_CORE,ROTATION,TOPOGRAPHY, &
ONE_CRUST,NCHUNKS, &
- SIMULATION_TYPE,MOVIE_VOLUME,SAVE_FORWARD
+ SIMULATION_TYPE,SAVE_FORWARD, &
+ MOVIE_VOLUME,MOVIE_VOLUME_TYPE
implicit none
@@ -434,6 +436,12 @@
endif
+ ! div_displ_outer_core
+ if( MOVIE_VOLUME .and. MOVIE_VOLUME_TYPE == 4 ) then
+ static_memory_size = static_memory_size + &
+ dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_OUTER_CORE)*dble(CUSTOM_REAL)
+ endif
+
! add arrays used for adjoint runs only (LQY: not very accurate)
! b_R_memory_crust_mantle
@@ -443,10 +451,9 @@
static_memory_size = static_memory_size + (5.d0*dble(N_SLS) + 9.d0)* &
dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_ADJOINT*dble(CUSTOM_REAL)
- ! b_div_displ_outer_core
! rho_kl_outer_core,alpha_kl_outer_core
static_memory_size = static_memory_size + &
- 3.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_OUTER_CORE_ADJOINT*dble(CUSTOM_REAL)
+ 2.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_OUTER_CORE_ADJOINT*dble(CUSTOM_REAL)
! b_R_memory_inner_core
! b_epsilondev_inner_core
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/parallel.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/parallel.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/parallel.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -500,12 +500,9 @@
subroutine bcast_all_r(buffer, count)
use mpi
- use constants
implicit none
- include "precision.h"
-
integer :: count
real, dimension(count) :: buffer
@@ -519,6 +516,25 @@
!-------------------------------------------------------------------------------------------------
!
+ subroutine bcast_all_singler(buffer)
+
+ use mpi
+
+ implicit none
+
+ real :: buffer
+
+ integer :: ier
+
+ call MPI_BCAST(buffer,1,MPI_REAL,0,MPI_COMM_WORLD,ier)
+
+ end subroutine bcast_all_singler
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine bcast_all_dp(buffer, count)
use :: mpi
@@ -534,7 +550,25 @@
end subroutine bcast_all_dp
+!
+!-------------------------------------------------------------------------------------------------
+!
+ subroutine bcast_all_singledp(buffer)
+
+ use :: mpi
+
+ implicit none
+
+ double precision :: buffer
+
+ integer :: ier
+
+ call MPI_BCAST(buffer,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+
+ end subroutine bcast_all_singledp
+
+
!
!-------------------------------------------------------------------------------------------------
!
@@ -910,6 +944,28 @@
!-------------------------------------------------------------------------------------------------
!
+ subroutine gather_all_singlei(sendbuf, recvbuf, NPROC)
+
+ use mpi
+
+ implicit none
+
+ integer :: NPROC
+ integer :: sendbuf
+ integer, dimension(0:NPROC-1) :: recvbuf
+
+ integer :: ier
+
+ call MPI_GATHER(sendbuf,1,MPI_INTEGER, &
+ recvbuf,1,MPI_INTEGER, &
+ 0,MPI_COMM_WORLD,ier)
+
+ end subroutine gather_all_singlei
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine gather_all_cr(sendbuf, sendcnt, recvbuf, recvcount, NPROC)
use mpi
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/param_reader.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/param_reader.c 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/param_reader.c 2013-10-22 14:57:01 UTC (rev 22968)
@@ -1,31 +1,31 @@
-/*
- !=====================================================================
- !
- ! S p e c f e m 3 D G l o b e V e r s i o n 6 . 0
- ! --------------------------------------------------
- !
- ! Main authors: Dimitri Komatitsch and Jeroen Tromp
- ! Princeton University, USA
- ! and CNRS / INRIA / University of Pau, France
- ! (c) Princeton University and CNRS / INRIA / University of Pau
- ! August 2013
- !
- ! This program is free software; you can redistribute it and/or modify
- ! it under the terms of the GNU General Public License as published by
- ! the Free Software Foundation; either version 2 of the License, or
- ! (at your option) any later version.
- !
- ! This program is distributed in the hope that it will be useful,
- ! but WITHOUT ANY WARRANTY; without even the implied warranty of
- ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- ! GNU General Public License for more details.
- !
- ! You should have received a copy of the GNU General Public License along
- ! with this program; if not, write to the Free Software Foundation, Inc.,
- ! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
- !
- !=====================================================================
- */
+/*
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 6 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! August 2013
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+*/
/*
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/rules.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/rules.mk 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/rules.mk 2013-10-22 14:57:01 UTC (rev 22968)
@@ -91,7 +91,7 @@
$(EMPTY_MACRO)
adios_shared_STUBS = \
- $O/adios_method_stubs.shared.o \
+ $O/adios_method_stubs.cc.o \
$(EMPTY_MACRO)
ifeq ($(ADIOS),yes)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -78,7 +78,7 @@
ABSORBING_CONDITIONS,EXACT_MASS_MATRIX_FOR_ROTATION, &
ATTENUATION_1D_WITH_3D_STORAGE, &
ATT1,ATT2,ATT3,ATT4,ATT5, &
- MOVIE_VOLUME
+ MOVIE_VOLUME,MOVIE_VOLUME_TYPE
implicit none
@@ -537,6 +537,12 @@
endif
write(IOUT,*)
+ if( MOVIE_VOLUME .and. MOVIE_VOLUME_TYPE == 4 ) then
+ write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE_3DMOVIE = NSPEC_OUTER_CORE'
+ else
+ write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE_3DMOVIE = 1'
+ endif
+
if (SAVE_REGULAR_KL) then
write(IOUT,*) 'integer, parameter :: NM_KL_REG_PTS_VAL = NM_KL_REG_PTS'
else
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_coupling.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_coupling.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_coupling.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -25,24 +25,25 @@
!
!=====================================================================
- subroutine compute_coupling_fluid_CMB(displ_crust_mantle, &
- ibool_crust_mantle,ibelm_bottom_crust_mantle, &
- accel_outer_core, &
- normal_top_outer_core,jacobian2D_top_outer_core, &
- wgllwgll_xy,ibool_outer_core,ibelm_top_outer_core, &
- nspec2D_top)
+ subroutine compute_coupling_fluid_CMB(NGLOB_CM,displ_crust_mantle, &
+ ibool_crust_mantle,ibelm_bottom_crust_mantle, &
+ NGLOB_OC,accel_outer_core, &
+ normal_top_outer_core,jacobian2D_top_outer_core, &
+ wgllwgll_xy,ibool_outer_core,ibelm_top_outer_core, &
+ nspec2D_top)
use constants_solver
implicit none
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
- displ_crust_mantle
+ integer :: NGLOB_CM
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CM) :: displ_crust_mantle
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
integer, dimension(NSPEC2D_BOTTOM_CM) :: ibelm_bottom_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: accel_outer_core
+ integer :: NGLOB_OC
+ real(kind=CUSTOM_REAL), dimension(NGLOB_OC) :: accel_outer_core
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_TOP_OC) :: normal_top_outer_core
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_TOP_OC) :: jacobian2D_top_outer_core
@@ -104,9 +105,9 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine compute_coupling_fluid_ICB(displ_inner_core, &
+ subroutine compute_coupling_fluid_ICB(NGLOB_IC,displ_inner_core, &
ibool_inner_core,ibelm_top_inner_core, &
- accel_outer_core, &
+ NGLOB_OC,accel_outer_core, &
normal_bottom_outer_core,jacobian2D_bottom_outer_core, &
wgllwgll_xy,ibool_outer_core,ibelm_bottom_outer_core, &
nspec_bottom)
@@ -115,13 +116,15 @@
implicit none
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: &
+ integer :: NGLOB_IC
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_IC) :: &
displ_inner_core
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: ibool_inner_core
integer, dimension(NSPEC2D_TOP_IC) :: ibelm_top_inner_core
- real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: accel_outer_core
+ integer :: NGLOB_OC
+ real(kind=CUSTOM_REAL), dimension(NGLOB_OC) :: accel_outer_core
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM_OC) :: normal_bottom_outer_core
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_BOTTOM_OC) :: jacobian2D_bottom_outer_core
@@ -183,10 +186,9 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine compute_coupling_CMB_fluid(displ_crust_mantle, &
- accel_crust_mantle, &
+ subroutine compute_coupling_CMB_fluid(NGLOB_CM,displ_crust_mantle,accel_crust_mantle, &
ibool_crust_mantle,ibelm_bottom_crust_mantle, &
- accel_outer_core, &
+ NGLOB_OC,accel_outer_core, &
normal_top_outer_core,jacobian2D_top_outer_core, &
wgllwgll_xy,ibool_outer_core,ibelm_top_outer_core, &
RHO_TOP_OC,minus_g_cmb, &
@@ -196,13 +198,14 @@
implicit none
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
- displ_crust_mantle,accel_crust_mantle
+ integer :: NGLOB_CM
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CM) :: displ_crust_mantle,accel_crust_mantle
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
integer, dimension(NSPEC2D_BOTTOM_CM) :: ibelm_bottom_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: accel_outer_core
+ integer :: NGLOB_OC
+ real(kind=CUSTOM_REAL), dimension(NGLOB_OC) :: accel_outer_core
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_TOP_OC) :: normal_top_outer_core
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_TOP_OC) :: jacobian2D_top_outer_core
@@ -270,9 +273,9 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine compute_coupling_ICB_fluid(displ_inner_core,accel_inner_core, &
+ subroutine compute_coupling_ICB_fluid(NGLOB_IC,displ_inner_core,accel_inner_core, &
ibool_inner_core,ibelm_top_inner_core, &
- accel_outer_core, &
+ NGLOB_OC,accel_outer_core, &
normal_bottom_outer_core,jacobian2D_bottom_outer_core, &
wgllwgll_xy,ibool_outer_core,ibelm_bottom_outer_core, &
RHO_BOTTOM_OC,minus_g_icb, &
@@ -282,13 +285,14 @@
implicit none
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: &
- displ_inner_core,accel_inner_core
+ integer :: NGLOB_IC
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_IC) :: displ_inner_core,accel_inner_core
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: ibool_inner_core
integer, dimension(NSPEC2D_TOP_IC) :: ibelm_top_inner_core
- real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: accel_outer_core
+ integer :: NGLOB_OC
+ real(kind=CUSTOM_REAL), dimension(NGLOB_OC) :: accel_outer_core
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM_OC) :: normal_bottom_outer_core
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_BOTTOM_OC) :: jacobian2D_bottom_outer_core
@@ -356,7 +360,7 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine compute_coupling_ocean(accel_crust_mantle, &
+ subroutine compute_coupling_ocean(NGLOB,accel_crust_mantle, &
rmassx_crust_mantle, rmassy_crust_mantle, rmassz_crust_mantle, &
rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
@@ -367,7 +371,8 @@
implicit none
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: accel_crust_mantle
+ integer :: NGLOB
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: accel_crust_mantle
! mass matrices
!
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -32,7 +32,10 @@
!
!--------------------------------------------------------------------------------------------
+#ifndef FORCE_VECTORIZATION
+ ! default, without vectorization
+
subroutine compute_element_iso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
xstore,ystore,zstore, &
@@ -119,7 +122,6 @@
double precision factor,sx_l,sy_l,sz_l,gxl,gyl,gzl
double precision Hxxl,Hyyl,Hzzl,Hxyl,Hxzl,Hyzl
- integer :: ispec_strain
integer :: i,j,k
integer :: int_radius
integer :: iglob
@@ -171,15 +173,21 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
+!ZN from the expression in which we use the strain later in the code.
+
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
- ispec_strain = 1
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ if( ispec == 1) then
+ epsilon_trace_over_3(i,j,k,1) = templ
+ endif
else
- ispec_strain = ispec
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ epsilon_trace_over_3(i,j,k,ispec) = templ
endif
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
@@ -372,13 +380,347 @@
end subroutine compute_element_iso
+#else
+
+! FORCE_VECTORIZATION
+! vectorized routine
+
+ subroutine compute_element_iso(ispec, &
+ minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ xstore,ystore,zstore, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ wgll_cube, &
+ kappavstore,muvstore, &
+ ibool, &
+ R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilon_trace_over_3, &
+ one_minus_sum_beta,vnspec, &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
+ dummyx_loc,dummyy_loc,dummyz_loc, &
+ epsilondev_loc,rho_s_H,is_backward_field)
+
+ use constants_solver
+ use specfem_par,only: COMPUTE_AND_STORE_STRAIN
+
+ implicit none
+
+
+ ! element id
+ integer :: ispec
+
+ ! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool
+
+ ! x y and z contain r theta and phi
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: xstore,ystore,zstore
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ ! array with derivatives of Lagrange polynomials and precalculated products
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+
+ ! store anisotropic properties only where needed to save memory
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: &
+ kappavstore,muvstore
+
+ ! variable sized array variables
+ integer :: vnspec
+
+ ! attenuation
+ ! memory variables for attenuation
+ ! memory variables R_ij are stored at the local rather than global level
+ ! to allow for optimization of cache access by compiler
+ real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUATION) :: R_xx,R_yy,R_xy,R_xz,R_yz
+
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
+
+ real(kind=CUSTOM_REAL), dimension(ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: one_minus_sum_beta
+
+ ! gravity
+ double precision, dimension(NRAD_GRAVITY) :: minus_gravity_table,density_table,minus_deriv_gravity_table
+
+ ! element info
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+
+ logical :: is_backward_field
+
+! local parameters
+ real(kind=CUSTOM_REAL) one_minus_sum_beta_use
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+ real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+ real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
+
+ real(kind=CUSTOM_REAL) templ
+ real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
+ real(kind=CUSTOM_REAL) kappal
+
+ ! for gravity
+ double precision dphi,dtheta
+ double precision radius,rho,minus_g,minus_dg
+ double precision minus_g_over_radius,minus_dg_plus_g_over_radius
+ double precision cos_theta,sin_theta,cos_phi,sin_phi
+ double precision cos_theta_sq,sin_theta_sq,cos_phi_sq,sin_phi_sq
+ double precision factor,sx_l,sy_l,sz_l,gxl,gyl,gzl
+ double precision Hxxl,Hyyl,Hzzl,Hxyl,Hxzl,Hyzl
+
+ integer :: int_radius
+ integer :: iglob
+
+ logical :: dummyl
+
+ ! force vectorization
+ integer :: ijk
+ real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
+
+! in this vectorized version we have to assume that N_SLS == 3 in order to be able to unroll and thus suppress
+! an inner loop that would otherwise prevent vectorization; this is safe in practice in all cases because N_SLS == 3
+! in all known applications, and in the main program we check that N_SLS == 3 if FORCE_VECTORIZATION is used and we stop otherwise
+
+!daniel debug
+ ! dummy to avoid compiler warning
+ dummyl = is_backward_field
+
+ ! isotropic element
+
+ do ijk=1,NGLLCUBE
+
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(ijk,1,1,ispec)
+ xiyl = xiy(ijk,1,1,ispec)
+ xizl = xiz(ijk,1,1,ispec)
+ etaxl = etax(ijk,1,1,ispec)
+ etayl = etay(ijk,1,1,ispec)
+ etazl = etaz(ijk,1,1,ispec)
+ gammaxl = gammax(ijk,1,1,ispec)
+ gammayl = gammay(ijk,1,1,ispec)
+ gammazl = gammaz(ijk,1,1,ispec)
+
+ ! compute the jacobian
+ jacobianl = 1.0_CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1(ijk,1,1) + etaxl*tempx2(ijk,1,1) + gammaxl*tempx3(ijk,1,1)
+ duxdyl = xiyl*tempx1(ijk,1,1) + etayl*tempx2(ijk,1,1) + gammayl*tempx3(ijk,1,1)
+ duxdzl = xizl*tempx1(ijk,1,1) + etazl*tempx2(ijk,1,1) + gammazl*tempx3(ijk,1,1)
+
+ duydxl = xixl*tempy1(ijk,1,1) + etaxl*tempy2(ijk,1,1) + gammaxl*tempy3(ijk,1,1)
+ duydyl = xiyl*tempy1(ijk,1,1) + etayl*tempy2(ijk,1,1) + gammayl*tempy3(ijk,1,1)
+ duydzl = xizl*tempy1(ijk,1,1) + etazl*tempy2(ijk,1,1) + gammazl*tempy3(ijk,1,1)
+
+ duzdxl = xixl*tempz1(ijk,1,1) + etaxl*tempz2(ijk,1,1) + gammaxl*tempz3(ijk,1,1)
+ duzdyl = xiyl*tempz1(ijk,1,1) + etayl*tempz2(ijk,1,1) + gammayl*tempz3(ijk,1,1)
+ duzdzl = xizl*tempz1(ijk,1,1) + etazl*tempz2(ijk,1,1) + gammazl*tempz3(ijk,1,1)
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
+!ZN from the expression in which we use the strain later in the code.
+
+ ! compute deviatoric strain
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if( NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1 ) then
+ if( ispec == 1) then
+ epsilon_trace_over_3(ijk,1,1,1) = templ
+ endif
+ else
+ epsilon_trace_over_3(ijk,1,1,ispec) = templ
+ endif
+ epsilondev_loc(1,ijk,1,1) = duxdxl - templ
+ epsilondev_loc(2,ijk,1,1) = duydyl - templ
+ epsilondev_loc(3,ijk,1,1) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
+ epsilondev_loc(4,ijk,1,1) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
+ epsilondev_loc(5,ijk,1,1) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
+ endif
+
+ !
+ ! compute isotropic elements
+ !
+
+ ! layer with no transverse isotropy, use kappav and muv
+ kappal = kappavstore(ijk,1,1,ispec)
+ mul = muvstore(ijk,1,1,ispec)
+
+ ! use unrelaxed parameters if attenuation
+ if(ATTENUATION_VAL) then
+ ! precompute terms for attenuation if needed
+ if( ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL ) then
+ one_minus_sum_beta_use = one_minus_sum_beta(ijk,1,1,ispec)
+ else
+ one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
+ endif
+ mul = mul * one_minus_sum_beta_use
+ endif
+
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.0_CUSTOM_REAL*mul
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+ sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+ sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+
+ sigma_xy = mul*duxdyl_plus_duydxl
+ sigma_xz = mul*duzdxl_plus_duxdzl
+ sigma_yz = mul*duzdyl_plus_duydzl
+
+ ! subtract memory variables if attenuation
+ if(ATTENUATION_VAL .and. .not. PARTIAL_PHYS_DISPERSION_ONLY_VAL) then
+ ! do NOT put this is a subroutine, otherwise the call to the subroutine prevents compilers from vectorizing the outer loop
+ ! here we assume that N_SLS == 3 in order to be able to unroll and suppress the loop in order to vectorize the outer loop
+ R_xx_val = R_xx(1,ijk,1,1,ispec)
+ R_yy_val = R_yy(1,ijk,1,1,ispec)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(1,ijk,1,1,ispec)
+ sigma_xz = sigma_xz - R_xz(1,ijk,1,1,ispec)
+ sigma_yz = sigma_yz - R_yz(1,ijk,1,1,ispec)
+
+ R_xx_val = R_xx(2,ijk,1,1,ispec)
+ R_yy_val = R_yy(2,ijk,1,1,ispec)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(2,ijk,1,1,ispec)
+ sigma_xz = sigma_xz - R_xz(2,ijk,1,1,ispec)
+ sigma_yz = sigma_yz - R_yz(2,ijk,1,1,ispec)
+
+ R_xx_val = R_xx(3,ijk,1,1,ispec)
+ R_yy_val = R_yy(3,ijk,1,1,ispec)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(3,ijk,1,1,ispec)
+ sigma_xz = sigma_xz - R_xz(3,ijk,1,1,ispec)
+ sigma_yz = sigma_yz - R_yz(3,ijk,1,1,ispec)
+ endif ! ATTENUATION_VAL
+
+ ! define symmetric components of sigma for gravity
+ sigma_yx = sigma_xy
+ sigma_zx = sigma_xz
+ sigma_zy = sigma_yz
+
+ ! compute non-symmetric terms for gravity
+ if(GRAVITY_VAL) then
+ ! use mesh coordinates to get theta and phi
+ ! x y and z contain r theta and phi
+ iglob = ibool(ijk,1,1,ispec)
+
+ dtheta = dble(ystore(iglob))
+ dphi = dble(zstore(iglob))
+
+ cos_theta = dcos(dtheta)
+ sin_theta = dsin(dtheta)
+ cos_phi = dcos(dphi)
+ sin_phi = dsin(dphi)
+
+ cos_theta_sq = cos_theta*cos_theta
+ sin_theta_sq = sin_theta*sin_theta
+ cos_phi_sq = cos_phi*cos_phi
+ sin_phi_sq = sin_phi*sin_phi
+
+ ! get g, rho and dg/dr=dg
+ ! spherical components of the gravitational acceleration
+ ! for efficiency replace with lookup table every 100 m in radial direction
+ radius = dble(xstore(iglob))
+
+ int_radius = nint(10.d0 * radius * R_EARTH_KM )
+ minus_g = minus_gravity_table(int_radius)
+ minus_dg = minus_deriv_gravity_table(int_radius)
+ rho = density_table(int_radius)
+
+ ! Cartesian components of the gravitational acceleration
+ gxl = minus_g*sin_theta*cos_phi
+ gyl = minus_g*sin_theta*sin_phi
+ gzl = minus_g*cos_theta
+
+ ! Cartesian components of gradient of gravitational acceleration
+ ! obtained from spherical components
+ minus_g_over_radius = minus_g / radius
+ minus_dg_plus_g_over_radius = minus_dg - minus_g_over_radius
+
+ Hxxl = minus_g_over_radius*(cos_phi_sq*cos_theta_sq + sin_phi_sq) + cos_phi_sq*minus_dg*sin_theta_sq
+ Hyyl = minus_g_over_radius*(cos_phi_sq + cos_theta_sq*sin_phi_sq) + minus_dg*sin_phi_sq*sin_theta_sq
+ Hzzl = cos_theta_sq*minus_dg + minus_g_over_radius*sin_theta_sq
+ Hxyl = cos_phi*minus_dg_plus_g_over_radius*sin_phi*sin_theta_sq
+ Hxzl = cos_phi*cos_theta*minus_dg_plus_g_over_radius*sin_theta
+ Hyzl = cos_theta*minus_dg_plus_g_over_radius*sin_phi*sin_theta
+
+ ! get displacement and multiply by density to compute G tensor
+ sx_l = rho * dummyx_loc(ijk,1,1)
+ sy_l = rho * dummyy_loc(ijk,1,1)
+ sz_l = rho * dummyz_loc(ijk,1,1)
+
+ ! compute G tensor from s . g and add to sigma (not symmetric)
+ sigma_xx = sigma_xx + sy_l*gyl + sz_l*gzl
+ sigma_yy = sigma_yy + sx_l*gxl + sz_l*gzl
+ sigma_zz = sigma_zz + sx_l*gxl + sy_l*gyl
+
+ sigma_xy = sigma_xy - sx_l * gyl
+ sigma_yx = sigma_yx - sy_l * gxl
+
+ sigma_xz = sigma_xz - sx_l * gzl
+ sigma_zx = sigma_zx - sz_l * gxl
+
+ sigma_yz = sigma_yz - sy_l * gzl
+ sigma_zy = sigma_zy - sz_l * gyl
+
+ ! precompute vector
+ factor = jacobianl * wgll_cube(ijk,1,1)
+ rho_s_H(1,ijk,1,1) = factor * (sx_l * Hxxl + sy_l * Hxyl + sz_l * Hxzl)
+ rho_s_H(2,ijk,1,1) = factor * (sx_l * Hxyl + sy_l * Hyyl + sz_l * Hyzl)
+ rho_s_H(3,ijk,1,1) = factor * (sx_l * Hxzl + sy_l * Hyzl + sz_l * Hzzl)
+ endif ! end of section with gravity terms
+
+ ! form dot product with test vector, non-symmetric form
+ tempx1(ijk,1,1) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(ijk,1,1) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(ijk,1,1) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
+
+ tempx2(ijk,1,1) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(ijk,1,1) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(ijk,1,1) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+
+ tempx3(ijk,1,1) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(ijk,1,1) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(ijk,1,1) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
+ enddo
+
+ end subroutine compute_element_iso
+
+! end of FORCE_VECTORIZATION
+#endif
+
!--------------------------------------------------------------------------------------------
!
! transversely isotropic element
!
!--------------------------------------------------------------------------------------------
+#ifndef FORCE_VECTORIZATION
+ ! default, without vectorization
+
subroutine compute_element_tiso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
xstore,ystore,zstore, &
@@ -483,7 +825,6 @@
double precision Hxxl,Hyyl,Hzzl,Hxyl,Hxzl,Hyzl
real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
- integer :: ispec_strain
integer :: i,j,k
integer :: int_radius
integer :: iglob
@@ -535,15 +876,21 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
+!ZN from the expression in which we use the strain later in the code.
+
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
- ispec_strain = 1
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ if( ispec == 1) then
+ epsilon_trace_over_3(i,j,k,1) = templ
+ endif
else
- ispec_strain = ispec
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ epsilon_trace_over_3(i,j,k,ispec) = templ
endif
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
@@ -911,14 +1258,545 @@
end subroutine compute_element_tiso
+#else
+ ! FORCE_VECTORIZATION
+
+ subroutine compute_element_tiso(ispec, &
+ minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ xstore,ystore,zstore, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ wgll_cube, &
+ kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
+ ibool, &
+ R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilon_trace_over_3, &
+ one_minus_sum_beta,vnspec, &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
+ dummyx_loc,dummyy_loc,dummyz_loc, &
+ epsilondev_loc,rho_s_H,is_backward_field)
+
+! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
+
+ use constants_solver
+ use specfem_par,only: COMPUTE_AND_STORE_STRAIN
+
+ implicit none
+
+ ! element id
+ integer :: ispec
+
+ ! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool
+
+ ! x y and z contain r theta and phi
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: xstore,ystore,zstore
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ ! array with derivatives of Lagrange polynomials and precalculated products
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+
+ ! store anisotropic properties only where needed to save memory
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_TISO_MANTLE) :: &
+ kappahstore,muhstore,eta_anisostore
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: &
+ kappavstore,muvstore
+
+ ! attenuation
+ ! memory variables for attenuation
+ ! memory variables R_ij are stored at the local rather than global level
+ ! to allow for optimization of cache access by compiler
+ real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUATION) :: R_xx,R_yy,R_xy,R_xz,R_yz
+
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
+
+ ! variable sized array variables
+ integer :: vnspec
+
+ ! [alpha,beta,gamma]val reduced to N_SLS to N_SLS*NUM_NODES
+ real(kind=CUSTOM_REAL), dimension(ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: one_minus_sum_beta
+
+ ! gravity
+ double precision, dimension(NRAD_GRAVITY) :: minus_gravity_table,density_table,minus_deriv_gravity_table
+
+ ! element info
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+
+ logical :: is_backward_field
+
+! local parameters
+ real(kind=CUSTOM_REAL) one_minus_sum_beta_use
+ ! the 21 coefficients for an anisotropic medium in reduced notation
+ real(kind=CUSTOM_REAL) c11,c22,c33,c44,c55,c66,c12,c13,c23,c14,c24,c34,c15,c25,c35,c45,c16,c26,c36,c46,c56
+
+ real(kind=CUSTOM_REAL) rhovphsq,sinphifour,cosphisq,sinphisq,costhetasq,rhovsvsq,sinthetasq, &
+ cosphifour,costhetafour,rhovpvsq,sinthetafour,rhovshsq,cosfourphi, &
+ costwotheta,cosfourtheta,sintwophisq,costheta,sinphi,sintheta,cosphi, &
+ sintwotheta,costwophi,sintwophi,costwothetasq,costwophisq,phi,theta
+
+ real(kind=CUSTOM_REAL) two_rhovsvsq,two_rhovshsq ! two_rhovpvsq,two_rhovphsq
+ real(kind=CUSTOM_REAL) four_rhovsvsq,four_rhovshsq ! four_rhovpvsq,four_rhovphsq
+
+ real(kind=CUSTOM_REAL) twoetaminone,etaminone,eta_aniso
+ real(kind=CUSTOM_REAL) two_eta_aniso,four_eta_aniso,six_eta_aniso
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+ real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+
+ real(kind=CUSTOM_REAL) templ
+ real(kind=CUSTOM_REAL) templ1,templ1_cos,templ2,templ2_cos,templ3,templ3_two,templ3_cos
+ real(kind=CUSTOM_REAL) kappavl,kappahl,muvl,muhl
+
+ ! for gravity
+ double precision dphi,dtheta
+ double precision radius,rho,minus_g,minus_dg
+ double precision minus_g_over_radius,minus_dg_plus_g_over_radius
+ double precision cos_theta,sin_theta,cos_phi,sin_phi
+ double precision cos_theta_sq,sin_theta_sq,cos_phi_sq,sin_phi_sq
+ double precision factor,sx_l,sy_l,sz_l,gxl,gyl,gzl
+ double precision Hxxl,Hyyl,Hzzl,Hxyl,Hxzl,Hyzl
+ real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
+
+ integer :: int_radius
+ integer :: iglob
+
+ logical :: dummyl
+
+ ! force vectorization
+ integer :: ijk
+ real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
+
+! in this vectorized version we have to assume that N_SLS == 3 in order to be able to unroll and thus suppress
+! an inner loop that would otherwise prevent vectorization; this is safe in practice in all cases because N_SLS == 3
+! in all known applications, and in the main program we check that N_SLS == 3 if FORCE_VECTORIZATION is used and we stop
+
+!daniel debug
+ ! dummy to avoid compiler warning
+ dummyl = is_backward_field
+
+ do ijk=1,NGLLCUBE
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(ijk,1,1,ispec)
+ xiyl = xiy(ijk,1,1,ispec)
+ xizl = xiz(ijk,1,1,ispec)
+ etaxl = etax(ijk,1,1,ispec)
+ etayl = etay(ijk,1,1,ispec)
+ etazl = etaz(ijk,1,1,ispec)
+ gammaxl = gammax(ijk,1,1,ispec)
+ gammayl = gammay(ijk,1,1,ispec)
+ gammazl = gammaz(ijk,1,1,ispec)
+
+ ! compute the jacobian
+ jacobianl = 1.0_CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1(ijk,1,1) + etaxl*tempx2(ijk,1,1) + gammaxl*tempx3(ijk,1,1)
+ duxdyl = xiyl*tempx1(ijk,1,1) + etayl*tempx2(ijk,1,1) + gammayl*tempx3(ijk,1,1)
+ duxdzl = xizl*tempx1(ijk,1,1) + etazl*tempx2(ijk,1,1) + gammazl*tempx3(ijk,1,1)
+
+ duydxl = xixl*tempy1(ijk,1,1) + etaxl*tempy2(ijk,1,1) + gammaxl*tempy3(ijk,1,1)
+ duydyl = xiyl*tempy1(ijk,1,1) + etayl*tempy2(ijk,1,1) + gammayl*tempy3(ijk,1,1)
+ duydzl = xizl*tempy1(ijk,1,1) + etazl*tempy2(ijk,1,1) + gammazl*tempy3(ijk,1,1)
+
+ duzdxl = xixl*tempz1(ijk,1,1) + etaxl*tempz2(ijk,1,1) + gammaxl*tempz3(ijk,1,1)
+ duzdyl = xiyl*tempz1(ijk,1,1) + etayl*tempz2(ijk,1,1) + gammayl*tempz3(ijk,1,1)
+ duzdzl = xizl*tempz1(ijk,1,1) + etazl*tempz2(ijk,1,1) + gammazl*tempz3(ijk,1,1)
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
+!ZN from the expression in which we use the strain later in the code.
+
+ ! compute deviatoric strain
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+ if( ispec == 1) then
+ epsilon_trace_over_3(ijk,1,1,1) = templ
+ endif
+ else
+ epsilon_trace_over_3(ijk,1,1,ispec) = templ
+ endif
+ epsilondev_loc(1,ijk,1,1) = duxdxl - templ
+ epsilondev_loc(2,ijk,1,1) = duydyl - templ
+ epsilondev_loc(3,ijk,1,1) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
+ epsilondev_loc(4,ijk,1,1) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
+ epsilondev_loc(5,ijk,1,1) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
+ endif
+
+ !
+ ! compute either isotropic or anisotropic elements
+ !
+
+! note: the mesh is built such that anisotropic elements are created first in anisotropic layers,
+! thus they are listed first ( see in create_regions_mesh.f90: perm_layer() ordering )
+! this is therefore still in bounds of 1:NSPECMAX_TISO_MANTLE even if NSPECMAX_TISO is less than NSPEC
+
+ ! use kappa and mu from transversely isotropic model
+ kappavl = kappavstore(ijk,1,1,ispec)
+ muvl = muvstore(ijk,1,1,ispec)
+
+ kappahl = kappahstore(ijk,1,1,ispec)
+ muhl = muhstore(ijk,1,1,ispec)
+
+ ! use unrelaxed parameters if attenuation
+ ! eta does not need to be shifted since it is a ratio
+ if(ATTENUATION_VAL) then
+ ! precompute terms for attenuation if needed
+ if( ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL ) then
+ one_minus_sum_beta_use = one_minus_sum_beta(ijk,1,1,ispec)
+ else
+ one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,ispec)
+ endif
+ muvl = muvl * one_minus_sum_beta_use
+ muhl = muhl * one_minus_sum_beta_use
+ endif
+
+ rhovpvsq = kappavl + FOUR_THIRDS * muvl !!! that is C
+ rhovphsq = kappahl + FOUR_THIRDS * muhl !!! that is A
+
+ rhovsvsq = muvl !!! that is L
+ rhovshsq = muhl !!! that is N
+
+ eta_aniso = eta_anisostore(ijk,1,1,ispec) !!! that is F / (A - 2 L)
+
+ ! use mesh coordinates to get theta and phi
+ ! ystore and zstore contain theta and phi
+ iglob = ibool(ijk,1,1,ispec)
+
+ theta = ystore(iglob)
+ phi = zstore(iglob)
+
+ ! precompute some products to reduce the CPU time
+
+ costheta = cos(theta)
+ sintheta = sin(theta)
+ cosphi = cos(phi)
+ sinphi = sin(phi)
+
+ costhetasq = costheta * costheta
+ sinthetasq = sintheta * sintheta
+ cosphisq = cosphi * cosphi
+ sinphisq = sinphi * sinphi
+
+ costhetafour = costhetasq * costhetasq
+ sinthetafour = sinthetasq * sinthetasq
+ cosphifour = cosphisq * cosphisq
+ sinphifour = sinphisq * sinphisq
+
+ costwotheta = cos(2.0_CUSTOM_REAL*theta)
+ sintwotheta = sin(2.0_CUSTOM_REAL*theta)
+ costwophi = cos(2.0_CUSTOM_REAL*phi)
+ sintwophi = sin(2.0_CUSTOM_REAL*phi)
+
+ cosfourtheta = cos(4.0_CUSTOM_REAL*theta)
+ cosfourphi = cos(4.0_CUSTOM_REAL*phi)
+
+ costwothetasq = costwotheta * costwotheta
+
+ costwophisq = costwophi * costwophi
+ sintwophisq = sintwophi * sintwophi
+
+ etaminone = eta_aniso - 1.0_CUSTOM_REAL
+ twoetaminone = 2.0_CUSTOM_REAL * eta_aniso - 1.0_CUSTOM_REAL
+
+ ! precompute some products to reduce the CPU time
+ two_eta_aniso = 2.0_CUSTOM_REAL*eta_aniso
+ four_eta_aniso = 4.0_CUSTOM_REAL*eta_aniso
+ six_eta_aniso = 6.0_CUSTOM_REAL*eta_aniso
+
+ two_rhovsvsq = 2.0_CUSTOM_REAL*rhovsvsq
+ two_rhovshsq = 2.0_CUSTOM_REAL*rhovshsq
+ four_rhovsvsq = 4.0_CUSTOM_REAL*rhovsvsq
+ four_rhovshsq = 4.0_CUSTOM_REAL*rhovshsq
+
+ ! pre-compute temporary values
+ templ1 = four_rhovsvsq - rhovpvsq + twoetaminone*rhovphsq - four_eta_aniso*rhovsvsq
+ templ1_cos = rhovphsq - rhovpvsq + costwotheta*templ1
+ templ2 = four_rhovsvsq - rhovpvsq - rhovphsq + two_eta_aniso*rhovphsq - four_eta_aniso*rhovsvsq
+ templ2_cos = rhovpvsq - rhovphsq + costwotheta*templ2
+ templ3 = rhovphsq + rhovpvsq - two_eta_aniso*rhovphsq + four_eta_aniso*rhovsvsq
+ templ3_two = templ3 - two_rhovshsq - two_rhovsvsq
+ templ3_cos = templ3_two + costwotheta*templ2
+
+ ! reordering operations to facilitate compilation, avoiding divisions, using locality for temporary values
+ c11 = rhovphsq*sinphifour &
+ + 2.0_CUSTOM_REAL*cosphisq*sinphisq* &
+ ( rhovphsq*costhetasq + sinthetasq*(eta_aniso*rhovphsq + two_rhovsvsq - two_eta_aniso*rhovsvsq) ) &
+ + cosphifour*(rhovphsq*costhetafour &
+ + 2.0_CUSTOM_REAL*costhetasq*sinthetasq*(eta_aniso*rhovphsq + two_rhovsvsq - two_eta_aniso*rhovsvsq) &
+ + rhovpvsq*sinthetafour)
+
+ c12 = 0.25_CUSTOM_REAL*costhetasq*(rhovphsq - two_rhovshsq)*(3.0_CUSTOM_REAL + cosfourphi) &
+ - four_rhovshsq*cosphisq*costhetasq*sinphisq &
+ + 0.03125_CUSTOM_REAL*rhovphsq*sintwophisq*(11.0_CUSTOM_REAL + cosfourtheta + 4.0*costwotheta) &
+ + eta_aniso*sinthetasq*(rhovphsq - two_rhovsvsq) &
+ *(cosphifour + sinphifour + 2.0_CUSTOM_REAL*cosphisq*costhetasq*sinphisq) &
+ + rhovpvsq*cosphisq*sinphisq*sinthetafour &
+ - rhovsvsq*sintwophisq*sinthetafour
+
+ c13 = 0.125_CUSTOM_REAL*cosphisq*(rhovphsq + six_eta_aniso*rhovphsq + rhovpvsq - four_rhovsvsq &
+ - 12.0_CUSTOM_REAL*eta_aniso*rhovsvsq + cosfourtheta*templ1) &
+ + sinphisq*(eta_aniso*costhetasq*(rhovphsq - two_rhovsvsq) + sinthetasq*(rhovphsq - two_rhovshsq))
+
+ ! uses temporary templ1 from c13
+ c15 = cosphi*costheta*sintheta* &
+ ( 0.5_CUSTOM_REAL*cosphisq* (rhovpvsq - rhovphsq + costwotheta*templ1) &
+ + etaminone*sinphisq*(rhovphsq - two_rhovsvsq))
+
+ c14 = costheta*sinphi*sintheta* &
+ ( 0.5_CUSTOM_REAL*cosphisq*(templ2_cos + four_rhovshsq - four_rhovsvsq) &
+ + sinphisq*(etaminone*rhovphsq + 2.0_CUSTOM_REAL*(rhovshsq - eta_aniso*rhovsvsq)) )
+
+ ! uses temporary templ2_cos from c14
+ c16 = 0.5_CUSTOM_REAL*cosphi*sinphi*sinthetasq* &
+ ( cosphisq*templ2_cos &
+ + 2.0_CUSTOM_REAL*etaminone*sinphisq*(rhovphsq - two_rhovsvsq) )
+
+ c22 = rhovphsq*cosphifour + 2.0_CUSTOM_REAL*cosphisq*sinphisq* &
+ (rhovphsq*costhetasq + sinthetasq*(eta_aniso*rhovphsq + two_rhovsvsq - two_eta_aniso*rhovsvsq)) &
+ + sinphifour* &
+ (rhovphsq*costhetafour + 2.0_CUSTOM_REAL*costhetasq*sinthetasq*(eta_aniso*rhovphsq &
+ + two_rhovsvsq - two_eta_aniso*rhovsvsq) + rhovpvsq*sinthetafour)
+
+ ! uses temporary templ1 from c13
+ c23 = 0.125_CUSTOM_REAL*sinphisq*(rhovphsq + six_eta_aniso*rhovphsq &
+ + rhovpvsq - four_rhovsvsq - 12.0_CUSTOM_REAL*eta_aniso*rhovsvsq + cosfourtheta*templ1) &
+ + cosphisq*(eta_aniso*costhetasq*(rhovphsq - two_rhovsvsq) + sinthetasq*(rhovphsq - two_rhovshsq))
+
+ ! uses temporary templ1 from c13
+ c24 = costheta*sinphi*sintheta* &
+ ( etaminone*cosphisq*(rhovphsq - two_rhovsvsq) &
+ + 0.5_CUSTOM_REAL*sinphisq*(rhovpvsq - rhovphsq + costwotheta*templ1) )
+
+ ! uses temporary templ2_cos from c14
+ c25 = cosphi*costheta*sintheta* &
+ ( cosphisq*(etaminone*rhovphsq + 2.0_CUSTOM_REAL*(rhovshsq - eta_aniso*rhovsvsq)) &
+ + 0.5_CUSTOM_REAL*sinphisq*(templ2_cos + four_rhovshsq - four_rhovsvsq) )
+
+ ! uses temporary templ2_cos from c14
+ c26 = 0.5_CUSTOM_REAL*cosphi*sinphi*sinthetasq* &
+ ( 2.0_CUSTOM_REAL*etaminone*cosphisq*(rhovphsq - two_rhovsvsq) &
+ + sinphisq*templ2_cos )
+
+ c33 = rhovpvsq*costhetafour &
+ + 2.0_CUSTOM_REAL*costhetasq*sinthetasq*(two_rhovsvsq + eta_aniso*(rhovphsq - two_rhovsvsq)) &
+ + rhovphsq*sinthetafour
+
+ ! uses temporary templ1_cos from c13
+ c34 = - 0.25_CUSTOM_REAL*sinphi*sintwotheta*templ1_cos
+
+ ! uses temporary templ1_cos from c34
+ c35 = - 0.25_CUSTOM_REAL*cosphi*sintwotheta*templ1_cos
+
+ ! uses temporary templ1_cos from c34
+ c36 = - 0.25_CUSTOM_REAL*sintwophi*sinthetasq*(templ1_cos - four_rhovshsq + four_rhovsvsq)
+
+ c44 = cosphisq*(rhovsvsq*costhetasq + rhovshsq*sinthetasq) &
+ + sinphisq*(rhovsvsq*costwothetasq + costhetasq*sinthetasq*templ3)
+
+ ! uses temporary templ3 from c44
+ c46 = - cosphi*costheta*sintheta* &
+ ( cosphisq*(rhovshsq - rhovsvsq) - 0.5_CUSTOM_REAL*sinphisq*templ3_cos )
+
+ ! uses templ3 from c46
+ c45 = 0.25_CUSTOM_REAL*sintwophi*sinthetasq* &
+ (templ3_two + costwotheta*(rhovphsq + rhovpvsq - two_eta_aniso*rhovphsq + 4.0_CUSTOM_REAL*etaminone*rhovsvsq))
+
+ c55 = sinphisq*(rhovsvsq*costhetasq + rhovshsq*sinthetasq) &
+ + cosphisq*(rhovsvsq*costwothetasq &
+ + costhetasq*sinthetasq*(rhovphsq - two_eta_aniso*rhovphsq + rhovpvsq + four_eta_aniso*rhovsvsq) )
+
+ ! uses temporary templ3_cos from c46
+ c56 = costheta*sinphi*sintheta* &
+ ( 0.5_CUSTOM_REAL*cosphisq*templ3_cos + sinphisq*(rhovsvsq - rhovshsq) )
+
+ c66 = rhovshsq*costwophisq*costhetasq &
+ - 2.0_CUSTOM_REAL*cosphisq*costhetasq*sinphisq*(rhovphsq - two_rhovshsq) &
+ + 0.03125_CUSTOM_REAL*rhovphsq*sintwophisq*(11.0_CUSTOM_REAL + 4.0_CUSTOM_REAL*costwotheta + cosfourtheta) &
+ - 0.125_CUSTOM_REAL*rhovsvsq*sinthetasq* &
+ ( -6.0_CUSTOM_REAL - 2.0_CUSTOM_REAL*costwotheta - 2.0_CUSTOM_REAL*cosfourphi &
+ + cos(4.0_CUSTOM_REAL*phi - 2.0_CUSTOM_REAL*theta) &
+ + cos(2.0_CUSTOM_REAL*(2.0_CUSTOM_REAL*phi + theta)) ) &
+ + rhovpvsq*cosphisq*sinphisq*sinthetafour &
+ - 0.5_CUSTOM_REAL*eta_aniso*sintwophisq*sinthetafour*(rhovphsq - two_rhovsvsq)
+
+ ! general expression of stress tensor for full Cijkl with 21 coefficients
+ sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
+ c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
+
+ sigma_yy = c12*duxdxl + c26*duxdyl_plus_duydxl + c22*duydyl + &
+ c25*duzdxl_plus_duxdzl + c24*duzdyl_plus_duydzl + c23*duzdzl
+
+ sigma_zz = c13*duxdxl + c36*duxdyl_plus_duydxl + c23*duydyl + &
+ c35*duzdxl_plus_duxdzl + c34*duzdyl_plus_duydzl + c33*duzdzl
+
+ sigma_xy = c16*duxdxl + c66*duxdyl_plus_duydxl + c26*duydyl + &
+ c56*duzdxl_plus_duxdzl + c46*duzdyl_plus_duydzl + c36*duzdzl
+
+ sigma_xz = c15*duxdxl + c56*duxdyl_plus_duydxl + c25*duydyl + &
+ c55*duzdxl_plus_duxdzl + c45*duzdyl_plus_duydzl + c35*duzdzl
+
+ sigma_yz = c14*duxdxl + c46*duxdyl_plus_duydxl + c24*duydyl + &
+ c45*duzdxl_plus_duxdzl + c44*duzdyl_plus_duydzl + c34*duzdzl
+
+ ! subtract memory variables if attenuation
+ if(ATTENUATION_VAL .and. .not. PARTIAL_PHYS_DISPERSION_ONLY_VAL) then
+! do NOT put this is a subroutine, otherwise the call to the subroutine prevents compilers from vectorizing the outer loop
+! here we assume that N_SLS == 3 in order to be able to unroll and suppress the loop in order to vectorize the outer loop
+ R_xx_val = R_xx(1,ijk,1,1,ispec)
+ R_yy_val = R_yy(1,ijk,1,1,ispec)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(1,ijk,1,1,ispec)
+ sigma_xz = sigma_xz - R_xz(1,ijk,1,1,ispec)
+ sigma_yz = sigma_yz - R_yz(1,ijk,1,1,ispec)
+
+ R_xx_val = R_xx(2,ijk,1,1,ispec)
+ R_yy_val = R_yy(2,ijk,1,1,ispec)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(2,ijk,1,1,ispec)
+ sigma_xz = sigma_xz - R_xz(2,ijk,1,1,ispec)
+ sigma_yz = sigma_yz - R_yz(2,ijk,1,1,ispec)
+
+ R_xx_val = R_xx(3,ijk,1,1,ispec)
+ R_yy_val = R_yy(3,ijk,1,1,ispec)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(3,ijk,1,1,ispec)
+ sigma_xz = sigma_xz - R_xz(3,ijk,1,1,ispec)
+ sigma_yz = sigma_yz - R_yz(3,ijk,1,1,ispec)
+ endif ! ATTENUATION_VAL
+
+ ! define symmetric components of sigma for gravity
+ sigma_yx = sigma_xy
+ sigma_zx = sigma_xz
+ sigma_zy = sigma_yz
+
+ ! compute non-symmetric terms for gravity
+ if(GRAVITY_VAL) then
+ ! use mesh coordinates to get theta and phi
+ ! x y and z contain r theta and phi
+ iglob = ibool(ijk,1,1,ispec)
+
+ dtheta = dble(ystore(iglob))
+ dphi = dble(zstore(iglob))
+ radius = dble(xstore(iglob))
+
+ cos_theta = dcos(dtheta)
+ sin_theta = dsin(dtheta)
+ cos_phi = dcos(dphi)
+ sin_phi = dsin(dphi)
+
+ cos_theta_sq = cos_theta*cos_theta
+ sin_theta_sq = sin_theta*sin_theta
+ cos_phi_sq = cos_phi*cos_phi
+ sin_phi_sq = sin_phi*sin_phi
+
+ ! get g, rho and dg/dr=dg
+ ! spherical components of the gravitational acceleration
+ ! for efficiency replace with lookup table every 100 m in radial direction
+ int_radius = nint(10.d0 * radius * R_EARTH_KM )
+ minus_g = minus_gravity_table(int_radius)
+ minus_dg = minus_deriv_gravity_table(int_radius)
+ rho = density_table(int_radius)
+
+ ! Cartesian components of the gravitational acceleration
+ gxl = minus_g*sin_theta*cos_phi
+ gyl = minus_g*sin_theta*sin_phi
+ gzl = minus_g*cos_theta
+
+ ! Cartesian components of gradient of gravitational acceleration
+ ! obtained from spherical components
+ minus_g_over_radius = minus_g / radius
+ minus_dg_plus_g_over_radius = minus_dg - minus_g_over_radius
+
+ Hxxl = minus_g_over_radius*(cos_phi_sq*cos_theta_sq + sin_phi_sq) + cos_phi_sq*minus_dg*sin_theta_sq
+ Hyyl = minus_g_over_radius*(cos_phi_sq + cos_theta_sq*sin_phi_sq) + minus_dg*sin_phi_sq*sin_theta_sq
+ Hzzl = cos_theta_sq*minus_dg + minus_g_over_radius*sin_theta_sq
+ Hxyl = cos_phi*minus_dg_plus_g_over_radius*sin_phi*sin_theta_sq
+ Hxzl = cos_phi*cos_theta*minus_dg_plus_g_over_radius*sin_theta
+ Hyzl = cos_theta*minus_dg_plus_g_over_radius*sin_phi*sin_theta
+
+ ! get displacement and multiply by density to compute G tensor
+ sx_l = rho * dummyx_loc(ijk,1,1)
+ sy_l = rho * dummyy_loc(ijk,1,1)
+ sz_l = rho * dummyz_loc(ijk,1,1)
+
+ ! compute G tensor from s . g and add to sigma (not symmetric)
+ sigma_xx = sigma_xx + sy_l*gyl + sz_l*gzl
+ sigma_yy = sigma_yy + sx_l*gxl + sz_l*gzl
+ sigma_zz = sigma_zz + sx_l*gxl + sy_l*gyl
+
+ sigma_xy = sigma_xy - sx_l * gyl
+ sigma_yx = sigma_yx - sy_l * gxl
+
+ sigma_xz = sigma_xz - sx_l * gzl
+ sigma_zx = sigma_zx - sz_l * gxl
+
+ sigma_yz = sigma_yz - sy_l * gzl
+ sigma_zy = sigma_zy - sz_l * gyl
+
+ ! precompute vector
+ factor = jacobianl * wgll_cube(ijk,1,1)
+ rho_s_H(1,ijk,1,1) = factor * (sx_l * Hxxl + sy_l * Hxyl + sz_l * Hxzl)
+ rho_s_H(2,ijk,1,1) = factor * (sx_l * Hxyl + sy_l * Hyyl + sz_l * Hyzl)
+ rho_s_H(3,ijk,1,1) = factor * (sx_l * Hxzl + sy_l * Hyzl + sz_l * Hzzl)
+ endif ! end of section with gravity terms
+
+ ! form dot product with test vector, non-symmetric form
+ tempx1(ijk,1,1) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(ijk,1,1) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(ijk,1,1) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
+
+ tempx2(ijk,1,1) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(ijk,1,1) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(ijk,1,1) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+
+ tempx3(ijk,1,1) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(ijk,1,1) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(ijk,1,1) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
+ enddo
+
+ end subroutine compute_element_tiso
+
+! end of FORCE_VECTORIZATION
+#endif
+
+
!--------------------------------------------------------------------------------------------
!
! anisotropic element
!
!--------------------------------------------------------------------------------------------
+#ifndef FORCE_VECTORIZATION
+ ! default, without vectorization
+
subroutine compute_element_aniso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
xstore,ystore,zstore, &
@@ -1012,7 +1890,6 @@
double precision Hxxl,Hyyl,Hzzl,Hxyl,Hxzl,Hyzl
real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
- integer :: ispec_strain
integer :: i,j,k
integer :: int_radius
integer :: iglob
@@ -1064,15 +1941,21 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
+!ZN from the expression in which we use the strain later in the code.
+
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
- ispec_strain = 1
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ if( ispec == 1) then
+ epsilon_trace_over_3(i,j,k,1) = templ
+ endif
else
- ispec_strain = ispec
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+ epsilon_trace_over_3(i,j,k,ispec) = templ
endif
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
@@ -1287,6 +2170,374 @@
end subroutine compute_element_aniso
+
+#else
+
+ ! FORCE_VECTORIZATION
+
+ subroutine compute_element_aniso(ispec, &
+ minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ xstore,ystore,zstore, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ wgll_cube, &
+ c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
+ c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
+ c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
+ ibool, &
+ R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilon_trace_over_3, &
+ one_minus_sum_beta,vnspec, &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
+ dummyx_loc,dummyy_loc,dummyz_loc, &
+ epsilondev_loc,rho_s_H,is_backward_field)
+
+ use constants_solver
+ use specfem_par,only: COMPUTE_AND_STORE_STRAIN
+
+ implicit none
+
+ ! element id
+ integer :: ispec
+
+ ! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool
+
+ ! x y and z contain r theta and phi
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: xstore,ystore,zstore
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ ! array with derivatives of Lagrange polynomials and precalculated products
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+
+ ! arrays for full anisotropy only when needed
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ANISO_MANTLE) :: &
+ c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
+ c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
+ c36store,c44store,c45store,c46store,c55store,c56store,c66store
+
+ ! attenuation
+ ! memory variables for attenuation
+ ! memory variables R_ij are stored at the local rather than global level
+ ! to allow for optimization of cache access by compiler
+ real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUATION) :: R_xx,R_yy,R_xy,R_xz,R_yz
+
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
+
+ ! variable sized array variables
+ integer :: vnspec
+
+ ! [alpha,beta,gamma]val reduced to N_SLS to N_SLS*NUM_NODES
+ real(kind=CUSTOM_REAL), dimension(ATT1_VAL,ATT2_VAL,ATT3_VAL,vnspec) :: one_minus_sum_beta
+
+ ! gravity
+ double precision, dimension(NRAD_GRAVITY) :: minus_gravity_table,density_table,minus_deriv_gravity_table
+
+ ! element info
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+
+ logical :: is_backward_field
+
+! local parameters
+ ! real(kind=CUSTOM_REAL) one_minus_sum_beta_use
+ real(kind=CUSTOM_REAL) minus_sum_beta,mul
+ ! the 21 coefficients for an anisotropic medium in reduced notation
+ real(kind=CUSTOM_REAL) c11,c22,c33,c44,c55,c66,c12,c13,c23,c14,c24,c34,c15,c25,c35,c45,c16,c26,c36,c46,c56
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+ real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+
+ real(kind=CUSTOM_REAL) templ
+
+ ! for gravity
+ double precision dphi,dtheta
+ double precision radius,rho,minus_g,minus_dg
+ double precision minus_g_over_radius,minus_dg_plus_g_over_radius
+ double precision cos_theta,sin_theta,cos_phi,sin_phi
+ double precision cos_theta_sq,sin_theta_sq,cos_phi_sq,sin_phi_sq
+ double precision factor,sx_l,sy_l,sz_l,gxl,gyl,gzl
+ double precision Hxxl,Hyyl,Hzzl,Hxyl,Hxzl,Hyzl
+ real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
+
+ integer :: int_radius
+ integer :: iglob
+
+ logical :: dummyl
+
+ ! force vectorization
+ integer :: ijk
+ real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
+
+! in this vectorized version we have to assume that N_SLS == 3 in order to be able to unroll and thus suppress
+! an inner loop that would otherwise prevent vectorization; this is safe in practice in all cases because N_SLS == 3
+! in all known applications, and in the main program we check that N_SLS == 3 if FORCE_VECTORIZATION is used and we stop
+
+!daniel debug
+ ! dummy to avoid compiler warning
+ dummyl = is_backward_field
+
+ ! anisotropic elements
+
+ do ijk=1,NGLLCUBE
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(ijk,1,1,ispec)
+ xiyl = xiy(ijk,1,1,ispec)
+ xizl = xiz(ijk,1,1,ispec)
+ etaxl = etax(ijk,1,1,ispec)
+ etayl = etay(ijk,1,1,ispec)
+ etazl = etaz(ijk,1,1,ispec)
+ gammaxl = gammax(ijk,1,1,ispec)
+ gammayl = gammay(ijk,1,1,ispec)
+ gammazl = gammaz(ijk,1,1,ispec)
+
+ ! compute the jacobian
+ jacobianl = 1.0_CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1(ijk,1,1) + etaxl*tempx2(ijk,1,1) + gammaxl*tempx3(ijk,1,1)
+ duxdyl = xiyl*tempx1(ijk,1,1) + etayl*tempx2(ijk,1,1) + gammayl*tempx3(ijk,1,1)
+ duxdzl = xizl*tempx1(ijk,1,1) + etazl*tempx2(ijk,1,1) + gammazl*tempx3(ijk,1,1)
+
+ duydxl = xixl*tempy1(ijk,1,1) + etaxl*tempy2(ijk,1,1) + gammaxl*tempy3(ijk,1,1)
+ duydyl = xiyl*tempy1(ijk,1,1) + etayl*tempy2(ijk,1,1) + gammayl*tempy3(ijk,1,1)
+ duydzl = xizl*tempy1(ijk,1,1) + etazl*tempy2(ijk,1,1) + gammazl*tempy3(ijk,1,1)
+
+ duzdxl = xixl*tempz1(ijk,1,1) + etaxl*tempz2(ijk,1,1) + gammaxl*tempz3(ijk,1,1)
+ duzdyl = xiyl*tempz1(ijk,1,1) + etayl*tempz2(ijk,1,1) + gammayl*tempz3(ijk,1,1)
+ duzdzl = xizl*tempz1(ijk,1,1) + etazl*tempz2(ijk,1,1) + gammazl*tempz3(ijk,1,1)
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
+!ZN from the expression in which we use the strain later in the code.
+
+ ! compute deviatoric strain
+ if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+ if( ispec == 1) then
+ epsilon_trace_over_3(ijk,1,1,1) = templ
+ endif
+ else
+ epsilon_trace_over_3(ijk,1,1,ispec) = templ
+ endif
+ epsilondev_loc(1,ijk,1,1) = duxdxl - templ
+ epsilondev_loc(2,ijk,1,1) = duydyl - templ
+ epsilondev_loc(3,ijk,1,1) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
+ epsilondev_loc(4,ijk,1,1) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
+ epsilondev_loc(5,ijk,1,1) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
+ endif
+
+ !
+ ! compute anisotropic elements
+ !
+
+ c11 = c11store(ijk,1,1,ispec)
+ c12 = c12store(ijk,1,1,ispec)
+ c13 = c13store(ijk,1,1,ispec)
+ c14 = c14store(ijk,1,1,ispec)
+ c15 = c15store(ijk,1,1,ispec)
+ c16 = c16store(ijk,1,1,ispec)
+ c22 = c22store(ijk,1,1,ispec)
+ c23 = c23store(ijk,1,1,ispec)
+ c24 = c24store(ijk,1,1,ispec)
+ c25 = c25store(ijk,1,1,ispec)
+ c26 = c26store(ijk,1,1,ispec)
+ c33 = c33store(ijk,1,1,ispec)
+ c34 = c34store(ijk,1,1,ispec)
+ c35 = c35store(ijk,1,1,ispec)
+ c36 = c36store(ijk,1,1,ispec)
+ c44 = c44store(ijk,1,1,ispec)
+ c45 = c45store(ijk,1,1,ispec)
+ c46 = c46store(ijk,1,1,ispec)
+ c55 = c55store(ijk,1,1,ispec)
+ c56 = c56store(ijk,1,1,ispec)
+ c66 = c66store(ijk,1,1,ispec)
+
+ if(ATTENUATION_VAL) then
+ ! precompute terms for attenuation if needed
+ if( ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL ) then
+ minus_sum_beta = one_minus_sum_beta(ijk,1,1,ispec) - 1.0_CUSTOM_REAL
+ else
+ minus_sum_beta = one_minus_sum_beta(1,1,1,ispec) - 1.0_CUSTOM_REAL
+ endif
+ mul = c44 * minus_sum_beta
+ c11 = c11 + FOUR_THIRDS * mul ! * minus_sum_beta * mul
+ c12 = c12 - TWO_THIRDS * mul
+ c13 = c13 - TWO_THIRDS * mul
+ c22 = c22 + FOUR_THIRDS * mul
+ c23 = c23 - TWO_THIRDS * mul
+ c33 = c33 + FOUR_THIRDS * mul
+ c44 = c44 + mul
+ c55 = c55 + mul
+ c66 = c66 + mul
+ endif
+
+ sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
+ c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
+
+ sigma_yy = c12*duxdxl + c26*duxdyl_plus_duydxl + c22*duydyl + &
+ c25*duzdxl_plus_duxdzl + c24*duzdyl_plus_duydzl + c23*duzdzl
+
+ sigma_zz = c13*duxdxl + c36*duxdyl_plus_duydxl + c23*duydyl + &
+ c35*duzdxl_plus_duxdzl + c34*duzdyl_plus_duydzl + c33*duzdzl
+
+ sigma_xy = c16*duxdxl + c66*duxdyl_plus_duydxl + c26*duydyl + &
+ c56*duzdxl_plus_duxdzl + c46*duzdyl_plus_duydzl + c36*duzdzl
+
+ sigma_xz = c15*duxdxl + c56*duxdyl_plus_duydxl + c25*duydyl + &
+ c55*duzdxl_plus_duxdzl + c45*duzdyl_plus_duydzl + c35*duzdzl
+
+ sigma_yz = c14*duxdxl + c46*duxdyl_plus_duydxl + c24*duydyl + &
+ c45*duzdxl_plus_duxdzl + c44*duzdyl_plus_duydzl + c34*duzdzl
+
+ ! subtract memory variables if attenuation
+ if(ATTENUATION_VAL .and. .not. PARTIAL_PHYS_DISPERSION_ONLY_VAL) then
+! do NOT put this is a subroutine, otherwise the call to the subroutine prevents compilers from vectorizing the outer loop
+! here we assume that N_SLS == 3 in order to be able to unroll and suppress the loop in order to vectorize the outer loop
+ R_xx_val = R_xx(1,ijk,1,1,ispec)
+ R_yy_val = R_yy(1,ijk,1,1,ispec)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(1,ijk,1,1,ispec)
+ sigma_xz = sigma_xz - R_xz(1,ijk,1,1,ispec)
+ sigma_yz = sigma_yz - R_yz(1,ijk,1,1,ispec)
+
+ R_xx_val = R_xx(2,ijk,1,1,ispec)
+ R_yy_val = R_yy(2,ijk,1,1,ispec)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(2,ijk,1,1,ispec)
+ sigma_xz = sigma_xz - R_xz(2,ijk,1,1,ispec)
+ sigma_yz = sigma_yz - R_yz(2,ijk,1,1,ispec)
+
+ R_xx_val = R_xx(3,ijk,1,1,ispec)
+ R_yy_val = R_yy(3,ijk,1,1,ispec)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(3,ijk,1,1,ispec)
+ sigma_xz = sigma_xz - R_xz(3,ijk,1,1,ispec)
+ sigma_yz = sigma_yz - R_yz(3,ijk,1,1,ispec)
+ endif ! ATTENUATION_VAL
+
+ ! define symmetric components of sigma for gravity
+ sigma_yx = sigma_xy
+ sigma_zx = sigma_xz
+ sigma_zy = sigma_yz
+
+ ! compute non-symmetric terms for gravity
+ if(GRAVITY_VAL) then
+ ! use mesh coordinates to get theta and phi
+ ! x y and z contain r theta and phi
+ iglob = ibool(ijk,1,1,ispec)
+
+ dtheta = dble(ystore(iglob))
+ dphi = dble(zstore(iglob))
+ radius = dble(xstore(iglob))
+
+ cos_theta = dcos(dtheta)
+ sin_theta = dsin(dtheta)
+ cos_phi = dcos(dphi)
+ sin_phi = dsin(dphi)
+
+ cos_theta_sq = cos_theta*cos_theta
+ sin_theta_sq = sin_theta*sin_theta
+ cos_phi_sq = cos_phi*cos_phi
+ sin_phi_sq = sin_phi*sin_phi
+
+ ! get g, rho and dg/dr=dg
+ ! spherical components of the gravitational acceleration
+ ! for efficiency replace with lookup table every 100 m in radial direction
+ int_radius = nint(10.d0 * radius * R_EARTH_KM )
+ minus_g = minus_gravity_table(int_radius)
+ minus_dg = minus_deriv_gravity_table(int_radius)
+ rho = density_table(int_radius)
+
+ ! Cartesian components of the gravitational acceleration
+ gxl = minus_g*sin_theta*cos_phi
+ gyl = minus_g*sin_theta*sin_phi
+ gzl = minus_g*cos_theta
+
+ ! Cartesian components of gradient of gravitational acceleration
+ ! obtained from spherical components
+ minus_g_over_radius = minus_g / radius
+ minus_dg_plus_g_over_radius = minus_dg - minus_g_over_radius
+
+ Hxxl = minus_g_over_radius*(cos_phi_sq*cos_theta_sq + sin_phi_sq) + cos_phi_sq*minus_dg*sin_theta_sq
+ Hyyl = minus_g_over_radius*(cos_phi_sq + cos_theta_sq*sin_phi_sq) + minus_dg*sin_phi_sq*sin_theta_sq
+ Hzzl = cos_theta_sq*minus_dg + minus_g_over_radius*sin_theta_sq
+ Hxyl = cos_phi*minus_dg_plus_g_over_radius*sin_phi*sin_theta_sq
+ Hxzl = cos_phi*cos_theta*minus_dg_plus_g_over_radius*sin_theta
+ Hyzl = cos_theta*minus_dg_plus_g_over_radius*sin_phi*sin_theta
+
+ ! get displacement and multiply by density to compute G tensor
+ sx_l = rho * dummyx_loc(ijk,1,1)
+ sy_l = rho * dummyy_loc(ijk,1,1)
+ sz_l = rho * dummyz_loc(ijk,1,1)
+
+ ! compute G tensor from s . g and add to sigma (not symmetric)
+ sigma_xx = sigma_xx + sy_l*gyl + sz_l*gzl
+ sigma_yy = sigma_yy + sx_l*gxl + sz_l*gzl
+ sigma_zz = sigma_zz + sx_l*gxl + sy_l*gyl
+
+ sigma_xy = sigma_xy - sx_l * gyl
+ sigma_yx = sigma_yx - sy_l * gxl
+
+ sigma_xz = sigma_xz - sx_l * gzl
+ sigma_zx = sigma_zx - sz_l * gxl
+
+ sigma_yz = sigma_yz - sy_l * gzl
+ sigma_zy = sigma_zy - sz_l * gyl
+
+ ! precompute vector
+ factor = jacobianl * wgll_cube(ijk,1,1)
+ rho_s_H(1,ijk,1,1) = factor * (sx_l * Hxxl + sy_l * Hxyl + sz_l * Hxzl)
+ rho_s_H(2,ijk,1,1) = factor * (sx_l * Hxyl + sy_l * Hyyl + sz_l * Hyzl)
+ rho_s_H(3,ijk,1,1) = factor * (sx_l * Hxzl + sy_l * Hyzl + sz_l * Hzzl)
+ endif ! end of section with gravity terms
+
+ ! form dot product with test vector, non-symmetric form
+ tempx1(ijk,1,1) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(ijk,1,1) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(ijk,1,1) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
+
+ tempx2(ijk,1,1) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(ijk,1,1) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(ijk,1,1) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+
+ tempx3(ijk,1,1) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(ijk,1,1) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(ijk,1,1) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
+ enddo
+
+ end subroutine compute_element_aniso
+
+! end of FORCE_VECTORIZATION
+#endif
+
+
!--------------------------------------------------------------------------------------------
!
! helper functions
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element_att_memory.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element_att_memory.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element_att_memory.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -89,14 +89,63 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: factor_common_c44_muv
integer :: i_SLS
+#ifdef FORCE_VECTORIZATION
+ integer :: ijk
+#endif
+
! use Runge-Kutta scheme to march in time
! get coefficients for that standard linear solid
! IMPROVE we use mu_v here even if there is some anisotropy
! IMPROVE we should probably use an average value instead
+#ifdef FORCE_VECTORIZATION
+
do i_SLS = 1,N_SLS
+ ! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
+ if( ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL ) then
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ do ijk=1,NGLLCUBE
+ factor_common_c44_muv(ijk,1,1) = factor_common(i_SLS,ijk,1,1,ispec) * c44store(ijk,1,1,ispec)
+ enddo
+ else
+ do ijk=1,NGLLCUBE
+ factor_common_c44_muv(ijk,1,1) = factor_common(i_SLS,ijk,1,1,ispec) * muvstore(ijk,1,1,ispec)
+ enddo
+ endif
+ else
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ do ijk=1,NGLLCUBE
+ factor_common_c44_muv(ijk,1,1) = factor_common(i_SLS,1,1,1,ispec) * c44store(ijk,1,1,ispec)
+ enddo
+ else
+ do ijk=1,NGLLCUBE
+ factor_common_c44_muv(ijk,1,1) = factor_common(i_SLS,1,1,1,ispec) * muvstore(ijk,1,1,ispec)
+ enddo
+ endif
+ endif
+ do ijk=1,NGLLCUBE
+ R_xx(i_SLS,ijk,1,1,ispec) = alphaval(i_SLS) * R_xx(i_SLS,ijk,1,1,ispec) + factor_common_c44_muv(ijk,1,1) * &
+ (betaval(i_SLS) * epsilondev_xx(ijk,1,1,ispec) + gammaval(i_SLS) * epsilondev_loc(1,ijk,1,1))
+
+ R_yy(i_SLS,ijk,1,1,ispec) = alphaval(i_SLS) * R_yy(i_SLS,ijk,1,1,ispec) + factor_common_c44_muv(ijk,1,1) * &
+ (betaval(i_SLS) * epsilondev_yy(ijk,1,1,ispec) + gammaval(i_SLS) * epsilondev_loc(2,ijk,1,1))
+
+ R_xy(i_SLS,ijk,1,1,ispec) = alphaval(i_SLS) * R_xy(i_SLS,ijk,1,1,ispec) + factor_common_c44_muv(ijk,1,1) * &
+ (betaval(i_SLS) * epsilondev_xy(ijk,1,1,ispec) + gammaval(i_SLS) * epsilondev_loc(3,ijk,1,1))
+
+ R_xz(i_SLS,ijk,1,1,ispec) = alphaval(i_SLS) * R_xz(i_SLS,ijk,1,1,ispec) + factor_common_c44_muv(ijk,1,1) * &
+ (betaval(i_SLS) * epsilondev_xz(ijk,1,1,ispec) + gammaval(i_SLS) * epsilondev_loc(4,ijk,1,1))
+
+ R_yz(i_SLS,ijk,1,1,ispec) = alphaval(i_SLS) * R_yz(i_SLS,ijk,1,1,ispec) + factor_common_c44_muv(ijk,1,1) * &
+ (betaval(i_SLS) * epsilondev_yz(ijk,1,1,ispec) + gammaval(i_SLS) * epsilondev_loc(5,ijk,1,1))
+ enddo
+ enddo ! i_SLS
+
+#else
+
+ do i_SLS = 1,N_SLS
! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
if( ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL ) then
if(ANISOTROPIC_3D_MANTLE_VAL) then
@@ -129,6 +178,8 @@
enddo ! i_SLS
+#endif
+
end subroutine compute_element_att_memory_cm
@@ -312,6 +363,10 @@
integer :: i_SLS
+#ifdef FORCE_VECTORIZATION
+ integer :: ijk
+#endif
+
! use Runge-Kutta scheme to march in time
! get coefficients for that standard linear solid
@@ -321,8 +376,41 @@
! note: epsilondev_loc is calculated based on displ( n + 1 ), thus corresponds to strain at time (n + 1)
! epsilondev_xx,.. are stored from previous step, thus corresponds now to strain at time n
+#ifdef FORCE_VECTORIZATION
+
do i_SLS = 1,N_SLS
+ ! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
+ if( ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL ) then
+ do ijk=1,NGLLCUBE
+ factor_common_use(ijk,1,1) = factor_common(i_SLS,ijk,1,1,ispec) * muvstore(ijk,1,1,ispec)
+ enddo
+ else
+ do ijk=1,NGLLCUBE
+ factor_common_use(ijk,1,1) = factor_common(i_SLS,1,1,1,ispec) * muvstore(ijk,1,1,ispec)
+ enddo
+ endif
+ do ijk=1,NGLLCUBE
+ R_xx(i_SLS,ijk,1,1,ispec) = alphaval(i_SLS) * R_xx(i_SLS,ijk,1,1,ispec) + factor_common_use(ijk,1,1) * &
+ (betaval(i_SLS) * epsilondev_xx(ijk,1,1,ispec) + gammaval(i_SLS) * epsilondev_loc(1,ijk,1,1))
+
+ R_yy(i_SLS,ijk,1,1,ispec) = alphaval(i_SLS) * R_yy(i_SLS,ijk,1,1,ispec) + factor_common_use(ijk,1,1) * &
+ (betaval(i_SLS) * epsilondev_yy(ijk,1,1,ispec) + gammaval(i_SLS) * epsilondev_loc(2,ijk,1,1))
+
+ R_xy(i_SLS,ijk,1,1,ispec) = alphaval(i_SLS) * R_xy(i_SLS,ijk,1,1,ispec) + factor_common_use(ijk,1,1) * &
+ (betaval(i_SLS) * epsilondev_xy(ijk,1,1,ispec) + gammaval(i_SLS) * epsilondev_loc(3,ijk,1,1))
+
+ R_xz(i_SLS,ijk,1,1,ispec) = alphaval(i_SLS) * R_xz(i_SLS,ijk,1,1,ispec) + factor_common_use(ijk,1,1) * &
+ (betaval(i_SLS) * epsilondev_xz(ijk,1,1,ispec) + gammaval(i_SLS) * epsilondev_loc(4,ijk,1,1))
+
+ R_yz(i_SLS,ijk,1,1,ispec) = alphaval(i_SLS) * R_yz(i_SLS,ijk,1,1,ispec) + factor_common_use(ijk,1,1) * &
+ (betaval(i_SLS) * epsilondev_yz(ijk,1,1,ispec) + gammaval(i_SLS) * epsilondev_loc(5,ijk,1,1))
+ enddo
+ enddo
+
+#else
+
+ do i_SLS = 1,N_SLS
! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
if( ATTENUATION_3D_VAL .or. ATTENUATION_1D_WITH_3D_STORAGE_VAL ) then
factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
@@ -330,12 +418,6 @@
factor_common_use(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
endif
-! do i_memory = 1,5
-! R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
-! + muvstore(:,:,:,ispec) * factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
-! enddo
-
R_xx(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
(betaval(i_SLS) * epsilondev_xx(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(1,:,:,:))
@@ -350,9 +432,10 @@
R_yz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
(betaval(i_SLS) * epsilondev_yz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(5,:,:,:))
-
enddo
+#endif
+
end subroutine compute_element_att_memory_ic
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element_strain.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element_strain.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element_strain.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -432,8 +432,8 @@
! strain separated into single xx,yy,xy,xz,yz-component arrays
!
!--------------------------------------------------------------------------------------------
-!
+
subroutine compute_element_strain_att_Dev(ispec,nglob,nspec, &
displ,veloc,deltat, &
ibool, &
@@ -444,7 +444,7 @@
epsilondev_xy_loc_nplus1, &
epsilondev_xz_loc_nplus1, &
epsilondev_yz_loc_nplus1, &
- eps_trace_over_3_loc_nplus1)
+ nspec_strain_only,eps_trace_over_3_loc_nplus1)
use constants
@@ -460,14 +460,15 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
- !real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_yy_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xy_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xz_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_yz_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_nplus1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: epsilondev_xx_loc_nplus1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: epsilondev_yy_loc_nplus1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: epsilondev_xy_loc_nplus1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: epsilondev_xz_loc_nplus1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: epsilondev_yz_loc_nplus1
+ integer :: nspec_strain_only
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_strain_only) :: eps_trace_over_3_loc_nplus1
+
! local variable
integer :: i,j,k,iglob
real(kind=CUSTOM_REAL) :: templ
@@ -630,12 +631,18 @@
duzdyl_plus_duydzl = duzdyl + duydzl
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- eps_trace_over_3_loc_nplus1 = templ
- epsilondev_xx_loc_nplus1(ijk,1,1) = duxdxl - templ
- epsilondev_yy_loc_nplus1(ijk,1,1) = duydyl - templ
- epsilondev_xy_loc_nplus1(ijk,1,1) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
- epsilondev_xz_loc_nplus1(ijk,1,1) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
- epsilondev_yz_loc_nplus1(ijk,1,1) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
+ if( nspec_strain_only == 1 ) then
+ if( ispec == 1 ) then
+ eps_trace_over_3_loc_nplus1(ijk,1,1,1) = templ
+ endif
+ else
+ eps_trace_over_3_loc_nplus1(ijk,1,1,ispec) = templ
+ endif
+ epsilondev_xx_loc_nplus1(ijk,1,1,ispec) = duxdxl - templ
+ epsilondev_yy_loc_nplus1(ijk,1,1,ispec) = duydyl - templ
+ epsilondev_xy_loc_nplus1(ijk,1,1,ispec) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
+ epsilondev_xz_loc_nplus1(ijk,1,1,ispec) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
+ epsilondev_yz_loc_nplus1(ijk,1,1,ispec) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
enddo
#else
do k=1,NGLLZ
@@ -678,12 +685,18 @@
duzdyl_plus_duydzl = duzdyl + duydzl
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- eps_trace_over_3_loc_nplus1 = templ
- epsilondev_xx_loc_nplus1(i,j,k) = duxdxl - templ
- epsilondev_yy_loc_nplus1(i,j,k) = duydyl - templ
- epsilondev_xy_loc_nplus1(i,j,k) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
- epsilondev_xz_loc_nplus1(i,j,k) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
- epsilondev_yz_loc_nplus1(i,j,k) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
+ if( nspec_strain_only == 1 ) then
+ if( ispec == 1 ) then
+ eps_trace_over_3_loc_nplus1(i,j,k,1) = templ
+ endif
+ else
+ eps_trace_over_3_loc_nplus1(i,j,k,ispec) = templ
+ endif
+ epsilondev_xx_loc_nplus1(i,j,k,ispec) = duxdxl - templ
+ epsilondev_yy_loc_nplus1(i,j,k,ispec) = duydyl - templ
+ epsilondev_xy_loc_nplus1(i,j,k,ispec) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
+ epsilondev_xz_loc_nplus1(i,j,k,ispec) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
+ epsilondev_yz_loc_nplus1(i,j,k,ispec) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
enddo
enddo
enddo
@@ -706,7 +719,7 @@
epsilondev_xy_loc_nplus1, &
epsilondev_xz_loc_nplus1, &
epsilondev_yz_loc_nplus1, &
- eps_trace_over_3_loc_nplus1)
+ nspec_strain_only,eps_trace_over_3_loc_nplus1)
use constants
@@ -728,13 +741,14 @@
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
!real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_yy_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xy_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xz_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_yz_loc_nplus1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: epsilondev_xx_loc_nplus1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: epsilondev_yy_loc_nplus1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: epsilondev_xy_loc_nplus1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: epsilondev_xz_loc_nplus1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: epsilondev_yz_loc_nplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_nplus1
+ integer :: nspec_strain_only
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_strain_only) :: eps_trace_over_3_loc_nplus1
! local variable
integer :: i,j,k,l,iglob
@@ -837,12 +851,18 @@
duzdyl_plus_duydzl = duzdyl + duydzl
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- eps_trace_over_3_loc_nplus1 = templ
- epsilondev_xx_loc_nplus1(i,j,k) = duxdxl - templ
- epsilondev_yy_loc_nplus1(i,j,k) = duydyl - templ
- epsilondev_xy_loc_nplus1(i,j,k) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
- epsilondev_xz_loc_nplus1(i,j,k) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
- epsilondev_yz_loc_nplus1(i,j,k) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
+ if( nspec_strain_only == 1 ) then
+ if( ispec == 1 ) then
+ eps_trace_over_3_loc_nplus1(i,j,k,1) = templ
+ endif
+ else
+ eps_trace_over_3_loc_nplus1(i,j,k,ispec) = templ
+ endif
+ epsilondev_xx_loc_nplus1(i,j,k,ispec) = duxdxl - templ
+ epsilondev_yy_loc_nplus1(i,j,k,ispec) = duydyl - templ
+ epsilondev_xy_loc_nplus1(i,j,k,ispec) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
+ epsilondev_xz_loc_nplus1(i,j,k,ispec) = 0.5_CUSTOM_REAL * duzdxl_plus_duxdzl
+ epsilondev_yz_loc_nplus1(i,j,k,ispec) = 0.5_CUSTOM_REAL * duzdyl_plus_duydzl
enddo
enddo
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_acoustic_calling_routine.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -35,6 +35,7 @@
use specfem_par_innercore,only: displ_inner_core, &
ibool_inner_core,ibelm_top_inner_core
use specfem_par_outercore
+ use specfem_par_movie,only: div_displ_outer_core
implicit none
! local parameters
@@ -122,9 +123,9 @@
!--- couple with mantle at the top of the outer core
!---
if(ACTUALLY_COUPLE_FLUID_CMB) &
- call compute_coupling_fluid_CMB(displ_crust_mantle, &
+ call compute_coupling_fluid_CMB(NGLOB_CRUST_MANTLE,displ_crust_mantle, &
ibool_crust_mantle,ibelm_bottom_crust_mantle, &
- accel_outer_core, &
+ NGLOB_OUTER_CORE,accel_outer_core, &
normal_top_outer_core,jacobian2D_top_outer_core, &
wgllwgll_xy,ibool_outer_core,ibelm_top_outer_core, &
NSPEC2D_TOP(IREGION_OUTER_CORE))
@@ -133,9 +134,9 @@
!--- couple with inner core at the bottom of the outer core
!---
if(ACTUALLY_COUPLE_FLUID_ICB) &
- call compute_coupling_fluid_ICB(displ_inner_core, &
+ call compute_coupling_fluid_ICB(NGLOB_INNER_CORE,displ_inner_core, &
ibool_inner_core,ibelm_top_inner_core, &
- accel_outer_core, &
+ NGLOB_OUTER_CORE,accel_outer_core, &
normal_bottom_outer_core,jacobian2D_bottom_outer_core, &
wgllwgll_xy,ibool_outer_core,ibelm_bottom_outer_core, &
NSPEC2D_BOTTOM(IREGION_OUTER_CORE))
@@ -241,6 +242,7 @@
use specfem_par_innercore,only: b_displ_inner_core, &
ibool_inner_core,ibelm_top_inner_core
use specfem_par_outercore
+ use specfem_par_movie,only: div_displ_outer_core
implicit none
! local parameters
@@ -313,14 +315,14 @@
b_A_array_rotation,b_B_array_rotation, &
b_A_array_rotation_lddrk,b_B_array_rotation_lddrk, &
b_displ_outer_core,b_accel_outer_core, &
- b_div_displ_outer_core,phase_is_inner)
+ div_displ_outer_core,phase_is_inner)
else
call compute_forces_outer_core(b_time,b_deltat,b_two_omega_earth, &
NSPEC_OUTER_CORE_ROT_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
b_A_array_rotation,b_B_array_rotation, &
b_A_array_rotation_lddrk,b_B_array_rotation_lddrk, &
b_displ_outer_core,b_accel_outer_core, &
- b_div_displ_outer_core,phase_is_inner)
+ div_displ_outer_core,phase_is_inner)
endif
else
! on GPU
@@ -351,9 +353,9 @@
!--- couple with mantle at the top of the outer core
!---
if(ACTUALLY_COUPLE_FLUID_CMB) &
- call compute_coupling_fluid_CMB(b_displ_crust_mantle, &
+ call compute_coupling_fluid_CMB(NGLOB_CRUST_MANTLE_ADJOINT,b_displ_crust_mantle, &
ibool_crust_mantle,ibelm_bottom_crust_mantle, &
- b_accel_outer_core, &
+ NGLOB_OUTER_CORE_ADJOINT,b_accel_outer_core, &
normal_top_outer_core,jacobian2D_top_outer_core, &
wgllwgll_xy,ibool_outer_core,ibelm_top_outer_core, &
NSPEC2D_TOP(IREGION_OUTER_CORE))
@@ -362,9 +364,9 @@
!--- couple with inner core at the bottom of the outer core
!---
if(ACTUALLY_COUPLE_FLUID_ICB) &
- call compute_coupling_fluid_ICB(b_displ_inner_core, &
+ call compute_coupling_fluid_ICB(NGLOB_INNER_CORE_ADJOINT,b_displ_inner_core, &
ibool_inner_core,ibelm_top_inner_core, &
- b_accel_outer_core, &
+ NGLOB_OUTER_CORE_ADJOINT,b_accel_outer_core, &
normal_bottom_outer_core,jacobian2D_bottom_outer_core, &
wgllwgll_xy,ibool_outer_core,ibelm_bottom_outer_core, &
NSPEC2D_BOTTOM(IREGION_OUTER_CORE))
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -115,7 +115,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
- integer ispec,iglob,ispec_strain
+ integer ispec,iglob
integer i,j,k,l
integer i_SLS
@@ -150,6 +150,7 @@
real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
+ real(kind=CUSTOM_REAL) templ
! for gravity
integer int_radius
@@ -264,14 +265,16 @@
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
- ispec_strain = 1
+ if( ispec == 1 ) then
+ epsilon_trace_over_3(i,j,k,1) = templ
+ endif
else
- ispec_strain = ispec
+ epsilon_trace_over_3(i,j,k,ispec) = templ
endif
- epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilondev_loc(1,i,j,k) = duxdxl - epsilon_trace_over_3(i,j,k,ispec_strain)
- epsilondev_loc(2,i,j,k) = duydyl - epsilon_trace_over_3(i,j,k,ispec_strain)
+ epsilondev_loc(1,i,j,k) = duxdxl - templ
+ epsilondev_loc(2,i,j,k) = duydyl - templ
epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -175,7 +175,7 @@
real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
integer :: int_radius
- integer :: ispec,ispec_strain
+ integer :: ispec
integer :: i,j,k
integer :: iglob
! integer :: computed_elements
@@ -336,12 +336,13 @@
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
if(NSPEC_INNER_CORE_STRAIN_ONLY == 1) then
- ispec_strain = 1
+ if( ispec == 1) then
+ epsilon_trace_over_3(i,j,k,1) = templ
+ endif
else
- ispec_strain = ispec
+ epsilon_trace_over_3(i,j,k,ispec) = templ
endif
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
epsilondev_loc(3,i,j,k) = 0.5_CUSTOM_REAL * duxdyl_plus_duydxl
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -103,7 +103,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
- integer :: ispec,iglob,ispec_strain
+ integer :: ispec,iglob
integer :: i,j,k,l
real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
@@ -126,6 +126,7 @@
real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
+ real(kind=CUSTOM_REAL) templ
real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
integer :: i_SLS
@@ -245,14 +246,16 @@
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
if(NSPEC_INNER_CORE_STRAIN_ONLY == 1) then
- ispec_strain = 1
+ if( ispec == 1 ) then
+ epsilon_trace_over_3(i,j,k,1) = templ
+ endif
else
- ispec_strain = ispec
+ epsilon_trace_over_3(i,j,k,ispec) = templ
endif
- epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilondev_loc(1,i,j,k) = duxdxl - epsilon_trace_over_3(i,j,k,ispec_strain)
- epsilondev_loc(2,i,j,k) = duydyl - epsilon_trace_over_3(i,j,k,ispec_strain)
+ epsilondev_loc(1,i,j,k) = duxdxl - templ
+ epsilondev_loc(2,i,j,k) = duydyl - templ
epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_Dev.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_Dev.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -75,7 +75,7 @@
real(kind=CUSTOM_REAL), dimension(NGLOB) :: displfluid,accelfluid
! divergence of displacement
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ADJOINT) :: div_displfluid
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_3DMOVIE) :: div_displfluid
! inner/outer element run flag
logical :: phase_is_inner
@@ -137,7 +137,9 @@
! big loop over all spectral elements in the fluid
! ****************************************************
- if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. (.not. phase_is_inner) ) div_displfluid(:,:,:,:) = 0._CUSTOM_REAL
+ if( MOVIE_VOLUME .and. NSPEC_OUTER_CORE_3DMOVIE /= 1 .and. (.not. phase_is_inner) ) then
+ div_displfluid(:,:,:,:) = 0._CUSTOM_REAL
+ endif
! computed_elements = 0
if( .not. phase_is_inner ) then
@@ -365,7 +367,8 @@
dpotentialdy_with_rot = dpotentialdy_with_rot + displ_times_grad_y_ln_rho(i,j,k)
dpotentialdzl = dpotentialdzl + displ_times_grad_z_ln_rho(i,j,k)
- else ! if gravity is turned on
+ else
+ ! if gravity is turned on
! compute divergence of displacment
gxl = temp_gxl(i,j,k)
@@ -392,7 +395,7 @@
! note: these calculations are only considered for SIMULATION_TYPE == 1 .and. SAVE_FORWARD
! and one has set MOVIE_VOLUME_TYPE == 4 when MOVIE_VOLUME is .true.;
! in case of SIMULATION_TYPE == 3, it gets overwritten by compute_kernels_outer_core()
- if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. MOVIE_VOLUME) then
+ if( MOVIE_VOLUME .and. NSPEC_OUTER_CORE_3DMOVIE /= 1 ) then
div_displfluid(i,j,k,ispec) = &
minus_rho_g_over_kappa_fluid(int_radius) &
* (dpotentialdx_with_rot * gxl &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_noDev.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_outer_core_noDev.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -67,7 +67,7 @@
real(kind=CUSTOM_REAL), dimension(NGLOB) :: displfluid,accelfluid
! divergence of displacement
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ADJOINT) :: div_displfluid
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_3DMOVIE) :: div_displfluid
! inner/outer element run flag
logical :: phase_is_inner
@@ -102,7 +102,9 @@
! big loop over all spectral elements in the fluid
! ****************************************************
- if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. ( .not. phase_is_inner )) div_displfluid(:,:,:,:) = 0._CUSTOM_REAL
+ if (MOVIE_VOLUME .and. NSPEC_OUTER_CORE_3DMOVIE /= 1 .and. ( .not. phase_is_inner )) then
+ div_displfluid(:,:,:,:) = 0._CUSTOM_REAL
+ endif
! computed_elements = 0
if( .not. phase_is_inner ) then
@@ -244,7 +246,8 @@
dpotentialdzl = dpotentialdzl + displfluid(iglob) * grad_z_ln_rho
- else ! if gravity is turned on
+ else
+ ! if gravity is turned on
! compute divergence of displacment
! precompute and store gravity term
@@ -290,7 +293,7 @@
! note: these calculations are only considered for SIMULATION_TYPE == 1 .and. SAVE_FORWARD
! and one has set MOVIE_VOLUME_TYPE == 4 when MOVIE_VOLUME is .true.;
! in case of SIMULATION_TYPE == 3, it gets overwritten by compute_kernels_outer_core()
- if (NSPEC_OUTER_CORE_ADJOINT /= 1 .and. MOVIE_VOLUME ) then
+ if( MOVIE_VOLUME .and. NSPEC_OUTER_CORE_3DMOVIE /= 1 ) then
div_displfluid(i,j,k,ispec) = &
minus_rho_g_over_kappa_fluid(int_radius) * (dpotentialdx_with_rot * gxl + &
dpotentialdy_with_rot * gyl + dpotentialdzl * gzl)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -212,9 +212,9 @@
!--- couple with outer core at the bottom of the mantle
!---
if(ACTUALLY_COUPLE_FLUID_CMB) &
- call compute_coupling_CMB_fluid(displ_crust_mantle,accel_crust_mantle, &
+ call compute_coupling_CMB_fluid(NGLOB_CRUST_MANTLE,displ_crust_mantle,accel_crust_mantle, &
ibool_crust_mantle,ibelm_bottom_crust_mantle, &
- accel_outer_core, &
+ NGLOB_OUTER_CORE,accel_outer_core, &
normal_top_outer_core,jacobian2D_top_outer_core, &
wgllwgll_xy,ibool_outer_core,ibelm_top_outer_core, &
RHO_TOP_OC,minus_g_cmb, &
@@ -224,9 +224,9 @@
!--- couple with outer core at the top of the inner core
!---
if(ACTUALLY_COUPLE_FLUID_ICB) &
- call compute_coupling_ICB_fluid(displ_inner_core,accel_inner_core, &
+ call compute_coupling_ICB_fluid(NGLOB_INNER_CORE,displ_inner_core,accel_inner_core, &
ibool_inner_core,ibelm_top_inner_core, &
- accel_outer_core, &
+ NGLOB_OUTER_CORE,accel_outer_core, &
normal_bottom_outer_core,jacobian2D_bottom_outer_core, &
wgllwgll_xy,ibool_outer_core,ibelm_bottom_outer_core, &
RHO_BOTTOM_OC,minus_g_icb, &
@@ -352,7 +352,7 @@
if( OCEANS_VAL ) then
if(.NOT. GPU_MODE) then
! on CPU
- call compute_coupling_ocean(accel_crust_mantle, &
+ call compute_coupling_ocean(NGLOB_CRUST_MANTLE,accel_crust_mantle, &
rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
rmass_ocean_load,normal_top_crust_mantle, &
ibool_crust_mantle,ibelm_top_crust_mantle, &
@@ -588,9 +588,9 @@
!--- couple with outer core at the bottom of the mantle
!---
if(ACTUALLY_COUPLE_FLUID_CMB) &
- call compute_coupling_CMB_fluid(b_displ_crust_mantle,b_accel_crust_mantle, &
+ call compute_coupling_CMB_fluid(NGLOB_CRUST_MANTLE_ADJOINT,b_displ_crust_mantle,b_accel_crust_mantle, &
ibool_crust_mantle,ibelm_bottom_crust_mantle, &
- b_accel_outer_core, &
+ NGLOB_OUTER_CORE_ADJOINT,b_accel_outer_core, &
normal_top_outer_core,jacobian2D_top_outer_core, &
wgllwgll_xy,ibool_outer_core,ibelm_top_outer_core, &
RHO_TOP_OC,minus_g_cmb, &
@@ -600,9 +600,9 @@
!--- couple with outer core at the top of the inner core
!---
if(ACTUALLY_COUPLE_FLUID_ICB) &
- call compute_coupling_ICB_fluid(b_displ_inner_core,b_accel_inner_core, &
+ call compute_coupling_ICB_fluid(NGLOB_INNER_CORE_ADJOINT,b_displ_inner_core,b_accel_inner_core, &
ibool_inner_core,ibelm_top_inner_core, &
- b_accel_outer_core, &
+ NGLOB_OUTER_CORE_ADJOINT,b_accel_outer_core, &
normal_bottom_outer_core,jacobian2D_bottom_outer_core, &
wgllwgll_xy,ibool_outer_core,ibelm_bottom_outer_core, &
RHO_BOTTOM_OC,minus_g_icb, &
@@ -732,7 +732,7 @@
if( OCEANS_VAL ) then
if(.NOT. GPU_MODE) then
! on CPU
- call compute_coupling_ocean(b_accel_crust_mantle, &
+ call compute_coupling_ocean(NGLOB_CRUST_MANTLE_ADJOINT,b_accel_crust_mantle, &
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, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -252,7 +252,7 @@
real(kind=CUSTOM_REAL), dimension(5) :: b_epsilondev_loc
real(kind=CUSTOM_REAL), dimension(5) :: epsilondev_loc
real(kind=CUSTOM_REAL), dimension(3) :: vector_accel
-
+ real(kind=CUSTOM_REAL) :: div_displ,b_div_displ
integer :: i,j,k,l,ispec,iglob
! outer_core -- compute the actual displacement and acceleration (NDIM,NGLOBMAX_OUTER_CORE)
@@ -323,11 +323,11 @@
! bulk modulus kernel
kappal = rhostore_outer_core(i,j,k,ispec)/kappavstore_outer_core(i,j,k,ispec)
- div_displ_outer_core(i,j,k,ispec) = kappal * accel_outer_core(iglob)
- b_div_displ_outer_core(i,j,k,ispec) = kappal * b_accel_outer_core(iglob)
+ div_displ = kappal * accel_outer_core(iglob)
+ b_div_displ = kappal * b_accel_outer_core(iglob)
alpha_kl_outer_core(i,j,k,ispec) = alpha_kl_outer_core(i,j,k,ispec) &
- + deltat * div_displ_outer_core(i,j,k,ispec) * b_div_displ_outer_core(i,j,k,ispec)
+ + deltat * div_displ * b_div_displ
! calculates gradient grad(displ) (also needed for boundary kernels)
if(SAVE_BOUNDARY_MESH .or. deviatoric_outercore) then
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -76,26 +76,26 @@
endif
! broadcast the information read on the master to the nodes
- call bcast_all_i(yr,1)
- call bcast_all_i(jda,1)
- call bcast_all_i(ho,1)
- call bcast_all_i(mi,1)
+ call bcast_all_singlei(yr)
+ call bcast_all_singlei(jda)
+ call bcast_all_singlei(ho)
+ call bcast_all_singlei(mi)
- call bcast_all_dp(sec,1)
+ call bcast_all_singledp(sec)
- call bcast_all_dp(tshift_cmt,1)
- call bcast_all_dp(t_shift,1)
+ call bcast_all_singledp(tshift_cmt)
+ call bcast_all_singledp(t_shift)
! event location given on first, PDE line
- call bcast_all_dp(elat,1)
- call bcast_all_dp(elon,1)
- call bcast_all_dp(depth,1)
+ call bcast_all_singledp(elat)
+ call bcast_all_singledp(elon)
+ call bcast_all_singledp(depth)
! cmt location given in CMT file
- call bcast_all_dp(cmt_lat,1)
- call bcast_all_dp(cmt_lon,1)
- call bcast_all_dp(cmt_depth,1)
- call bcast_all_dp(cmt_hdur,1)
+ call bcast_all_singledp(cmt_lat)
+ call bcast_all_singledp(cmt_lon)
+ call bcast_all_singledp(cmt_depth)
+ call bcast_all_singledp(cmt_hdur)
call bcast_all_ch(event_name,20)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/iterate_time_undoatt.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/iterate_time_undoatt.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/iterate_time_undoatt.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -343,12 +343,12 @@
xix_inner_core,xiy_inner_core,xiz_inner_core, &
etax_inner_core,etay_inner_core,etaz_inner_core, &
gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
- epsilondev_xx_inner_core(1,1,1,ispec), &
- epsilondev_yy_inner_core(1,1,1,ispec), &
- epsilondev_xy_inner_core(1,1,1,ispec), &
- epsilondev_xz_inner_core(1,1,1,ispec), &
- epsilondev_yz_inner_core(1,1,1,ispec), &
- eps_trace_over_3_inner_core(1,1,1,ispec))
+ epsilondev_xx_inner_core, &
+ epsilondev_yy_inner_core, &
+ epsilondev_xy_inner_core, &
+ epsilondev_xz_inner_core, &
+ epsilondev_yz_inner_core, &
+ NSPEC_INNER_CORE_STRAIN_ONLY,eps_trace_over_3_inner_core)
enddo
! crust mantle
do ispec = 1, NSPEC_crust_mantle
@@ -359,12 +359,12 @@
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, &
- epsilondev_xx_crust_mantle(1,1,1,ispec), &
- epsilondev_yy_crust_mantle(1,1,1,ispec), &
- epsilondev_xy_crust_mantle(1,1,1,ispec), &
- epsilondev_xz_crust_mantle(1,1,1,ispec), &
- epsilondev_yz_crust_mantle(1,1,1,ispec), &
- eps_trace_over_3_crust_mantle(1,1,1,ispec))
+ epsilondev_xx_crust_mantle, &
+ epsilondev_yy_crust_mantle, &
+ epsilondev_xy_crust_mantle, &
+ epsilondev_xz_crust_mantle, &
+ epsilondev_yz_crust_mantle, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,eps_trace_over_3_crust_mantle)
enddo
else
@@ -378,12 +378,12 @@
xix_inner_core,xiy_inner_core,xiz_inner_core, &
etax_inner_core,etay_inner_core,etaz_inner_core, &
gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
- epsilondev_xx_inner_core(1,1,1,ispec), &
- epsilondev_yy_inner_core(1,1,1,ispec), &
- epsilondev_xy_inner_core(1,1,1,ispec), &
- epsilondev_xz_inner_core(1,1,1,ispec), &
- epsilondev_yz_inner_core(1,1,1,ispec), &
- eps_trace_over_3_inner_core(1,1,1,ispec))
+ epsilondev_xx_inner_core, &
+ epsilondev_yy_inner_core, &
+ epsilondev_xy_inner_core, &
+ epsilondev_xz_inner_core, &
+ epsilondev_yz_inner_core, &
+ NSPEC_INNER_CORE_STRAIN_ONLY,eps_trace_over_3_inner_core)
enddo
! crust mantle
do ispec = 1, NSPEC_CRUST_MANTLE
@@ -394,12 +394,12 @@
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, &
- epsilondev_xx_crust_mantle(1,1,1,ispec), &
- epsilondev_yy_crust_mantle(1,1,1,ispec), &
- epsilondev_xy_crust_mantle(1,1,1,ispec), &
- epsilondev_xz_crust_mantle(1,1,1,ispec), &
- epsilondev_yz_crust_mantle(1,1,1,ispec), &
- eps_trace_over_3_crust_mantle(1,1,1,ispec))
+ epsilondev_xx_crust_mantle, &
+ epsilondev_yy_crust_mantle, &
+ epsilondev_xy_crust_mantle, &
+ epsilondev_xz_crust_mantle, &
+ epsilondev_yz_crust_mantle, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,eps_trace_over_3_crust_mantle)
enddo
endif
@@ -431,15 +431,16 @@
xix_inner_core,xiy_inner_core,xiz_inner_core, &
etax_inner_core,etay_inner_core,etaz_inner_core, &
gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
- b_epsilondev_xx_inner_core(1,1,1,ispec), &
- b_epsilondev_yy_inner_core(1,1,1,ispec), &
- b_epsilondev_xy_inner_core(1,1,1,ispec), &
- b_epsilondev_xz_inner_core(1,1,1,ispec), &
- b_epsilondev_yz_inner_core(1,1,1,ispec), &
- b_eps_trace_over_3_inner_core(1,1,1,ispec))
+ b_epsilondev_xx_inner_core, &
+ b_epsilondev_yy_inner_core, &
+ b_epsilondev_xy_inner_core, &
+ b_epsilondev_xz_inner_core, &
+ b_epsilondev_yz_inner_core, &
+ NSPEC_INNER_CORE_STRAIN_ONLY,b_eps_trace_over_3_inner_core)
enddo
+
! crust mantle
- do ispec = 1, NSPEC_crust_mantle
+ do ispec = 1, NSPEC_CRUST_MANTLE
call compute_element_strain_att_Dev(ispec,NGLOB_CRUST_MANTLE_ADJOINT,NSPEC_CRUST_MANTLE_ADJOINT, &
b_displ_crust_mantle,b_veloc_crust_mantle,0._CUSTOM_REAL, &
ibool_crust_mantle, &
@@ -447,12 +448,12 @@
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, &
- b_epsilondev_xx_crust_mantle(1,1,1,ispec), &
- b_epsilondev_yy_crust_mantle(1,1,1,ispec), &
- b_epsilondev_xy_crust_mantle(1,1,1,ispec), &
- b_epsilondev_xz_crust_mantle(1,1,1,ispec), &
- b_epsilondev_yz_crust_mantle(1,1,1,ispec), &
- b_eps_trace_over_3_crust_mantle(1,1,1,ispec))
+ b_epsilondev_xx_crust_mantle, &
+ b_epsilondev_yy_crust_mantle, &
+ b_epsilondev_xy_crust_mantle, &
+ b_epsilondev_xz_crust_mantle, &
+ b_epsilondev_yz_crust_mantle, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,b_eps_trace_over_3_crust_mantle)
enddo
else
@@ -466,12 +467,12 @@
xix_inner_core,xiy_inner_core,xiz_inner_core, &
etax_inner_core,etay_inner_core,etaz_inner_core, &
gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
- b_epsilondev_xx_inner_core(1,1,1,ispec), &
- b_epsilondev_yy_inner_core(1,1,1,ispec), &
- b_epsilondev_xy_inner_core(1,1,1,ispec), &
- b_epsilondev_xz_inner_core(1,1,1,ispec), &
- b_epsilondev_yz_inner_core(1,1,1,ispec), &
- b_eps_trace_over_3_inner_core(1,1,1,ispec))
+ b_epsilondev_xx_inner_core, &
+ b_epsilondev_yy_inner_core, &
+ b_epsilondev_xy_inner_core, &
+ b_epsilondev_xz_inner_core, &
+ b_epsilondev_yz_inner_core, &
+ NSPEC_INNER_CORE_STRAIN_ONLY,b_eps_trace_over_3_inner_core)
enddo
! crust mantle
do ispec = 1, NSPEC_crust_mantle
@@ -482,12 +483,12 @@
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, &
- b_epsilondev_xx_crust_mantle(1,1,1,ispec), &
- b_epsilondev_yy_crust_mantle(1,1,1,ispec), &
- b_epsilondev_xy_crust_mantle(1,1,1,ispec), &
- b_epsilondev_xz_crust_mantle(1,1,1,ispec), &
- b_epsilondev_yz_crust_mantle(1,1,1,ispec), &
- b_eps_trace_over_3_crust_mantle(1,1,1,ispec))
+ b_epsilondev_xx_crust_mantle, &
+ b_epsilondev_yy_crust_mantle, &
+ b_epsilondev_xy_crust_mantle, &
+ b_epsilondev_xz_crust_mantle, &
+ b_epsilondev_yz_crust_mantle, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,b_eps_trace_over_3_crust_mantle)
enddo
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -673,7 +673,6 @@
if(myrank == 0) then
! MPI loop on all the results to determine the best slice
- islice_selected_rec(:) = -1
do irec_in_this_subset = 1,nrec_SUBSET_current_size
! mapping from station number in current subset to real station number in all the subsets
@@ -854,7 +853,7 @@
deallocate(final_distance)
! main process broadcasts the results to all the slices
- call bcast_all_i(nrec,1)
+ call bcast_all_singlei(nrec)
call bcast_all_i(islice_selected_rec,nrec)
call bcast_all_i(ispec_selected_rec,nrec)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -156,7 +156,7 @@
call bcast_all_dp(depth,NSOURCES)
call bcast_all_dp(moment_tensor,6*NSOURCES)
- call bcast_all_dp(min_tshift_cmt_original,1)
+ call bcast_all_singledp(min_tshift_cmt_original)
! define topology of the control element
call hex_nodes(iaddx,iaddy,iaddr)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -1122,7 +1122,6 @@
alpha_kl_inner_core(:,:,:,:) = 0._CUSTOM_REAL
div_displ_outer_core(:,:,:,:) = 0._CUSTOM_REAL
- b_div_displ_outer_core(:,:,:,:) = 0._CUSTOM_REAL
! deviatoric kernel check
if( deviatoric_outercore) then
@@ -2426,7 +2425,7 @@
if( ier /= 0 ) stop 'error allocating arrays'
free_points_all(:) = 0
- call gather_all_i(free_np,1,free_points_all,1,NPROC)
+ call gather_all_singlei(free_np,free_points_all,NPROC)
! array offsets
allocate(free_offset_all(NPROC),stat=ier)
@@ -2446,7 +2445,7 @@
if( ier /= 0 ) stop 'error allocating arrays'
free_conn_nspec_all(:) = 0
- call gather_all_i(4*free_nspec,1,free_conn_nspec_all,1,NPROC)
+ call gather_all_singlei(4*free_nspec,free_conn_nspec_all,NPROC)
! array offsets
allocate(free_conn_offset_all(NPROC),stat=ier)
@@ -2694,7 +2693,7 @@
if( ier /= 0 ) stop 'error allocating arrays'
vtkdata_points_all(:) = 0
- call gather_all_i(vol_np,1,vtkdata_points_all,1,NPROC)
+ call gather_all_singlei(vol_np,vtkdata_points_all,NPROC)
! array offsets
allocate(vtkdata_offset_all(NPROC),stat=ier)
@@ -2714,7 +2713,7 @@
if( ier /= 0 ) stop 'error allocating arrays'
vol_conn_nspec_all(:) = 0
- call gather_all_i(8*vol_nspec,1,vol_conn_nspec_all,1,NPROC)
+ call gather_all_singlei(8*vol_nspec,vol_conn_nspec_all,NPROC)
! array offsets
allocate(vol_conn_offset_all(NPROC),stat=ier)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -37,7 +37,8 @@
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
ibool,idoubling,ispec_is_tiso, &
- rmassx,rmassy,rmassz,rmass_ocean_load, &
+ rmassx,rmassy,rmassz, &
+ nglob_oceans,rmass_ocean_load, &
READ_KAPPA_MU,READ_TISO, &
b_rmassx,b_rmassy)
@@ -85,8 +86,10 @@
real(kind=CUSTOM_REAL), dimension(nglob_xy) :: b_rmassx,b_rmassy
real(kind=CUSTOM_REAL), dimension(nglob) :: rmassz
- real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
+ integer :: nglob_oceans
+ real(kind=CUSTOM_REAL), dimension(nglob_oceans) :: rmass_ocean_load
+
! flags to know if we should read Vs and anisotropy arrays
logical :: READ_KAPPA_MU,READ_TISO
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver_adios.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver_adios.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_arrays_solver_adios.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -37,7 +37,8 @@
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
ibool,idoubling,ispec_is_tiso, &
- rmassx,rmassy,rmassz,rmass_ocean_load, &
+ rmassx,rmassy,rmassz, &
+ nglob_oceans,rmass_ocean_load, &
READ_KAPPA_MU,READ_TISO, &
b_rmassx,b_rmassy)
@@ -88,8 +89,10 @@
real(kind=CUSTOM_REAL), dimension(nglob_xy) :: b_rmassx,b_rmassy
real(kind=CUSTOM_REAL), dimension(nglob) :: rmassz
- real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
+ integer :: nglob_oceans
+ real(kind=CUSTOM_REAL), dimension(nglob_oceans) :: rmass_ocean_load
+
! flags to know if we should read Vs and anisotropy arrays
logical :: READ_KAPPA_MU,READ_TISO
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -238,7 +238,8 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,dummy_idoubling,ispec_is_tiso_crust_mantle, &
- rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle,rmass_ocean_load, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ NGLOB_CRUST_MANTLE_OCEANS,rmass_ocean_load, &
READ_KAPPA_MU,READ_TISO, &
b_rmassx_crust_mantle,b_rmassy_crust_mantle)
else
@@ -260,7 +261,8 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,dummy_idoubling,ispec_is_tiso_crust_mantle, &
- rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle,rmass_ocean_load, &
+ rmassx_crust_mantle,rmassy_crust_mantle,rmassz_crust_mantle, &
+ NGLOB_CRUST_MANTLE_OCEANS,rmass_ocean_load, &
READ_KAPPA_MU,READ_TISO, &
b_rmassx_crust_mantle,b_rmassy_crust_mantle)
endif
@@ -372,7 +374,8 @@
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,dummy_array, &
+ dummy_rmass,dummy_rmass,rmass_outer_core, &
+ 1,dummy_array, &
READ_KAPPA_MU,READ_TISO, &
dummy_rmass,dummy_rmass)
else
@@ -394,7 +397,8 @@
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,dummy_array, &
+ dummy_rmass,dummy_rmass,rmass_outer_core, &
+ 1, dummy_array, &
READ_KAPPA_MU,READ_TISO, &
dummy_rmass,dummy_rmass)
endif
@@ -494,7 +498,8 @@
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,dummy_array, &
+ rmassx_inner_core,rmassy_inner_core,rmassz_inner_core, &
+ 1,dummy_array, &
READ_KAPPA_MU,READ_TISO, &
b_rmassx_inner_core,b_rmassy_inner_core)
else
@@ -516,7 +521,8 @@
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,dummy_array, &
+ rmassx_inner_core,rmassy_inner_core,rmassz_inner_core, &
+ 1,dummy_array, &
READ_KAPPA_MU,READ_TISO, &
b_rmassx_inner_core,b_rmassy_inner_core)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/rules.mk 2013-10-22 14:57:01 UTC (rev 22968)
@@ -214,7 +214,7 @@
$(EMPTY_MACRO)
adios_specfem3D_SHARED_STUBS = \
- $O/adios_method_stubs.shared.o \
+ $O/adios_method_stubs.cc.o \
$(EMPTY_MACRO)
# conditional adios linking
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -365,6 +365,8 @@
! count number of receivers located in this slice
nrec_local = 0
if (SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then
+ ! note: for 1-chunk simulations, nrec is now the actual number of receivers found in this chunk
+ ! (excludes stations located outside of chunk)
nrec_simulation = nrec
do irec = 1,nrec
if(myrank == islice_selected_rec(irec)) nrec_local = nrec_local + 1
@@ -445,16 +447,22 @@
if(myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) 'found a total of ',nrec_tot_found,' receivers in all slices'
- if(nrec_tot_found <= 1) then
- call exit_MPI(myrank,'error: no receiver found in the mesh')
-!! DK DK we do not check if the total found is equal to nrec_simulation here because some stations in the STATIONS file
-!! DK DK can be located outside of the mesh, thus the only thing we can check is that the total found is smaller or equal
- else if(nrec_tot_found > nrec_simulation) then
- call exit_MPI(myrank,'error: too many receivers found, problem when dispatching the receivers')
+ ! checks for number of receivers
+ ! note: for 1-chunk simulations, nrec_simulations is the number of receivers/sources found in this chunk
+ if(nrec_tot_found /= nrec_simulation) then
+ call exit_MPI(myrank,'problem when dispatching the receivers')
+ else
+ write(IMAIN,*) 'this total is okay'
endif
call flush_IMAIN()
endif
+ ! check that the sum of the number of receivers in each slice is nrec
+ if( SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3 ) then
+ if(myrank == 0 .and. nrec_tot_found /= nrec) &
+ call exit_MPI(myrank,'total number of receivers is incorrect')
+ endif
+
end subroutine setup_receivers
!
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D_par.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D_par.F90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D_par.F90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -680,12 +680,6 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ADJOINT) :: &
rho_kl_outer_core,alpha_kl_outer_core
- ! kernel runs
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ADJOINT) :: &
- div_displ_outer_core
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ADJOINT) :: &
- b_div_displ_outer_core
-
! check for deviatoric kernel for outer core region
real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: beta_kl_outer_core
integer :: nspec_beta_kl_outer_core
@@ -747,6 +741,10 @@
logical, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: mask_3dmovie
logical, dimension(NGLOB_CRUST_MANTLE_3DMOVIE) :: mask_ibool
+ ! outer core
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_3DMOVIE) :: &
+ div_displ_outer_core
+
! vtk run-time visualization
#ifdef HAVE_VTK
! vtk window
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_output.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_output.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_output.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -39,10 +39,12 @@
character(len=256) :: filename
integer,dimension(:),allocatable :: dummy_i
+ !-----------------------------------------------------------------------------
+ ! user parameters
+ ! outputs volume snapshot vtk-files of displacement in crust/mantle region for debugging
logical, parameter :: DEBUG_SNAPSHOT = .false.
+ !-----------------------------------------------------------------------------
- logical, parameter :: RUN_EXTERNAL_SCRIPT = .true.
- character(len=256) :: script_name = "tar_databases_file.sh"
character(len=256) :: system_command
! save movie on surface
@@ -62,6 +64,18 @@
! save velocity here to avoid static offset on displacement for movies
call write_movie_surface()
+ ! executes an external script on the node
+ if( RUN_EXTERNAL_MOVIE_SCRIPT ) then
+ ! synchronizes outputs
+ call synchronize_all()
+ ! calls shell external command
+ if( myrank == 0 ) then
+ write(system_command,"(a,1x,i6.6,' >& out.',i6.6,'.log &')") trim(MOVIE_SCRIPT_NAME),it,it
+ !print*,trim(system_command)
+ call system(system_command)
+ endif
+ endif
+
endif
endif
@@ -232,14 +246,17 @@
end select ! MOVIE_VOLUME_TYPE
! executes an external script on the node
- if( RUN_EXTERNAL_SCRIPT ) then
+ if( RUN_EXTERNAL_MOVIE_SCRIPT ) then
+ ! synchronizes outputs
call synchronize_all()
+ ! calls shell external command
if( myrank == 0 ) then
- write(system_command,"('./',a,1x,i6.6,' >& out.',i6.6,'.log &')") trim(script_name),it,it
+ write(system_command,"(a,1x,i6.6,' >& out.',i6.6,'.log &')") trim(MOVIE_SCRIPT_NAME),it,it
!print*,trim(system_command)
call system(system_command)
endif
endif
+
endif
endif ! MOVIE_VOLUME
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90 2013-10-22 09:55:50 UTC (rev 22967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90 2013-10-22 14:57:01 UTC (rev 22968)
@@ -317,7 +317,7 @@
! input
integer :: myrank,npoints_3dmovie,MOVIE_VOLUME_TYPE,it
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: eps_trace_over_3_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: eps_trace_over_3_crust_mantle
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: &
epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
@@ -573,14 +573,14 @@
!
subroutine write_movie_volume_divcurl(myrank,it,eps_trace_over_3_crust_mantle,&
- div_displ_outer_core, &
- accel_outer_core,kappavstore_outer_core,rhostore_outer_core,ibool_outer_core, &
- eps_trace_over_3_inner_core, &
- epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
- epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
- epsilondev_xx_inner_core,epsilondev_yy_inner_core,epsilondev_xy_inner_core, &
- epsilondev_xz_inner_core,epsilondev_yz_inner_core, &
- LOCAL_TMP_PATH)
+ div_displ_outer_core, &
+ accel_outer_core,kappavstore_outer_core,rhostore_outer_core,ibool_outer_core, &
+ eps_trace_over_3_inner_core, &
+ epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle,epsilondev_xy_crust_mantle, &
+ epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle, &
+ epsilondev_xx_inner_core,epsilondev_yy_inner_core,epsilondev_xy_inner_core, &
+ epsilondev_xz_inner_core,epsilondev_yz_inner_core, &
+ LOCAL_TMP_PATH)
! outputs divergence and curl: MOVIE_VOLUME_TYPE == 4
@@ -590,7 +590,7 @@
integer :: myrank,it
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: eps_trace_over_3_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ADJOINT) :: div_displ_outer_core
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_3DMOVIE) :: div_displ_outer_core
real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: accel_outer_core
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE) :: &
@@ -632,7 +632,7 @@
close(27)
! outer core
- if (NSPEC_OUTER_CORE_ADJOINT /= 1 ) then
+ if (NSPEC_OUTER_CORE_3DMOVIE /= 1 ) then
write(outputname,"('proc',i6.6,'_reg2_div_displ_it',i6.6,'.bin')") myrank,it
open(unit=27,file=trim(LOCAL_TMP_PATH)//'/'//trim(outputname),status='unknown',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error opening file '//trim(outputname))
More information about the CIG-COMMITS
mailing list