[cig-commits] r20011 - in seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER: . src/create_header_file src/cuda src/meshfem3D src/shared src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Tue May 1 07:58:17 PDT 2012
Author: danielpeter
Date: 2012-05-01 07:58:16 -0700 (Tue, 01 May 2012)
New Revision: 20011
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/create_header_file/create_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_crust_mantle_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_inner_core_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_MPI_interfaces.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_central_cube_buffers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_chunk_buffers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_meshes.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/finalize_mesher.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_interfaces.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_models.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_aniso_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_color_perm.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_model.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/test_MPI_interfaces.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/define_all_layers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/get_timestep_and_layers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/memory_eval.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/model_topo_bathy.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/write_VTK_file.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_boundary_kernel.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_elastic.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_kernels.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_topography_bathymetry.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_output.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90
Log:
updates memory evaluation and source location routine; updates tabs and whitespaces
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/Makefile.in 2012-05-01 14:58:16 UTC (rev 20011)
@@ -122,7 +122,12 @@
clean: required
- (rm -rf bin lib obj src/meshfem3D/*.mod src/specfem3D/*.mod)
+ (rm -rf bin lib obj )
+ (cd src/meshfem3D; make clean)
+ (cd src/specfem3D; make clean)
+ (cd src/create_header_file; make clean)
+ (cd src/auxiliaries; make clean)
+ (rm -rf src/meshfem3D/*.mod src/specfem3D/*.mod src/auxiliaries/*.mod)
help:
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/create_header_file/create_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/create_header_file/create_header_file.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/create_header_file/create_header_file.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -155,10 +155,11 @@
endif
! evaluate the amount of static memory needed by the solver
- call memory_eval(OCEANS,ABSORBING_CONDITIONS,ATTENUATION,ANISOTROPIC_3D_MANTLE,&
- TRANSVERSE_ISOTROPY,ANISOTROPIC_INNER_CORE,ROTATION,&
- ONE_CRUST,doubling_index,this_region_has_a_doubling,&
- ner,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_sampling_array,&
+ call memory_eval(OCEANS,ABSORBING_CONDITIONS,ATTENUATION,ANISOTROPIC_3D_MANTLE, &
+ TRANSVERSE_ISOTROPY,ANISOTROPIC_INNER_CORE,ROTATION,TOPOGRAPHY, &
+ ONE_CRUST,doubling_index,this_region_has_a_doubling, &
+ ner,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ ratio_sampling_array,NPROCTOT, &
NSPEC,nglob,SIMULATION_TYPE,MOVIE_VOLUME,SAVE_FORWARD, &
NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
@@ -171,7 +172,9 @@
NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
- NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION,static_memory_size)
+ NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ static_memory_size)
NGLOB1D_RADIAL_TEMP(:) = &
(/maxval(NGLOB1D_RADIAL_CORNER(1,:)),maxval(NGLOB1D_RADIAL_CORNER(2,:)),maxval(NGLOB1D_RADIAL_CORNER(3,:))/)
@@ -179,7 +182,8 @@
! create include file for the solver
call save_header_file(NSPEC,nglob,NEX_XI,NEX_ETA,NPROC,NPROCTOT, &
TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- ELLIPTICITY,GRAVITY,ROTATION,OCEANS,ATTENUATION,ATTENUATION_3D, &
+ ELLIPTICITY,GRAVITY,ROTATION, &
+ TOPOGRAPHY,OCEANS,ATTENUATION,ATTENUATION_3D, &
ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,NCHUNKS, &
INCLUDE_CENTRAL_CUBE,CENTER_LONGITUDE_IN_DEGREES, &
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,NSOURCES,NSTEP,&
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu 2012-05-01 14:58:16 UTC (rev 20011)
@@ -606,8 +606,8 @@
// make updated component of right-hand side
// we divide by rmass() which is 1 / M
// we use the total force which includes the Coriolis term above
- force_normal_comp = ( accel_crust_mantle[iglob*3]*nx
- + accel_crust_mantle[iglob*3+1]*ny
+ force_normal_comp = ( accel_crust_mantle[iglob*3]*nx
+ + accel_crust_mantle[iglob*3+1]*ny
+ accel_crust_mantle[iglob*3+2]*nz ) / rmass_crust_mantle[iglob];
additional_term = (rmass_ocean_load[iglob] - rmass_crust_mantle[iglob]) * force_normal_comp;
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_crust_mantle_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_crust_mantle_cuda.cu 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_crust_mantle_cuda.cu 2012-05-01 14:58:16 UTC (rev 20011)
@@ -896,7 +896,7 @@
int* d_phase_ispec_inner,
int num_phase_ispec,
int d_iphase,
- realw d_deltat,
+ realw d_deltat,
int use_mesh_coloring_gpu,
realw* d_displ,
realw* d_veloc,
@@ -1031,16 +1031,16 @@
#endif
if( ATTENUATION){
- // use first order Taylor expansion of displacement for local storage of stresses
- // at this current time step, to fix attenuation in a consistent way
+ // use first order Taylor expansion of displacement for local storage of stresses
+ // at this current time step, to fix attenuation in a consistent way
#ifdef USE_TEXTURES
- s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob);
- s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + NGLOB);
- s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + 2*NGLOB);
+ s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob);
+ s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + NGLOB);
+ s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + 2*NGLOB);
#else
- s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * d_veloc[iglob*3];
- s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * d_veloc[iglob*3 + 1];
- s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * d_veloc[iglob*3 + 2];
+ s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * d_veloc[iglob*3];
+ s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * d_veloc[iglob*3 + 1];
+ s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * d_veloc[iglob*3 + 2];
#endif
}
@@ -1092,21 +1092,21 @@
if( ATTENUATION){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
- tempx1l_att = 0.f;
- tempx2l_att = 0.f;
- tempx3l_att = 0.f;
+ tempx1l_att = 0.f;
+ tempx2l_att = 0.f;
+ tempx3l_att = 0.f;
- tempy1l_att = 0.f;
- tempy2l_att = 0.f;
- tempy3l_att = 0.f;
+ tempy1l_att = 0.f;
+ tempy2l_att = 0.f;
+ tempy3l_att = 0.f;
- tempz1l_att = 0.f;
- tempz2l_att = 0.f;
- tempz3l_att = 0.f;
+ tempz1l_att = 0.f;
+ tempz2l_att = 0.f;
+ tempz3l_att = 0.f;
- for (l=0;l<NGLLX;l++) {
+ for (l=0;l<NGLLX;l++) {
hp1 = d_hprime_xx[l*NGLLX+I];
offset = K*NGLL2+J*NGLLX+l;
tempx1l_att += s_dummyx_loc_att[offset]*hp1;
@@ -1125,7 +1125,7 @@
tempy3l_att += s_dummyy_loc_att[offset]*hp3;
tempz3l_att += s_dummyz_loc_att[offset]*hp3;
- }
+ }
}
#else
@@ -1184,61 +1184,61 @@
+ s_dummyz_loc[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
if( ATTENUATION){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
- tempx1l_att = s_dummyx_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempx1l_att = s_dummyx_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempy1l_att = s_dummyy_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempy1l_att = s_dummyy_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempz1l_att = s_dummyz_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempz1l_att = s_dummyz_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempx2l_att = s_dummyx_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyx_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempx2l_att = s_dummyx_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyx_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempy2l_att = s_dummyy_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyy_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempy2l_att = s_dummyy_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyy_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempz2l_att = s_dummyz_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyz_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempz2l_att = s_dummyz_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyz_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempx3l_att = s_dummyx_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyx_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyx_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyx_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyx_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
+ tempx3l_att = s_dummyx_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyx_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyx_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyx_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyx_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
- tempy3l_att = s_dummyy_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyy_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyy_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyy_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyy_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
+ tempy3l_att = s_dummyy_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyy_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyy_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyy_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyy_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
- tempz3l_att = s_dummyz_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyz_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyz_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyz_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyz_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
+ tempz3l_att = s_dummyz_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyz_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyz_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyz_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyz_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
}
#endif
@@ -1277,60 +1277,60 @@
duzdyl_plus_duydzl = duzdyl + duydzl;
if( ATTENUATION){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
- duxdxl_att = xixl*tempx1l_att + etaxl*tempx2l_att + gammaxl*tempx3l_att;
- duxdyl_att = xiyl*tempx1l_att + etayl*tempx2l_att + gammayl*tempx3l_att;
- duxdzl_att = xizl*tempx1l_att + etazl*tempx2l_att + gammazl*tempx3l_att;
+ duxdxl_att = xixl*tempx1l_att + etaxl*tempx2l_att + gammaxl*tempx3l_att;
+ duxdyl_att = xiyl*tempx1l_att + etayl*tempx2l_att + gammayl*tempx3l_att;
+ duxdzl_att = xizl*tempx1l_att + etazl*tempx2l_att + gammazl*tempx3l_att;
- duydxl_att = xixl*tempy1l_att + etaxl*tempy2l_att + gammaxl*tempy3l_att;
- duydyl_att = xiyl*tempy1l_att + etayl*tempy2l_att + gammayl*tempy3l_att;
- duydzl_att = xizl*tempy1l_att + etazl*tempy2l_att + gammazl*tempy3l_att;
+ duydxl_att = xixl*tempy1l_att + etaxl*tempy2l_att + gammaxl*tempy3l_att;
+ duydyl_att = xiyl*tempy1l_att + etayl*tempy2l_att + gammayl*tempy3l_att;
+ duydzl_att = xizl*tempy1l_att + etazl*tempy2l_att + gammazl*tempy3l_att;
- duzdxl_att = xixl*tempz1l_att + etaxl*tempz2l_att + gammaxl*tempz3l_att;
- duzdyl_att = xiyl*tempz1l_att + etayl*tempz2l_att + gammayl*tempz3l_att;
- duzdzl_att = xizl*tempz1l_att + etazl*tempz2l_att + gammazl*tempz3l_att;
+ duzdxl_att = xixl*tempz1l_att + etaxl*tempz2l_att + gammaxl*tempz3l_att;
+ duzdyl_att = xiyl*tempz1l_att + etayl*tempz2l_att + gammayl*tempz3l_att;
+ duzdzl_att = xizl*tempz1l_att + etazl*tempz2l_att + gammazl*tempz3l_att;
- // precompute some sums to save CPU time
- duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att;
- duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att;
- duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att;
+ // precompute some sums to save CPU time
+ duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att;
+ duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att;
+ duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att;
- // computes deviatoric strain attenuation and/or for kernel calculations
- if(COMPUTE_AND_STORE_STRAIN) {
- realw templ = 0.33333333333333333333f * (duxdxl_att + duydyl_att + duzdzl_att); // 1./3. = 0.33333
+ // computes deviatoric strain attenuation and/or for kernel calculations
+ if(COMPUTE_AND_STORE_STRAIN) {
+ realw templ = 0.33333333333333333333f * (duxdxl_att + duydyl_att + duzdzl_att); // 1./3. = 0.33333
- // local storage: stresses at this current time step
- epsilondev_xx_loc = duxdxl_att - templ;
- epsilondev_yy_loc = duydyl_att - templ;
- epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl_att;
- epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl_att;
- epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl_att;
+ // local storage: stresses at this current time step
+ epsilondev_xx_loc = duxdxl_att - templ;
+ epsilondev_yy_loc = duydyl_att - templ;
+ epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl_att;
+ epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl_att;
+ epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl_att;
- if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) {
- epsilon_trace_over_3[tx] = templ;
- }else{
- epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
- }
- }
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) {
+ epsilon_trace_over_3[tx] = templ;
+ }else{
+ epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
+ }
+ }
}else{
- // computes deviatoric strain attenuation and/or for kernel calculations
- if(COMPUTE_AND_STORE_STRAIN) {
- realw templ = 0.33333333333333333333f * (duxdxl + duydyl + duzdzl); // 1./3. = 0.33333
+ // computes deviatoric strain attenuation and/or for kernel calculations
+ if(COMPUTE_AND_STORE_STRAIN) {
+ realw templ = 0.33333333333333333333f * (duxdxl + duydyl + duzdzl); // 1./3. = 0.33333
- // local storage: stresses at this current time step
- epsilondev_xx_loc = duxdxl - templ;
- epsilondev_yy_loc = duydyl - templ;
- epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl;
- epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl;
- epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl;
+ // local storage: stresses at this current time step
+ epsilondev_xx_loc = duxdxl - templ;
+ epsilondev_yy_loc = duydyl - templ;
+ epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl;
+ epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl;
+ epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl;
- if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) {
- epsilon_trace_over_3[tx] = templ;
- }else{
- epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
- }
- }
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) {
+ epsilon_trace_over_3[tx] = templ;
+ }else{
+ epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
+ }
+ }
}
// attenuation
@@ -1619,7 +1619,7 @@
/* ----------------------------------------------------------------------------------------------- */
void Kernel_2_crust_mantle(int nb_blocks_to_compute,Mesh* mp,
- realw d_deltat,
+ realw d_deltat,
int d_iphase,
int* d_ibool,
int* d_ispec_is_tiso,
@@ -1695,7 +1695,7 @@
mp->d_phase_ispec_inner_crust_mantle,
mp->num_phase_ispec_crust_mantle,
d_iphase,
- d_deltat,
+ d_deltat,
mp->use_mesh_coloring_gpu,
mp->d_displ_crust_mantle,
mp->d_veloc_crust_mantle,
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_inner_core_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_inner_core_cuda.cu 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_forces_inner_core_cuda.cu 2012-05-01 14:58:16 UTC (rev 20011)
@@ -298,10 +298,10 @@
int* d_phase_ispec_inner,
int num_phase_ispec,
int d_iphase,
- realw d_deltat,
+ realw d_deltat,
int use_mesh_coloring_gpu,
realw* d_displ,
- realw* d_veloc,
+ realw* d_veloc,
realw* d_accel,
realw* d_xix, realw* d_xiy, realw* d_xiz,
realw* d_etax, realw* d_etay, realw* d_etaz,
@@ -433,19 +433,19 @@
s_dummyz_loc[tx] = d_displ[iglob*3 + 2];
#endif
- if( ATTENUATION ){
- // use first order Taylor expansion of displacement for local storage of stresses
- // at this current time step, to fix attenuation in a consistent way
+ if( ATTENUATION ){
+ // use first order Taylor expansion of displacement for local storage of stresses
+ // at this current time step, to fix attenuation in a consistent way
#ifdef USE_TEXTURES
- s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob);
- s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + NGLOB);
- s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + 2*NGLOB);
+ s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob);
+ s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + NGLOB);
+ s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * tex1Dfetch(tex_veloc, iglob + 2*NGLOB);
#else
- s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * d_veloc[iglob*3];
- s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * d_veloc[iglob*3 + 1];
- s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * d_veloc[iglob*3 + 2];
+ s_dummyx_loc_att[tx] = s_dummyx_loc[tx] + d_deltat * d_veloc[iglob*3];
+ s_dummyy_loc_att[tx] = s_dummyy_loc[tx] + d_deltat * d_veloc[iglob*3 + 1];
+ s_dummyz_loc_att[tx] = s_dummyz_loc[tx] + d_deltat * d_veloc[iglob*3 + 2];
#endif
- }
+ }
}
}
@@ -495,21 +495,21 @@
}
if( ATTENUATION ){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
- tempx1l_att = 0.f;
- tempx2l_att = 0.f;
- tempx3l_att = 0.f;
+ tempx1l_att = 0.f;
+ tempx2l_att = 0.f;
+ tempx3l_att = 0.f;
- tempy1l_att = 0.f;
- tempy2l_att = 0.f;
- tempy3l_att = 0.f;
+ tempy1l_att = 0.f;
+ tempy2l_att = 0.f;
+ tempy3l_att = 0.f;
- tempz1l_att = 0.f;
- tempz2l_att = 0.f;
- tempz3l_att = 0.f;
+ tempz1l_att = 0.f;
+ tempz2l_att = 0.f;
+ tempz3l_att = 0.f;
- for (l=0;l<NGLLX;l++) {
+ for (l=0;l<NGLLX;l++) {
hp1 = d_hprime_xx[l*NGLLX+I];
offset = K*NGLL2+J*NGLLX+l;
tempx1l_att += s_dummyx_loc_att[offset]*hp1;
@@ -528,7 +528,7 @@
tempy3l_att += s_dummyy_loc_att[offset]*hp3;
tempz3l_att += s_dummyz_loc_att[offset]*hp3;
- }
+ }
}
#else
@@ -587,61 +587,61 @@
+ s_dummyz_loc[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
if( ATTENUATION ){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
- tempx1l_att = s_dummyx_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyx_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempx1l_att = s_dummyx_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyx_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempy1l_att = s_dummyy_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyy_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempy1l_att = s_dummyy_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyy_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempz1l_att = s_dummyz_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
- + s_dummyz_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
+ tempz1l_att = s_dummyz_loc_att[K*NGLL2+J*NGLLX]*d_hprime_xx[I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+1]*d_hprime_xx[NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+2]*d_hprime_xx[2*NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+3]*d_hprime_xx[3*NGLLX+I]
+ + s_dummyz_loc_att[K*NGLL2+J*NGLLX+4]*d_hprime_xx[4*NGLLX+I];
- tempx2l_att = s_dummyx_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyx_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyx_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempx2l_att = s_dummyx_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyx_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyx_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempy2l_att = s_dummyy_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyy_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyy_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempy2l_att = s_dummyy_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyy_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyy_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempz2l_att = s_dummyz_loc_att[K*NGLL2+I]*d_hprime_xx[J]
- + s_dummyz_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
- + s_dummyz_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
+ tempz2l_att = s_dummyz_loc_att[K*NGLL2+I]*d_hprime_xx[J]
+ + s_dummyz_loc_att[K*NGLL2+NGLLX+I]*d_hprime_xx[NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+2*NGLLX+I]*d_hprime_xx[2*NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+3*NGLLX+I]*d_hprime_xx[3*NGLLX+J]
+ + s_dummyz_loc_att[K*NGLL2+4*NGLLX+I]*d_hprime_xx[4*NGLLX+J];
- tempx3l_att = s_dummyx_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyx_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyx_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyx_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyx_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
+ tempx3l_att = s_dummyx_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyx_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyx_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyx_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyx_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
- tempy3l_att = s_dummyy_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyy_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyy_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyy_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyy_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
+ tempy3l_att = s_dummyy_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyy_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyy_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyy_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyy_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
- tempz3l_att = s_dummyz_loc_att[J*NGLLX+I]*d_hprime_xx[K]
- + s_dummyz_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
- + s_dummyz_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
- + s_dummyz_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
- + s_dummyz_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
+ tempz3l_att = s_dummyz_loc_att[J*NGLLX+I]*d_hprime_xx[K]
+ + s_dummyz_loc_att[NGLL2+J*NGLLX+I]*d_hprime_xx[NGLLX+K]
+ + s_dummyz_loc_att[2*NGLL2+J*NGLLX+I]*d_hprime_xx[2*NGLLX+K]
+ + s_dummyz_loc_att[3*NGLL2+J*NGLLX+I]*d_hprime_xx[3*NGLLX+K]
+ + s_dummyz_loc_att[4*NGLL2+J*NGLLX+I]*d_hprime_xx[4*NGLLX+K];
}
#endif
@@ -680,56 +680,56 @@
duzdyl_plus_duydzl = duzdyl + duydzl;
if( ATTENUATION ){
- // temporary variables used for fixing attenuation in a consistent way
+ // temporary variables used for fixing attenuation in a consistent way
- duxdxl_att = xixl*tempx1l_att + etaxl*tempx2l_att + gammaxl*tempx3l_att;
- duxdyl_att = xiyl*tempx1l_att + etayl*tempx2l_att + gammayl*tempx3l_att;
- duxdzl_att = xizl*tempx1l_att + etazl*tempx2l_att + gammazl*tempx3l_att;
+ duxdxl_att = xixl*tempx1l_att + etaxl*tempx2l_att + gammaxl*tempx3l_att;
+ duxdyl_att = xiyl*tempx1l_att + etayl*tempx2l_att + gammayl*tempx3l_att;
+ duxdzl_att = xizl*tempx1l_att + etazl*tempx2l_att + gammazl*tempx3l_att;
- duydxl_att = xixl*tempy1l_att + etaxl*tempy2l_att + gammaxl*tempy3l_att;
- duydyl_att = xiyl*tempy1l_att + etayl*tempy2l_att + gammayl*tempy3l_att;
- duydzl_att = xizl*tempy1l_att + etazl*tempy2l_att + gammazl*tempy3l_att;
+ duydxl_att = xixl*tempy1l_att + etaxl*tempy2l_att + gammaxl*tempy3l_att;
+ duydyl_att = xiyl*tempy1l_att + etayl*tempy2l_att + gammayl*tempy3l_att;
+ duydzl_att = xizl*tempy1l_att + etazl*tempy2l_att + gammazl*tempy3l_att;
- duzdxl_att = xixl*tempz1l_att + etaxl*tempz2l_att + gammaxl*tempz3l_att;
- duzdyl_att = xiyl*tempz1l_att + etayl*tempz2l_att + gammayl*tempz3l_att;
- duzdzl_att = xizl*tempz1l_att + etazl*tempz2l_att + gammazl*tempz3l_att;
+ duzdxl_att = xixl*tempz1l_att + etaxl*tempz2l_att + gammaxl*tempz3l_att;
+ duzdyl_att = xiyl*tempz1l_att + etayl*tempz2l_att + gammayl*tempz3l_att;
+ duzdzl_att = xizl*tempz1l_att + etazl*tempz2l_att + gammazl*tempz3l_att;
- // precompute some sums to save CPU time
- duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att;
- duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att;
- duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att;
+ // precompute some sums to save CPU time
+ duxdyl_plus_duydxl_att = duxdyl_att + duydxl_att;
+ duzdxl_plus_duxdzl_att = duzdxl_att + duxdzl_att;
+ duzdyl_plus_duydzl_att = duzdyl_att + duydzl_att;
- // computes deviatoric strain attenuation and/or for kernel calculations
- if(COMPUTE_AND_STORE_STRAIN) {
- realw templ = 0.33333333333333333333f * (duxdxl_att + duydyl_att + duzdzl_att); // 1./3. = 0.33333
+ // computes deviatoric strain attenuation and/or for kernel calculations
+ if(COMPUTE_AND_STORE_STRAIN) {
+ realw templ = 0.33333333333333333333f * (duxdxl_att + duydyl_att + duzdzl_att); // 1./3. = 0.33333
- // local storage: stresses at this current time step
- epsilondev_xx_loc = duxdxl_att - templ;
- epsilondev_yy_loc = duydyl_att - templ;
- epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl_att;
- epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl_att;
- epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl_att;
+ // local storage: stresses at this current time step
+ epsilondev_xx_loc = duxdxl_att - templ;
+ epsilondev_yy_loc = duydyl_att - templ;
+ epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl_att;
+ epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl_att;
+ epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl_att;
- if(SIMULATION_TYPE == 3) {
- epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
- }
- }
+ if(SIMULATION_TYPE == 3) {
+ epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
+ }
+ }
}else{
- // computes deviatoric strain attenuation and/or for kernel calculations
- if(COMPUTE_AND_STORE_STRAIN) {
- realw templ = 0.33333333333333333333f * (duxdxl + duydyl + duzdzl); // 1./3. = 0.33333
+ // computes deviatoric strain attenuation and/or for kernel calculations
+ if(COMPUTE_AND_STORE_STRAIN) {
+ realw templ = 0.33333333333333333333f * (duxdxl + duydyl + duzdzl); // 1./3. = 0.33333
- // local storage: stresses at this current time step
- epsilondev_xx_loc = duxdxl - templ;
- epsilondev_yy_loc = duydyl - templ;
- epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl;
- epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl;
- epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl;
+ // local storage: stresses at this current time step
+ epsilondev_xx_loc = duxdxl - templ;
+ epsilondev_yy_loc = duydyl - templ;
+ epsilondev_xy_loc = 0.5f * duxdyl_plus_duydxl;
+ epsilondev_xz_loc = 0.5f * duzdxl_plus_duxdzl;
+ epsilondev_yz_loc = 0.5f * duzdyl_plus_duydzl;
- if(SIMULATION_TYPE == 3) {
- epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
- }
- }
+ if(SIMULATION_TYPE == 3) {
+ epsilon_trace_over_3[tx + working_element*NGLL3] = templ;
+ }
+ }
}
// compute elements with an elastic isotropic rheology
@@ -1035,7 +1035,7 @@
/* ----------------------------------------------------------------------------------------------- */
void Kernel_2_inner_core(int nb_blocks_to_compute,Mesh* mp,
- realw d_deltat,
+ realw d_deltat,
int d_iphase,
int* d_ibool,
int* d_idoubling,
@@ -1105,10 +1105,10 @@
mp->d_phase_ispec_inner_inner_core,
mp->num_phase_ispec_inner_core,
d_iphase,
- d_deltat,
+ d_deltat,
mp->use_mesh_coloring_gpu,
mp->d_displ_inner_core,
- mp->d_veloc_inner_core,
+ mp->d_veloc_inner_core,
mp->d_accel_inner_core,
d_xix, d_xiy, d_xiz,
d_etax, d_etay, d_etaz,
@@ -1150,7 +1150,7 @@
mp->d_phase_ispec_inner_inner_core,
mp->num_phase_ispec_inner_core,
d_iphase,
- d_deltat,
+ d_deltat,
mp->use_mesh_coloring_gpu,
mp->d_b_displ_inner_core,
mp->d_b_veloc_inner_core,
@@ -1207,8 +1207,8 @@
extern "C"
void FC_FUNC_(compute_forces_inner_core_cuda,
COMPUTE_FORCES_INNER_CORE_CUDA)(long* Mesh_pointer_f,
- realw* deltat,
- int* iphase) {
+ realw* deltat,
+ int* iphase) {
TRACE("compute_forces_inner_core_cuda");
@@ -1275,7 +1275,7 @@
//}
Kernel_2_inner_core(nb_blocks_to_compute,mp,
- *deltat,
+ *deltat,
*iphase,
mp->d_ibool_inner_core + color_offset_nonpadded,
mp->d_idoubling_inner_core + color_offset_ispec,
@@ -1336,7 +1336,7 @@
// no mesh coloring: uses atomic updates
Kernel_2_inner_core(num_elements,mp,
- *deltat,
+ *deltat,
*iphase,
mp->d_ibool_inner_core,
mp->d_idoubling_inner_core,
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/specfem3D_gpu_cuda_method_stubs.c 2012-05-01 14:58:16 UTC (rev 20011)
@@ -161,15 +161,15 @@
void FC_FUNC_(compute_coupling_cmb_fluid_cuda,
COMPUTE_COUPLING_CMB_FLUID_CUDA)(long* Mesh_pointer_f,
- double RHO_TOP_OC,
- realw minus_g_cmb,
- int GRAVITY_VAL) {}
+ double RHO_TOP_OC,
+ realw minus_g_cmb,
+ int GRAVITY_VAL) {}
void FC_FUNC_(compute_coupling_icb_fluid_cuda,
COMPUTE_COUPLING_ICB_FLUID_CUDA)(long* Mesh_pointer_f,
- double RHO_BOTTOM_OC,
- realw minus_g_icb,
- int GRAVITY_VAL) {}
+ double RHO_BOTTOM_OC,
+ realw minus_g_icb,
+ int GRAVITY_VAL) {}
void FC_FUNC_(compute_coupling_ocean_cuda,
COMPUTE_COUPLING_OCEAN_CUDA)(long* Mesh_pointer_f) {}
@@ -181,6 +181,7 @@
void FC_FUNC_(compute_forces_crust_mantle_cuda,
COMPUTE_FORCES_CRUST_MANTLE_CUDA)(long* Mesh_pointer_f,
+ realw* deltat,
int* iphase) {}
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_MPI_interfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_MPI_interfaces.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_MPI_interfaces.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -28,7 +28,7 @@
module create_MPI_interfaces_par
use constants,only: CUSTOM_REAL,NUMFACES_SHARED,NB_SQUARE_EDGES_ONEDIR,NDIM,IMAIN
-
+
! indirect addressing for each message for faces and corners of the chunks
! a given slice can belong to at most one corner and at most two faces
integer :: NGLOB2DMAX_XY
@@ -42,7 +42,7 @@
! number of message types
integer :: NUM_MSG_TYPES
- integer :: NGLOB1D_RADIAL_CM
+ integer :: NGLOB1D_RADIAL_CM
integer :: NGLOB1D_RADIAL_OC
integer :: NGLOB1D_RADIAL_IC
@@ -76,7 +76,7 @@
integer :: NGLOB_CRUST_MANTLE
integer :: NGLOB_INNER_CORE
integer :: NGLOB_OUTER_CORE
-
+
!-----------------------------------------------------------------
! assembly
!-----------------------------------------------------------------
@@ -179,7 +179,7 @@
xstore_outer_core,ystore_outer_core,zstore_outer_core
integer, dimension(:),allocatable :: idoubling_outer_core
integer, dimension(:,:,:,:),allocatable :: ibool_outer_core
-
+
! assembly
integer :: npoin2D_faces_outer_core(NUMFACES_SHARED)
integer, dimension(NB_SQUARE_EDGES_ONEDIR) :: npoin2D_xi_outer_core,npoin2D_eta_outer_core
@@ -248,7 +248,7 @@
! mesh coloring
integer :: num_colors_outer_inner_core,num_colors_inner_inner_core
integer,dimension(:),allocatable :: num_elem_colors_inner_core
-
+
end module create_MPI_interfaces_par
!
@@ -260,7 +260,7 @@
use meshfem3D_par
use create_MPI_interfaces_par
implicit none
-
+
! sets up arrays
call cmi_read_addressing()
@@ -275,15 +275,15 @@
! sets up mesh coloring
call cmi_setup_color_perm()
-
+
! saves interface infos
call cmi_save_interfaces()
-
+
! frees memory
call cmi_free_arrays()
end subroutine create_MPI_interfaces
-
+
!
!-------------------------------------------------------------------------------------------------
!
@@ -297,7 +297,7 @@
! local parameters
integer :: NUM_FACES,NPROC_ONE_DIRECTION
integer :: ier
-
+
! define maximum size for message buffers
! use number of elements found in the mantle since it is the largest region
NGLOB2DMAX_XY = max(NGLOB2DMAX_XMIN_XMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_YMIN_YMAX(IREGION_CRUST_MANTLE))
@@ -328,14 +328,14 @@
iprocto_faces(NUMMSGS_FACES), &
imsg_type(NUMMSGS_FACES),stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating iproc faces arrays')
-
+
! communication pattern for corners between chunks
allocate(iproc_master_corners(NCORNERSCHUNKS), &
iproc_worker1_corners(NCORNERSCHUNKS), &
iproc_worker2_corners(NCORNERSCHUNKS),stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating iproc corner arrays')
-
+
! parameters from header file
NGLOB1D_RADIAL_CM = NGLOB1D_RADIAL(IREGION_CRUST_MANTLE)
NGLOB1D_RADIAL_OC = NGLOB1D_RADIAL(IREGION_OUTER_CORE)
@@ -383,7 +383,7 @@
! crust mantle
allocate(iboolcorner_crust_mantle(NGLOB1D_RADIAL_CM,NUMCORNERS_SHARED))
allocate(iboolleft_xi_crust_mantle(NGLOB2DMAX_XMIN_XMAX_CM), &
- iboolright_xi_crust_mantle(NGLOB2DMAX_XMIN_XMAX_CM))
+ iboolright_xi_crust_mantle(NGLOB2DMAX_XMIN_XMAX_CM))
allocate(iboolleft_eta_crust_mantle(NGLOB2DMAX_YMIN_YMAX_CM), &
iboolright_eta_crust_mantle(NGLOB2DMAX_YMIN_YMAX_CM))
allocate(iboolfaces_crust_mantle(NGLOB2DMAX_XY,NUMFACES_SHARED))
@@ -447,8 +447,8 @@
is_on_a_slice_edge_outer_core(NSPEC_OUTER_CORE), &
stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating temporary is_on_a_slice_edge arrays')
-
-
+
+
! read coordinates of the mesh
! crust mantle
ibool_crust_mantle(:,:,:,:) = -1
@@ -465,7 +465,7 @@
call exit_MPI(myrank,'incorrect global numbering: iboolmax does not equal nglob in crust and mantle')
! outer core
- ibool_outer_core(:,:,:,:) = -1
+ ibool_outer_core(:,:,:,:) = -1
call cmi_read_solver_data(myrank,IREGION_OUTER_CORE, &
NSPEC_OUTER_CORE,NGLOB_OUTER_CORE, &
xstore_outer_core,ystore_outer_core,zstore_outer_core,&
@@ -479,7 +479,7 @@
call exit_MPI(myrank,'incorrect global numbering: iboolmax does not equal nglob in outer core')
! inner core
- ibool_inner_core(:,:,:,:) = -1
+ ibool_inner_core(:,:,:,:) = -1
call cmi_read_solver_data(myrank,IREGION_INNER_CORE, &
NSPEC_INNER_CORE,NGLOB_INNER_CORE, &
xstore_inner_core,ystore_inner_core,zstore_inner_core,&
@@ -493,9 +493,9 @@
! synchronize processes
call sync_all()
-
- end subroutine cmi_read_addressing
+ end subroutine cmi_read_addressing
+
!
!-------------------------------------------------------------------------------------------------
!
@@ -518,7 +518,7 @@
! mantle and crust
if(myrank == 0) then
- write(IMAIN,*)
+ write(IMAIN,*)
write(IMAIN,*) 'crust/mantle region:'
endif
@@ -572,7 +572,7 @@
open(unit=IIN,file=prname(1:len_trim(prname))//'boundary.bin', &
status='old',form='unformatted',action='read',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening boundary.bin file')
-
+
read(IIN) nspec2D_xmin_inner_core
read(IIN) nspec2D_xmax_inner_core
read(IIN) nspec2D_ymin_inner_core
@@ -607,7 +607,7 @@
write(IMAIN,*) 'including central cube'
endif
call sync_all()
-
+
! compute number of messages to expect in cube as well as their size
call comp_central_cube_buffer_size(iproc_xi,iproc_eta,ichunk, &
NPROC_XI,NPROC_ETA,NSPEC2D_BOTTOM(IREGION_INNER_CORE), &
@@ -748,11 +748,13 @@
logical,parameter :: DEBUG_INTERFACES = .false.
! estimates a maximum size of needed arrays
- MAX_NEIGHBOURS = 8 + NCORNERSCHUNKS
+ MAX_NEIGHBOURS = 8 + NCORNERSCHUNKS
+ if( INCLUDE_CENTRAL_CUBE ) MAX_NEIGHBOURS = MAX_NEIGHBOURS + NUMMSGS_FACES
+
allocate(my_neighbours(MAX_NEIGHBOURS), &
nibool_neighbours(MAX_NEIGHBOURS),stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating my_neighbours array')
-
+
! estimates initial maximum ibool array
max_nibool = npoin2D_max_all_CM_IC * NUMFACES_SHARED &
+ non_zero_nb_msgs_theor_in_cube*npoin2D_cube_from_slices
@@ -960,7 +962,7 @@
enddo
call sync_all()
endif
-
+
! checks addressing
call test_MPI_neighbours(IREGION_OUTER_CORE, &
num_interfaces_outer_core,max_nibool_interfaces_outer_core, &
@@ -1023,7 +1025,7 @@
!call write_VTK_glob_points(NGLOB_INNER_CORE, &
! xstore_inner_core,ystore_inner_core,zstore_inner_core, &
! test_flag,filename)
-
+
! debug: idoubling inner core
if( DEBUG_INTERFACES ) then
write(filename,'(a,i6.6)') trim(OUTPUT_FILES)//'/MPI_idoubling_inner_core_proc',myrank
@@ -1033,7 +1035,7 @@
idoubling_inner_core,filename)
call sync_all()
endif
-
+
! including central cube
if(INCLUDE_CENTRAL_CUBE) then
! user output
@@ -1126,7 +1128,7 @@
enddo
call sync_all()
endif
-
+
! checks addressing
call test_MPI_neighbours(IREGION_INNER_CORE, &
num_interfaces_inner_core,max_nibool_interfaces_inner_core, &
@@ -1159,7 +1161,7 @@
!
subroutine cmi_setup_InnerOuter()
-
+
use meshfem3D_par
use create_MPI_interfaces_par
implicit none
@@ -1167,10 +1169,10 @@
! local parameters
real :: percentage_edge
integer :: ier,ispec,iinner,iouter
- ! debug
+ ! debug
character(len=150) :: filename
logical,parameter :: DEBUG_INTERFACES = .false.
-
+
! stores inner / outer elements
!
! note: arrays is_on_a_slice_edge_.. have flags set for elements which need to
@@ -1293,7 +1295,7 @@
ibool_inner_core, &
is_on_a_slice_edge_inner_core,filename)
endif
-
+
end subroutine cmi_setup_InnerOuter
@@ -1302,7 +1304,7 @@
!
subroutine cmi_setup_color_perm()
-
+
use meshfem3D_par
use create_MPI_interfaces_par
implicit none
@@ -1363,11 +1365,11 @@
deallocate(perm)
else
- ! dummy array
+ ! dummy array
allocate(num_elem_colors_outer_core(num_colors_outer_outer_core+num_colors_inner_outer_core),stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating num_elem_colors_outer_core array')
endif
-
+
! inner core
! initializes
num_colors_outer_inner_core = 0
@@ -1393,15 +1395,15 @@
allocate(num_elem_colors_inner_core(num_colors_outer_inner_core+num_colors_inner_inner_core),stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating num_elem_colors_inner_core array')
endif
-
+
end subroutine cmi_setup_color_perm
-
+
!
!-------------------------------------------------------------------------------------------------
!
subroutine cmi_save_interfaces()
-
+
use meshfem3D_par
use create_MPI_interfaces_par
implicit none
@@ -1441,8 +1443,8 @@
end subroutine cmi_save_interfaces
-
-
+
+
!
!-------------------------------------------------------------------------------------------------
!
@@ -1464,14 +1466,14 @@
! crust mantle
deallocate(iboolcorner_crust_mantle)
deallocate(iboolleft_xi_crust_mantle, &
- iboolright_xi_crust_mantle)
+ iboolright_xi_crust_mantle)
deallocate(iboolleft_eta_crust_mantle, &
iboolright_eta_crust_mantle)
deallocate(iboolfaces_crust_mantle)
deallocate(phase_ispec_inner_crust_mantle)
deallocate(num_elem_colors_crust_mantle)
-
+
! outer core
deallocate(iboolcorner_outer_core)
deallocate(iboolleft_xi_outer_core, &
@@ -1511,7 +1513,7 @@
deallocate(num_elem_colors_inner_core)
deallocate(mask_ibool)
-
+
! frees temporary allocated arrays
deallocate(is_on_a_slice_edge_crust_mantle, &
is_on_a_slice_edge_outer_core, &
@@ -1535,25 +1537,25 @@
integer :: iregion_code,myrank
integer :: nspec,nglob
-
+
real(kind=CUSTOM_REAL), dimension(nglob) :: xstore,ystore,zstore
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
integer, dimension(nspec) :: idoubling
logical, dimension(nspec) :: is_on_a_slice_edge
character(len=150) :: LOCAL_PATH
-
+
! local parameters
character(len=150) prname
integer :: ier
-
+
! create the name for the database of the current slide and region
call create_name_database(prname,myrank,iregion_code,LOCAL_PATH)
-
+
open(unit=IIN,file=prname(1:len_trim(prname))//'solver_data_2.bin', &
status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening solver_data_2.bin')
-
+
read(IIN) xstore
read(IIN) ystore
read(IIN) zstore
@@ -1562,7 +1564,7 @@
read(IIN) is_on_a_slice_edge
close(IIN)
-
+
end subroutine cmi_read_solver_data
!
@@ -1601,18 +1603,18 @@
integer :: num_colors_outer,num_colors_inner
integer, dimension(num_colors_outer + num_colors_inner) :: &
num_elem_colors
-
+
! local parameters
character(len=150) prname
integer :: ier
-
+
! create the name for the database of the current slide and region
call create_name_database(prname,myrank,iregion_code,LOCAL_PATH)
-
+
open(unit=IOUT,file=prname(1:len_trim(prname))//'solver_data_mpi.bin', &
status='unknown',action='write',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening solver_data_mpi.bin')
-
+
! MPI interfaces
write(IOUT) num_interfaces
if( num_interfaces > 0 ) then
@@ -1636,5 +1638,5 @@
close(IOUT)
end subroutine cmi_save_solver_data
-
-
+
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_central_cube_buffers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_central_cube_buffers.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_central_cube_buffers.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -106,8 +106,8 @@
! check that the number of points in this slice is correct
if(minval(ibool_inner_core(:,:,:,:)) /= 1 .or. maxval(ibool_inner_core(:,:,:,:)) /= NGLOB_INNER_CORE) &
call exit_MPI(myrank,'incorrect global numbering: iboolmax does not equal nglob in inner core')
-
+
!--- processor to send information to in cube from slices
! four vertical sides first
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_chunk_buffers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_chunk_buffers.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_chunk_buffers.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -74,7 +74,7 @@
integer iproc_eta_slice(0:NPROCTOT-1)
integer NCHUNKS
-
+
! local parameters
integer NGLOB1D_RADIAL
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -129,7 +129,7 @@
! fluid in outer core
case( IREGION_OUTER_CORE )
-
+
! no anisotropy in the fluid, use kappav
! distinguish between single and double precision for reals
@@ -145,7 +145,7 @@
call exit_MPI(myrank,'wrong region code')
end select
-
+
enddo
enddo
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_meshes.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_meshes.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_meshes.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -30,7 +30,7 @@
use meshfem3D_par
implicit none
-
+
! local parameters
! parameters needed to store the radii of the grid points
! in the spherically symmetric Earth
@@ -39,7 +39,7 @@
! arrays with the mesh in double precision
double precision, dimension(:,:,:,:), allocatable :: xstore,ystore,zstore
integer :: ier
-
+
! get addressing for this process
ichunk = ichunk_slice(myrank)
iproc_xi = iproc_xi_slice(myrank)
@@ -87,13 +87,13 @@
zstore(NGLLX,NGLLY,NGLLZ,NSPEC(iregion_code)), &
stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating memory for arrays')
-
+
! this for non blocking MPI
allocate(is_on_a_slice_edge(NSPEC(iregion_code)), &
stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating is_on_a_slice_edge array')
-
+
! create all the regions of the mesh
! perform two passes in this part to be able to save memory
do ipass = 1,2
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -26,9 +26,9 @@
!=====================================================================
module create_regions_mesh_par
-
+
use constants,only: NGLLX,NGLLY,NGLLZ,NGNOD,NGNOD2D,NDIM,NDIM2D
-
+
! topology of the elements
integer, dimension(NGNOD) :: iaddx,iaddy,iaddz
@@ -47,8 +47,8 @@
double precision, dimension(NGNOD2D,NGLLX,NGLLY) :: shape2D_bottom,shape2D_top
double precision, dimension(NDIM2D,NGNOD2D,NGLLY,NGLLZ) :: dershape2D_x
double precision, dimension(NDIM2D,NGNOD2D,NGLLX,NGLLZ) :: dershape2D_y
- double precision, dimension(NDIM2D,NGNOD2D,NGLLX,NGLLY) :: dershape2D_bottom,dershape2D_top
-
+ double precision, dimension(NDIM2D,NGNOD2D,NGLLX,NGLLY) :: dershape2D_bottom,dershape2D_top
+
end module create_regions_mesh_par
@@ -80,7 +80,7 @@
! creates the different regions of the mesh
use meshfem3D_models_par
- use create_regions_mesh_par
+ use create_regions_mesh_par
implicit none
! code for the four regions of the mesh
@@ -102,7 +102,7 @@
double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
! meshing parameters
- double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: rmins,rmaxs
+ double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: rmins,rmaxs
integer :: iproc_xi,iproc_eta,ichunk
@@ -115,7 +115,7 @@
integer :: NEX_XI,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
integer :: NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP
integer :: NPROC_XI,NPROC_ETA
-
+
! this to cut the doubling brick
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_EDGES_ONEDIR) :: NSPEC2D_XI_FACE,NSPEC2D_ETA_FACE
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_CORNERS) :: NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER
@@ -131,16 +131,16 @@
double precision :: ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD
logical :: SAVE_MESH_FILES
-
+
integer :: NCHUNKS
logical :: INCLUDE_CENTRAL_CUBE,ABSORBING_CONDITIONS
double precision :: R_CENTRAL_CUBE,RICB
- double precision :: RHO_OCEANS
+ double precision :: RHO_OCEANS
double precision :: RCMB,R670,RMOHO,RMOHO_FICTITIOUS_IN_MESHER, &
RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN
-
+
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
@@ -252,7 +252,7 @@
logical, dimension(:), allocatable :: ispec_is_tiso
integer i,j,k,ispec
-
+
! name of the database file
character(len=150) :: prname
character(len=150) :: errmsg
@@ -685,9 +685,7 @@
call get_global(nspec,xp,yp,zp,ibool,locval,ifseg,nglob,npointot)
- deallocate(xp,stat=ier); if(ier /= 0) stop 'error in deallocate'
- deallocate(yp,stat=ier); if(ier /= 0) stop 'error in deallocate'
- deallocate(zp,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(xp,yp,zp)
! check that number of points found equals theoretical value
if(nglob /= nglob_theor) then
@@ -755,8 +753,7 @@
! rhostore,kappavstore,muvstore,Qmu_store,ATTENUATION)
endif
- deallocate(locval,stat=ier); if(ier /= 0) stop 'error in deallocate'
- deallocate(ifseg,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(locval,ifseg)
! only create mass matrix and save all the final arrays in the second pass
case( 2 )
@@ -839,8 +836,7 @@
size(tau_e_store,2),size(tau_e_store,3),size(tau_e_store,4),size(tau_e_store,5),&
ABSORBING_CONDITIONS,SAVE_MESH_FILES,ispec_is_tiso)
- deallocate(rmass,stat=ier); if(ier /= 0) stop 'error in deallocate'
- deallocate(rmass_ocean_load,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(rmass,rmass_ocean_load)
! boundary mesh
if (SAVE_BOUNDARY_MESH .and. iregion_code == IREGION_CRUST_MANTLE) then
@@ -887,10 +883,11 @@
deallocate(stretch_tab)
deallocate(perm_layer)
- ! deallocate these arrays after each pass because they have a different size in each pass to save memory
- deallocate(xixstore,xiystore,xizstore,stat=ier); if(ier /= 0) stop 'error in deallocate'
- deallocate(etaxstore,etaystore,etazstore,stat=ier); if(ier /= 0) stop 'error in deallocate'
- deallocate(gammaxstore,gammaystore,gammazstore,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ ! deallocate these arrays after each pass
+ ! because they have a different size in each pass to save memory
+ deallocate(xixstore,xiystore,xizstore)
+ deallocate(etaxstore,etaystore,etazstore)
+ deallocate(gammaxstore,gammaystore,gammazstore)
! deallocate arrays
deallocate(rhostore,dvpstore,kappavstore,kappahstore)
@@ -918,5 +915,8 @@
deallocate(normal_moho,normal_400,normal_670)
deallocate(jacobian2D_moho,jacobian2D_400,jacobian2D_670)
+ ! synchronizes processes
+ call sync_all()
+
end subroutine create_regions_mesh
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/finalize_mesher.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/finalize_mesher.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/finalize_mesher.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -95,10 +95,11 @@
write(IMAIN,*)
! evaluate the amount of static memory needed by the solver
- call memory_eval(OCEANS,ABSORBING_CONDITIONS,ATTENUATION,ANISOTROPIC_3D_MANTLE,&
- TRANSVERSE_ISOTROPY,ANISOTROPIC_INNER_CORE,ROTATION,&
- ONE_CRUST,doubling_index,this_region_has_a_doubling,&
- ner,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_sampling_array,&
+ call memory_eval(OCEANS,ABSORBING_CONDITIONS,ATTENUATION,ANISOTROPIC_3D_MANTLE, &
+ TRANSVERSE_ISOTROPY,ANISOTROPIC_INNER_CORE,ROTATION,TOPOGRAPHY, &
+ ONE_CRUST,doubling_index,this_region_has_a_doubling, &
+ ner,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ ratio_sampling_array,NPROCTOT, &
NSPEC,nglob,SIMULATION_TYPE,MOVIE_VOLUME,SAVE_FORWARD, &
NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
@@ -111,7 +112,9 @@
NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
- NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION,static_memory_size)
+ NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ static_memory_size)
NGLOB1D_RADIAL_TEMP(:) = &
(/maxval(NGLOB1D_RADIAL_CORNER(1,:)),maxval(NGLOB1D_RADIAL_CORNER(2,:)),maxval(NGLOB1D_RADIAL_CORNER(3,:))/)
@@ -119,7 +122,8 @@
! create include file for the solver
call save_header_file(NSPEC,nglob,NEX_XI,NEX_ETA,NPROC,NPROCTOT, &
TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- ELLIPTICITY,GRAVITY,ROTATION,OCEANS,ATTENUATION,ATTENUATION_3D, &
+ ELLIPTICITY,GRAVITY,ROTATION, &
+ TOPOGRAPHY,OCEANS,ATTENUATION,ATTENUATION_3D, &
ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,NCHUNKS, &
INCLUDE_CENTRAL_CUBE,CENTER_LONGITUDE_IN_DEGREES,&
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,NSOURCES,NSTEP, &
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_interfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_interfaces.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_interfaces.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -62,7 +62,7 @@
real(kind=CUSTOM_REAL),dimension(NGLOB),intent(in) :: xstore,ystore,zstore
integer :: NPROCTOT
-
+
! local parameters
integer :: ispec,iglob,j,k
integer :: iface,iedge,icorner
@@ -164,8 +164,10 @@
! updates interfaces array
if( .not. is_done ) then
iinterface = iinterface + 1
- if( iinterface > MAX_NEIGHBOURS ) &
+ if( iinterface > MAX_NEIGHBOURS ) then
+ print*,'error interfaces rank:',myrank,'iinterface = ',iinterface,MAX_NEIGHBOURS
call exit_mpi(myrank,'interface face exceeds MAX_NEIGHBOURS range')
+ endif
! adds as neighbor new interface
my_neighbours(iinterface) = rank
icurrent = iinterface
@@ -298,8 +300,10 @@
! updates interfaces array
if( .not. is_done ) then
iinterface = iinterface + 1
- if( iinterface > MAX_NEIGHBOURS ) &
+ if( iinterface > MAX_NEIGHBOURS ) then
+ print*,'error interfaces rank:',myrank,'iinterface = ',iinterface,MAX_NEIGHBOURS
call exit_mpi(myrank,'interface edge exceeds MAX_NEIGHBOURS range')
+ endif
! adds as neighbor new interface
my_neighbours(iinterface) = rank
icurrent = iinterface
@@ -449,8 +453,10 @@
! updates interfaces array
if( .not. is_done ) then
iinterface = iinterface + 1
- if( iinterface > MAX_NEIGHBOURS ) &
+ if( iinterface > MAX_NEIGHBOURS ) then
+ print*,'error interfaces rank:',myrank,'iinterface = ',iinterface,MAX_NEIGHBOURS
call exit_mpi(myrank,'interface corner exceed MAX_NEIGHBOURS range')
+ endif
! adds as neighbor new interface
my_neighbours(iinterface) = rank
icurrent = iinterface
@@ -567,33 +573,33 @@
is_on_a_slice_edge(:) = work_ispec_is_outer(:)
end subroutine get_MPI_interfaces
-
+
!
!-------------------------------------------------------------------------------------------------
!
-
+
subroutine sort_MPI_interface(myrank,npoin,ibool_n, &
NGLOB,xstore,ystore,zstore)
use constants,only: CUSTOM_REAL
-
+
implicit none
-
+
integer,intent(in) :: myrank,npoin
integer,dimension(npoin),intent(inout) :: ibool_n
integer,intent(in) :: NGLOB
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: xstore,ystore,zstore
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: xstore,ystore,zstore
! local parameters
! arrays for sorting routine
double precision, dimension(:), allocatable :: work
double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
- integer, dimension(:), allocatable :: ibool_selected
+ integer, dimension(:), allocatable :: ibool_selected
integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
logical, dimension(:), allocatable :: ifseg
integer :: nglob_selected,i,ipoin,ier
-
+
! allocate arrays for buffers with maximum size
allocate(ibool_selected(npoin), &
xstore_selected(npoin), &
@@ -607,17 +613,17 @@
iwork(npoin), &
work(npoin),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error sort MPI interface: allocating temporary sorting arrays')
-
+
! sets up working arrays
do i=1,npoin
ipoin = ibool_n(i)
-
+
ibool_selected(i) = ipoin
xstore_selected(i) = xstore(ipoin)
ystore_selected(i) = ystore(ipoin)
zstore_selected(i) = zstore(ipoin)
enddo
-
+
! sort buffer obtained to be conforming with neighbor in other chunk
! sort on x, y and z, the other arrays will be swapped as well
call sort_array_coordinates(npoin,xstore_selected,ystore_selected,zstore_selected, &
@@ -634,7 +640,7 @@
deallocate(ibool_selected,xstore_selected,ystore_selected,zstore_selected, &
ind,ninseg,iglob,locval,ifseg,iwork,work)
-
+
end subroutine sort_MPI_interface
!
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -350,7 +350,7 @@
! setup mpi communication interfaces
call create_MPI_interfaces()
-
+
! outputs mesh infos and saves new header file
call finalize_mesher()
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_models.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_models.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_models.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -44,7 +44,7 @@
!---
use constants
-
+
implicit none
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_aniso_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_aniso_mantle.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_aniso_mantle.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -368,7 +368,7 @@
include "constants.h"
-! model_aniso_mantle_variables
+ ! model_aniso_mantle_variables
type model_aniso_mantle_variables
sequence
double precision beta(14,34,37,73)
@@ -378,15 +378,21 @@
end type model_aniso_mantle_variables
type (model_aniso_mantle_variables) AMM_V
-! model_aniso_mantle_variables
+ ! model_aniso_mantle_variables
- integer nx,ny,np1,np2,ipar,ipa1,ipa,ilat,ilon,il,idep,nfin,nfi0,nf,nri
- double precision xinf,yinf,pxy,ppp,angle,A,A2L,AL,af
- double precision ra(47),pari(14,47)
- double precision bet2(14,34,37,73)
- double precision alph(73,37),ph(73,37)
- character(len=150) glob_prem3sm01, globpreman3sm01
+ ! local parameters
+ integer :: nx,ny,np1,np2,ipar,ipa1,ipa,ilat,ilon,il,idep,nfin,nfi0,nf,nri
+ double precision :: xinf,yinf,pxy,ppp,angle,A,A2L,AL,af
+ double precision :: ra(47),pari(14,47)
+ double precision,dimension(:,:,:,:),allocatable :: bet2 ! bet2(14,34,37,73)
+ double precision :: alph(73,37),ph(73,37)
+ integer :: ier
+ character(len=150) :: glob_prem3sm01, globpreman3sm01
+ ! dynamic allocation
+ allocate(bet2(14,34,37,73),stat=ier)
+ if( ier /= 0 ) stop 'error allocating bet2 array'
+
np1 = 1
np2 = 34
AMM_V%npar1 = (np2 - np1 + 1)
@@ -395,7 +401,8 @@
! glob-prem3sm01: model with rho,A,L,xi-1,1-phi,eta
!
call get_value_string(glob_prem3sm01, 'model.glob_prem3sm01', 'DATA/Montagner_model/glob-prem3sm01')
- open(19,file=glob_prem3sm01,status='old',action='read')
+ open(19,file=glob_prem3sm01,status='old',action='read',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file DATA/Montagner_model/glob-prem3sm01'
!
! read the models
@@ -456,7 +463,8 @@
! normalized, in percents: 100 G/L
!
call get_value_string(globpreman3sm01, 'model.globpreman3sm01', 'DATA/Montagner_model/globpreman3sm01')
- open(unit=15,file=globpreman3sm01,status='old',action='read')
+ open(unit=15,file=globpreman3sm01,status='old',action='read',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file DATA/Montagner_model/globpreman3sm01'
do nf = 7,nfin,2
ipa = nf
@@ -528,8 +536,10 @@
enddo
enddo
- end subroutine read_aniso_mantle_model
+ deallocate(bet2)
+ end subroutine read_aniso_mantle_model
+
!--------------------------------------------------------------------
subroutine lecmod(nri,pari,ra)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_color_perm.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_color_perm.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_color_perm.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -28,7 +28,7 @@
subroutine setup_color_perm(myrank,iregion_code,nspec,nglob, &
ibool,is_on_a_slice_edge,prname, &
npoin2D_xi,npoin2D_eta)
-
+
use constants
use meshfem3D_par,only: NSTEP,DT,NPROC_XI,NPROC_ETA
implicit none
@@ -49,7 +49,7 @@
character(len=150) :: prname
integer :: npoin2D_xi,npoin2D_eta
-
+
! local parameters
integer :: nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer
integer, dimension(:), allocatable :: perm
@@ -174,7 +174,7 @@
open(unit=99,file=prname(1:len_trim(prname))//'values_from_mesher_f90.h', &
status='unknown',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening file values_from_mesher_f90.h')
-
+
write(99,*) 'integer, parameter :: NSPEC = ',nspec
write(99,*) 'integer, parameter :: NGLOB = ',nglob
!!! DK DK use 1000 time steps only for the scaling tests
@@ -211,7 +211,7 @@
open(unit=99,file=prname(1:len_trim(prname))//'values_from_mesher_C.h', &
status='unknown',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening file values_from_mesher_C.h')
-
+
write(99,*) '#define NSPEC ',nspec
write(99,*) '#define NGLOB ',nglob
!! write(99,*) '#define NSTEP ',nstep
@@ -239,7 +239,7 @@
open(unit=99,file=prname(1:len_trim(prname))//'values_from_mesher_nspec_outer.h', &
status='unknown',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening values_from_mesher_nspec_outer.h file')
-
+
write(99,*) '#define NSPEC_OUTER ',nspec_outer_max_global
write(99,*) '// NSPEC_OUTER_min = ',nspec_outer_min_global
write(99,*) '// NSPEC_OUTER_max = ',nspec_outer_max_global
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_model.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_model.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -69,7 +69,7 @@
write(IMAIN,*)
endif
call sync_all()
-
+
end subroutine setup_model
!
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/test_MPI_interfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/test_MPI_interfaces.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/test_MPI_interfaces.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -33,7 +33,7 @@
use constants
use meshfem3D_par,only: NPROCTOT,myrank
- use create_MPI_interfaces_par,only: NGLOB_CRUST_MANTLE,NGLOB_OUTER_CORE,NGLOB_INNER_CORE
+ use create_MPI_interfaces_par,only: NGLOB_CRUST_MANTLE,NGLOB_OUTER_CORE,NGLOB_INNER_CORE
implicit none
include 'mpif.h'
@@ -42,7 +42,7 @@
integer,intent(in) :: num_interfaces,max_nibool_interfaces
integer,dimension(num_interfaces),intent(in) :: my_neighbours,nibool_interfaces
integer,dimension(max_nibool_interfaces,num_interfaces),intent(in):: ibool_interfaces
-
+
! local parameters
integer,dimension(:),allocatable :: dummy_i
integer,dimension(:,:),allocatable :: test_interfaces
@@ -51,10 +51,10 @@
integer :: ineighbour,iproc,inum,i,j,ier,ipoints,max_num,iglob
logical :: is_okay
logical,dimension(:),allocatable :: mask
-
+
! daniel: debug output
!do iproc=0,NPROCTOT-1
- ! if( myrank == iproc ) then
+ ! if( myrank == iproc ) then
! print*, 'mpi rank',myrank,'interfaces : ',num_interfaces,'region',iregion_code
! do j=1,num_interfaces
! print*, ' my_neighbours: ',my_neighbours(j),nibool_interfaces(j)
@@ -63,13 +63,13 @@
! endif
! call sync_all()
!enddo
-
+
! checks maximum number of interface points
if( max_nibool_interfaces == 0 .and. NPROCTOT > 1 ) then
print*,'test MPI: rank ',myrank,'max_nibool_interfaces is zero'
call exit_mpi(myrank,'error test max_nibool_interfaces zero')
endif
-
+
! allocates global mask
select case(iregion_code)
case( IREGION_CRUST_MANTLE )
@@ -81,7 +81,7 @@
case default
call exit_mpi(myrank,'error test MPI: iregion_code not recognized')
end select
-
+
! test ibool entries
! (must be non-zero and unique)
do i = 1,num_interfaces
@@ -90,19 +90,19 @@
print*,'error test MPI: rank',myrank,'nibool values:',nibool_interfaces(i),max_nibool_interfaces
call exit_mpi(myrank,'error test MPI: nibool exceeds max_nibool_interfaces')
endif
-
+
mask(:) = .false.
-
- ! ibool entries
+
+ ! ibool entries
do j = 1,nibool_interfaces(i)
iglob = ibool_interfaces(j,i)
-
+
! checks zero entry
if( iglob <= 0 ) then
print*,'error test MPI: rank ',myrank,'ibool value:',iglob,'interface:',i,'point:',j
call exit_mpi(myrank,'error test MPI: ibool values invalid')
endif
-
+
! checks duplicate
if( j < nibool_interfaces(i) ) then
if( iglob == ibool_interfaces(j+1,i) ) then
@@ -110,18 +110,18 @@
call exit_mpi(myrank,'error test MPI: ibool duplicates')
endif
endif
-
+
! checks if unique global value
if( .not. mask(iglob) ) then
mask(iglob) = .true.
else
print*,'error test MPI: rank',myrank,'ibool masked:',iglob,'interface:',i,'point:',j
call exit_mpi(myrank,'error test MPI: ibool masked already')
- endif
+ endif
enddo
enddo
deallocate(mask)
-
+
! checks neighbors
! gets maximum interfaces from all processes
call MPI_REDUCE(num_interfaces,max_num,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,ier)
@@ -249,7 +249,7 @@
integer :: inum,icount,ival
integer :: num_unique,num_max_valence
integer,dimension(:),allocatable :: valence
-
+
! crust mantle
allocate(test_flag_vector(NDIM,NGLOB_CRUST_MANTLE),stat=ier)
if( ier /= 0 ) stop 'error allocating array test_flag etc.'
@@ -258,10 +258,10 @@
! points defined by interfaces
valence(:) = 0
- test_flag_vector(:,:) = 0.0
+ test_flag_vector(:,:) = 0.0
do i=1,num_interfaces_crust_mantle
do j=1,nibool_interfaces_crust_mantle(i)
- iglob = ibool_interfaces_crust_mantle(j,i)
+ iglob = ibool_interfaces_crust_mantle(j,i)
! sets flag on
test_flag_vector(1,iglob) = 1.0_CUSTOM_REAL
! counts valence (occurrences)
@@ -269,19 +269,19 @@
enddo
enddo
! total number of interface points
- i = sum(nibool_interfaces_crust_mantle)
+ i = sum(nibool_interfaces_crust_mantle)
call MPI_REDUCE(i,inum,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ier)
-
+
! total number of unique points (some could be shared between different processes)
i = nint( sum(test_flag_vector) )
num_unique= i
call MPI_REDUCE(i,icount,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ier)
-
+
! maximum valence
i = maxval( valence(:) )
num_max_valence = i
call MPI_REDUCE(i,ival,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,ier)
-
+
! user output
if( myrank == 0 ) then
write(IMAIN,*) ' total MPI interface points : ',inum
@@ -307,7 +307,7 @@
do iglob=1,NGLOB_CRUST_MANTLE
! only counts flags with MPI contributions
if( test_flag_vector(1,iglob) > 0.0 ) i = i + 1
-
+
! checks valence
if( valence(iglob) /= nint(test_flag_vector(1,iglob)) .or. &
valence(iglob) /= nint(test_flag_vector(2,iglob)) .or. &
@@ -320,9 +320,9 @@
! checks within slice
if( i /= num_unique ) then
print*,'error test crust mantle : rank',myrank,'unique mpi points:',i,num_unique
- call exit_mpi(myrank,'error MPI assembly crust mantle')
+ call exit_mpi(myrank,'error MPI assembly crust mantle')
endif
-
+
! total number of assembly points
call MPI_REDUCE(i,inum,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ier)
@@ -333,10 +333,10 @@
print*,'error crust mantle : total mpi points:',myrank,'total: ',inum,icount
call exit_mpi(myrank,'error MPI assembly crust mantle')
endif
-
+
! user output
write(IMAIN,*) ' total unique MPI interface points:',inum
- write(IMAIN,*)
+ write(IMAIN,*)
endif
deallocate(test_flag_vector)
@@ -384,20 +384,20 @@
enddo
i = sum(nibool_interfaces_outer_core)
call MPI_REDUCE(i,inum,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ier)
-
+
i = nint( sum(test_flag) )
num_unique = i
call MPI_REDUCE(i,icount,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ier)
-
+
! maximum valence
i = maxval( valence(:) )
num_max_valence = i
call MPI_REDUCE(i,ival,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,ier)
-
+
if( myrank == 0 ) then
write(IMAIN,*) ' total MPI interface points : ',inum
write(IMAIN,*) ' unique MPI interface points: ',icount
- write(IMAIN,*) ' maximum valence : ',ival
+ write(IMAIN,*) ' maximum valence : ',ival
endif
! initialized for assembly
@@ -424,15 +424,15 @@
if( valence(iglob) /= nint(test_flag(iglob)) ) then
print*,'error test MPI: rank',myrank,'valence:',valence(iglob),'flag:',test_flag(iglob)
call exit_mpi(myrank,'error test outer core valence')
- endif
+ endif
enddo
! checks within slice
if( i /= num_unique ) then
print*,'error test outer core : rank',myrank,'unique mpi points:',i,num_unique
- call exit_mpi(myrank,'error MPI assembly outer core')
+ call exit_mpi(myrank,'error MPI assembly outer core')
endif
-
+
call MPI_REDUCE(i,inum,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ier)
! output
@@ -445,7 +445,7 @@
! user output
write(IMAIN,*) ' total assembled MPI interface points:',inum
- write(IMAIN,*)
+ write(IMAIN,*)
endif
deallocate(test_flag)
@@ -472,7 +472,7 @@
real(kind=CUSTOM_REAL),dimension(:,:),allocatable :: test_flag_vector
integer :: i,j,iglob,ier
integer :: inum,icount,ival
- integer :: num_unique,num_max_valence
+ integer :: num_unique,num_max_valence
integer,dimension(:),allocatable :: valence
! inner core
@@ -482,7 +482,7 @@
if( ier /= 0 ) stop 'error allocating array valence'
! points defined by interfaces
- valence(:) = 0
+ valence(:) = 0
test_flag_vector(:,:) = 0.0
do i=1,num_interfaces_inner_core
do j=1,nibool_interfaces_inner_core(i)
@@ -495,20 +495,20 @@
enddo
i = sum(nibool_interfaces_inner_core)
call MPI_REDUCE(i,inum,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ier)
-
+
i = nint( sum(test_flag_vector) )
- num_unique= i
+ num_unique= i
call MPI_REDUCE(i,icount,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ier)
-
+
! maximum valence
i = maxval( valence(:) )
num_max_valence = i
call MPI_REDUCE(i,ival,1,MPI_INTEGER,MPI_MAX,0,MPI_COMM_WORLD,ier)
-
+
if( myrank == 0 ) then
write(IMAIN,*) ' total MPI interface points : ',inum
write(IMAIN,*) ' unique MPI interface points: ',icount
- write(IMAIN,*) ' maximum valence : ',ival
+ write(IMAIN,*) ' maximum valence : ',ival
endif
! initializes for assembly
@@ -537,15 +537,15 @@
print*,'error test MPI: rank',myrank,'valence:',valence(iglob),'flag:',test_flag_vector(:,:)
call exit_mpi(myrank,'error test MPI inner core valence')
endif
-
+
enddo
! checks within slice
if( i /= num_unique ) then
print*,'error test inner core : rank',myrank,'unique mpi points:',i,num_unique
- call exit_mpi(myrank,'error MPI assembly inner core')
+ call exit_mpi(myrank,'error MPI assembly inner core')
endif
-
+
call MPI_REDUCE(i,inum,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ier)
if( myrank == 0 ) then
@@ -557,7 +557,7 @@
! user output
write(IMAIN,*) ' total assembled MPI interface points:',inum
- write(IMAIN,*)
+ write(IMAIN,*)
endif
deallocate(test_flag_vector)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/define_all_layers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/define_all_layers.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/define_all_layers.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -56,7 +56,7 @@
integer NER_CRUST,NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
NER_TOP_CENTRAL_CUBE_ICB
-
+
! radii
double precision RMIDDLE_CRUST,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
R_CENTRAL_CUBE,RMOHO_FICTITIOUS_IN_MESHER,R80_FICTITIOUS_IN_MESHER
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/get_timestep_and_layers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/get_timestep_and_layers.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/get_timestep_and_layers.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -417,7 +417,7 @@
if( .false. ) then
if( THREE_D_MODEL == THREE_D_MODEL_PPM ) DT = DT * (1.d0 - 0.2d0)
endif
-
+
! takes a 5% safety margin on the maximum stable time step
! which was obtained by trial and error
DT = DT * (1.d0 - 0.05d0)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/memory_eval.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/memory_eval.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/memory_eval.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -28,42 +28,47 @@
! compute the approximate amount of static memory needed to run the solver
subroutine memory_eval(OCEANS,ABSORBING_CONDITIONS,ATTENUATION,ANISOTROPIC_3D_MANTLE,&
- TRANSVERSE_ISOTROPY,ANISOTROPIC_INNER_CORE,ROTATION,&
- ONE_CRUST,doubling_index,this_region_has_a_doubling,&
- ner,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_sampling_array,&
- NSPEC,nglob,SIMULATION_TYPE,MOVIE_VOLUME,SAVE_FORWARD, &
- NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
- NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
- NSPEC_INNER_CORE_ATTENUATION, &
- NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
- NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
- NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
- NSPEC_CRUST_MANTLE_ADJOINT, &
- NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
- NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
- NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
- NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
- NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION,static_memory_size)
+ TRANSVERSE_ISOTROPY,ANISOTROPIC_INNER_CORE,ROTATION,TOPOGRAPHY,&
+ ONE_CRUST,doubling_index,this_region_has_a_doubling,&
+ ner,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ ratio_sampling_array,NPROCTOT, &
+ NSPEC,nglob,SIMULATION_TYPE,MOVIE_VOLUME,SAVE_FORWARD, &
+ NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
+ NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
+ NSPEC_INNER_CORE_ATTENUATION, &
+ NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, &
+ NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, &
+ NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, &
+ NSPEC_CRUST_MANTLE_ADJOINT, &
+ NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, &
+ NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, &
+ NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, &
+ NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
+ NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ static_memory_size)
implicit none
include "constants.h"
-! input
+ ! input
logical, intent(in) :: TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- ROTATION,ATTENUATION,ONE_CRUST,OCEANS,ABSORBING_CONDITIONS,MOVIE_VOLUME,SAVE_FORWARD
- integer, dimension(MAX_NUM_REGIONS), intent(in) :: NSPEC, nglob
+ ROTATION,TOPOGRAPHY, &
+ ATTENUATION,ONE_CRUST,OCEANS,ABSORBING_CONDITIONS, &
+ MOVIE_VOLUME,SAVE_FORWARD
+ integer, dimension(MAX_NUM_REGIONS), intent(in) :: NSPEC, nglob, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP
+
integer, intent(in) :: NEX_PER_PROC_XI,NEX_PER_PROC_ETA,SIMULATION_TYPE
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS), intent(in) :: doubling_index
logical, dimension(MAX_NUMBER_OF_MESH_LAYERS), intent(in) :: this_region_has_a_doubling
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS), intent(in) :: ner,ratio_sampling_array
+ integer, intent(in) :: NPROCTOT
-! output
+ ! output
double precision, intent(out) :: static_memory_size
-! variables
- integer :: ilayer,NUMBER_OF_MESH_LAYERS,ner_without_doubling,ispec_aniso
-
integer, intent(out) :: NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, &
NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUAT, &
NSPEC_INNER_CORE_ATTENUATION, &
@@ -77,7 +82,10 @@
NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION
-! generate the elements in all the regions of the mesh
+ ! local variables
+ integer :: ilayer,NUMBER_OF_MESH_LAYERS,ner_without_doubling,ispec_aniso
+
+ ! generate the elements in all the regions of the mesh
ispec_aniso = 0
if (ONE_CRUST) then
@@ -86,23 +94,22 @@
NUMBER_OF_MESH_LAYERS = MAX_NUMBER_OF_MESH_LAYERS
endif
-! count anisotropic elements
+ ! count anisotropic elements
do ilayer = 1, NUMBER_OF_MESH_LAYERS
- if(doubling_index(ilayer) == IFLAG_220_80 .or. doubling_index(ilayer) == IFLAG_80_MOHO) then
- ner_without_doubling = ner(ilayer)
- if(this_region_has_a_doubling(ilayer)) then
- ner_without_doubling = ner_without_doubling - 2
- ispec_aniso = ispec_aniso + &
- (NSPEC_DOUBLING_SUPERBRICK*(NEX_PER_PROC_XI/ratio_sampling_array(ilayer)/2)* &
- (NEX_PER_PROC_ETA/ratio_sampling_array(ilayer)/2))
- endif
- ispec_aniso = ispec_aniso + &
- ((NEX_PER_PROC_XI/ratio_sampling_array(ilayer))*(NEX_PER_PROC_ETA/ratio_sampling_array(ilayer))*ner_without_doubling)
+ if(doubling_index(ilayer) == IFLAG_220_80 .or. doubling_index(ilayer) == IFLAG_80_MOHO) then
+ ner_without_doubling = ner(ilayer)
+ if(this_region_has_a_doubling(ilayer)) then
+ ner_without_doubling = ner_without_doubling - 2
+ ispec_aniso = ispec_aniso + &
+ (NSPEC_DOUBLING_SUPERBRICK*(NEX_PER_PROC_XI/ratio_sampling_array(ilayer)/2)* &
+ (NEX_PER_PROC_ETA/ratio_sampling_array(ilayer)/2))
endif
+ ispec_aniso = ispec_aniso + &
+ ((NEX_PER_PROC_XI/ratio_sampling_array(ilayer))*(NEX_PER_PROC_ETA/ratio_sampling_array(ilayer))*ner_without_doubling)
+ endif
enddo
-! define static size of the arrays whose size depends on logical tests
-
+ ! define static size of the arrays whose size depends on logical tests
if(ANISOTROPIC_INNER_CORE) then
NSPECMAX_ANISO_IC = NSPEC(IREGION_INNER_CORE)
else
@@ -117,10 +124,10 @@
NSPECMAX_ISO_MANTLE = NSPEC(IREGION_CRUST_MANTLE)
if(TRANSVERSE_ISOTROPY) then
-! note: the number of transverse isotropic elements is ispec_aniso
-! however for transverse isotropic kernels, the arrays muhstore,kappahstore,eta_anisostore,
-! will be needed for the crust_mantle region everywhere still...
-! originally: NSPECMAX_TISO_MANTLE = ispec_aniso
+ ! note: the number of transverse isotropic elements is ispec_aniso
+ ! however for transverse isotropic kernels, the arrays muhstore,kappahstore,eta_anisostore,
+ ! will be needed for the crust_mantle region everywhere still...
+ ! originally: NSPECMAX_TISO_MANTLE = ispec_aniso
NSPECMAX_TISO_MANTLE = NSPEC(IREGION_CRUST_MANTLE)
else
NSPECMAX_TISO_MANTLE = 1
@@ -129,7 +136,7 @@
NSPECMAX_ANISO_MANTLE = 1
endif
-! if attenuation is off, set dummy size of arrays to one
+ ! if attenuation is off, set dummy size of arrays to one
if(ATTENUATION) then
NSPEC_CRUST_MANTLE_ATTENUAT = NSPEC(IREGION_CRUST_MANTLE)
NSPEC_INNER_CORE_ATTENUATION = NSPEC(IREGION_INNER_CORE)
@@ -190,7 +197,7 @@
NSPEC_OUTER_CORE_ROT_ADJOINT = 1
endif
-! if absorbing conditions are off, set dummy size of arrays to one
+ ! if absorbing conditions are off, set dummy size of arrays to one
if(ABSORBING_CONDITIONS) then
NSPEC_CRUST_MANTLE_STACEY = NSPEC(IREGION_CRUST_MANTLE)
NSPEC_OUTER_CORE_STACEY = NSPEC(IREGION_OUTER_CORE)
@@ -199,9 +206,9 @@
NSPEC_OUTER_CORE_STACEY = 1
endif
-! if oceans are off, set dummy size of arrays to one
+ ! if oceans are off, set dummy size of arrays to one
if(OCEANS) then
- NGLOB_CRUST_MANTLE_OCEANS = NGLOB(IREGION_CRUST_MANTLE)
+ NGLOB_CRUST_MANTLE_OCEANS = nglob(IREGION_CRUST_MANTLE)
else
NGLOB_CRUST_MANTLE_OCEANS = 1
endif
@@ -216,151 +223,247 @@
static_memory_size = 0.d0
-! R_memory_crust_mantle
- static_memory_size = static_memory_size + 5.d0*dble(N_SLS)*dble(NGLLX)* &
- dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_ATTENUAT*dble(CUSTOM_REAL)
+! crust/mantle
-! R_memory_inner_core
- static_memory_size = static_memory_size + 5.d0*dble(N_SLS)*dble(NGLLX)* &
- dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_ATTENUATION*dble(CUSTOM_REAL)
+ ! ibool_crust_mantle
+ static_memory_size = static_memory_size + &
+ dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_CRUST_MANTLE)*dble(SIZE_INTEGER)
-! 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
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_CRUST_MANTLE)*9.d0*dble(CUSTOM_REAL)
+ ! 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
+ static_memory_size = static_memory_size + &
+ 9.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
-! ibool_crust_mantle
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_CRUST_MANTLE)*dble(SIZE_INTEGER)
+ ! xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle,rmass_crust_mantle
+ static_memory_size = static_memory_size + &
+ 4.d0*nglob(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
-! xix_outer_core,xiy_outer_core,xiz_outer_core,
-! etax_outer_core,etay_outer_core,etaz_outer_core,
-! gammax_outer_core,gammay_outer_core,gammaz_outer_core
-! rhostore_outer_core,kappavstore_outer_core
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_OUTER_CORE)*11.d0*dble(CUSTOM_REAL)
+ ! rhostore_crust_mantle,kappavstore_crust_mantle,muvstore_crust_mantle
+ static_memory_size = static_memory_size + &
+ 3.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPECMAX_ISO_MANTLE*dble(CUSTOM_REAL)
-! ibool_outer_core
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_OUTER_CORE)*dble(SIZE_INTEGER)
+ ! kappahstore_crust_mantle,muhstore_crust_mantle,eta_anisostore_crust_mantle
+ static_memory_size = static_memory_size + &
+ 3.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPECMAX_TISO_MANTLE*dble(CUSTOM_REAL)
-! idoubling_inner_core
- static_memory_size = static_memory_size + NSPEC(IREGION_INNER_CORE)*dble(SIZE_INTEGER)
+ ! c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle,
+ ! c14store_crust_mantle,c15store_crust_mantle,c16store_crust_mantle,
+ ! c22store_crust_mantle,c23store_crust_mantle,c24store_crust_mantle,
+ ! c25store_crust_mantle,c26store_crust_mantle,c33store_crust_mantle,
+ ! c34store_crust_mantle,c35store_crust_mantle,c36store_crust_mantle,
+ ! c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle,
+ ! c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle
+ static_memory_size = static_memory_size + &
+ 21.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPECMAX_ANISO_MANTLE*dble(CUSTOM_REAL)
-! ispec_is_tiso_crust_mantle
+ ! ispec_is_tiso_crust_mantle
static_memory_size = static_memory_size + NSPEC(IREGION_CRUST_MANTLE)*dble(SIZE_LOGICAL)
-! ispec_is_tiso_outer_core
- static_memory_size = static_memory_size + NSPEC(IREGION_OUTER_CORE)*dble(SIZE_LOGICAL)
-! ispec_is_tiso_inner_core
- static_memory_size = static_memory_size + NSPEC(IREGION_INNER_CORE)*dble(SIZE_LOGICAL)
-! xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle,rmass_crust_mantle
- static_memory_size = static_memory_size + NGLOB(IREGION_CRUST_MANTLE)*4.d0*dble(CUSTOM_REAL)
+ ! displ_crust_mantle,veloc_crust_mantle,accel_crust_mantle
+ static_memory_size = static_memory_size + &
+ 3.d0*dble(NDIM)*nglob(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
-! rhostore_crust_mantle,kappavstore_crust_mantle,muvstore_crust_mantle
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPECMAX_ISO_MANTLE*3.d0*dble(CUSTOM_REAL)
+ ! one_minus_sum_beta_crust_mantle, factor_scale_crust_mantle
+ static_memory_size = static_memory_size + &
+ 2.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_ATTENUAT*dble(CUSTOM_REAL)
-! kappahstore_crust_mantle,muhstore_crust_mantle,eta_anisostore_crust_mantle
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPECMAX_TISO_MANTLE*3.d0*dble(CUSTOM_REAL)
+ ! factor_common_crust_mantle
+ static_memory_size = static_memory_size + &
+ dble(N_SLS)*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_ATTENUAT*dble(CUSTOM_REAL)
-! c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle,
-! c14store_crust_mantle,c15store_crust_mantle,c16store_crust_mantle,
-! c22store_crust_mantle,c23store_crust_mantle,c24store_crust_mantle,
-! c25store_crust_mantle,c26store_crust_mantle,c33store_crust_mantle,
-! c34store_crust_mantle,c35store_crust_mantle,c36store_crust_mantle,
-! c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle,
-! c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPECMAX_ANISO_MANTLE*21.d0*dble(CUSTOM_REAL)
+ ! R_memory_crust_mantle (R_xx, R_yy, ..)
+ static_memory_size = static_memory_size + &
+ 5.d0*dble(N_SLS)*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_ATTENUAT*dble(CUSTOM_REAL)
-! displ_crust_mantle,veloc_crust_mantle,accel_crust_mantle
- static_memory_size = static_memory_size + dble(NDIM)*NGLOB(IREGION_CRUST_MANTLE)*3.d0*dble(CUSTOM_REAL)
+ ! add arrays used to save strain for attenuation or for adjoint runs
+ ! epsilondev_crust_mantle
+ static_memory_size = static_memory_size + &
+ 5.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_STR_OR_ATT*dble(CUSTOM_REAL)
-! xstore_outer_core, ystore_outer_core, zstore_outer_core, rmass_outer_core, displ_outer_core, veloc_outer_core, accel_outer_core
- static_memory_size = static_memory_size + NGLOB(IREGION_OUTER_CORE)*7.d0*dble(CUSTOM_REAL)
+ ! eps_trace_over_3_crust_mantle
+ static_memory_size = static_memory_size + &
+ dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_STRAIN_ONLY*dble(CUSTOM_REAL)
-! ibool_inner_core
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_INNER_CORE)*dble(SIZE_INTEGER)
+ ! normal_bottom_crust_mantle
+ static_memory_size = static_memory_size + &
+ dble(NDIM)*dble(NGLLX)*dble(NGLLY)*NSPEC2D_BOTTOM(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
-! 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,
-! rhostore_inner_core,kappavstore_inner_core,muvstore_inner_core
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_INNER_CORE)*12.d0*dble(CUSTOM_REAL)
+ ! normal_top_crust_mantle
+ static_memory_size = static_memory_size + &
+ dble(NDIM)*dble(NGLLX)*dble(NGLLY)*NSPEC2D_TOP(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
-! xstore_inner_core,ystore_inner_core,zstore_inner_core,rmass_inner_core
- static_memory_size = static_memory_size + NGLOB(IREGION_INNER_CORE)*4.d0*dble(CUSTOM_REAL)
-! c11store_inner_core,c33store_inner_core,c12store_inner_core,c13store_inner_core,c44store_inner_core
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPECMAX_ANISO_IC*5.d0*dble(CUSTOM_REAL)
+ ! rho_vp_crust_mantle,rho_vs_crust_mantle
+ static_memory_size = static_memory_size + &
+ 2.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_STACEY*dble(CUSTOM_REAL)
-! displ_inner_core,veloc_inner_core,accel_inner_core
- static_memory_size = static_memory_size + dble(NDIM)*NGLOB(IREGION_INNER_CORE)*3.d0*dble(CUSTOM_REAL)
-! A_array_rotation,B_array_rotation
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_OUTER_CORE_ROTATION*2.d0*dble(CUSTOM_REAL)
+! inner core
- if(ABSORBING_CONDITIONS) then
+ ! ibool_inner_core
+ static_memory_size = static_memory_size + &
+ dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_INNER_CORE)*dble(SIZE_INTEGER)
-! rho_vp_crust_mantle,rho_vs_crust_mantle
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_CRUST_MANTLE)*2.d0*dble(CUSTOM_REAL)
+ ! 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,
+ ! rhostore_inner_core,kappavstore_inner_core,muvstore_inner_core
+ static_memory_size = static_memory_size + &
+ 12.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_INNER_CORE)*dble(CUSTOM_REAL)
-! vp_outer_core
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_OUTER_CORE)*dble(CUSTOM_REAL)
+ ! xstore_inner_core,ystore_inner_core,zstore_inner_core,rmass_inner_core
+ static_memory_size = static_memory_size + &
+ 4.d0*nglob(IREGION_INNER_CORE)*dble(CUSTOM_REAL)
- endif
+ ! c11store_inner_core,c33store_inner_core,c12store_inner_core,c13store_inner_core,c44store_inner_core
+ static_memory_size = static_memory_size + &
+ 5.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPECMAX_ANISO_IC*dble(CUSTOM_REAL)
- if(OCEANS) then
+ ! idoubling_inner_core
+ static_memory_size = static_memory_size + NSPEC(IREGION_INNER_CORE)*dble(SIZE_INTEGER)
-! rmass_ocean_load
- static_memory_size = static_memory_size + NGLOB(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
+ ! displ_inner_core,veloc_inner_core,accel_inner_core
+ static_memory_size = static_memory_size + &
+ 3.d0*dble(NDIM)*nglob(IREGION_INNER_CORE)*dble(CUSTOM_REAL)
-! updated_dof_ocean_load
- static_memory_size = static_memory_size + NGLOB(IREGION_CRUST_MANTLE)*dble(SIZE_LOGICAL)
+ ! one_minus_sum_beta_inner_core, factor_scale_inner_core
+ static_memory_size = static_memory_size + &
+ 2.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_ATTENUATION*dble(CUSTOM_REAL)
+ ! factor_common_inner_core
+ static_memory_size = static_memory_size + &
+ dble(N_SLS)*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_ATTENUATION*dble(CUSTOM_REAL)
+
+ ! R_memory_inner_core
+ static_memory_size = static_memory_size + &
+ 5.d0*dble(N_SLS)*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_ATTENUATION*dble(CUSTOM_REAL)
+
+ ! add arrays used to save strain for attenuation or for adjoint runs
+ ! epsilondev_inner_core
+ static_memory_size = static_memory_size + &
+ 5.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_STR_OR_ATT*dble(CUSTOM_REAL)
+
+ ! eps_trace_over_3_inner_core
+ static_memory_size = static_memory_size + &
+ dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_STRAIN_ONLY*dble(CUSTOM_REAL)
+
+! outer core
+
+ ! ibool_outer_core
+ static_memory_size = static_memory_size + &
+ dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_OUTER_CORE)*dble(SIZE_INTEGER)
+
+ ! xix_outer_core,xiy_outer_core,xiz_outer_core,
+ ! etax_outer_core,etay_outer_core,etaz_outer_core,
+ ! gammax_outer_core,gammay_outer_core,gammaz_outer_core
+ ! rhostore_outer_core,kappavstore_outer_core
+ static_memory_size = static_memory_size + &
+ 11.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_OUTER_CORE)*dble(CUSTOM_REAL)
+
+ ! xstore_outer_core, ystore_outer_core, zstore_outer_core, rmass_outer_core,
+ ! displ_outer_core, veloc_outer_core, accel_outer_core
+ static_memory_size = static_memory_size + &
+ 7.d0*nglob(IREGION_OUTER_CORE)*dble(CUSTOM_REAL)
+
+ ! vp_outer_core
+ static_memory_size = static_memory_size + &
+ dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_OUTER_CORE_STACEY*dble(CUSTOM_REAL)
+
+ ! ispec_is_tiso_outer_core
+ static_memory_size = static_memory_size + NSPEC(IREGION_OUTER_CORE)*dble(SIZE_LOGICAL)
+
+
+
+! additional arrays
+
+ ! A_array_rotation,B_array_rotation
+ static_memory_size = static_memory_size + &
+ 2.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_OUTER_CORE_ROTATION*dble(CUSTOM_REAL)
+
+ ! minus_gravity_table, &
+ ! minus_deriv_gravity_table,density_table,d_ln_density_dr_table,minus_rho_g_over_kappa_fluid
+ static_memory_size = static_memory_size + &
+ 5.d0*NRAD_GRAVITY*dble(SIZE_DOUBLE)
+
+ ! rspl,espl,espl2
+ static_memory_size = static_memory_size + &
+ 3.d0*NR*dble(SIZE_DOUBLE)
+
+ ! rmass_ocean_load
+ static_memory_size = static_memory_size + &
+ NGLOB_CRUST_MANTLE_OCEANS*dble(CUSTOM_REAL)
+
+ ! updated_dof_ocean_load
+ static_memory_size = static_memory_size + &
+ NGLOB_CRUST_MANTLE_OCEANS*dble(SIZE_LOGICAL)
+
+ ! ichunk_slice,iproc_xi_slice,iproc_eta_slice,addressing
+ static_memory_size = static_memory_size + &
+ 4.d0*NPROCTOT*dble(SIZE_INTEGER)
+
+ if( TOPOGRAPHY ) then
+ ! ibathy_topo
+ static_memory_size = static_memory_size + &
+ NX_BATHY*NY_BATHY*dble(SIZE_INTEGER)
endif
-! add arrays used to save strain for attenuation or for adjoint runs
+ if( MOVIE_VOLUME ) then
+ ! iepsilondev_.. crust_mantle
+ static_memory_size = static_memory_size + &
+ 5.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
-! epsilondev_crust_mantle
- static_memory_size = static_memory_size + 5.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_STR_OR_ATT*dble(CUSTOM_REAL)
+ ! muvstore_crust_mantle_3dmovie
+ static_memory_size = static_memory_size + &
+ dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
-! eps_trace_over_3_crust_mantle
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_STR_OR_ATT*dble(CUSTOM_REAL)
+ ! mask_ibool
+ static_memory_size = static_memory_size + &
+ nglob(IREGION_CRUST_MANTLE)*dble(SIZE_LOGICAL)
-! epsilondev_inner_core
- static_memory_size = static_memory_size + 5.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_STR_OR_ATT*dble(CUSTOM_REAL)
+ endif
-! eps_trace_over_3_inner_core
- static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_STR_OR_ATT*dble(CUSTOM_REAL)
-
! add arrays used for adjoint runs only (LQY: not very accurate)
-! b_R_memory_crust_mantle
-! b_epsilondev_crust_mantle
-! b_eps_trace_over_3_crust_mantle
-! rho_kl_crust_mantle,beta_kl_crust_mantle, alpha_kl_crust_mantle
+ ! b_R_memory_crust_mantle
+ ! b_epsilondev_crust_mantle
+ ! b_eps_trace_over_3_crust_mantle
+ ! rho_kl_crust_mantle,beta_kl_crust_mantle, alpha_kl_crust_mantle
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)
+ ! 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)
-! b_R_memory_inner_core
-! b_epsilondev_inner_core
-! b_eps_trace_over_3_inner_core
-! rho_kl_inner_core,beta_kl_inner_core, alpha_kl_inner_core
+ ! b_R_memory_inner_core
+ ! b_epsilondev_inner_core
+ ! b_eps_trace_over_3_inner_core
+ ! rho_kl_inner_core,beta_kl_inner_core, alpha_kl_inner_core
static_memory_size = static_memory_size + (5.d0*dble(N_SLS) + 9.d0)* &
dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_ADJOINT*dble(CUSTOM_REAL)
-! b_displ_crust_mantle,b_veloc_crust_mantle,b_accel_crust_mantle
- static_memory_size = static_memory_size + 3.d0*dble(NDIM)*NGLOB_CRUST_MANTLE_ADJOINT*dble(CUSTOM_REAL)
+ ! b_displ_crust_mantle,b_veloc_crust_mantle,b_accel_crust_mantle
+ static_memory_size = static_memory_size + &
+ 3.d0*dble(NDIM)*NGLOB_CRUST_MANTLE_ADJOINT*dble(CUSTOM_REAL)
-! b_displ_outer_core,b_veloc_outer_core,b_accel_outer_core
- static_memory_size = static_memory_size + 3.d0*NGLOB_OUTER_CORE_ADJOINT*dble(CUSTOM_REAL)
+ ! b_displ_outer_core,b_veloc_outer_core,b_accel_outer_core
+ static_memory_size = static_memory_size + &
+ 3.d0*NGLOB_OUTER_CORE_ADJOINT*dble(CUSTOM_REAL)
-! b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core
- static_memory_size = static_memory_size + 3.d0*dble(NDIM)*NGLOB_INNER_CORE_ADJOINT*dble(CUSTOM_REAL)
+ ! vector_accel_outer_core,vector_displ_outer_core,b_vector_displ_outer_core
+ static_memory_size = static_memory_size + &
+ 3.d0*dble(NDIM)*NGLOB_OUTER_CORE_ADJOINT*dble(CUSTOM_REAL)
-! b_A_array_rotation,b_B_array_rotation
- static_memory_size = static_memory_size + 2.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_OUTER_CORE_ROT_ADJOINT*dble(CUSTOM_REAL)
+ ! b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core
+ static_memory_size = static_memory_size + &
+ 3.d0*dble(NDIM)*NGLOB_INNER_CORE_ADJOINT*dble(CUSTOM_REAL)
+ ! b_A_array_rotation,b_B_array_rotation
+ static_memory_size = static_memory_size + &
+ 2.d0*dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_OUTER_CORE_ROT_ADJOINT*dble(CUSTOM_REAL)
+
+
end subroutine memory_eval
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/model_topo_bathy.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/model_topo_bathy.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/model_topo_bathy.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -58,7 +58,7 @@
! user output
write(IMAIN,*)
write(IMAIN,*) 'incorporating topography'
-
+
! read/save topo file on master
call read_topo_bathy_file(ibathy_topo)
call save_topo_bathy_database(ibathy_topo,LOCAL_PATH)
@@ -91,7 +91,7 @@
character(len=150) :: topo_bathy_file
integer,parameter :: TOPO_MINIMUM = - 10000 ! (depth in m )
integer,parameter :: TOPO_MAXIMUM = + 10000 ! (height in m )
-
+
call get_value_string(topo_bathy_file, 'model.topoBathy.PATHNAME_TOPO_FILE', PATHNAME_TOPO_FILE)
! reads in topography values from file
@@ -105,7 +105,7 @@
do itopo_y=1,NY_BATHY
do itopo_x=1,NX_BATHY
read(13,*,iostat=ier) val
-
+
! checks
if( ier /= 0 ) then
print*,'error read topo_bathy: ix,iy = ',itopo_x,itopo_y,val
@@ -179,7 +179,7 @@
! only master needs to save this
call create_name_database(prname,0,IREGION_CRUST_MANTLE,LOCAL_PATH)
- ! saves topography and bathymetry file for solver
+ ! saves topography and bathymetry file for solver
open(unit=27,file=prname(1:len_trim(prname))//'topo.bin', &
status='unknown',form='unformatted',action='write',iostat=ier)
if( ier /= 0 ) then
@@ -189,9 +189,9 @@
print*,'please check if path exists and rerun mesher'
call exit_mpi(0,'error opening file for database topo')
endif
-
+
write(27) ibathy_topo
- close(27)
+ close(27)
end subroutine save_topo_bathy_database
@@ -226,17 +226,17 @@
print*,'error opening file: ',prname(1:len_trim(prname))//'topo.bin'
!print*,'please check if file exists and rerun solver'
!call exit_mpi(0,'error opening file for database topo')
-
+
! read by original file
print*,'trying original topography file...'
call read_topo_bathy_file(ibathy_topo)
! saves database topo file for next time
- call save_topo_bathy_database(ibathy_topo,LOCAL_PATH)
+ call save_topo_bathy_database(ibathy_topo,LOCAL_PATH)
else
! database topo file exists
read(27) ibathy_topo
- close(27)
+ close(27)
! user output
write(IMAIN,*) " topography/bathymetry: min/max = ",minval(ibathy_topo),maxval(ibathy_topo)
@@ -263,7 +263,7 @@
! location latitude/longitude (in degree)
double precision,intent(in):: xlat,xlon
-
+
! returns elevation (in meters)
double precision,intent(out):: value
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -439,7 +439,7 @@
if(NEX_ETA < 48) &
stop 'NEX_ETA must be greater than 48 to cut the sphere into slices with positive Jacobian'
endif
-
+
! check that topology is correct if more than two chunks
if(NCHUNKS > 2 .and. NEX_XI /= NEX_ETA) &
stop 'must have NEX_XI = NEX_ETA for more than two chunks'
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_parameter_file.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -195,15 +195,15 @@
! initializes
LOCAL_TMP_PATH = LOCAL_PATH
-
+
! opens file Par_file
call open_parameter_file()
call read_value_string(LOCAL_TMP_PATH, 'LOCAL_TMP_PATH')
call read_value_clear_err()
-
+
! close parameter file
- call close_parameter_file()
+ call close_parameter_file()
end subroutine read_parameter_file
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -29,7 +29,8 @@
subroutine save_header_file(NSPEC,nglob,NEX_XI,NEX_ETA,NPROC,NPROCTOT, &
TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- ELLIPTICITY,GRAVITY,ROTATION,OCEANS,ATTENUATION,ATTENUATION_3D, &
+ ELLIPTICITY,GRAVITY,ROTATION, &
+ TOPOGRAPHY,OCEANS,ATTENUATION,ATTENUATION_3D, &
ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,NCHUNKS, &
INCLUDE_CENTRAL_CUBE,CENTER_LONGITUDE_IN_DEGREES, &
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,NSOURCES,NSTEP,&
@@ -60,31 +61,15 @@
integer NEX_XI,NEX_ETA,NPROC,NPROCTOT,NCHUNKS,NSOURCES,NSTEP
logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- ELLIPTICITY,GRAVITY,ROTATION,OCEANS,ATTENUATION,ATTENUATION_3D,INCLUDE_CENTRAL_CUBE
+ ELLIPTICITY,GRAVITY,ROTATION, &
+ TOPOGRAPHY,OCEANS,ATTENUATION,ATTENUATION_3D,INCLUDE_CENTRAL_CUBE
double precision ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES, &
CENTER_LONGITUDE_IN_DEGREES,CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH
- double precision :: subtract_central_cube_elems,subtract_central_cube_points
-
- character(len=150) HEADER_FILE
-
-! for regional code
- double precision x,y,gamma,rgt,xi,eta
- double precision x_top,y_top,z_top
- double precision ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD
-
-! rotation matrix from Euler angles
- integer i,j,ix,iy,icorner
- double precision rotation_matrix(3,3)
- double precision vector_ori(3),vector_rotated(3)
- double precision r_corner,theta_corner,phi_corner,lat,long,colat_corner
-
-! static memory size needed by the solver
+ ! static memory size needed by the solver
double precision :: static_memory_size
- integer :: att1,att2,att3,att4,att5,NCORNERSCHUNKS,NUM_FACES,NUM_MSG_TYPES
-
integer, dimension(MAX_NUM_REGIONS) :: NGLOB1D_RADIAL,NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
NSPEC2D_TOP,NSPEC2D_BOTTOM,NSPEC2DMAX_YMIN_YMAX,NSPEC2DMAX_XMIN_XMAX
integer :: NPROC_XI,NPROC_ETA
@@ -106,10 +91,26 @@
integer :: SIMULATION_TYPE
logical :: SAVE_FORWARD,MOVIE_VOLUME
+ ! local parameters
+ double precision :: subtract_central_cube_elems,subtract_central_cube_points
+ ! for regional code
+ double precision x,y,gamma,rgt,xi,eta
+ double precision x_top,y_top,z_top
+ double precision ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD
+ ! rotation matrix from Euler angles
+ integer i,j,ix,iy,icorner
+ double precision rotation_matrix(3,3)
+ double precision vector_ori(3),vector_rotated(3)
+ double precision r_corner,theta_corner,phi_corner,lat,long,colat_corner
+ integer :: att1,att2,att3,att4,att5,NCORNERSCHUNKS,NUM_FACES,NUM_MSG_TYPES
+ integer :: ier
+ character(len=150) HEADER_FILE
-! copy number of elements and points in an include file for the solver
+ ! copy number of elements and points in an include file for the solver
call get_value_string(HEADER_FILE, 'solver.HEADER_FILE', 'OUTPUT_FILES/values_from_mesher.h')
- open(unit=IOUT,file=HEADER_FILE,status='unknown')
+ open(unit=IOUT,file=HEADER_FILE,status='unknown',iostat=ier)
+ if( ier /= 0 ) stop 'error opening OUTPUT_FILES/values_from_mesher.h'
+
write(IOUT,*)
write(IOUT,*) '!'
@@ -387,6 +388,15 @@
endif
write(IOUT,*)
+ if(TOPOGRAPHY .or. OCEANS) then
+ write(IOUT,*) 'integer, parameter :: NX_BATHY_VAL = NX_BATHY'
+ write(IOUT,*) 'integer, parameter :: NY_BATHY_VAL = NY_BATHY'
+ else
+ write(IOUT,*) 'integer, parameter :: NX_BATHY_VAL = 1'
+ write(IOUT,*) 'integer, parameter :: NY_BATHY_VAL = 1'
+ endif
+ write(IOUT,*)
+
if(ROTATION) then
write(IOUT,*) 'logical, parameter :: ROTATION_VAL = .true.'
else
@@ -525,8 +535,14 @@
write(IOUT,*) 'logical, parameter :: COMPUTE_AND_STORE_STRAIN = .false.'
endif
+ if (MOVIE_VOLUME) then
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_3DMOVIE = NSPEC_CRUST_MANTLE'
+ write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_3DMOVIE = NGLOB_CRUST_MANTLE'
+ else
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_3DMOVIE = 1'
+ write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_3DMOVIE = 1'
+ endif
-
close(IOUT)
end subroutine save_header_file
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/write_VTK_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/write_VTK_file.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/write_VTK_file.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -431,7 +431,7 @@
integer :: myrank,NPROCTOT
integer ::nspec,nglob
-
+
integer, dimension(nspec):: idoubling
! global coordinates
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_boundary_kernel.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_boundary_kernel.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_boundary_kernel.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -381,7 +381,7 @@
real(kind=CUSTOM_REAL) :: b_kl(NGLLX,NGLLY,NSPEC2D_DISC)
logical :: fluid_solid_boundary
-! --- local variables ---
+ ! --- local variables ---
integer ispec2D,i,j,k,iglob,ispec
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: displl, accell, b_displl, Kdvect
real(kind=CUSTOM_REAL), dimension(NDIM) :: normal, temp1, temp2, temp3
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -1307,7 +1307,7 @@
integer :: i_SLS1,i_SLS2
#endif
-#ifdef __HANDOPT_ATT
+#ifdef _HANDOPT_ATT
! way 2:
! note: this should help compilers to pipeline the code and make better use of the cache;
! depending on compilers, it can further decrease the computation time by ~ 30%.
@@ -1332,11 +1332,11 @@
do i_SLS = imodulo_N_SLS+1,N_SLS,3
R_xx_val1 = R_xx_loc(i_SLS)
R_yy_val1 = R_yy_loc(i_SLS)
-
+
i_SLS1=i_SLS+1
R_xx_val2 = R_xx_loc(i_SLS1)
R_yy_val2 = R_yy_loc(i_SLS1)
-
+
i_SLS2 =i_SLS+2
R_xx_val3 = R_xx_loc(i_SLS2)
R_yy_val3 = R_yy_loc(i_SLS2)
@@ -1346,7 +1346,7 @@
sigma_zz = sigma_zz + R_xx_val1 + R_yy_val1 &
+ R_xx_val2 + R_yy_val2 &
+ R_xx_val3 + R_yy_val3
-
+
sigma_xy = sigma_xy - R_xy_loc(i_SLS)
sigma_xz = sigma_xz - R_xz_loc(i_SLS)
sigma_yz = sigma_yz - R_yz_loc(i_SLS)
@@ -1437,7 +1437,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX, NGLLY, NGLLZ) :: factor_common_c44_muv
integer :: i_SLS
-#ifdef __HANDOPT_ATT
+#ifdef _HANDOPT_ATT
real(kind=CUSTOM_REAL) :: alphal,betal,gammal
integer :: i,j,k
#endif
@@ -1448,7 +1448,7 @@
! IMPROVE we use mu_v here even if there is some anisotropy
! IMPROVE we should probably use an average value instead
-#ifdef __HANDOPT_ATT
+#ifdef _HANDOPT_ATT
! way 2:
do i_SLS = 1,N_SLS
@@ -1604,7 +1604,7 @@
integer :: i_SLS
-#ifdef __HANDOPT_ATT
+#ifdef _HANDOPT_ATT
real(kind=CUSTOM_REAL) :: alphal,betal,gammal
integer :: i,j,k
#endif
@@ -1615,7 +1615,7 @@
! IMPROVE we use mu_v here even if there is some anisotropy
! IMPROVE we should probably use an average value instead
-#ifdef __HANDOPT_ATT
+#ifdef _HANDOPT_ATT
! way 2:
do i_SLS = 1,N_SLS
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -246,7 +246,7 @@
tempz2l_att = 0._CUSTOM_REAL
tempz3l_att = 0._CUSTOM_REAL
- ! use first order Taylor expansion of displacement for local storage of stresses
+ ! use first order Taylor expansion of displacement for local storage of stresses
! at this current time step, to fix attenuation in a consistent way
do l=1,NGLLX
hp1 = hprime_xx(i,l)
@@ -352,7 +352,7 @@
epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl_att
endif
- else
+ else
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -76,7 +76,7 @@
phase_ispec_inner => phase_ispec_inner_crust_mantle, &
nspec_outer => nspec_outer_crust_mantle, &
nspec_inner => nspec_inner_crust_mantle
-
+
implicit none
integer :: NSPEC,NGLOB,NSPEC_ATT
@@ -247,9 +247,9 @@
! way 1:
do i=1,NGLLX
iglob1 = ibool(i,j,k,ispec)
- dummyx_loc(i,j,k) = displ_crust_mantle(1,iglob1)
- dummyy_loc(i,j,k) = displ_crust_mantle(2,iglob1)
- dummyz_loc(i,j,k) = displ_crust_mantle(3,iglob1)
+ dummyx_loc(i,j,k) = displ_crust_mantle(1,iglob1)
+ dummyy_loc(i,j,k) = displ_crust_mantle(2,iglob1)
+ dummyz_loc(i,j,k) = displ_crust_mantle(3,iglob1)
enddo
#endif
@@ -257,7 +257,7 @@
enddo
if( ATTENUATION_VAL ) then
- ! use first order Taylor expansion of displacement for local storage of stresses
+ ! use first order Taylor expansion of displacement for local storage of stresses
! at this current time step, to fix attenuation in a consistent way
do k=1,NGLLZ
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_elastic.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_elastic.F90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_elastic.F90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -242,14 +242,14 @@
! add sources for backward/reconstructed wavefield
if (SIMULATION_TYPE == 3 .and. nsources_local > 0) &
call compute_add_sources_backward()
-
+
case( 1 )
! the first step of noise tomography is to use |S(\omega)|^2 as a point force source at one of the receivers.
! hence, instead of a moment tensor 'sourcearrays', a 'noise_sourcearray' for a point force is needed.
! furthermore, the CMTSOLUTION needs to be zero, i.e., no earthquakes.
! now this must be manually set in DATA/CMTSOLUTION, by USERS.
call noise_add_source_master_rec()
-
+
case( 2 )
! second step of noise tomography, i.e., read the surface movie saved at every timestep
! use the movie to drive the ensemble forward wavefield
@@ -259,14 +259,14 @@
! note the ensemble forward sources are generally distributed on the surface of the earth
! that's to say, the ensemble forward source is kind of a surface force density, not a body force density
! therefore, we must add it here, before applying the inverse of mass matrix
-
+
case( 3 )
! third step of noise tomography, i.e., read the surface movie saved at every timestep
! use the movie to reconstruct the ensemble forward wavefield
! the ensemble adjoint wavefield is done as usual
! note instead of "NSTEP-it+1", now we us "it", since reconstruction is a reversal of reversal
call noise_read_add_surface_movie(NGLOB_CRUST_MANTLE_ADJOINT,b_accel_crust_mantle,it)
-
+
end select
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -59,7 +59,7 @@
phase_ispec_inner => phase_ispec_inner_inner_core, &
nspec_outer => nspec_outer_inner_core, &
nspec_inner => nspec_inner_inner_core
-
+
implicit none
integer :: NSPEC,NGLOB,NSPEC_ATT
@@ -227,7 +227,7 @@
tempz2l_att = 0._CUSTOM_REAL
tempz3l_att = 0._CUSTOM_REAL
- ! use first order Taylor expansion of displacement for local storage of stresses
+ ! use first order Taylor expansion of displacement for local storage of stresses
! at this current time step, to fix attenuation in a consistent way
do l=1,NGLLX
hp1 = hprime_xx(i,l)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_inner_core_Dev.F90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -275,9 +275,9 @@
#endif
enddo
enddo
-
+
if( ATTENUATION_VAL ) then
- ! use first order Taylor expansion of displacement for local storage of stresses
+ ! use first order Taylor expansion of displacement for local storage of stresses
! at this current time step, to fix attenuation in a consistent way
do k=1,NGLLZ
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -37,7 +37,7 @@
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
minus_rho_g_over_kappa_fluid,d_ln_density_dr_table, &
- MOVIE_VOLUME
+ MOVIE_VOLUME
use specfem_par_outercore,only: &
xstore => xstore_outer_core,ystore => ystore_outer_core,zstore => zstore_outer_core, &
@@ -48,7 +48,7 @@
phase_ispec_inner => phase_ispec_inner_outer_core, &
nspec_outer => nspec_outer_outer_core, &
nspec_inner => nspec_inner_outer_core
-
+
implicit none
integer :: NSPEC,NGLOB
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -49,7 +49,7 @@
phase_ispec_inner => phase_ispec_inner_outer_core, &
nspec_outer => nspec_outer_outer_core, &
nspec_inner => nspec_inner_outer_core
-
+
implicit none
integer :: NSPEC,NGLOB
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_kernels.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_kernels.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -207,6 +207,7 @@
real(kind=CUSTOM_REAL) :: tempz1l,tempz2l,tempz3l
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
integer :: i,j,k,l,ispec,iglob
@@ -260,13 +261,13 @@
do l=1,NGLLZ
tempx3l = tempx3l + accel_outer_core(ibool_outer_core(i,j,l,ispec)) * hprime_zz(k,l)
enddo
- vector_accel_outer_core(1,iglob) = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
- vector_accel_outer_core(2,iglob) = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
- vector_accel_outer_core(3,iglob) = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+ vector_accel(1) = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ vector_accel(2) = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ vector_accel(3) = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
! density kernel
rho_kl_outer_core(i,j,k,ispec) = rho_kl_outer_core(i,j,k,ispec) &
- + deltat * dot_product(vector_accel_outer_core(:,iglob), b_vector_displ_outer_core(:,iglob))
+ + deltat * dot_product(vector_accel(:), b_vector_displ_outer_core(:,iglob))
! bulk modulus kernel
kappal = rhostore_outer_core(i,j,k,ispec)/kappavstore_outer_core(i,j,k,ispec)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -147,7 +147,7 @@
request_send_scalar_outer_core,request_recv_scalar_outer_core)
deallocate(buffer_send_vector_inner_core,buffer_recv_vector_inner_core, &
request_send_vector_inner_core,request_recv_vector_inner_core)
-
+
if( SIMULATION_TYPE == 3 ) then
deallocate(b_buffer_send_vector_crust_mantle,b_buffer_recv_vector_crust_mantle, &
b_request_send_vector_crust_mantle,b_request_recv_vector_crust_mantle)
@@ -173,7 +173,7 @@
deallocate(num_elem_colors_crust_mantle)
deallocate(num_elem_colors_outer_core)
deallocate(num_elem_colors_inner_core)
-
+
! sources
deallocate(islice_selected_source, &
ispec_selected_source, &
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/iterate_time.F90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -67,12 +67,12 @@
imodulo_NGLOB_OUTER_CORE = mod(NGLOB_OUTER_CORE,3)
#endif
-! get MPI starting time
+ ! get MPI starting time
time_start = MPI_WTIME()
-! *********************************************************
-! ************* MAIN LOOP OVER THE TIME STEPS *************
-! *********************************************************
+ ! *********************************************************
+ ! ************* MAIN LOOP OVER THE TIME STEPS *************
+ ! *********************************************************
do it = it_begin,it_end
@@ -120,9 +120,9 @@
enddo ! end of main time loop
-!
-!---- end of time iteration loop
-!
+ !
+ !---- end of time iteration loop
+ !
call it_print_elapsed_time()
@@ -131,7 +131,9 @@
end subroutine iterate_time
-!=====================================================================
+!
+!-------------------------------------------------------------------------------------------------
+!
subroutine it_update_displacement_scheme()
@@ -168,7 +170,6 @@
use specfem_par_crustmantle
use specfem_par_innercore
use specfem_par_outercore
- use specfem_par_movie
implicit none
! local parameters
@@ -418,61 +419,38 @@
b_deltat, b_deltatsqover2, b_deltatover2)
endif
- ! integral of strain for adjoint movie volume
- if(MOVIE_VOLUME .and. (MOVIE_VOLUME_TYPE == 2 .or. MOVIE_VOLUME_TYPE == 3) ) then
- if( GPU_MODE ) then
- ! transfers strain arrays onto CPU
- call transfer_strain_cm_from_device(Mesh_pointer,eps_trace_over_3_crust_mantle, &
- epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle, &
- epsilondev_xy_crust_mantle,epsilondev_xz_crust_mantle, &
- epsilondev_yz_crust_mantle)
- endif
-
- ! updates integral values
- Iepsilondev_xx_crust_mantle(:,:,:,:) = Iepsilondev_xx_crust_mantle(:,:,:,:) &
- + deltat*epsilondev_xx_crust_mantle(:,:,:,:)
- Iepsilondev_yy_crust_mantle(:,:,:,:) = Iepsilondev_yy_crust_mantle(:,:,:,:) &
- + deltat*epsilondev_yy_crust_mantle(:,:,:,:)
- Iepsilondev_xy_crust_mantle(:,:,:,:) = Iepsilondev_xy_crust_mantle(:,:,:,:) &
- + deltat*epsilondev_xy_crust_mantle(:,:,:,:)
- Iepsilondev_xz_crust_mantle(:,:,:,:) = Iepsilondev_xz_crust_mantle(:,:,:,:) &
- + deltat*epsilondev_xz_crust_mantle(:,:,:,:)
- Iepsilondev_yz_crust_mantle(:,:,:,:) = Iepsilondev_yz_crust_mantle(:,:,:,:) &
- + deltat*epsilondev_yz_crust_mantle(:,:,:,:)
-
- Ieps_trace_over_3_crust_mantle(:,:,:,:) = Ieps_trace_over_3_crust_mantle(:,:,:,:) &
- + deltat*eps_trace_over_3_crust_mantle(:,:,:,:)
- endif
-
-
-
end subroutine it_update_displacement_scheme
-!=====================================================================
+!
+!-------------------------------------------------------------------------------------------------
+!
subroutine it_print_elapsed_time()
- use specfem_par
- implicit none
- include 'mpif.h'
+ use specfem_par
+ implicit none
- ! local parameters
- integer :: ihours,iminutes,iseconds,int_tCPU
+ include 'mpif.h'
- if(myrank == 0) then
- ! elapsed time since beginning of the simulation
- tCPU = MPI_WTIME() - time_start
- int_tCPU = int(tCPU)
- ihours = int_tCPU / 3600
- iminutes = (int_tCPU - 3600*ihours) / 60
- iseconds = int_tCPU - 3600*ihours - 60*iminutes
- write(IMAIN,*) 'Time-Loop Complete. Timing info:'
- write(IMAIN,*) 'Total elapsed time in seconds = ',tCPU
- write(IMAIN,"(' Total elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
- endif
+ ! local parameters
+ integer :: ihours,iminutes,iseconds,int_tCPU
+
+ if(myrank == 0) then
+ ! elapsed time since beginning of the simulation
+ tCPU = MPI_WTIME() - time_start
+ int_tCPU = int(tCPU)
+ ihours = int_tCPU / 3600
+ iminutes = (int_tCPU - 3600*ihours) / 60
+ iseconds = int_tCPU - 3600*ihours - 60*iminutes
+ write(IMAIN,*) 'Time-Loop Complete. Timing info:'
+ write(IMAIN,*) 'Total elapsed time in seconds = ',tCPU
+ write(IMAIN,"(' Total elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
+ endif
end subroutine it_print_elapsed_time
-!=====================================================================
+!
+!-------------------------------------------------------------------------------------------------
+!
subroutine it_check_stability()
@@ -519,7 +497,9 @@
end subroutine it_check_stability
-!=====================================================================
+!
+!-------------------------------------------------------------------------------------------------
+!
subroutine it_transfer_from_GPU()
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -114,7 +114,7 @@
double precision, dimension(:,:), allocatable :: final_distance_all
double precision distmin,final_distance_max
-
+
! receiver information
! timing information for the stations
! station information for writing the seismograms
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -66,7 +66,6 @@
integer i,j,k,ispec,iglob
integer ier
- double precision t0, hdur_gaussian(NSOURCES)
double precision ell
double precision elevation
@@ -109,11 +108,10 @@
double precision, dimension(:), allocatable :: xi_source_subset,eta_source_subset,gamma_source_subset
double precision, dimension(NSOURCES) :: lat,long,depth
- double precision scalar_moment
double precision moment_tensor(6,NSOURCES)
double precision radius
- character(len=150) OUTPUT_FILES,plot_file
+ character(len=150) OUTPUT_FILES
double precision, dimension(:), allocatable :: x_found_source,y_found_source,z_found_source
double precision r_found_source
@@ -127,20 +125,10 @@
logical located_target
- ! for calculation of source time function and spectrum
- integer it,iom
- double precision time_source,om
- double precision, external :: comp_source_time_function,comp_source_spectrum
- double precision, external :: comp_source_time_function_rickr
-
- ! number of points to plot the source time function and spectrum
- integer, parameter :: NSAMP_PLOT_SOURCE = 1000
-
integer iorientation
double precision stazi,stdip,thetan,phin,n(3)
integer imin,imax,jmin,jmax,kmin,kmax
double precision :: f0,t0_ricker
- double precision t_cmt_used(NSOURCES)
! mask source region (mask values are between 0 and 1, with 0 around sources)
real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable :: mask_source
@@ -373,7 +361,7 @@
kmin = 2
kmax = NGLLZ - 1
endif
-
+
do k = kmin,kmax
do j = jmin,jmax
do i = imin,imax
@@ -500,7 +488,7 @@
dx = - (x - x_target_source)
dy = - (y - y_target_source)
dz = - (z - z_target_source)
-
+
! compute increments
dxi = xix*dx + xiy*dy + xiz*dz
deta = etax*dx + etay*dy + etaz*dz
@@ -577,7 +565,8 @@
if(myrank == 0) then
! check that the gather operation went well
- if(minval(ispec_selected_source_all) <= 0) call exit_MPI(myrank,'gather operation failed for source')
+ if(minval(ispec_selected_source_all) <= 0) &
+ call exit_MPI(myrank,'gather operation failed for source')
! loop on all the sources within subsets
do isource_in_this_subset = 1,NSOURCES_SUBSET_current_size
@@ -663,7 +652,7 @@
! convert geocentric to geographic colatitude
colat_source = PI/2.0d0 &
- - datan(1.006760466d0*dcos(theta_source(isource))/dmax1(TINYVAL,dsin(theta_source(isource))))
+ - datan(1.006760466d0*dcos(theta_source(isource))/dmax1(TINYVAL,dsin(theta_source(isource))))
if(phi_source(isource)>PI) phi_source(isource)=phi_source(isource)-TWO_PI
write(IMAIN,*)
@@ -697,90 +686,9 @@
endif
! print source time function and spectrum
- if(PRINT_SOURCE_TIME_FUNCTION) then
+ if(PRINT_SOURCE_TIME_FUNCTION) call print_stf(NSOURCES,isource,moment_tensor,tshift_cmt,hdur, &
+ min_tshift_cmt_original,NSTEP,DT,OUTPUT_FILES)
- write(IMAIN,*)
- write(IMAIN,*) 'printing the source-time function'
-
- ! print the source-time function
- if(NSOURCES == 1) then
- plot_file = '/plot_source_time_function.txt'
- else
- if(isource < 10) then
- write(plot_file,"('/plot_source_time_function',i1,'.txt')") isource
- elseif(isource < 100) then
- write(plot_file,"('/plot_source_time_function',i2,'.txt')") isource
- else
- write(plot_file,"('/plot_source_time_function',i3,'.txt')") isource
- endif
- endif
- open(unit=27,file=trim(OUTPUT_FILES)//plot_file,status='unknown')
-
- scalar_moment = 0.
- do i = 1,6
- scalar_moment = scalar_moment + moment_tensor(i,isource)**2
- enddo
- scalar_moment = dsqrt(scalar_moment/2.)
-
- ! define t0 as the earliest start time
- ! note: this calculation here is only used for outputting the plot_source_time_function file
- ! (see setup_sources_receivers.f90)
- t0 = - 1.5d0*minval( tshift_cmt(:) - hdur(:) )
- if( USE_FORCE_POINT_SOURCE ) t0 = - 1.2d0 * minval(tshift_cmt(:) - 1.0d0/hdur(:))
- t_cmt_used(:) = t_cmt_used(:)
- if( USER_T0 > 0.d0 ) then
- if( t0 <= USER_T0 + min_tshift_cmt_original ) then
- t_cmt_used(:) = tshift_cmt(:) + min_tshift_cmt_original
- t0 = USER_T0
- endif
- endif
- ! convert the half duration for triangle STF to the one for gaussian STF
- ! note: this calculation here is only used for outputting the plot_source_time_function file
- ! (see setup_sources_receivers.f90)
- hdur_gaussian(:) = hdur(:)/SOURCE_DECAY_MIMIC_TRIANGLE
-
- ! writes out source time function to file
- do it=1,NSTEP
- time_source = dble(it-1)*DT-t0-t_cmt_used(isource)
- if( USE_FORCE_POINT_SOURCE ) then
- ! Ricker source time function
- f0 = hdur(isource)
- write(27,*) sngl(dble(it-1)*DT-t0), &
- sngl(FACTOR_FORCE_SOURCE*comp_source_time_function_rickr(time_source,f0))
- else
- ! Gaussian source time function
- write(27,*) sngl(dble(it-1)*DT-t0), &
- sngl(scalar_moment*comp_source_time_function(time_source,hdur_gaussian(isource)))
- endif
- enddo
- close(27)
-
- write(IMAIN,*)
- write(IMAIN,*) 'printing the source spectrum'
-
- ! print the spectrum of the derivative of the source from 0 to 1/8 Hz
- if(NSOURCES == 1) then
- plot_file = '/plot_source_spectrum.txt'
- else
- if(isource < 10) then
- write(plot_file,"('/plot_source_spectrum',i1,'.txt')") isource
- elseif(isource < 100) then
- write(plot_file,"('/plot_source_spectrum',i2,'.txt')") isource
- else
- write(plot_file,"('/plot_source_spectrum',i3,'.txt')") isource
- endif
- endif
- open(unit=27,file=trim(OUTPUT_FILES)//plot_file,status='unknown')
-
- do iom=1,NSAMP_PLOT_SOURCE
- om=TWO_PI*(1.0d0/8.0d0)*(iom-1)/dble(NSAMP_PLOT_SOURCE-1)
- write(27,*) sngl(om/TWO_PI), &
- sngl(scalar_moment*om*comp_source_spectrum(om,hdur(isource)))
- enddo
- close(27)
-
- endif !PRINT_SOURCE_TIME_FUNCTION
-
enddo ! end of loop on all the sources within current source subset
endif ! end of section executed by main process only
@@ -803,7 +711,6 @@
write(IMAIN,*)
endif
-
! main process broadcasts the results to all the slices
call MPI_BCAST(islice_selected_source,NSOURCES,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
call MPI_BCAST(ispec_selected_source,NSOURCES,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
@@ -827,7 +734,7 @@
write(IMAIN,*)
endif
call sync_all()
-
+
end subroutine locate_sources
!
@@ -913,8 +820,139 @@
open(unit=27,file=trim(prname)//'mask_source.bin', &
status='unknown',form='unformatted',action='write',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening mask_source.bin file')
-
+
write(27) mask_source
close(27)
end subroutine save_mask_source
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine print_stf(NSOURCES,isource,moment_tensor,tshift_cmt,hdur, &
+ min_tshift_cmt_original,NSTEP,DT,OUTPUT_FILES)
+
+! prints source time function
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: NSOURCES,isource
+
+ double precision :: moment_tensor(6,NSOURCES)
+ double precision :: tshift_cmt(NSOURCES)
+ double precision :: hdur(NSOURCES)
+
+ double precision :: min_tshift_cmt_original
+ integer :: NSTEP
+ double precision :: DT
+
+ character(len=150) OUTPUT_FILES
+
+ ! local parameters
+ integer :: i,it,iom,ier
+ double precision :: scalar_moment
+ double precision :: t0, hdur_gaussian(NSOURCES)
+ double precision :: t_cmt_used(NSOURCES)
+ double precision time_source,om
+ double precision :: f0
+
+ double precision, external :: comp_source_time_function,comp_source_spectrum
+ double precision, external :: comp_source_time_function_rickr
+
+ character(len=150) plot_file
+
+ ! number of points to plot the source time function and spectrum
+ integer, parameter :: NSAMP_PLOT_SOURCE = 1000
+
+ ! user output
+ write(IMAIN,*)
+ write(IMAIN,*) 'printing the source-time function'
+
+ ! print the source-time function
+ if(NSOURCES == 1) then
+ plot_file = '/plot_source_time_function.txt'
+ else
+ if(isource < 10) then
+ write(plot_file,"('/plot_source_time_function',i1,'.txt')") isource
+ elseif(isource < 100) then
+ write(plot_file,"('/plot_source_time_function',i2,'.txt')") isource
+ else
+ write(plot_file,"('/plot_source_time_function',i3,'.txt')") isource
+ endif
+ endif
+
+ open(unit=27,file=trim(OUTPUT_FILES)//plot_file, &
+ status='unknown',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(0,'error opening plot_source_time_function file')
+
+ scalar_moment = 0.
+ do i = 1,6
+ scalar_moment = scalar_moment + moment_tensor(i,isource)**2
+ enddo
+ scalar_moment = dsqrt(scalar_moment/2.0d0)
+
+ ! define t0 as the earliest start time
+ ! note: this calculation here is only used for outputting the plot_source_time_function file
+ ! (see setup_sources_receivers.f90)
+ t0 = - 1.5d0*minval( tshift_cmt(:) - hdur(:) )
+ if( USE_FORCE_POINT_SOURCE ) t0 = - 1.2d0 * minval(tshift_cmt(:) - 1.0d0/hdur(:))
+ t_cmt_used(:) = t_cmt_used(:)
+ if( USER_T0 > 0.d0 ) then
+ if( t0 <= USER_T0 + min_tshift_cmt_original ) then
+ t_cmt_used(:) = tshift_cmt(:) + min_tshift_cmt_original
+ t0 = USER_T0
+ endif
+ endif
+ ! convert the half duration for triangle STF to the one for gaussian STF
+ ! note: this calculation here is only used for outputting the plot_source_time_function file
+ ! (see setup_sources_receivers.f90)
+ hdur_gaussian(:) = hdur(:)/SOURCE_DECAY_MIMIC_TRIANGLE
+
+ ! writes out source time function to file
+ do it=1,NSTEP
+ time_source = dble(it-1)*DT-t0-t_cmt_used(isource)
+ if( USE_FORCE_POINT_SOURCE ) then
+ ! Ricker source time function
+ f0 = hdur(isource)
+ write(27,*) sngl(dble(it-1)*DT-t0), &
+ sngl(FACTOR_FORCE_SOURCE*comp_source_time_function_rickr(time_source,f0))
+ else
+ ! Gaussian source time function
+ write(27,*) sngl(dble(it-1)*DT-t0), &
+ sngl(scalar_moment*comp_source_time_function(time_source,hdur_gaussian(isource)))
+ endif
+ enddo
+ close(27)
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'printing the source spectrum'
+
+ ! print the spectrum of the derivative of the source from 0 to 1/8 Hz
+ if(NSOURCES == 1) then
+ plot_file = '/plot_source_spectrum.txt'
+ else
+ if(isource < 10) then
+ write(plot_file,"('/plot_source_spectrum',i1,'.txt')") isource
+ elseif(isource < 100) then
+ write(plot_file,"('/plot_source_spectrum',i2,'.txt')") isource
+ else
+ write(plot_file,"('/plot_source_spectrum',i3,'.txt')") isource
+ endif
+ endif
+
+ open(unit=27,file=trim(OUTPUT_FILES)//plot_file, &
+ status='unknown',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(0,'error opening plot_source_spectrum file')
+
+ do iom=1,NSAMP_PLOT_SOURCE
+ om=TWO_PI*(1.0d0/8.0d0)*(iom-1)/dble(NSAMP_PLOT_SOURCE-1)
+ write(27,*) sngl(om/TWO_PI), &
+ sngl(scalar_moment*om*comp_source_spectrum(om,hdur(isource)))
+ enddo
+ close(27)
+
+ end subroutine print_stf
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -405,18 +405,12 @@
if (NSPEC_CRUST_MANTLE_STRAIN_ONLY /= NSPEC_CRUST_MANTLE) &
stop 'NSPEC_CRUST_MANTLE_STRAIN_ONLY /= NSPEC_CRUST_MANTLE'
- call count_points_movie_volume()
-! LOCAL_TMP_PATH,ibool_crust_mantle, xstore_crust_mantle,ystore_crust_mantle, &
-! zstore_crust_mantle,MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH, &
-! MOVIE_COARSE,npoints_3dmovie,nspecel_3dmovie,num_ibool_3dmovie,mask_ibool,mask_3dmovie)
+ call movie_volume_count_points()
allocate(nu_3dmovie(3,3,npoints_3dmovie),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating nu for 3d movie')
call write_movie_volume_mesh()
-! npoints_3dmovie,LOCAL_TMP_PATH,ibool_crust_mantle,xstore_crust_mantle, &
-! ystore_crust_mantle,zstore_crust_mantle, muvstore_crust_mantle_3dmovie, &
-! mask_3dmovie,mask_ibool,num_ibool_3dmovie,nu_3dmovie,MOVIE_COARSE)
if(myrank == 0) then
write(IMAIN,*)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_arrays_solver.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_arrays_solver.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -98,7 +98,7 @@
logical,dimension(nspec) :: dummy_l
! processor identification
character(len=150) prname
-
+
! create the name for the database of the current slide and region
call create_name_database(prname,myrank,iregion_code,LOCAL_PATH)
@@ -185,7 +185,7 @@
open(unit=IIN,file=prname(1:len_trim(prname))//'solver_data_2.bin', &
status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening file solver_data_2.bin')
-
+
read(IIN) xstore
read(IIN) ystore
read(IIN) zstore
@@ -195,7 +195,7 @@
read(IIN) idoubling
read(IIN) dummy_l ! is_on_a_slice_edge
-
+
read(IIN) ispec_is_tiso
close(IIN)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -92,10 +92,14 @@
! local parameters
integer :: nspec_iso,nspec_tiso,nspec_ani
logical :: READ_KAPPA_MU,READ_TISO
-
! dummy array that does not need to be actually read
- integer, dimension(NSPEC_CRUST_MANTLE) :: dummy_i
+ integer, dimension(:),allocatable :: dummy_i
+ integer :: ier
+ ! allocates dummy array
+ allocate(dummy_i(NSPEC_CRUST_MANTLE),stat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error allocating dummy_i array in read_mesh_databases_CM')
+
! crust and mantle
if(ANISOTROPIC_3D_MANTLE_VAL) then
READ_KAPPA_MU = .false.
@@ -147,6 +151,7 @@
maxval(ibool_crust_mantle(:,:,:,:)) /= NGLOB_CRUST_MANTLE) &
call exit_MPI(myrank,'incorrect global numbering: iboolmax does not equal nglob in crust and mantle')
+ deallocate(dummy_i)
end subroutine read_mesh_databases_CM
@@ -301,7 +306,7 @@
! local parameters
integer :: njunk1,njunk2,njunk3
integer :: ier
-
+
! crust and mantle
! create name of database
call create_name_database(prname,myrank,IREGION_CRUST_MANTLE,LOCAL_PATH)
@@ -310,7 +315,7 @@
open(unit=27,file=prname(1:len_trim(prname))//'boundary.bin', &
status='old',form='unformatted',action='read',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening crust_mantle boundary.bin file')
-
+
read(27) nspec2D_xmin_crust_mantle
read(27) nspec2D_xmax_crust_mantle
read(27) nspec2D_ymin_crust_mantle
@@ -355,7 +360,7 @@
open(unit=27,file=prname(1:len_trim(prname))//'boundary.bin', &
status='old',form='unformatted',action='read',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening outer_core boundary.bin file')
-
+
read(27) nspec2D_xmin_outer_core
read(27) nspec2D_xmax_outer_core
read(27) nspec2D_ymin_outer_core
@@ -399,7 +404,7 @@
open(unit=27,file=prname(1:len_trim(prname))//'boundary.bin', &
status='old',form='unformatted',action='read',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening inner_core boundary.bin file')
-
+
read(27) nspec2D_xmin_inner_core
read(27) nspec2D_xmax_inner_core
read(27) nspec2D_ymin_inner_core
@@ -425,7 +430,7 @@
open(unit=27,file=prname(1:len_trim(prname))//'boundary_disc.bin', &
status='old',form='unformatted',action='read',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening boundary_disc.bin file')
-
+
read(27) njunk1,njunk2,njunk3
if (njunk1 /= NSPEC2D_MOHO .and. njunk2 /= NSPEC2D_400 .and. njunk3 /= NSPEC2D_670) &
call exit_mpi(myrank, 'Error reading ibelm_disc.bin file')
@@ -472,12 +477,12 @@
if(myrank == 0) then
open(unit=IIN,file=trim(OUTPUT_FILES)//'/addressing.txt',status='old',action='read',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening addressing.txt')
-
+
do iproc = 0,NPROCTOT_VAL-1
read(IIN,*) iproc_read,ichunk,iproc_xi,iproc_eta
-
+
if(iproc_read /= iproc) call exit_MPI(myrank,'incorrect slice number read')
-
+
addressing(ichunk,iproc_xi,iproc_eta) = iproc
ichunk_slice(iproc) = ichunk
iproc_xi_slice(iproc) = iproc_xi
@@ -557,11 +562,11 @@
!
subroutine read_mesh_databases_MPI()
-
+
use specfem_par
use specfem_par_crustmantle
use specfem_par_innercore
- use specfem_par_outercore
+ use specfem_par_outercore
implicit none
! local parameters
@@ -590,7 +595,7 @@
! outer core
call read_mesh_databases_MPI_OC()
-
+
allocate(buffer_send_scalar_outer_core(max_nibool_interfaces_outer_core,num_interfaces_outer_core), &
buffer_recv_scalar_outer_core(max_nibool_interfaces_outer_core,num_interfaces_outer_core), &
request_send_scalar_outer_core(num_interfaces_outer_core), &
@@ -647,15 +652,15 @@
endif
! synchronizes MPI processes
call sync_all()
-
+
end subroutine read_mesh_databases_MPI
-
+
!
!-------------------------------------------------------------------------------------------------
!
subroutine read_mesh_databases_MPI_CM()
-
+
use specfem_par
use specfem_par_crustmantle
implicit none
@@ -664,10 +669,10 @@
integer :: ier
! crust mantle region
-
+
! create the name for the database of the current slide and region
call create_name_database(prname,myrank,IREGION_CRUST_MANTLE,LOCAL_PATH)
-
+
open(unit=IIN,file=prname(1:len_trim(prname))//'solver_data_mpi.bin', &
status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening solver_data_mpi.bin')
@@ -679,21 +684,21 @@
stat=ier)
if( ier /= 0 ) &
call exit_mpi(myrank,'error allocating array my_neighbours_crust_mantle etc.')
-
+
if( num_interfaces_crust_mantle > 0 ) then
read(IIN) max_nibool_interfaces_crust_mantle
allocate(ibool_interfaces_crust_mantle(max_nibool_interfaces_crust_mantle,num_interfaces_crust_mantle), &
stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating array ibool_interfaces_crust_mantle')
-
+
read(IIN) my_neighbours_crust_mantle
read(IIN) nibool_interfaces_crust_mantle
read(IIN) ibool_interfaces_crust_mantle
else
! dummy array
max_nibool_interfaces_crust_mantle = 0
- allocate(ibool_interfaces_crust_mantle(0,0),stat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error allocating array dummy ibool_interfaces_crust_mantle')
+ allocate(ibool_interfaces_crust_mantle(0,0),stat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error allocating array dummy ibool_interfaces_crust_mantle')
endif
! inner / outer elements
@@ -706,7 +711,7 @@
stat=ier)
if( ier /= 0 ) &
call exit_mpi(myrank,'error allocating array phase_ispec_inner_crust_mantle')
-
+
if(num_phase_ispec_crust_mantle > 0 ) read(IIN) phase_ispec_inner_crust_mantle
! mesh coloring for GPUs
@@ -729,9 +734,9 @@
if( ier /= 0 ) &
call exit_mpi(myrank,'error allocating num_elem_colors_crust_mantle array')
endif
-
+
close(IIN)
-
+
end subroutine read_mesh_databases_MPI_CM
@@ -740,19 +745,19 @@
!
subroutine read_mesh_databases_MPI_OC()
-
+
use specfem_par
- use specfem_par_outercore
+ use specfem_par_outercore
implicit none
! local parameters
integer :: ier
! crust mantle region
-
+
! create the name for the database of the current slide and region
call create_name_database(prname,myrank,IREGION_OUTER_CORE,LOCAL_PATH)
-
+
open(unit=IIN,file=prname(1:len_trim(prname))//'solver_data_mpi.bin', &
status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening solver_data_mpi.bin')
@@ -764,21 +769,21 @@
stat=ier)
if( ier /= 0 ) &
call exit_mpi(myrank,'error allocating array my_neighbours_outer_core etc.')
-
+
if( num_interfaces_outer_core > 0 ) then
read(IIN) max_nibool_interfaces_outer_core
allocate(ibool_interfaces_outer_core(max_nibool_interfaces_outer_core,num_interfaces_outer_core), &
stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating array ibool_interfaces_outer_core')
-
+
read(IIN) my_neighbours_outer_core
read(IIN) nibool_interfaces_outer_core
read(IIN) ibool_interfaces_outer_core
else
! dummy array
max_nibool_interfaces_outer_core = 0
- allocate(ibool_interfaces_outer_core(0,0),stat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error allocating array dummy ibool_interfaces_outer_core')
+ allocate(ibool_interfaces_outer_core(0,0),stat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error allocating array dummy ibool_interfaces_outer_core')
endif
! inner / outer elements
@@ -791,7 +796,7 @@
stat=ier)
if( ier /= 0 ) &
call exit_mpi(myrank,'error allocating array phase_ispec_inner_outer_core')
-
+
if(num_phase_ispec_outer_core > 0 ) read(IIN) phase_ispec_inner_outer_core
! mesh coloring for GPUs
@@ -814,9 +819,9 @@
if( ier /= 0 ) &
call exit_mpi(myrank,'error allocating num_elem_colors_outer_core array')
endif
-
+
close(IIN)
-
+
end subroutine read_mesh_databases_MPI_OC
!
@@ -824,7 +829,7 @@
!
subroutine read_mesh_databases_MPI_IC()
-
+
use specfem_par
use specfem_par_innercore
implicit none
@@ -833,10 +838,10 @@
integer :: ier
! crust mantle region
-
+
! create the name for the database of the current slide and region
call create_name_database(prname,myrank,IREGION_INNER_CORE,LOCAL_PATH)
-
+
open(unit=IIN,file=prname(1:len_trim(prname))//'solver_data_mpi.bin', &
status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening solver_data_mpi.bin')
@@ -848,21 +853,21 @@
stat=ier)
if( ier /= 0 ) &
call exit_mpi(myrank,'error allocating array my_neighbours_inner_core etc.')
-
+
if( num_interfaces_inner_core > 0 ) then
read(IIN) max_nibool_interfaces_inner_core
allocate(ibool_interfaces_inner_core(max_nibool_interfaces_inner_core,num_interfaces_inner_core), &
stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating array ibool_interfaces_inner_core')
-
+
read(IIN) my_neighbours_inner_core
read(IIN) nibool_interfaces_inner_core
read(IIN) ibool_interfaces_inner_core
else
! dummy array
max_nibool_interfaces_inner_core = 0
- allocate(ibool_interfaces_inner_core(0,0),stat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error allocating array dummy ibool_interfaces_inner_core')
+ allocate(ibool_interfaces_inner_core(0,0),stat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error allocating array dummy ibool_interfaces_inner_core')
endif
! inner / outer elements
@@ -875,7 +880,7 @@
stat=ier)
if( ier /= 0 ) &
call exit_mpi(myrank,'error allocating array phase_ispec_inner_inner_core')
-
+
if(num_phase_ispec_inner_core > 0 ) read(IIN) phase_ispec_inner_inner_core
! mesh coloring for GPUs
@@ -898,9 +903,9 @@
if( ier /= 0 ) &
call exit_mpi(myrank,'error allocating num_elem_colors_inner_core array')
endif
-
+
close(IIN)
-
+
end subroutine read_mesh_databases_MPI_IC
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_topography_bathymetry.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_topography_bathymetry.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_topography_bathymetry.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -43,25 +43,27 @@
! make ellipticity
if( ELLIPTICITY_VAL ) then
+ ! splines used for locating exact source/receivers positions
+ ! in locate_sources() and locate_receivers() routines
call make_ellipticity(nspl,rspl,espl,espl2,ONE_CRUST)
endif
! read topography and bathymetry file
- if( TOPOGRAPHY ) then
+ if( TOPOGRAPHY .or. OCEANS_VAL ) then
! initializes
ibathy_topo(:,:) = 0
! master reads file
if(myrank == 0 ) then
! user output
- write(IMAIN,*) 'topography:'
+ write(IMAIN,*) 'topography:'
! reads topo file
call read_topo_bathy_database(ibathy_topo,LOCAL_PATH)
endif
-
+
! broadcast the information read on the master to the nodes
- call MPI_BCAST(ibathy_topo,NX_BATHY*NY_BATHY,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(ibathy_topo,NX_BATHY_VAL*NY_BATHY_VAL,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
endif
! user output
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -287,7 +287,7 @@
write(IMAIN,*)
endif
call sync_all()
-
+
! locate receivers in the crust in the mesh
call locate_receivers(NSPEC_CRUST_MANTLE,NGLOB_CRUST_MANTLE,ibool_crust_mantle, &
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
@@ -395,7 +395,7 @@
deallocate(theta_source,phi_source)
call sync_all()
-
+
end subroutine setup_receivers
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -61,7 +61,7 @@
!-----------------------------------------------------------------
! use integer array to store values
- integer, dimension(NX_BATHY,NY_BATHY) :: ibathy_topo
+ integer, dimension(NX_BATHY_VAL,NY_BATHY_VAL) :: ibathy_topo
! additional mass matrix for ocean load
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
@@ -75,7 +75,7 @@
! for ellipticity
integer :: nspl
- double precision :: rspl(NR),espl(NR),espl2(NR)
+ double precision,dimension(NR) :: rspl,espl,espl2
!-----------------------------------------------------------------
! rotation
@@ -310,10 +310,6 @@
integer, dimension(:), allocatable :: request_send_scalar_outer_core,request_recv_scalar_outer_core
integer, dimension(:), allocatable :: b_request_send_scalar_outer_core,b_request_recv_scalar_outer_core
-
- ! temporary array
- logical, dimension(NGLOB_CRUST_MANTLE) :: mask_ibool
-
!-----------------------------------------------------------------
! gpu
!-----------------------------------------------------------------
@@ -373,6 +369,7 @@
! displacement, velocity, acceleration
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
displ_crust_mantle,veloc_crust_mantle,accel_crust_mantle
+
! ADJOINT
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE_ADJOINT) :: &
b_displ_crust_mantle,b_veloc_crust_mantle,b_accel_crust_mantle
@@ -535,6 +532,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STRAIN_ONLY) :: &
eps_trace_over_3_inner_core
+
! ADJOINT
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_AND_ATT) :: &
b_R_xx_inner_core,b_R_yy_inner_core,b_R_xy_inner_core,b_R_xz_inner_core,b_R_yz_inner_core
@@ -604,11 +602,11 @@
! velocity potential
real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: &
displ_outer_core,veloc_outer_core,accel_outer_core
+
! ADJOINT
real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE_ADJOINT) :: &
b_displ_outer_core,b_veloc_outer_core,b_accel_outer_core
-
! Stacey
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_STACEY) :: vp_outer_core
integer :: nspec2D_xmin_outer_core,nspec2D_xmax_outer_core, &
@@ -625,8 +623,8 @@
reclen_ymin_outer_core, reclen_ymax_outer_core
integer :: reclen_zmin
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_OUTER_CORE) :: vector_accel_outer_core,&
- vector_displ_outer_core, b_vector_displ_outer_core
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_OUTER_CORE_ADJOINT) :: &
+ vector_accel_outer_core,vector_displ_outer_core,b_vector_displ_outer_core
! arrays to couple with the fluid regions by pointwise matching
integer, dimension(NSPEC2DMAX_XMIN_XMAX_OC) :: ibelm_xmin_outer_core,ibelm_xmax_outer_core
@@ -683,30 +681,34 @@
! to save movie frames
integer :: nmovie_points,NIT
+
real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
store_val_x,store_val_y,store_val_z, &
store_val_ux,store_val_uy,store_val_uz
+
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
store_val_x_all,store_val_y_all,store_val_z_all, &
store_val_ux_all,store_val_uy_all,store_val_uz_all
! to save movie volume
- integer :: npoints_3dmovie,nspecel_3dmovie
- integer, dimension(NGLOB_CRUST_MANTLE) :: num_ibool_3dmovie
double precision :: scalingval
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: &
muvstore_crust_mantle_3dmovie
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: nu_3dmovie
- logical, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: mask_3dmovie
-! real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
-! Iepsilondev_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: &
Iepsilondev_xx_crust_mantle,Iepsilondev_yy_crust_mantle,Iepsilondev_xy_crust_mantle, &
Iepsilondev_xz_crust_mantle,Iepsilondev_yz_crust_mantle
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: &
Ieps_trace_over_3_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: nu_3dmovie
+
+ integer :: npoints_3dmovie,nspecel_3dmovie
+ integer, dimension(NGLOB_CRUST_MANTLE_3DMOVIE) :: num_ibool_3dmovie
+
+ logical, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: mask_3dmovie
+ logical, dimension(NGLOB_CRUST_MANTLE_3DMOVIE) :: mask_ibool
+
end module specfem_par_movie
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_output.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_output.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_output.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -60,6 +60,30 @@
! save movie in full 3D mesh
if(MOVIE_VOLUME ) then
+
+ ! updates integral of strain for adjoint movie volume
+ if( MOVIE_VOLUME_TYPE == 2 .or. MOVIE_VOLUME_TYPE == 3 ) then
+ ! transfers strain arrays onto CPU
+ if( GPU_MODE ) then
+ call transfer_strain_cm_from_device(Mesh_pointer,eps_trace_over_3_crust_mantle, &
+ epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle, &
+ epsilondev_xy_crust_mantle,epsilondev_xz_crust_mantle, &
+ epsilondev_yz_crust_mantle)
+ endif
+
+ ! integrates strain
+ call movie_volume_integrate_strain(deltat,size(Ieps_trace_over_3_crust_mantle,4), &
+ eps_trace_over_3_crust_mantle, &
+ epsilondev_xx_crust_mantle,epsilondev_yy_crust_mantle, &
+ epsilondev_xy_crust_mantle,epsilondev_xz_crust_mantle, &
+ epsilondev_yz_crust_mantle, &
+ Ieps_trace_over_3_crust_mantle, &
+ Iepsilondev_xx_crust_mantle,Iepsilondev_yy_crust_mantle, &
+ Iepsilondev_xy_crust_mantle,Iepsilondev_xz_crust_mantle, &
+ Iepsilondev_yz_crust_mantle)
+ endif
+
+ ! file output
if( mod(it-MOVIE_START,NTSTEP_BETWEEN_FRAMES) == 0 &
.and. it >= MOVIE_START .and. it <= MOVIE_STOP) then
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90 2012-05-01 02:19:33 UTC (rev 20010)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90 2012-05-01 14:58:16 UTC (rev 20011)
@@ -26,7 +26,7 @@
!=====================================================================
- subroutine count_points_movie_volume()
+ subroutine movie_volume_count_points()
! this subroutine counts the number of points and elements within the movie volume
! in this processor slice, and returns arrays that keep track of them, both in global and local indexing schemes
@@ -39,12 +39,14 @@
! local parameters
integer :: ipoints_3dmovie,ispecel_3dmovie,ispec,iglob,i,j,k,iNIT
real(kind=custom_real) :: rval,thetaval,phival
+ integer :: ier
if(MOVIE_COARSE) then
iNIT = NGLLX-1
else
iNIT = 1
endif
+
ipoints_3dmovie=0
num_ibool_3dmovie(:) = -99
ispecel_3dmovie = 0
@@ -54,48 +56,94 @@
! create name of database
write(prname,'(a,i6.6,a)') trim(LOCAL_TMP_PATH)//'/'//'proc',myrank,'_'
- open(unit=IOUT,file=trim(prname)//'movie3D_info.txt',status='unknown')
+ open(unit=IOUT,file=trim(prname)//'movie3D_info.txt', &
+ status='unknown',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening file movie3D_info.txt')
!find and count points within given region for storing movie
- do ispec = 1,NSPEC_CRUST_MANTLE
- !output element if center of element is in the given region
- iglob = ibool_crust_mantle((NGLLX+1)/2,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
- rval = xstore_crust_mantle(iglob)
- thetaval = ystore_crust_mantle(iglob)
- phival = zstore_crust_mantle(iglob)
- ! we alread changed xyz back to rthetaphi
- if( (rval < MOVIE_TOP .and. rval > MOVIE_BOTTOM) .and. &
- (thetaval > MOVIE_NORTH .and. thetaval < MOVIE_SOUTH) .and. &
- ( (phival < MOVIE_EAST .and. phival > MOVIE_WEST) .or. &
- ( (MOVIE_EAST < MOVIE_WEST) .and. (phival >MOVIE_EAST .or. phival < MOVIE_WEST) ) ) ) then
- ispecel_3dmovie=ispecel_3dmovie+1
- do k=1,NGLLZ,iNIT
- do j=1,NGLLY,iNIT
- do i=1,NGLLX,iNIT
- iglob = ibool_crust_mantle(i,j,k,ispec)
- if(.not. mask_ibool(iglob)) then
- ipoints_3dmovie = ipoints_3dmovie + 1
- mask_ibool(iglob)=.true.
- mask_3dmovie(i,j,k,ispec)=.true.
- num_ibool_3dmovie(iglob)= ipoints_3dmovie
- endif
- enddo !i
- enddo !j
- enddo !k
- endif !in region
- enddo !ispec
- npoints_3dmovie=ipoints_3dmovie
- nspecel_3dmovie=ispecel_3dmovie
+ do ispec = 1,NSPEC_CRUST_MANTLE
+ !output element if center of element is in the given region
+ iglob = ibool_crust_mantle((NGLLX+1)/2,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
+ rval = xstore_crust_mantle(iglob)
+ thetaval = ystore_crust_mantle(iglob)
+ phival = zstore_crust_mantle(iglob)
- write(IOUT,*) npoints_3dmovie, nspecel_3dmovie
- close(IOUT)
+ ! we alread changed xyz back to rthetaphi
+ if( (rval < MOVIE_TOP .and. rval > MOVIE_BOTTOM) .and. &
+ (thetaval > MOVIE_NORTH .and. thetaval < MOVIE_SOUTH) .and. &
+ ( (phival < MOVIE_EAST .and. phival > MOVIE_WEST) .or. &
+ ( (MOVIE_EAST < MOVIE_WEST) .and. (phival >MOVIE_EAST .or. phival < MOVIE_WEST) ) ) ) then
+ ispecel_3dmovie=ispecel_3dmovie+1
+ do k=1,NGLLZ,iNIT
+ do j=1,NGLLY,iNIT
+ do i=1,NGLLX,iNIT
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+ if(.not. mask_ibool(iglob)) then
+ ipoints_3dmovie = ipoints_3dmovie + 1
+ mask_ibool(iglob)=.true.
+ mask_3dmovie(i,j,k,ispec)=.true.
+ num_ibool_3dmovie(iglob)= ipoints_3dmovie
+ endif
+ enddo !i
+ enddo !j
+ enddo !k
+ endif !in region
+ enddo !ispec
+ npoints_3dmovie=ipoints_3dmovie
+ nspecel_3dmovie=ispecel_3dmovie
- end subroutine count_points_movie_volume
+ write(IOUT,*) npoints_3dmovie, nspecel_3dmovie
+ close(IOUT)
+ end subroutine movie_volume_count_points
+
!
!-------------------------------------------------------------------------------------------------
!
+ subroutine movie_volume_integrate_strain(deltat,vnspec, &
+ eps_trace_over_3, &
+ epsilondev_xx,epsilondev_yy, &
+ epsilondev_xy,epsilondev_xz, &
+ epsilondev_yz, &
+ Ieps_trace_over_3, &
+ Iepsilondev_xx,Iepsilondev_yy, &
+ Iepsilondev_xy,Iepsilondev_xz, &
+ Iepsilondev_yz)
+
+ use constants_solver
+ implicit none
+
+ real(kind=CUSTOM_REAL) :: deltat
+
+ integer :: vnspec
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,vnspec) :: &
+ eps_trace_over_3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,vnspec) :: &
+ epsilondev_xx,epsilondev_yy,epsilondev_xy, &
+ epsilondev_xz,epsilondev_yz
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,vnspec) :: &
+ Ieps_trace_over_3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,vnspec) :: &
+ Iepsilondev_xx,Iepsilondev_yy,Iepsilondev_xy, &
+ Iepsilondev_xz,Iepsilondev_yz
+
+ ! updates integral values
+ Iepsilondev_xx(:,:,:,:) = Iepsilondev_xx(:,:,:,:) + deltat*epsilondev_xx(:,:,:,:)
+ Iepsilondev_yy(:,:,:,:) = Iepsilondev_yy(:,:,:,:) + deltat*epsilondev_yy(:,:,:,:)
+ Iepsilondev_xy(:,:,:,:) = Iepsilondev_xy(:,:,:,:) + deltat*epsilondev_xy(:,:,:,:)
+ Iepsilondev_xz(:,:,:,:) = Iepsilondev_xz(:,:,:,:) + deltat*epsilondev_xz(:,:,:,:)
+ Iepsilondev_yz(:,:,:,:) = Iepsilondev_yz(:,:,:,:) + deltat*epsilondev_yz(:,:,:,:)
+
+ Ieps_trace_over_3(:,:,:,:) = Ieps_trace_over_3(:,:,:,:) + deltat*eps_trace_over_3(:,:,:,:)
+
+ end subroutine movie_volume_integrate_strain
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine write_movie_volume_mesh()
! writes meshfiles to merge with solver snapshots for 3D volume movies. Also computes and outputs
@@ -251,44 +299,52 @@
! input
integer :: myrank,npoints_3dmovie,MOVIE_VOLUME_TYPE,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(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: &
+ 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, &
epsilondev_xz_crust_mantle,epsilondev_yz_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: muvstore_crust_mantle_3dmovie
- logical, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: mask_3dmovie
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: &
+ muvstore_crust_mantle_3dmovie
+
+ logical, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: mask_3dmovie
+
logical :: MOVIE_COARSE
real(kind=CUSTOM_REAL), dimension(3,3,npoints_3dmovie) :: nu_3dmovie
character(len=150) LOCAL_TMP_PATH,outputname
! variables
!character(len=150) prname
- integer :: ipoints_3dmovie,i,j,k,ispec,iNIT
real(kind=CUSTOM_REAL) :: muv_3dmovie
real(kind=CUSTOM_REAL),dimension(3,3) :: eps_loc,eps_loc_new
real(kind=CUSTOM_REAL),dimension(:),allocatable :: store_val3d_NN,store_val3d_EE,store_val3d_ZZ,&
store_val3d_NE,store_val3d_NZ,store_val3d_EZ
+ integer :: ipoints_3dmovie,i,j,k,ispec,iNIT,ier
character(len=1) movie_prefix
- allocate(store_val3d_NN(npoints_3dmovie))
- allocate(store_val3d_EE(npoints_3dmovie))
- allocate(store_val3d_ZZ(npoints_3dmovie))
- allocate(store_val3d_NE(npoints_3dmovie))
- allocate(store_val3d_NZ(npoints_3dmovie))
- allocate(store_val3d_EZ(npoints_3dmovie))
-
+ ! check
if(NDIM /= 3) call exit_MPI(myrank, 'write_movie_volume requires NDIM = 3')
+ ! allocates arrays
+ allocate(store_val3d_NN(npoints_3dmovie), &
+ store_val3d_EE(npoints_3dmovie), &
+ store_val3d_ZZ(npoints_3dmovie), &
+ store_val3d_NE(npoints_3dmovie), &
+ store_val3d_NZ(npoints_3dmovie), &
+ store_val3d_EZ(npoints_3dmovie), &
+ stat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error allocating store_val3d_ .. arrays')
+
if(MOVIE_VOLUME_TYPE == 1) then
- movie_prefix='E' ! strain
+ movie_prefix='E' ! strain
else if(MOVIE_VOLUME_TYPE == 2) then
- movie_prefix='S' ! time integral of strain
+ movie_prefix='S' ! time integral of strain
else if(MOVIE_VOLUME_TYPE == 3) then
- movie_prefix='P' ! potency, or integral of strain x \mu
+ movie_prefix='P' ! potency, or integral of strain x \mu
endif
if(MOVIE_COARSE) then
iNIT = NGLLX-1
@@ -306,20 +362,28 @@
if(mask_3dmovie(i,j,k,ispec)) then
ipoints_3dmovie=ipoints_3dmovie+1
muv_3dmovie=muvstore_crust_mantle_3dmovie(i,j,k,ispec)
- eps_loc(1,1)=eps_trace_over_3_crust_mantle(i,j,k,ispec) + epsilondev_xx_crust_mantle(i,j,k,ispec)
- eps_loc(2,2)=eps_trace_over_3_crust_mantle(i,j,k,ispec) + epsilondev_yy_crust_mantle(i,j,k,ispec)
- eps_loc(3,3)=eps_trace_over_3_crust_mantle(i,j,k,ispec)- &
- epsilondev_xx_crust_mantle(i,j,k,ispec) - epsilondev_yy_crust_mantle(i,j,k,ispec)
+
+ eps_loc(1,1)=eps_trace_over_3_crust_mantle(i,j,k,ispec) + &
+ epsilondev_xx_crust_mantle(i,j,k,ispec)
+ eps_loc(2,2)=eps_trace_over_3_crust_mantle(i,j,k,ispec) + &
+ epsilondev_yy_crust_mantle(i,j,k,ispec)
+ eps_loc(3,3)=eps_trace_over_3_crust_mantle(i,j,k,ispec) - &
+ epsilondev_xx_crust_mantle(i,j,k,ispec) - &
+ epsilondev_yy_crust_mantle(i,j,k,ispec)
+
eps_loc(1,2)=epsilondev_xy_crust_mantle(i,j,k,ispec)
eps_loc(1,3)=epsilondev_xz_crust_mantle(i,j,k,ispec)
eps_loc(2,3)=epsilondev_yz_crust_mantle(i,j,k,ispec)
+
eps_loc(2,1)=eps_loc(1,2)
eps_loc(3,1)=eps_loc(1,3)
eps_loc(3,2)=eps_loc(2,3)
- ! rotate eps_loc to spherical coordinates
- eps_loc_new(:,:) = matmul(matmul(nu_3dmovie(:,:,ipoints_3dmovie),eps_loc(:,:)), transpose(nu_3dmovie(:,:,ipoints_3dmovie)))
+ ! rotate eps_loc to spherical coordinates
+ eps_loc_new(:,:) = matmul(matmul(nu_3dmovie(:,:,ipoints_3dmovie),eps_loc(:,:)), &
+ transpose(nu_3dmovie(:,:,ipoints_3dmovie)))
if(MOVIE_VOLUME_TYPE == 3) eps_loc_new(:,:) = eps_loc(:,:)*muv_3dmovie
+
store_val3d_NN(ipoints_3dmovie)=eps_loc_new(1,1)
store_val3d_EE(ipoints_3dmovie)=eps_loc_new(2,2)
store_val3d_ZZ(ipoints_3dmovie)=eps_loc_new(3,3)
@@ -364,6 +428,9 @@
write(27) store_val3d_EZ(1:npoints_3dmovie)
close(27)
+ deallocate(store_val3d_NN,store_val3d_EE,store_val3d_ZZ, &
+ store_val3d_NE,store_val3d_NZ,store_val3d_EZ)
+
end subroutine write_movie_volume_strains
!
@@ -379,33 +446,47 @@
include "OUTPUT_FILES/values_from_mesher.h"
! input
- integer :: myrank,npoints_3dmovie,MOVIE_VOLUME_TYPE,it
+ integer :: myrank,it
+ integer :: npoints_3dmovie
+ integer :: MOVIE_VOLUME_TYPE
+ logical :: MOVIE_COARSE
+
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
- real(kind=CUSTOM_REAL), dimension(3,NGLOB_CRUST_MANTLE) :: vector_crust_mantle,vector_scaled
+ real(kind=CUSTOM_REAL), dimension(3,NGLOB_CRUST_MANTLE) :: vector_crust_mantle
+
+ double precision :: scalingval
real(kind=CUSTOM_REAL), dimension(3,3,npoints_3dmovie) :: nu_3dmovie
- double precision :: scalingval
- real(kind=CUSTOM_REAL), dimension(3) :: vector_local,vector_local_new
+
logical, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: mask_3dmovie
- logical :: MOVIE_COARSE
+
character(len=150) LOCAL_TMP_PATH
- ! variables
- integer :: ipoints_3dmovie,i,j,k,ispec,iNIT,iglob
- real(kind=CUSTOM_REAL),dimension(:),allocatable :: store_val3d_N,store_val3d_E,store_val3d_Z
+ ! local parameters
+ real(kind=CUSTOM_REAL), dimension(3) :: vector_local,vector_local_new
+ real(kind=CUSTOM_REAL), dimension(:,:),allocatable :: vector_scaled
+ real(kind=CUSTOM_REAL), dimension(:),allocatable :: store_val3d_N,store_val3d_E,store_val3d_Z
+
+ integer :: ipoints_3dmovie,i,j,k,ispec,iNIT,iglob,ier
character(len=150) outputname
character(len=2) movie_prefix
+ ! check
if(NDIM /= 3) call exit_MPI(myrank,'write_movie_volume requires NDIM = 3')
- allocate(store_val3d_N(npoints_3dmovie))
- allocate(store_val3d_E(npoints_3dmovie))
- allocate(store_val3d_Z(npoints_3dmovie))
+ ! allocates arrays
+ allocate(store_val3d_N(npoints_3dmovie), &
+ store_val3d_E(npoints_3dmovie), &
+ store_val3d_Z(npoints_3dmovie), &
+ vector_scaled(3,NGLOB_CRUST_MANTLE), &
+ stat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error allocating store_val3d_N,.. movie arrays')
if(MOVIE_VOLUME_TYPE == 5) then
- movie_prefix='DI' ! displacement
+ movie_prefix='DI' ! displacement
else if(MOVIE_VOLUME_TYPE == 6) then
- movie_prefix='VE' ! velocity
+ movie_prefix='VE' ! velocity
endif
+
if(MOVIE_COARSE) then
iNIT = NGLLX-1
else
@@ -428,7 +509,7 @@
iglob = ibool_crust_mantle(i,j,k,ispec)
vector_local(:) = vector_scaled(:,iglob)
- ! rotate eps_loc to spherical coordinates
+ ! rotate eps_loc to spherical coordinates
vector_local_new(:) = matmul(nu_3dmovie(:,:,ipoints_3dmovie), vector_local(:))
store_val3d_N(ipoints_3dmovie)=vector_local_new(1)
store_val3d_E(ipoints_3dmovie)=vector_local_new(2)
@@ -456,6 +537,8 @@
write(27) store_val3d_Z(1:npoints_3dmovie)
close(27)
+ deallocate(store_val3d_N,store_val3d_E,store_val3d_Z, &
+ vector_scaled)
end subroutine write_movie_volume_vector
More information about the CIG-COMMITS
mailing list