[cig-commits] r20567 - in seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER: EXAMPLES/regional_Greece_small/DATA setup src/cuda src/meshfem3D src/shared src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Aug 13 20:23:42 PDT 2012
Author: danielpeter
Date: 2012-08-13 20:23:41 -0700 (Mon, 13 Aug 2012)
New Revision: 20567
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/regional_Greece_small/DATA/Par_file
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/it_update_displacement_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_MPI_interfaces.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_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_absorb.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_jacobian_boundaries.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_models.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_par.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_atten3D_QRFSI12.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_s362ani.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.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/shared/exit_mpi.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_simulation_stability.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_stacey_outer_core.f90
Log:
updates module usage for attenuation models; updates coupling terms; adds flush statement for main file output
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/regional_Greece_small/DATA/Par_file
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/regional_Greece_small/DATA/Par_file 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/regional_Greece_small/DATA/Par_file 2012-08-14 03:23:41 UTC (rev 20567)
@@ -37,15 +37,15 @@
# parameters describing the Earth model
OCEANS = .true.
-ELLIPTICITY = .true.
+ELLIPTICITY = .false.
TOPOGRAPHY = .true.
-GRAVITY = .true.
-ROTATION = .true.
-ATTENUATION = .true.
-ATTENUATION_NEW = .true.
+GRAVITY = .false.
+ROTATION = .false.
+ATTENUATION = .false.
+ATTENUATION_NEW = .false.
# absorbing boundary conditions for a regional simulation
-ABSORBING_CONDITIONS = .true.
+ABSORBING_CONDITIONS = .false.
# record length in minutes
RECORD_LENGTH_IN_MINUTES = 15.0d0
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in 2012-08-14 03:23:41 UTC (rev 20567)
@@ -560,10 +560,6 @@
!!
!!-----------------------------------------------------------
-! QRFSI12 constants
- integer,parameter :: NKQ=8,MAXL_Q=12
- integer,parameter :: NSQ=(MAXL_Q+1)**2,NDEPTHS_REFQ=913
-
! The meaningful range of Zhao et al. (1994) model is as follows:
! latitude : 32 - 45 N
! longitude: 130-145 E
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-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/compute_coupling_cuda.cu 2012-08-14 03:23:41 UTC (rev 20567)
@@ -68,16 +68,15 @@
// "-1" from index values to convert from Fortran-> C indexing
ispec = ibelm_top_outer_core[iface] - 1;
+ ispec_selected = ibelm_bottom_crust_mantle[iface] - 1;
// only for DOFs exactly on the CMB (top of these elements)
k = NGLLX - 1;
-
// get displacement on the solid side using pointwise matching
- ispec_selected = ibelm_bottom_crust_mantle[iface] - 1;
+ k_corresp = 0;
// get global point number
// corresponding points are located at the bottom of the mantle
- k_corresp = 0;
iglob_cm = ibool_crust_mantle[INDEX4(NGLLX,NGLLX,NGLLX,i,j,k_corresp,ispec_selected)] - 1;
// elastic displacement on global point
@@ -134,16 +133,15 @@
// "-1" from index values to convert from Fortran-> C indexing
ispec = ibelm_bottom_outer_core[iface] - 1;
+ ispec_selected = ibelm_top_inner_core[iface] - 1;
// only for DOFs exactly on the ICB (bottom of these elements)
k = 0;
-
// get displacement on the solid side using pointwise matching
- ispec_selected = ibelm_top_inner_core[iface] - 1;
+ k_corresp = NGLLX - 1;
// get global point number
// corresponding points are located at the bottom of the mantle
- k_corresp = NGLLX - 1;
iglob_ic = ibool_inner_core[INDEX4(NGLLX,NGLLX,NGLLX,i,j,k_corresp,ispec_selected)] - 1;
// elastic displacement on global point
@@ -315,12 +313,12 @@
// "-1" from index values to convert from Fortran-> C indexing
ispec = ibelm_bottom_crust_mantle[iface] - 1;
+ ispec_selected = ibelm_top_outer_core[iface] - 1;
// only for DOFs exactly on the CMB (bottom of these elements)
k = 0;
-
// get velocity potential on the fluid side using pointwise matching
- ispec_selected = ibelm_top_outer_core[iface] - 1;
+ k_corresp = NGLLX - 1;
// get normal on the CMB
nx = normal_top_outer_core[INDEX4(NDIM,NGLLX,NGLLX,0,i,j,iface)]; // (1,i,j,iface)
@@ -329,7 +327,6 @@
// get global point number
// corresponding points are located at the top of the outer core
- k_corresp = NGLLX - 1;
iglob_oc = ibool_outer_core[INDEX4(NGLLX,NGLLX,NGLLX,i,j,k_corresp,ispec_selected)] - 1;
iglob_cm = ibool_crust_mantle[INDEX4(NGLLX,NGLLX,NGLLX,i,j,k,ispec)] - 1;
@@ -355,6 +352,8 @@
}
}
+/* ----------------------------------------------------------------------------------------------- */
+
__global__ void compute_coupling_ICB_fluid_kernel(realw* displ_inner_core,
realw* accel_inner_core,
realw* accel_outer_core,
@@ -384,12 +383,12 @@
// "-1" from index values to convert from Fortran-> C indexing
ispec = ibelm_top_inner_core[iface] - 1;
+ ispec_selected = ibelm_bottom_outer_core[iface] - 1;
// only for DOFs exactly on the ICB (top of these elements)
k = NGLLX - 1;
-
// get velocity potential on the fluid side using pointwise matching
- ispec_selected = ibelm_bottom_outer_core[iface] - 1;
+ k_corresp = 0;
// get normal on the ICB
nx = normal_bottom_outer_core[INDEX4(NDIM,NGLLX,NGLLX,0,i,j,iface)]; // (1,i,j,iface)
@@ -398,7 +397,6 @@
// get global point number
// corresponding points are located at the bottom of the outer core
- k_corresp = 0;
iglob_oc = ibool_outer_core[INDEX4(NGLLX,NGLLX,NGLLX,i,j,k_corresp,ispec_selected)] - 1;
iglob_ic = ibool_inner_core[INDEX4(NGLLX,NGLLX,NGLLX,i,j,k,ispec)] - 1;
@@ -650,12 +648,8 @@
// initializes temporary array to zero
print_CUDA_error_if_any(cudaMemset(mp->d_updated_dof_ocean_load,0,
- sizeof(int)*mp->NGLOB_CRUST_MANTLE_OCEANS),88501);
+ sizeof(int)*(mp->NGLOB_CRUST_MANTLE_OCEANS)),88501);
-#ifdef ENABLE_VERY_SLOW_ERROR_CHECKING
- exit_on_cuda_error("before kernel compute_coupling_ocean_cuda");
-#endif
-
if( *NCHUNKS_VAL != 6 && mp->absorbing_conditions){
compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_accel_crust_mantle,
mp->d_rmassx_crust_mantle,
@@ -671,8 +665,8 @@
// for backward/reconstructed potentials
if( mp->simulation_type == 3 ) {
// re-initializes array
- print_CUDA_error_if_any(cudaMemset(mp->d_updated_dof_ocean_load,0,
- sizeof(int)*mp->NGLOB_CRUST_MANTLE_OCEANS),88502);
+ print_CUDA_error_if_any(cudaMemset(mp->d_b_updated_dof_ocean_load,0,
+ sizeof(int)*(mp->NGLOB_CRUST_MANTLE_OCEANS)),88502);
compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_b_accel_crust_mantle,
mp->d_rmassx_crust_mantle,
@@ -682,7 +676,7 @@
mp->d_normal_top_crust_mantle,
mp->d_ibool_crust_mantle,
mp->d_ibelm_top_crust_mantle,
- mp->d_updated_dof_ocean_load,
+ mp->d_b_updated_dof_ocean_load,
mp->nspec2D_top_crust_mantle);
}
}else{
@@ -700,8 +694,8 @@
// for backward/reconstructed potentials
if( mp->simulation_type == 3 ) {
// re-initializes array
- print_CUDA_error_if_any(cudaMemset(mp->d_updated_dof_ocean_load,0,
- sizeof(int)*mp->NGLOB_CRUST_MANTLE_OCEANS),88502);
+ print_CUDA_error_if_any(cudaMemset(mp->d_b_updated_dof_ocean_load,0,
+ sizeof(int)*(mp->NGLOB_CRUST_MANTLE_OCEANS)),88502);
compute_coupling_ocean_cuda_kernel<<<grid,threads>>>(mp->d_b_accel_crust_mantle,
mp->d_rmassz_crust_mantle,
@@ -711,7 +705,7 @@
mp->d_normal_top_crust_mantle,
mp->d_ibool_crust_mantle,
mp->d_ibelm_top_crust_mantle,
- mp->d_updated_dof_ocean_load,
+ mp->d_b_updated_dof_ocean_load,
mp->nspec2D_top_crust_mantle);
}
}
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/it_update_displacement_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/it_update_displacement_cuda.cu 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/it_update_displacement_cuda.cu 2012-08-14 03:23:41 UTC (rev 20567)
@@ -29,17 +29,10 @@
#include <stdio.h>
#include <cuda.h>
-//#include <cublas.h>
-
#include "config.h"
#include "mesh_constants_cuda.h"
-//#define CUBLAS_ERROR(s,n) if (s != CUBLAS_STATUS_SUCCESS) { \
-//fprintf (stderr, "CUBLAS Memory Write Error @ %d\n",n); \
-//exit(EXIT_FAILURE); }
-
-
/* ----------------------------------------------------------------------------------------------- */
// elastic wavefield
@@ -298,8 +291,8 @@
realw* accel, int size,
realw deltatover2,
realw* rmassx,
- realw* rmassy,
- realw* rmassz) {
+ realw* rmassy,
+ realw* rmassz) {
int id = threadIdx.x + blockIdx.x*blockDim.x + blockIdx.y*gridDim.x*blockDim.x;
/* because of block and grid sizing problems, there is a small */
@@ -320,8 +313,8 @@
__global__ void kernel_3_accel_cuda_device(realw* accel,
int size,
realw* rmassx,
- realw* rmassy,
- realw* rmassz) {
+ realw* rmassy,
+ realw* rmassz) {
int id = threadIdx.x + blockIdx.x*blockDim.x + blockIdx.y*gridDim.x*blockDim.x;
/* because of block and grid sizing problems, there is a small */
@@ -359,7 +352,7 @@
int* SIMULATION_TYPE_f,
realw* b_deltatover2_F,
int* OCEANS,
- int* NCHUNKS_VAL) {
+ int* NCHUNKS_VAL) {
TRACE("kernel_3_a_cuda");
Mesh* mp = (Mesh*)(*Mesh_pointer); // get Mesh from fortran integer wrapper
@@ -527,18 +520,18 @@
mp->d_accel_inner_core,
mp->NGLOB_INNER_CORE,
deltatover2,
- mp->d_rmass_inner_core,
- mp->d_rmass_inner_core,
- mp->d_rmass_inner_core);
+ mp->d_rmass_inner_core,
+ mp->d_rmass_inner_core,
+ mp->d_rmass_inner_core);
if(SIMULATION_TYPE == 3) {
kernel_3_cuda_device<<< grid, threads>>>(mp->d_b_veloc_inner_core,
mp->d_b_accel_inner_core,
mp->NGLOB_INNER_CORE,
b_deltatover2,
- mp->d_rmass_inner_core,
- mp->d_rmass_inner_core,
- mp->d_rmass_inner_core);
+ mp->d_rmass_inner_core,
+ mp->d_rmass_inner_core,
+ mp->d_rmass_inner_core);
}
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/mesh_constants_cuda.h 2012-08-14 03:23:41 UTC (rev 20567)
@@ -489,6 +489,7 @@
// temporary global array: used to synchronize updates on global accel array
int* d_updated_dof_ocean_load;
+ int* d_b_updated_dof_ocean_load;
// ------------------------------------------------------------------ //
// attenuation
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2012-08-14 03:23:41 UTC (rev 20567)
@@ -1107,7 +1107,7 @@
cudaMemcpyHostToDevice),1204);
// allocates mpi buffer for exchange with cpu
print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_send_accel_buffer_crust_mantle),
- 3*(mp->max_nibool_interfaces_crust_mantle)*(mp->num_interfaces_crust_mantle)*sizeof(realw)),4004);
+ NDIM*(mp->max_nibool_interfaces_crust_mantle)*(mp->num_interfaces_crust_mantle)*sizeof(realw)),4004);
}
// inner core mesh
@@ -1127,7 +1127,7 @@
cudaMemcpyHostToDevice),1204);
// allocates mpi buffer for exchange with cpu
print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_send_accel_buffer_inner_core),
- 3*(mp->max_nibool_interfaces_inner_core)*(mp->num_interfaces_inner_core)*sizeof(realw)),4004);
+ NDIM*(mp->max_nibool_interfaces_inner_core)*(mp->num_interfaces_inner_core)*sizeof(realw)),4004);
}
// outer core mesh
@@ -1346,7 +1346,8 @@
// no anisotropy
// transverse isotropy flag
- print_CUDA_error_if_any(cudaMalloc((void**) &mp->d_ispec_is_tiso_crust_mantle, (mp->NSPEC_CRUST_MANTLE)*sizeof(int)),1025);
+ print_CUDA_error_if_any(cudaMalloc((void**) &mp->d_ispec_is_tiso_crust_mantle,
+ (mp->NSPEC_CRUST_MANTLE)*sizeof(int)),1025);
print_CUDA_error_if_any(cudaMemcpy(mp->d_ispec_is_tiso_crust_mantle, h_ispec_is_tiso,
(mp->NSPEC_CRUST_MANTLE)*sizeof(int),cudaMemcpyHostToDevice),1025);
@@ -1510,13 +1511,19 @@
mp->nspec2D_bottom_crust_mantle = *NSPEC2D_BOTTOM_CM;
int size_tcm = NGLL2*(mp->nspec2D_top_crust_mantle);
- print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_normal_top_crust_mantle),sizeof(realw)*NDIM*size_tcm),40020);
- print_CUDA_error_if_any(cudaMemcpy(mp->d_normal_top_crust_mantle,h_normal_top_crust_mantle,sizeof(realw)*NDIM*size_tcm,cudaMemcpyHostToDevice),40030);
+ print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_normal_top_crust_mantle),
+ sizeof(realw)*NDIM*size_tcm),40020);
+ print_CUDA_error_if_any(cudaMemcpy(mp->d_normal_top_crust_mantle,h_normal_top_crust_mantle,
+ sizeof(realw)*NDIM*size_tcm,cudaMemcpyHostToDevice),40030);
- print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_ibelm_top_crust_mantle),sizeof(int)*(mp->nspec2D_top_crust_mantle)),40021);
- print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_ibelm_bottom_crust_mantle),sizeof(int)*(mp->nspec2D_bottom_crust_mantle)),40021);
- print_CUDA_error_if_any(cudaMemcpy(mp->d_ibelm_top_crust_mantle,h_ibelm_top_crust_mantle,sizeof(int)*(mp->nspec2D_top_crust_mantle),cudaMemcpyHostToDevice),40031);
- print_CUDA_error_if_any(cudaMemcpy(mp->d_ibelm_bottom_crust_mantle,h_ibelm_bottom_crust_mantle,sizeof(int)*(mp->nspec2D_bottom_crust_mantle),cudaMemcpyHostToDevice),40031);
+ print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_ibelm_top_crust_mantle),
+ sizeof(int)*(mp->nspec2D_top_crust_mantle)),40021);
+ print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_ibelm_bottom_crust_mantle),
+ sizeof(int)*(mp->nspec2D_bottom_crust_mantle)),40021);
+ print_CUDA_error_if_any(cudaMemcpy(mp->d_ibelm_top_crust_mantle,h_ibelm_top_crust_mantle,
+ sizeof(int)*(mp->nspec2D_top_crust_mantle),cudaMemcpyHostToDevice),40031);
+ print_CUDA_error_if_any(cudaMemcpy(mp->d_ibelm_bottom_crust_mantle,h_ibelm_bottom_crust_mantle,
+ sizeof(int)*(mp->nspec2D_bottom_crust_mantle),cudaMemcpyHostToDevice),40031);
// wavefield
int size = NDIM * mp->NGLOB_CRUST_MANTLE;
@@ -2030,7 +2037,7 @@
extern "C"
void FC_FUNC_(prepare_oceans_device,
PREPARE_OCEANS_DEVICE)(long* Mesh_pointer_f,
- realw* h_rmass_ocean_load) {
+ realw* h_rmass_ocean_load) {
TRACE("prepare_oceans_device");
@@ -2044,9 +2051,12 @@
// temporary global array: used to synchronize updates on global accel array
print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_updated_dof_ocean_load),
- sizeof(int)*mp->NGLOB_CRUST_MANTLE_OCEANS),4502);
+ sizeof(int)*mp->NGLOB_CRUST_MANTLE_OCEANS),4503);
+ if( mp->simulation_type == 3 ){
+ print_CUDA_error_if_any(cudaMalloc((void**)&(mp->d_b_updated_dof_ocean_load),
+ sizeof(int)*mp->NGLOB_CRUST_MANTLE_OCEANS),4504);
+ }
-
#ifdef ENABLE_VERY_SLOW_ERROR_CHECKING
exit_on_cuda_error("prepare_oceans_device");
#endif
@@ -2464,6 +2474,9 @@
if( mp->oceans ){
cudaFree(mp->d_rmass_ocean_load);
cudaFree(mp->d_updated_dof_ocean_load);
+ if( mp->simulation_type == 3 ){
+ cudaFree(mp->d_b_updated_dof_ocean_load);
+ }
}
if( mp->gravity ){
@@ -2509,23 +2522,6 @@
}
cudaFree(mp->d_rmass_inner_core);
-/*
-
- if( *OCEANS ){
- if( mp->num_free_surface_faces > 0 ){
- cudaFree(mp->d_rmass_ocean_load);
- cudaFree(mp->d_free_surface_normal);
- cudaFree(mp->d_updated_dof_ocean_load);
- if( *NOISE_TOMOGRAPHY == 0){
- cudaFree(mp->d_free_surface_ispec);
- cudaFree(mp->d_free_surface_ijk);
- }
- }
- }
- } // ELASTIC_SIMULATION
-
-*/
-
// releases previous contexts
cudaThreadExit();
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-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_MPI_interfaces.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -167,14 +167,6 @@
case( IREGION_INNER_CORE )
! inner core
- allocate(ibelm_xmin_inner_core(NSPEC2DMAX_XMIN_XMAX_IC), &
- ibelm_xmax_inner_core(NSPEC2DMAX_XMIN_XMAX_IC))
- allocate(ibelm_ymin_inner_core(NSPEC2DMAX_YMIN_YMAX_IC), &
- ibelm_ymax_inner_core(NSPEC2DMAX_YMIN_YMAX_IC))
- allocate(ibelm_bottom_inner_core(NSPEC2D_BOTTOM_IC))
- allocate(ibelm_top_inner_core(NSPEC2D_TOP_IC))
-
-
allocate(iboolcorner_inner_core(NGLOB1D_RADIAL_IC,NUMCORNERS_SHARED))
allocate(iboolleft_xi_inner_core(NGLOB2DMAX_XMIN_XMAX_IC), &
iboolright_xi_inner_core(NGLOB2DMAX_XMIN_XMAX_IC))
@@ -290,6 +282,12 @@
! local parameters
integer :: ier
+ ! for central cube buffers
+ integer :: nspec2D_xmin_inner_core,nspec2D_xmax_inner_core, &
+ nspec2D_ymin_inner_core,nspec2D_ymax_inner_core
+ integer, dimension(:),allocatable :: ibelm_xmin_inner_core,ibelm_xmax_inner_core
+ integer, dimension(:),allocatable :: ibelm_ymin_inner_core,ibelm_ymax_inner_core
+ integer, dimension(:),allocatable :: ibelm_top_inner_core
! debug file output
character(len=150) :: filename
@@ -402,19 +400,6 @@
iboolfaces_inner_core,npoin2D_faces_inner_core, &
iboolcorner_inner_core)
- ! gets coupling arrays for inner core
- nspec2D_xmin_inner_core = nspec2D_xmin
- nspec2D_xmax_inner_core = nspec2D_xmax
- nspec2D_ymin_inner_core = nspec2D_ymin
- nspec2D_ymax_inner_core = nspec2D_ymax
-
- ibelm_xmin_inner_core(:) = ibelm_xmin(:)
- ibelm_xmax_inner_core(:) = ibelm_xmax(:)
- ibelm_ymin_inner_core(:) = ibelm_ymin(:)
- ibelm_ymax_inner_core(:) = ibelm_ymax(:)
- ibelm_bottom_inner_core(:) = ibelm_bottom(:)
- ibelm_top_inner_core(:) = ibelm_top(:)
-
! central cube buffers
if(INCLUDE_CENTRAL_CUBE) then
@@ -424,6 +409,29 @@
endif
call sync_all()
+ ! allocates boundary indexing arrays for central cube
+ allocate(ibelm_xmin_inner_core(NSPEC2DMAX_XMIN_XMAX_IC), &
+ ibelm_xmax_inner_core(NSPEC2DMAX_XMIN_XMAX_IC), &
+ ibelm_ymin_inner_core(NSPEC2DMAX_YMIN_YMAX_IC), &
+ ibelm_ymax_inner_core(NSPEC2DMAX_YMIN_YMAX_IC), &
+ ibelm_top_inner_core(NSPEC2D_TOP_IC), &
+ ibelm_bottom_inner_core(NSPEC2D_BOTTOM_IC), &
+ stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating central cube index arrays')
+
+ ! gets coupling arrays for inner core
+ nspec2D_xmin_inner_core = nspec2D_xmin
+ nspec2D_xmax_inner_core = nspec2D_xmax
+ nspec2D_ymin_inner_core = nspec2D_ymin
+ nspec2D_ymax_inner_core = nspec2D_ymax
+
+ ibelm_xmin_inner_core(:) = ibelm_xmin(:)
+ ibelm_xmax_inner_core(:) = ibelm_xmax(:)
+ ibelm_ymin_inner_core(:) = ibelm_ymin(:)
+ ibelm_ymax_inner_core(:) = ibelm_ymax(:)
+ ibelm_bottom_inner_core(:) = ibelm_bottom(:)
+ ibelm_top_inner_core(:) = ibelm_top(:)
+
! 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), &
@@ -463,6 +471,11 @@
if(myrank == 0) write(IMAIN,*) ''
+ ! frees memory
+ deallocate(ibelm_xmin_inner_core,ibelm_xmax_inner_core)
+ deallocate(ibelm_ymin_inner_core,ibelm_ymax_inner_core)
+ deallocate(ibelm_top_inner_core)
+
else
! allocate fictitious buffers for cube and slices with a dummy size
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-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -70,8 +70,7 @@
double precision :: elevation,height_oceans
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
- integer :: ispec,i,j,k,iglob
- integer :: ix_oceans,iy_oceans,iz_oceans,ispec_oceans,ispec2D_top_crust
+ integer :: ispec,i,j,k,iglob,ispec2D
! initializes matrices
!
@@ -89,10 +88,6 @@
do k = 1,NGLLZ
do j = 1,NGLLY
do i = 1,NGLLX
-
- weight = wxgll(i)*wygll(j)*wzgll(k)
- iglob = ibool(i,j,k,ispec)
-
! compute the jacobian
xixl = xixstore(i,j,k,ispec)
xiyl = xiystore(i,j,k,ispec)
@@ -108,9 +103,12 @@
- xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ iglob = ibool(i,j,k,ispec)
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+
! definition depends if region is fluid or solid
select case( iregion_code)
-
case( IREGION_CRUST_MANTLE, IREGION_INNER_CORE )
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
@@ -152,24 +150,25 @@
! add contribution of the oceans
! for surface elements exactly at the top of the crust (ocean bottom)
- do ispec2D_top_crust = 1,NSPEC2D_TOP
+ do ispec2D = 1,NSPEC2D_TOP
- ispec_oceans = ibelm_top(ispec2D_top_crust)
+ ! gets spectral element index
+ ispec = ibelm_top(ispec2D)
- iz_oceans = NGLLZ
+ ! assumes elements are order such that k == NGLLZ is top surface
+ k = NGLLZ
- do ix_oceans = 1,NGLLX
- do iy_oceans = 1,NGLLY
-
- iglob=ibool(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
-
+ ! loops over surface points
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
! if 3D Earth with topography, compute local height of oceans
if( TOPOGRAPHY ) then
! get coordinates of current point
- xval = xstore(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
- yval = ystore(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
- zval = zstore(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
+ xval = xstore(i,j,k,ispec)
+ yval = ystore(i,j,k,ispec)
+ zval = zstore(i,j,k,ispec)
! map to latitude and longitude for bathymetry routine
! slightly move points to avoid roundoff problem when exactly on the polar axis
@@ -205,10 +204,13 @@
endif
! take into account inertia of water column
- weight = wxgll(ix_oceans) * wygll(iy_oceans) &
- * dble(jacobian2D_top(ix_oceans,iy_oceans,ispec2D_top_crust)) &
+ weight = wxgll(i) * wygll(j) &
+ * dble(jacobian2D_top(i,j,ispec2D)) &
* dble(RHO_OCEANS) * height_oceans
+ ! gets global point index
+ iglob = ibool(i,j,k,ispec)
+
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
rmass_ocean_load(iglob) = rmass_ocean_load(iglob) + sngl(weight)
@@ -321,19 +323,7 @@
enddo
enddo
-
-! ! read arrays for Stacey conditions
-! open(unit=27,file=prname(1:len_trim(prname))//'stacey.bin', &
-! status='old',form='unformatted',action='read',iostat=ier)
-! if( ier /= 0 ) call exit_mpi(myrank,'error opening stacey.bin in create_mass_matrices')
-! read(27) nimin
-! read(27) nimax
-! read(27) njmin
-! read(27) njmax
-! read(27) nkmin_xi
-! read(27) nkmin_eta
-! close(27)
-
+ ! adds contributions to mass matrix to stabilize stacey conditions
select case(iregion_code)
case(IREGION_CRUST_MANTLE)
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-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_regions_mesh.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -104,6 +104,7 @@
case default
call exit_MPI(myrank,'error ipass value in create_regions_mesh')
end select
+ call flush_IMAIN()
endif
! create the name for the database of the current slide and region
@@ -114,6 +115,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...allocating arrays '
+ call flush_IMAIN()
endif
call crm_allocate_arrays(iregion_code,nspec,ipass, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
@@ -125,6 +127,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...setting up layers '
+ call flush_IMAIN()
endif
call crm_setup_layers(iregion_code,nspec,ipass,NEX_PER_PROC_ETA)
@@ -133,6 +136,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...creating mesh elements '
+ call flush_IMAIN()
endif
call crm_create_elements(iregion_code,nspec,ipass, &
NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
@@ -147,6 +151,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...creating global addressing'
+ call flush_IMAIN()
endif
call crm_setup_indexing(nspec,nglob_theor,npointot)
@@ -156,6 +161,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...creating MPI buffers'
+ call flush_IMAIN()
endif
call crm_setup_mpi_buffers(npointot,nspec,iregion_code)
@@ -173,6 +179,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...precomputing jacobian'
+ call flush_IMAIN()
endif
call get_jacobian_boundaries(myrank,iboun,nspec,xstore,ystore,zstore, &
dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
@@ -193,6 +200,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...creating chunk buffers'
+ call flush_IMAIN()
endif
call create_chunk_buffers(iregion_code,nspec,ibool,idoubling, &
xstore,ystore,zstore,nglob_theor, &
@@ -210,6 +218,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...preparing MPI interfaces'
+ call flush_IMAIN()
endif
! creates MPI interface arrays
call create_MPI_interfaces(iregion_code)
@@ -227,6 +236,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...element inner/outer separation '
+ call flush_IMAIN()
endif
call setup_inner_outer(iregion_code)
@@ -235,6 +245,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...element mesh coloring '
+ call flush_IMAIN()
endif
call setup_color_perm(iregion_code)
@@ -266,6 +277,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...creating mass matrix'
+ call flush_IMAIN()
endif
! allocates mass matrices in this slice (will be fully assembled in the solver)
@@ -319,6 +331,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...saving binary files'
+ call flush_IMAIN()
endif
! saves mesh and model parameters
call save_arrays_solver(myrank,nspec,nglob,idoubling,ibool, &
@@ -344,8 +357,11 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...saving boundary mesh files'
+ call flush_IMAIN()
endif
+ ! saves boundary file
call save_arrays_solver_boundary()
+
endif
! compute volume, bottom and top area of that part of the slice
@@ -367,6 +383,7 @@
if( myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' ...saving AVS mesh files'
+ call flush_IMAIN()
endif
call crm_save_mesh_files(nspec,npointot,iregion_code)
endif
@@ -897,13 +914,20 @@
if(myrank == 0 ) then
! time estimate
tCPU = wtime() - time_start
- ! remaining
+
+ ! outputs remaining time (poor estimation)
tCPU = (1.0-(ilayer_loop-ifirst_region+1.0)/(ilast_region-ifirst_region+1.0)) &
/(ilayer_loop-ifirst_region+1.0)/(ilast_region-ifirst_region+1.0)*tCPU*10.0
+
write(IMAIN,*) " ",(ilayer_loop-ifirst_region+1.0)/(ilast_region-ifirst_region+1.0) * 100.0," %", &
" time remaining:", tCPU,"s"
+
+ ! outputs current time on system
call date_and_time(VALUES=tval)
write(IMAIN,*) " ",tval(5),"h",tval(6),"min",tval(7),".",tval(8),"sec"
+
+ ! flushes I/O buffer
+ call flush_IMAIN()
endif
enddo !ilayer_loop
@@ -1215,55 +1239,19 @@
integer,intent(in):: iregion_code
! free memory
- deallocate(iprocfrom_faces,iprocto_faces,imsg_type)
- deallocate(iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners)
- deallocate(buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar)
- deallocate(buffer_send_chunkcorn_vector,buffer_recv_chunkcorn_vector)
-
select case( iregion_code )
case( IREGION_CRUST_MANTLE )
! crust mantle
- deallocate(iboolcorner_crust_mantle)
- deallocate(iboolleft_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)
-
case( IREGION_OUTER_CORE )
! outer core
- deallocate(iboolcorner_outer_core)
- deallocate(iboolleft_xi_outer_core, &
- iboolright_xi_outer_core)
- deallocate(iboolleft_eta_outer_core, &
- iboolright_eta_outer_core)
- deallocate(iboolfaces_outer_core)
-
deallocate(phase_ispec_inner_outer_core)
deallocate(num_elem_colors_outer_core)
-
case( IREGION_INNER_CORE )
! inner core
- deallocate(ibelm_xmin_inner_core, &
- ibelm_xmax_inner_core)
- deallocate(ibelm_ymin_inner_core, &
- ibelm_ymax_inner_core)
- deallocate(ibelm_bottom_inner_core)
- deallocate(ibelm_top_inner_core)
-
- deallocate(iboolcorner_inner_core)
- deallocate(iboolleft_xi_inner_core, &
- iboolright_xi_inner_core)
- deallocate(iboolleft_eta_inner_core, &
- iboolright_eta_inner_core)
- deallocate(iboolfaces_inner_core)
-
deallocate(phase_ispec_inner_inner_core)
deallocate(num_elem_colors_inner_core)
-
end select
end subroutine crm_free_MPI_arrays
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_absorb.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_absorb.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_absorb.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -47,7 +47,7 @@
logical :: iboun(6,nspec)
! global element numbering
- integer :: ispecg
+ integer :: ispec
! counters to keep track of the number of elements on each of the
! five absorbing boundaries
@@ -64,11 +64,11 @@
ispecb4=0
ispecb5=0
- do ispecg=1,nspec
+ do ispec=1,nspec
! determine if the element falls on an absorbing boundary
- if(iboun(1,ispecg)) then
+ if(iboun(1,ispec)) then
! on boundary 1: xmin
ispecb1=ispecb1+1
@@ -79,10 +79,10 @@
! check for ovelap with other boundaries
nkmin_xi(1,ispecb1)=1
- if(iboun(5,ispecg)) nkmin_xi(1,ispecb1)=2
+ if(iboun(5,ispec)) nkmin_xi(1,ispecb1)=2
endif
- if(iboun(2,ispecg)) then
+ if(iboun(2,ispec)) then
! on boundary 2: xmax
ispecb2=ispecb2+1
@@ -93,39 +93,39 @@
! check for ovelap with other boundaries
nkmin_xi(2,ispecb2)=1
- if(iboun(5,ispecg)) nkmin_xi(2,ispecb2)=2
+ if(iboun(5,ispec)) nkmin_xi(2,ispecb2)=2
endif
- if(iboun(3,ispecg)) then
+ if(iboun(3,ispec)) then
! on boundary 3: ymin
ispecb3=ispecb3+1
! check for ovelap with other boundaries
nimin(1,ispecb3)=1
- if(iboun(1,ispecg)) nimin(1,ispecb3)=2
+ if(iboun(1,ispec)) nimin(1,ispecb3)=2
nimax(1,ispecb3)=NGLLX
- if(iboun(2,ispecg)) nimax(1,ispecb3)=NGLLX-1
+ if(iboun(2,ispec)) nimax(1,ispecb3)=NGLLX-1
nkmin_eta(1,ispecb3)=1
- if(iboun(5,ispecg)) nkmin_eta(1,ispecb3)=2
+ if(iboun(5,ispec)) nkmin_eta(1,ispecb3)=2
endif
- if(iboun(4,ispecg)) then
+ if(iboun(4,ispec)) then
! on boundary 4: ymax
ispecb4=ispecb4+1
! check for ovelap with other boundaries
nimin(2,ispecb4)=1
- if(iboun(1,ispecg)) nimin(2,ispecb4)=2
+ if(iboun(1,ispec)) nimin(2,ispecb4)=2
nimax(2,ispecb4)=NGLLX
- if(iboun(2,ispecg)) nimax(2,ispecb4)=NGLLX-1
+ if(iboun(2,ispec)) nimax(2,ispecb4)=NGLLX-1
nkmin_eta(2,ispecb4)=1
- if(iboun(5,ispecg)) nkmin_eta(2,ispecb4)=2
+ if(iboun(5,ispec)) nkmin_eta(2,ispecb4)=2
endif
! on boundary 5: bottom
- if(iboun(5,ispecg)) ispecb5=ispecb5+1
+ if(iboun(5,ispec)) ispecb5=ispecb5+1
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_jacobian_boundaries.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_jacobian_boundaries.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -42,46 +42,45 @@
include "constants.h"
- integer nspec,myrank
- integer NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX
+ integer :: nspec,myrank
+ integer :: NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX
- integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax
- integer ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),ibelm_xmax(NSPEC2DMAX_XMIN_XMAX)
- integer ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),ibelm_ymax(NSPEC2DMAX_YMIN_YMAX)
- integer ibelm_bottom(NSPEC2D_BOTTOM),ibelm_top(NSPEC2D_TOP)
+ integer :: nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax
+ integer :: ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),ibelm_xmax(NSPEC2DMAX_XMIN_XMAX)
+ integer :: ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),ibelm_ymax(NSPEC2DMAX_YMIN_YMAX)
+ integer :: ibelm_bottom(NSPEC2D_BOTTOM)
+ integer :: ibelm_top(NSPEC2D_TOP)
- logical iboun(6,nspec)
+ logical :: iboun(6,nspec)
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision,dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
- real(kind=CUSTOM_REAL) jacobian2D_xmin(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- real(kind=CUSTOM_REAL) jacobian2D_xmax(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- real(kind=CUSTOM_REAL) jacobian2D_ymin(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- real(kind=CUSTOM_REAL) jacobian2D_ymax(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- real(kind=CUSTOM_REAL) jacobian2D_bottom(NGLLX,NGLLY,NSPEC2D_BOTTOM)
- real(kind=CUSTOM_REAL) jacobian2D_top(NGLLX,NGLLY,NSPEC2D_TOP)
+ real(kind=CUSTOM_REAL) :: jacobian2D_xmin(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+ real(kind=CUSTOM_REAL) :: jacobian2D_xmax(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+ real(kind=CUSTOM_REAL) :: jacobian2D_ymin(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ real(kind=CUSTOM_REAL) :: jacobian2D_ymax(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ real(kind=CUSTOM_REAL) :: jacobian2D_bottom(NGLLX,NGLLY,NSPEC2D_BOTTOM)
+ real(kind=CUSTOM_REAL) :: jacobian2D_top(NGLLX,NGLLY,NSPEC2D_TOP)
- real(kind=CUSTOM_REAL) normal_xmin(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- real(kind=CUSTOM_REAL) normal_xmax(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- real(kind=CUSTOM_REAL) normal_ymin(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- real(kind=CUSTOM_REAL) normal_ymax(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- real(kind=CUSTOM_REAL) normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM)
- real(kind=CUSTOM_REAL) normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP)
+ real(kind=CUSTOM_REAL) :: normal_xmin(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+ real(kind=CUSTOM_REAL) :: normal_xmax(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+ real(kind=CUSTOM_REAL) :: normal_ymin(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ real(kind=CUSTOM_REAL) :: normal_ymax(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ real(kind=CUSTOM_REAL) :: normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM)
+ real(kind=CUSTOM_REAL) :: normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP)
- double precision dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ)
- double precision dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ)
- double precision dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY)
- double precision dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY)
+ double precision :: dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ)
+ double precision :: dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ)
+ double precision :: dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY)
+ double precision :: dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY)
! global element numbering
- integer ispec
+ integer :: ispec
! counters to keep track of number of elements on each of the boundaries
- integer ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6
+ integer :: ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6
- double precision xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D)
+ double precision :: xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D)
! Parameters used to calculate 2D Jacobian based upon 25 GLL points
integer:: i,j,k
@@ -421,7 +420,7 @@
zelm(9)=zstore((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec)
call compute_jacobian_2D(myrank,ispecb6,xelm,yelm,zelm,dershape2D_top, &
- jacobian2D_top,normal_top,NGLLX,NGLLY,NSPEC2D_TOP)
+ jacobian2D_top,normal_top,NGLLX,NGLLY,NSPEC2D_TOP)
else
! get 25 GLL points for zmax
do j = 1,NGLLY
@@ -433,8 +432,8 @@
end do
! recalcuate jacobian according to 2D gll points
call calc_jacobian_gll2D(myrank,ispecb6,xelm2D,yelm2D,zelm2D,&
- xigll,yigll,jacobian2D_top,normal_top,&
- NGLLX,NGLLY,NSPEC2D_TOP)
+ xigll,yigll,jacobian2D_top,normal_top,&
+ NGLLX,NGLLY,NSPEC2D_TOP)
end if
@@ -461,7 +460,8 @@
! -------------------------------------------------------
- subroutine compute_jacobian_2D(myrank,ispecb,xelm,yelm,zelm,dershape2D,jacobian2D,normal,NGLLA,NGLLB,NSPEC2DMAX_AB)
+ subroutine compute_jacobian_2D(myrank,ispecb,xelm,yelm,zelm,dershape2D, &
+ jacobian2D,normal,NGLLA,NGLLB,NSPEC2DMAX_AB)
implicit none
@@ -506,7 +506,7 @@
unz=xxi*yeta-xeta*yxi
jacobian=dsqrt(unx**2+uny**2+unz**2)
- if(jacobian == ZERO) call exit_MPI(myrank,'2D Jacobian undefined')
+ if(jacobian <= ZERO) call exit_MPI(myrank,'2D Jacobian undefined')
! normalize normal vector and store surface jacobian
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-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_models.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -160,12 +160,12 @@
! 3D attenuation
if( ATTENUATION_3D) then
! Colleen's model defined originally between 24.4km and 650km
- call model_atten3D_QRFSI12_broadcast(myrank,QRFSI12_Q)
+ call model_atten3D_QRFSI12_broadcast(myrank)
else
! sets up attenuation coefficients according to the chosen, "pure" 1D model
! (including their 1D-crustal profiles)
- call model_attenuation_setup(REFERENCE_1D_MODEL, RICB, RCMB, &
- R670, R220, R80,AM_V,AM_S,AS_V)
+ call model_attenuation_setup(myrank,REFERENCE_1D_MODEL,RICB,RCMB, &
+ R670,R220,R80,AM_V,AM_S,AS_V,CRUSTAL)
endif
endif
@@ -605,7 +605,7 @@
! ! Get the value of Qmu (Attenuation) dependedent on
! ! the radius (r_prem) and idoubling flag
! !call model_attenuation_1D_PREM(r_prem, Qmu, idoubling)
-! call model_atten3D_QRFSI12(r_prem*R_EARTH_KM,theta_degrees,phi_degrees,Qmu,QRFSI12_Q,idoubling)
+! call model_atten3D_QRFSI12(r_prem*R_EARTH_KM,theta_degrees,phi_degrees,Qmu,idoubling)
! ! Get tau_e from tau_s and Qmu
! call model_attenuation_getstored_tau(Qmu, T_c_source, tau_s, tau_e, AM_V, AM_S, AS_V)
! endif
@@ -826,16 +826,16 @@
double precision moho
! attenuation values
- double precision Qkappa,Qmu
+ double precision :: Qkappa,Qmu
double precision, dimension(N_SLS) :: tau_s, tau_e
- double precision T_c_source
+ double precision :: T_c_source
logical elem_in_crust
! local parameters
double precision r_dummy,theta,phi,theta_degrees,phi_degrees
- double precision, parameter :: rmoho_prem = 6371.0-24.4
double precision r_used
+ double precision, parameter :: rmoho_prem = 6371.d0 - 24.4d0
! initializes
tau_e(:) = 0.0d0
@@ -872,7 +872,7 @@
endif ! CRUSTAL
! gets value according to radius/theta/phi location and idoubling flag
- call model_atten3D_QRFSI12(r_used,theta_degrees,phi_degrees,Qmu,QRFSI12_Q,idoubling)
+ call model_atten3D_QRFSI12(r_used,theta_degrees,phi_degrees,Qmu,idoubling)
else
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_par.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_par.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -71,17 +71,6 @@
type (model_attenuation_variables) AM_V
! model_attenuation_variables
-! model_atten3D_QRFSI12_variables
- type model_atten3D_QRFSI12_variables
- sequence
- double precision dqmu(NKQ,NSQ)
- double precision spknt(NKQ)
- double precision refdepth(NDEPTHS_REFQ)
- double precision refqmu(NDEPTHS_REFQ)
- end type model_atten3D_QRFSI12_variables
- type (model_atten3D_QRFSI12_variables) QRFSI12_Q
-! model_atten3D_QRFSI12_variables
-
! model_attenuation_storage_var
type model_attenuation_storage_var
sequence
@@ -394,6 +383,7 @@
ibelm_670_top,ibelm_670_bot
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_moho,normal_400,normal_670
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: jacobian2D_moho,jacobian2D_400,jacobian2D_670
+
integer :: ispec2D_moho_top,ispec2D_moho_bot,ispec2D_400_top,ispec2D_400_bot, &
ispec2D_670_top,ispec2D_670_bot
double precision :: r_moho,r_400,r_670
@@ -584,13 +574,8 @@
integer :: nb_msgs_theor_in_cube,non_zero_nb_msgs_theor_in_cube, &
npoin2D_cube_from_slices,receiver_cube_from_slices
- integer :: nspec2D_xmin_inner_core,nspec2D_xmax_inner_core, &
- nspec2D_ymin_inner_core,nspec2D_ymax_inner_core
-
- integer, dimension(:),allocatable :: ibelm_xmin_inner_core,ibelm_xmax_inner_core
- integer, dimension(:),allocatable :: ibelm_ymin_inner_core,ibelm_ymax_inner_core
+ ! bottom inner core / top central cube
integer, dimension(:),allocatable :: ibelm_bottom_inner_core
- integer, dimension(:),allocatable :: ibelm_top_inner_core
integer, dimension(NUMFACES_SHARED) :: npoin2D_faces_inner_core
integer, dimension(NB_SQUARE_EDGES_ONEDIR) :: npoin2D_xi_inner_core,npoin2D_eta_inner_core
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_atten3D_QRFSI12.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_atten3D_QRFSI12.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_atten3D_QRFSI12.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -37,83 +37,99 @@
! Last edit: Colleen Dalton, March 25, 2008
!
! Q1: what are theta and phi?
+! A1: input theta is colatitude in degrees, phi is longitude in degrees
+!
! Q2: units for radius?
+! A2: radius is given in km
+!
! Q3: what to do about core?
+! A3: more research :)
+!
!--------------------------------------------------------------------------------------------------
+ module model_atten3D_QRFSI12_par
- subroutine model_atten3D_QRFSI12_broadcast(myrank,QRFSI12_Q)
+ ! QRFSI12 constants
+ integer,parameter :: NKQ=8,MAXL_Q=12
+ integer,parameter :: NSQ=(MAXL_Q+1)**2,NDEPTHS_REFQ=913
+
+ ! model_atten3D_QRFSI12_variables
+ double precision,dimension(:,:),allocatable :: QRFSI12_Q_dqmu
+ double precision,dimension(:),allocatable :: QRFSI12_Q_spknt
+ double precision,dimension(:),allocatable :: QRFSI12_Q_refdepth,QRFSI12_Q_refqmu
+
+ end module model_atten3D_QRFSI12_par
+!
+!--------------------------------------------------------------------------------------------------
+!
+
+ subroutine model_atten3D_QRFSI12_broadcast(myrank)
+
! standard routine to setup model
+ use model_atten3D_QRFSI12_par
+
implicit none
include "constants.h"
! standard include of the MPI library
include 'mpif.h'
- ! model_atten3D_QRFSI12_variables
- type model_atten3D_QRFSI12_variables
- sequence
- double precision dqmu(NKQ,NSQ)
- double precision spknt(NKQ)
- double precision refdepth(NDEPTHS_REFQ)
- double precision refqmu(NDEPTHS_REFQ)
- end type model_atten3D_QRFSI12_variables
-
- type (model_atten3D_QRFSI12_variables) QRFSI12_Q
- ! model_atten3D_QRFSI12_variables
-
integer :: myrank
+
+ ! local parameters
integer :: ier
- if(myrank == 0) call read_atten_model_3D_QRFSI12(QRFSI12_Q)
+ ! allocates model arrays
+ allocate(QRFSI12_Q_dqmu(NKQ,NSQ), &
+ QRFSI12_Q_spknt(NKQ), &
+ QRFSI12_Q_refdepth(NDEPTHS_REFQ), &
+ QRFSI12_Q_refqmu(NDEPTHS_REFQ), &
+ stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating QRFSI12_Q model arrays')
- call MPI_BCAST(QRFSI12_Q%dqmu, NKQ*NSQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
- call MPI_BCAST(QRFSI12_Q%spknt, NKQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
- call MPI_BCAST(QRFSI12_Q%refdepth, NDEPTHS_REFQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
- call MPI_BCAST(QRFSI12_Q%refqmu, NDEPTHS_REFQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+ ! master process reads in file values
+ if(myrank == 0) call read_atten_model_3D_QRFSI12()
+ ! broadcasts to all processes
+ call MPI_BCAST(QRFSI12_Q_dqmu, NKQ*NSQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+ call MPI_BCAST(QRFSI12_Q_spknt, NKQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+ call MPI_BCAST(QRFSI12_Q_refdepth, NDEPTHS_REFQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+ call MPI_BCAST(QRFSI12_Q_refqmu, NDEPTHS_REFQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+
if(myrank == 0) write(IMAIN,*) 'read 3D attenuation model'
-
end subroutine
!
!-------------------------------------------------------------------------------------------------
!
- subroutine read_atten_model_3D_QRFSI12(QRFSI12_Q)
+ subroutine read_atten_model_3D_QRFSI12()
+ use model_atten3D_QRFSI12_par
+
implicit none
include "constants.h"
-! three_d_model_atten3D_QRFSI12_variables
- type model_atten3D_QRFSI12_variables
- sequence
- double precision dqmu(NKQ,NSQ)
- double precision spknt(NKQ)
- double precision refdepth(NDEPTHS_REFQ)
- double precision refqmu(NDEPTHS_REFQ)
- end type model_atten3D_QRFSI12_variables
+ ! local parameters
+ integer :: j,k,l,m,ier
+ integer :: index,ll,mm
+ double precision :: v1,v2
- type (model_atten3D_QRFSI12_variables) QRFSI12_Q
-! three_d_model_atten3D_QRFSI12_variables
+ character(len=150) :: QRFSI12,QRFSI12_ref
- integer j,k,l,m
- integer index,ll,mm
- double precision v1,v2
-
- character(len=150) QRFSI12,QRFSI12_ref
-
! read in QRFSI12
! hard-wire for now
QRFSI12='DATA/QRFSI12/QRFSI12.dat'
QRFSI12_ref='DATA/QRFSI12/ref_QRFSI12'
! get the dq model coefficients
- open(unit=10,file=QRFSI12,status='old',action='read')
+ open(unit=10,file=QRFSI12,status='old',action='read',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(0,'error opening model file QRFSI12.dat')
+
do k=1,NKQ
read(10,*)index
j=0
@@ -122,13 +138,13 @@
if(m.eq.0)then
j=j+1
read(10,*)ll,mm,v1
- QRFSI12_Q%dqmu(k,j)=v1
+ QRFSI12_Q_dqmu(k,j)=v1
else
j=j+2
read(10,*)ll,mm,v1,v2
! write(*,*) 'k,l,m,ll,mm:',k,l,m,ll,mm,v1
- QRFSI12_Q%dqmu(k,j-1)=2.*v1
- QRFSI12_Q%dqmu(k,j)=-2.*v2
+ QRFSI12_Q_dqmu(k,j-1)=2.*v1
+ QRFSI12_Q_dqmu(k,j)=-2.*v2
endif
enddo
enddo
@@ -136,64 +152,65 @@
close(10)
! get the depths (km) of the spline knots
- QRFSI12_Q%spknt(1) = 24.4
- QRFSI12_Q%spknt(2) = 75.0
- QRFSI12_Q%spknt(3) = 150.0
- QRFSI12_Q%spknt(4) = 225.0
- QRFSI12_Q%spknt(5) = 300.0
- QRFSI12_Q%spknt(6) = 410.0
- QRFSI12_Q%spknt(7) = 530.0
- QRFSI12_Q%spknt(8) = 650.0
+ QRFSI12_Q_spknt(1) = 24.4d0
+ QRFSI12_Q_spknt(2) = 75.d0
+ QRFSI12_Q_spknt(3) = 150.d0
+ QRFSI12_Q_spknt(4) = 225.d0
+ QRFSI12_Q_spknt(5) = 300.d0
+ QRFSI12_Q_spknt(6) = 410.d0
+ QRFSI12_Q_spknt(7) = 530.d0
+ QRFSI12_Q_spknt(8) = 650.d0
! get the depths and 1/Q values of the reference model
- open(11,file=QRFSI12_ref,status='old',action='read')
+ open(11,file=QRFSI12_ref,status='old',action='read',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(0,'error opening model file ref_QRFSI12')
+
do j=1,NDEPTHS_REFQ
- read(11,*)QRFSI12_Q%refdepth(j),QRFSI12_Q%refqmu(j)
+ read(11,*)QRFSI12_Q_refdepth(j),QRFSI12_Q_refqmu(j)
enddo
close(11)
end subroutine read_atten_model_3D_QRFSI12
+!
!----------------------------------
-!----------------------------------
+!
- subroutine model_atten3D_QRFSI12(radius,theta,phi,Qmu,QRFSI12_Q,idoubling)
+ subroutine model_atten3D_QRFSI12(radius,theta,phi,Qmu,idoubling)
+ use model_atten3D_QRFSI12_par
+
implicit none
include "constants.h"
-! model_atten3D_QRFSI12_variables
- type model_atten3D_QRFSI12_variables
- sequence
- double precision dqmu(NKQ,NSQ)
- double precision spknt(NKQ)
- double precision refdepth(NDEPTHS_REFQ)
- double precision refqmu(NDEPTHS_REFQ)
- end type model_atten3D_QRFSI12_variables
+ double precision :: radius,theta,phi,Qmu
+ integer :: idoubling
- type (model_atten3D_QRFSI12_variables) QRFSI12_Q
-! model_atten3D_QRFSI12_variables
-
- integer i,j,k,n,idoubling
- integer ifnd
- double precision radius,theta,phi,Qmu,smallq,dqmu,smallq_ref
- real(kind=4) splpts(NKQ),splcon(NKQ),splcond(NKQ)
- real(kind=4) depth,ylat,xlon
- real(kind=4) shdep(NSQ)
- real(kind=4) xlmvec(NSQ),wk1(NSQ),wk2(NSQ),wk3(NSQ)
+ ! local parameters
+ integer :: i,j,k,n
+ integer :: ifnd
+ double precision :: smallq,dqmu,smallq_ref
+ real(kind=4) :: splpts(NKQ),splcon(NKQ),splcond(NKQ)
+ real(kind=4) :: depth,ylat,xlon
+ real(kind=4) :: shdep(NSQ)
+ real(kind=4) :: xlmvec(NSQ),wk1(NSQ),wk2(NSQ),wk3(NSQ)
double precision, parameter :: rmoho_prem = 6371.0-24.4
double precision, parameter :: rcmb = 3480.0
- !in Colleen's original code theta refers to the latitude. Here we have redefined theta to be colatitude
- !to agree with the rest of specfem
-! print *,'entering QRFSI12 subroutine'
+ ! in Colleen's original code theta refers to the latitude. Here we have redefined theta to be colatitude
+ ! to agree with the rest of specfem
+
+ ! debug
+ ! print *,'entering QRFSI12 subroutine'
ylat=90.0d0-theta
xlon=phi
-! only checks radius for crust, idoubling is missleading for oceanic crust when we want to expand mantle up to surface...
+ ! only checks radius for crust, idoubling is missleading for oceanic crust
+ ! when we want to expand mantle up to surface...
+
! !if(idoubling == IFLAG_CRUST .or. radius >= rmoho) then
if( radius >= rmoho_prem ) then
! print *,'QRFSI12: we are in the crust'
@@ -201,17 +218,34 @@
else if(idoubling == IFLAG_INNER_CORE_NORMAL .or. idoubling == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
idoubling == IFLAG_BOTTOM_CENTRAL_CUBE .or. idoubling == IFLAG_TOP_CENTRAL_CUBE .or. &
idoubling == IFLAG_IN_FICTITIOUS_CUBE) then
- ! print *,'QRFSI12: we are in the inner core'
- Qmu = 84.6d0
+ ! we are in the inner core
+
+ !debug
+ ! print *,'QRFSI12: we are in the inner core'
+
+ Qmu = 84.6d0
+
else if(idoubling == IFLAG_OUTER_CORE_NORMAL) then
- ! print *,'QRFSI12: we are in the outer core'
- Qmu = 0.0d0
- else !we are in the mantle
+
+ ! we are in the outer core
+
+ !debug
+ ! print *,'QRFSI12: we are in the outer core'
+
+ Qmu = 0.0d0
+
+ else
+
+ ! we are in the mantle
+
depth = 6371.-radius
-! print *,'QRFSI12: we are in the mantle at depth',depth
+
+ !debug
+ ! print *,'QRFSI12: we are in the mantle at depth',depth
+
ifnd=0
do i=2,NDEPTHS_REFQ
- if(depth >= QRFSI12_Q%refdepth(i-1) .and. depth < QRFSI12_Q%refdepth(i))then
+ if(depth >= QRFSI12_Q_refdepth(i-1) .and. depth < QRFSI12_Q_refdepth(i))then
ifnd=i
endif
enddo
@@ -219,7 +253,7 @@
write(6,"('problem finding reference Q value at depth: ',f8.3)") depth
stop
endif
- smallq_ref=QRFSI12_Q%refqmu(ifnd)
+ smallq_ref=QRFSI12_Q_refqmu(ifnd)
smallq = smallq_ref
if(depth < 650.d0) then !Colleen's model is only defined between depths of 24.4 and 650km
@@ -227,12 +261,12 @@
shdep(j)=0.
enddo
do n=1,NKQ
- splpts(n)=QRFSI12_Q%spknt(n)
+ splpts(n)=QRFSI12_Q_spknt(n)
enddo
call vbspl(depth,NKQ,splpts,splcon,splcond)
do n=1,NKQ
do j=1,NSQ
- shdep(j)=shdep(j)+(splcon(n)*QRFSI12_Q%dqmu(n,j))
+ shdep(j)=shdep(j)+(splcon(n)*QRFSI12_Q_dqmu(n,j))
enddo
enddo
call ylm(ylat,xlon,MAXL_Q,xlmvec,wk1,wk2,wk3)
@@ -242,19 +276,23 @@
enddo
smallq = smallq_ref + dqmu
endif
- ! if smallq is small and negative (due to numerical error), Qmu is very large:
+
+ ! if smallq is small and negative (due to numerical error), Qmu is very large:
if(smallq < 0.0d0) smallq = 1.0d0/ATTENUATION_COMP_MAXIMUM
+
Qmu = 1/smallq
- ! Qmu is larger than MAX_ATTENUATION_VALUE, set it to ATTENUATION_COMP_MAXIMUM. This assumes that this
- ! value is high enough that at this point there is almost no attenuation at all.
+
+ ! if Qmu is larger than MAX_ATTENUATION_VALUE, set it to ATTENUATION_COMP_MAXIMUM.
+ ! This assumes that this value is high enough that at this point there is almost no attenuation at all.
if(Qmu >= ATTENUATION_COMP_MAXIMUM) Qmu = 0.99d0*ATTENUATION_COMP_MAXIMUM
endif
end subroutine model_atten3D_QRFSI12
+!
!----------------------------------
-!----------------------------------
+!
!!$ subroutine vbspl(x,np,xarr,splcon,splcond)
!!$!
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_attenuation.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_attenuation.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -88,18 +88,19 @@
integer :: myrank
integer :: ier
+ allocate(AM_V%Qtau_s(N_SLS))
+
+ ! master process determines period ranges
if(myrank == 0) call read_attenuation_model(MIN_ATTENUATION_PERIOD, MAX_ATTENUATION_PERIOD, AM_V)
- if(myrank /= 0) allocate(AM_V%Qtau_s(N_SLS))
+ ! broadcasts to all others
call MPI_BCAST(AM_V%min_period, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
call MPI_BCAST(AM_V%max_period, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
call MPI_BCAST(AM_V%QT_c_source, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
- call MPI_BCAST(AM_V%Qtau_s(1), 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
- call MPI_BCAST(AM_V%Qtau_s(2), 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
- call MPI_BCAST(AM_V%Qtau_s(3), 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+ call MPI_BCAST(AM_V%Qtau_s, N_SLS, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
- end subroutine
+ end subroutine model_attenuation_broadcast
!
!-------------------------------------------------------------------------------------------------
@@ -134,13 +135,11 @@
type (model_attenuation_variables) AM_V
! model_attenuation_variables
- integer min_att_period, max_att_period
+ integer :: min_att_period, max_att_period
AM_V%min_period = min_att_period * 1.0d0
AM_V%max_period = max_att_period * 1.0d0
- allocate(AM_V%Qtau_s(N_SLS))
-
call attenuation_tau_sigma(AM_V%Qtau_s, N_SLS, AM_V%min_period, AM_V%max_period)
call attenuation_source_frequency(AM_V%QT_c_source, AM_V%min_period, AM_V%max_period)
@@ -157,8 +156,8 @@
!
! All this subroutine does is define the Attenuation vs Radius and then Compute the Attenuation
! Variables (tau_sigma and tau_epslion ( or tau_mu) )
- subroutine model_attenuation_setup(REFERENCE_1D_MODEL,RICB,RCMB,R670, &
- R220,R80,AM_V,AM_S,AS_V)
+ subroutine model_attenuation_setup(myrank,REFERENCE_1D_MODEL,RICB,RCMB, &
+ R670,R220,R80,AM_V,AM_S,AS_V,CRUSTAL)
use model_1dref_par, only: &
NR_REF,Mref_V_radius_ref,Mref_V_Qmu_ref
@@ -228,51 +227,55 @@
type(attenuation_simplex_variables) AS_V
! attenuation_simplex_variables
- integer myrank
- integer REFERENCE_1D_MODEL
- double precision RICB, RCMB, R670, R220, R80
- double precision tau_e(N_SLS)
+ integer :: myrank,REFERENCE_1D_MODEL
+ double precision :: RICB, RCMB, R670, R220, R80
+ logical :: CRUSTAL
+
+ ! local parameters
+ double precision :: tau_e(N_SLS)
+ double precision :: Qb
+ double precision :: R120
+ integer :: i,ier
- integer i
- double precision Qb
- double precision R120
-
+ ! parameter definitions
Qb = 57287.0d0
R120 = 6251.d3 ! as defined by IASP91
- call world_rank(myrank)
- if(myrank > 0) return
-
-
! uses "pure" 1D models including their 1D-crust profiles
! (uses USE_EXTERNAL_CRUSTAL_MODEL set to false)
if(REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM) then
- AM_V%Qn = 12
+ AM_V%Qn = 12
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_IASP91) then
- AM_V%Qn = 12
+ AM_V%Qn = 12
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_AK135) then
- call define_model_ak135(.FALSE.)
- AM_V%Qn = NR_AK135
+ ! redefines "pure" 1D model without crustal modification
+ call define_model_ak135(.FALSE.)
+ AM_V%Qn = NR_AK135
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) then
- call define_model_1066a(.FALSE.)
- AM_V%Qn = NR_1066A
+ ! redefines "pure" 1D model without crustal modification
+ call define_model_1066a(.FALSE.)
+ AM_V%Qn = NR_1066A
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_1DREF) then
- call define_model_1dref(.FALSE.)
- AM_V%Qn = NR_REF
+ ! redefines "pure" 1D model without crustal modification
+ call define_model_1dref(.FALSE.)
+ AM_V%Qn = NR_REF
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_JP1D) then
- AM_V%Qn = 12
+ AM_V%Qn = 12
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_SEA1D) then
- call define_model_sea1d(.FALSE.)
- AM_V%Qn = NR_SEA1D
+ ! redefines "pure" 1D model without crustal modification
+ call define_model_sea1d(.FALSE.)
+ AM_V%Qn = NR_SEA1D
else
- call exit_MPI(myrank, 'Reference 1D Model Not recognized')
+ call exit_MPI(myrank, 'Reference 1D Model Not recognized')
endif
! sets up attenuation storage (for all possible Qmu values defined in the 1D models)
- allocate(AM_V%Qr(AM_V%Qn))
- allocate(AM_V%Qmu(AM_V%Qn))
- allocate(AM_V%interval_Q(AM_V%Qn))
- allocate(AM_V%Qtau_e(N_SLS,AM_V%Qn))
+ allocate(AM_V%Qr(AM_V%Qn), &
+ AM_V%Qmu(AM_V%Qn), &
+ AM_V%interval_Q(AM_V%Qn), &
+ AM_V%Qtau_e(N_SLS,AM_V%Qn), &
+ stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating AM_V arrays')
if(REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM) then
AM_V%Qr(:) = (/ 0.0d0, RICB, RICB, RCMB, RCMB, R670, R670, R220, R220, R80, R80, R_EARTH /)
@@ -298,10 +301,27 @@
end if
do i = 1, AM_V%Qn
- call model_attenuation_getstored_tau(AM_V%Qmu(i), AM_V%QT_c_source, AM_V%Qtau_s, tau_e, AM_V, AM_S,AS_V)
- AM_V%Qtau_e(:,i) = tau_e(:)
+ call model_attenuation_getstored_tau(AM_V%Qmu(i), AM_V%QT_c_source, AM_V%Qtau_s, tau_e, AM_V, AM_S,AS_V)
+ AM_V%Qtau_e(:,i) = tau_e(:)
end do
+ ! re-defines 1D models with crustal modification if necessary
+ if( CRUSTAL ) then
+ if(REFERENCE_1D_MODEL == REFERENCE_MODEL_AK135) then
+ ! redefines 1D model with crustal modification
+ call define_model_ak135(CRUSTAL)
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) then
+ ! redefines 1D model with crustal modification
+ call define_model_1066a(CRUSTAL)
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_1DREF) then
+ ! redefines 1D model with crustal modification
+ call define_model_1dref(CRUSTAL)
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_SEA1D) then
+ ! redefines 1D model with crustal modification
+ call define_model_sea1d(CRUSTAL)
+ endif
+ endif
+
end subroutine model_attenuation_setup
!
@@ -368,10 +388,11 @@
type(attenuation_simplex_variables) AS_V
! attenuation_simplex_variables
- double precision Qmu_in, T_c_source
+ double precision :: Qmu_in, T_c_source
double precision, dimension(N_SLS) :: tau_s, tau_e
- integer rw
+ ! local parameters
+ integer :: rw
! READ
rw = 1
@@ -408,12 +429,14 @@
type (model_attenuation_storage_var) AM_S
! model_attenuation_storage_var
- integer myrank
- double precision Qmu, Qmu_new
+ double precision :: Qmu
double precision, dimension(N_SLS) :: tau_e
- integer rw
+ integer :: rw
- integer Qtmp
+ ! local parameters
+ double precision :: Qmu_new
+ integer :: myrank
+ integer :: Qtmp
integer, save :: first_time_called = 1
if(first_time_called == 1) then
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_s362ani.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_s362ani.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_s362ani.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -57,24 +57,21 @@
real(kind=4),dimension(:,:),allocatable :: ylmcof
real(kind=4),dimension(:),allocatable :: wk1,wk2,wk3
- integer lmxhpa(maxhpa)
- integer itypehpa(maxhpa)
- integer ihpakern(maxker)
- integer numcoe(maxhpa)
- integer ivarkern(maxker)
- integer itpspl(maxcoe,maxhpa)
+ integer, dimension(maxhpa) :: lmxhpa,itypehpa,numcoe,nconpt
+
+ integer,dimension(:,:),allocatable :: itpspl,iconpt
+ integer,dimension(:),allocatable :: ihpakern,ivarkern
- integer nconpt(maxhpa),iver
- integer iconpt(maxver,maxhpa)
- integer numker
- integer numhpa,numcof
- integer ihpa,lmax,nylm
+ integer :: iver
+ integer :: numker
+ integer :: numhpa,numcof
+ integer :: ihpa,lmax,nylm
- character(len=80) kerstr
- character(len=80) refmdl
- character(len=40) varstr(maxker)
- character(len=80) hsplfl(maxhpa)
- character(len=40) dskker(maxker)
+ character(len=80) :: kerstr
+ character(len=80) :: refmdl
+ character(len=40) :: varstr(maxker)
+ character(len=80) :: hsplfl(maxhpa)
+ character(len=40) :: dskker(maxker)
end module model_s362ani_par
@@ -110,6 +107,10 @@
wk1(maxl+1), &
wk2(maxl+1), &
wk3(maxl+1), &
+ itpspl(maxcoe,maxhpa), &
+ iconpt(maxver,maxhpa), &
+ ihpakern(maxker), &
+ ivarkern(maxker), &
stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating s362ani arrays')
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -33,7 +33,7 @@
use meshfem3D_models_par,only: &
OCEANS,TRANSVERSE_ISOTROPY,HETEROGEN_3D_MANTLE,ANISOTROPIC_3D_MANTLE, &
- ANISOTROPIC_INNER_CORE,ATTENUATION
+ ANISOTROPIC_INNER_CORE,ATTENUATION,ATTENUATION_3D
use meshfem3D_par,only: &
NCHUNKS,ABSORBING_CONDITIONS,SAVE_MESH_FILES
@@ -73,19 +73,15 @@
integer :: NSPEC2D_TOP,NSPEC2D_BOTTOM
! local parameters
- integer i,j,k,ispec,iglob,nspec1,nglob1,ier
- real(kind=CUSTOM_REAL) :: scaleval1,scaleval2
+ integer :: i,j,k,ispec,iglob,ier
! save nspec and nglob, to be used in combine_paraview_data
open(unit=27,file=prname(1:len_trim(prname))//'array_dims.txt', &
status='unknown',action='write',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening array_dims file')
-
- nspec1 = nspec
- nglob1 = nglob
-
- write(27,*) nspec1
- write(27,*) nglob1
+
+ write(27,*) nspec
+ write(27,*) nglob
close(27)
! mesh parameters
@@ -106,15 +102,6 @@
write(27) rhostore
write(27) kappavstore
- if(HETEROGEN_3D_MANTLE .and. iregion_code == IREGION_CRUST_MANTLE) then
- open(unit=29,file=prname(1:len_trim(prname))//'dvp.bin', &
- status='unknown',form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening dvp.bin file')
-
- write(29) dvpstore
- close(29)
- endif
-
! other terms needed in the solid regions only
if(iregion_code /= IREGION_OUTER_CORE) then
@@ -309,37 +296,120 @@
close(27)
if(ATTENUATION) then
- open(unit=27, file=prname(1:len_trim(prname))//'attenuation.bin', &
+ open(unit=27, file=prname(1:len_trim(prname))//'attenuation.bin', &
status='unknown', form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening attenuation.bin file')
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening attenuation.bin file')
- write(27) tau_s
- write(27) tau_e_store
- write(27) Qmu_store
- write(27) T_c_source
- close(27)
+ write(27) tau_s
+ write(27) tau_e_store
+ write(27) Qmu_store
+ write(27) T_c_source
+ close(27)
endif
+ if(HETEROGEN_3D_MANTLE .and. iregion_code == IREGION_CRUST_MANTLE) then
+ open(unit=27,file=prname(1:len_trim(prname))//'dvp.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening dvp.bin file')
+
+ write(27) dvpstore
+ close(27)
+ endif
+
+
! uncomment for vp & vs model storage
if( SAVE_MESH_FILES ) then
- scaleval1 = sngl( sqrt(PI*GRAV*RHOAV)*(R_EARTH/1000.0d0) )
- scaleval2 = sngl( RHOAV/1000.0d0 )
+ ! outputs model files in binary format
+ call save_arrays_solver_meshfiles(myrank,nspec)
+ endif
- ! isotropic model
- ! vp
- open(unit=27,file=prname(1:len_trim(prname))//'vp.bin', &
- status='unknown',form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening vp.bin file')
+ end subroutine save_arrays_solver
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine save_arrays_solver_meshfiles(myrank,nspec)
+
+! outputs model files in binary format
+
+ use constants
+
+ use meshfem3D_models_par,only: &
+ TRANSVERSE_ISOTROPY,ATTENUATION,ATTENUATION_3D
+
+ use create_regions_mesh_par2,only: &
+ rhostore,kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
+ Qmu_store, &
+ prname
+
+ implicit none
+
+ integer :: myrank
+ integer :: nspec
+
+ ! local parameters
+ integer :: i,j,k,ispec,ier
+ real(kind=CUSTOM_REAL) :: scaleval1,scaleval2
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable :: temp_store
+
+ ! scaling factors to re-dimensionalize units
+ scaleval1 = sngl( sqrt(PI*GRAV*RHOAV)*(R_EARTH/1000.0d0) )
+ scaleval2 = sngl( RHOAV/1000.0d0 )
+
+ ! isotropic model
+ ! vp
+ open(unit=27,file=prname(1:len_trim(prname))//'vp.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening vp.bin file')
+
+ write(27) sqrt( (kappavstore+4.*muvstore/3.)/rhostore )*scaleval1
+ close(27)
+ ! vs
+ open(unit=27,file=prname(1:len_trim(prname))//'vs.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening vs.bin file')
+
+ write(27) sqrt( muvstore/rhostore )*scaleval1
+ close(27)
+ ! rho
+ open(unit=27,file=prname(1:len_trim(prname))//'rho.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening rho.bin file')
+
+ write(27) rhostore*scaleval2
+ close(27)
+
+ ! transverse isotropic model
+ if( TRANSVERSE_ISOTROPY ) then
+ ! vpv
+ open(unit=27,file=prname(1:len_trim(prname))//'vpv.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening vpv.bin file')
+
write(27) sqrt( (kappavstore+4.*muvstore/3.)/rhostore )*scaleval1
close(27)
- ! vs
- open(unit=27,file=prname(1:len_trim(prname))//'vs.bin', &
+ ! vph
+ open(unit=27,file=prname(1:len_trim(prname))//'vph.bin', &
status='unknown',form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening vs.bin file')
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening vph.bin file')
+ write(27) sqrt( (kappahstore+4.*muhstore/3.)/rhostore )*scaleval1
+ close(27)
+ ! vsv
+ open(unit=27,file=prname(1:len_trim(prname))//'vsv.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening vsv.bin file')
+
write(27) sqrt( muvstore/rhostore )*scaleval1
close(27)
+ ! vsh
+ open(unit=27,file=prname(1:len_trim(prname))//'vsh.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening vsh.bin file')
+
+ write(27) sqrt( muhstore/rhostore )*scaleval1
+ close(27)
! rho
open(unit=27,file=prname(1:len_trim(prname))//'rho.bin', &
status='unknown',form='unformatted',action='write',iostat=ier)
@@ -347,56 +417,58 @@
write(27) rhostore*scaleval2
close(27)
+ ! eta
+ open(unit=27,file=prname(1:len_trim(prname))//'eta.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening eta.bin file')
- ! transverse isotropic model
- if( TRANSVERSE_ISOTROPY ) then
- ! vpv
- open(unit=27,file=prname(1:len_trim(prname))//'vpv.bin', &
- status='unknown',form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening vpv.bin file')
+ write(27) eta_anisostore
+ close(27)
+ endif ! TRANSVERSE_ISOTROPY
+
+ ! shear attenuation
+ if( ATTENUATION ) then
+ ! saves Qmu_store to full custom_real array
+ ! uses temporary array
+ allocate(temp_store(NGLLX,NGLLY,NGLLZ,nspec))
+ if (ATTENUATION_3D) then
+ ! attenuation arrays are fully 3D
+ if(CUSTOM_REAL == SIZE_REAL) then
+ temp_store(:,:,:,:) = sngl(Qmu_store(:,:,:,:))
+ else
+ temp_store(:,:,:,:) = Qmu_store(:,:,:,:)
+ endif
+ else
+ ! attenuation array dimensions: Q_mustore(1,1,1,nspec)
+ do ispec = 1,nspec
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ temp_store(i,j,k,ispec) = sngl(Qmu_store(1,1,1,ispec))
+ else
+ temp_store(i,j,k,ispec) = Qmu_store(1,1,1,ispec)
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+ endif
- write(27) sqrt( (kappavstore+4.*muvstore/3.)/rhostore )*scaleval1
- close(27)
- ! vph
- open(unit=27,file=prname(1:len_trim(prname))//'vph.bin', &
- status='unknown',form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening vph.bin file')
+ ! Qmu
+ open(unit=27,file=prname(1:len_trim(prname))//'qmu.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening qmu.bin file')
- write(27) sqrt( (kappahstore+4.*muhstore/3.)/rhostore )*scaleval1
- close(27)
- ! vsv
- open(unit=27,file=prname(1:len_trim(prname))//'vsv.bin', &
- status='unknown',form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening vsv.bin file')
+ write(27) temp_store
+ close(27)
- write(27) sqrt( muvstore/rhostore )*scaleval1
- close(27)
- ! vsh
- open(unit=27,file=prname(1:len_trim(prname))//'vsh.bin', &
- status='unknown',form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening vsh.bin file')
+ ! frees temporary memory
+ deallocate(temp_store)
+ endif ! ATTENUATION
- write(27) sqrt( muhstore/rhostore )*scaleval1
- close(27)
- ! rho
- open(unit=27,file=prname(1:len_trim(prname))//'rho.bin', &
- status='unknown',form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening rho.bin file')
-
- write(27) rhostore*scaleval2
- close(27)
- ! eta
- open(unit=27,file=prname(1:len_trim(prname))//'eta.bin', &
- status='unknown',form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening eta.bin file')
-
- write(27) eta_anisostore
- close(27)
- endif
- endif ! SAVE_MESH_FILES
-
- end subroutine save_arrays_solver
-
+ end subroutine save_arrays_solver_meshfiles
!
!-------------------------------------------------------------------------------------------------
!
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -32,8 +32,11 @@
INCLUDE_CENTRAL_CUBE,myrank,NUMFACES_SHARED
use create_MPI_interfaces_par
- use MPI_inner_core_par,only: &
- non_zero_nb_msgs_theor_in_cube,npoin2D_cube_from_slices
+
+ use MPI_crust_mantle_par
+ use MPI_outer_core_par
+ use MPI_inner_core_par
+
implicit none
integer,intent(in):: iregion_code
@@ -84,6 +87,32 @@
deallocate(ibool_neighbours)
deallocate(my_neighbours,nibool_neighbours)
+ ! frees arrays not needed any further
+ deallocate(iprocfrom_faces,iprocto_faces,imsg_type)
+ deallocate(iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners)
+ deallocate(buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar)
+ deallocate(buffer_send_chunkcorn_vector,buffer_recv_chunkcorn_vector)
+ select case( iregion_code )
+ case( IREGION_CRUST_MANTLE )
+ ! crust mantle
+ deallocate(iboolcorner_crust_mantle)
+ deallocate(iboolleft_xi_crust_mantle,iboolright_xi_crust_mantle)
+ deallocate(iboolleft_eta_crust_mantle,iboolright_eta_crust_mantle)
+ deallocate(iboolfaces_crust_mantle)
+ case( IREGION_OUTER_CORE )
+ ! outer core
+ deallocate(iboolcorner_outer_core)
+ deallocate(iboolleft_xi_outer_core,iboolright_xi_outer_core)
+ deallocate(iboolleft_eta_outer_core,iboolright_eta_outer_core)
+ deallocate(iboolfaces_outer_core)
+ case( IREGION_INNER_CORE )
+ ! inner core
+ deallocate(iboolcorner_inner_core)
+ deallocate(iboolleft_xi_inner_core,iboolright_xi_inner_core)
+ deallocate(iboolleft_eta_inner_core,iboolright_eta_inner_core)
+ deallocate(iboolfaces_inner_core)
+ end select
+
! synchronizes MPI processes
call sync_all()
@@ -364,6 +393,7 @@
use create_MPI_interfaces_par
use MPI_inner_core_par
+
implicit none
integer :: MAX_NEIGHBOURS,max_nibool
@@ -448,6 +478,10 @@
NGLOB_INNER_CORE, &
test_flag,ndim_assemble, &
iproc_eta,addressing,NCHUNKS,NPROC_XI,NPROC_ETA)
+
+ ! frees array not needed anymore
+ deallocate(ibelm_bottom_inner_core)
+
endif
! removes own myrank id (+1)
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-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_color_perm.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -607,10 +607,10 @@
implicit none
- integer :: myrank,nspec,nglob
+ integer,intent(in) :: myrank,nspec,nglob
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
- integer :: idomain
+ integer,intent(in) :: idomain
integer, dimension(nspec),intent(inout) :: perm
integer :: num_colors_outer,num_colors_inner
@@ -782,28 +782,12 @@
call permute_elements_real(gammaxstore,temp_array_real,perm,nspec)
call permute_elements_real(gammaystore,temp_array_real,perm,nspec)
call permute_elements_real(gammazstore,temp_array_real,perm,nspec)
+
! material parameters
call permute_elements_real(rhostore,temp_array_real,perm,nspec)
call permute_elements_real(kappavstore,temp_array_real,perm,nspec)
deallocate(temp_array_real)
- ! attenuation arrays
- if (ATTENUATION) then
- if (ATTENUATION_3D) then
- allocate(temp_array_dble(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(temp_array_dble_sls(N_SLS,NGLLX,NGLLY,NGLLZ,nspec))
- call permute_elements_dble(Qmu_store,temp_array_dble,perm,nspec)
- call permute_elements_dble_sls(tau_e_store,temp_array_dble_sls,perm,nspec)
- deallocate(temp_array_dble,temp_array_dble_sls)
- else
- allocate(temp_array_dble1(1,1,1,nspec))
- allocate(temp_array_dble_sls1(N_SLS,1,1,1,nspec))
- call permute_elements_dble1(Qmu_store,temp_array_dble1,perm,nspec)
- call permute_elements_dble_sls1(tau_e_store,temp_array_dble_sls1,perm,nspec)
- deallocate(temp_array_dble1,temp_array_dble_sls1)
- endif
- endif
-
! boundary surfaces
! note: only arrays pointing to ispec will have to be permutated since value of ispec will be different
!
@@ -844,11 +828,28 @@
ibelm_top(iface) = new_ispec
enddo
+ ! attenuation arrays
+ if (ATTENUATION) then
+ if (ATTENUATION_3D) then
+ allocate(temp_array_dble(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(temp_array_dble_sls(N_SLS,NGLLX,NGLLY,NGLLZ,nspec))
+ call permute_elements_dble(Qmu_store,temp_array_dble,perm,nspec)
+ call permute_elements_dble_sls(tau_e_store,temp_array_dble_sls,perm,nspec)
+ deallocate(temp_array_dble,temp_array_dble_sls)
+ else
+ allocate(temp_array_dble1(1,1,1,nspec))
+ allocate(temp_array_dble_sls1(N_SLS,1,1,1,nspec))
+ call permute_elements_dble1(Qmu_store,temp_array_dble1,perm,nspec)
+ call permute_elements_dble_sls1(tau_e_store,temp_array_dble_sls1,perm,nspec)
+ deallocate(temp_array_dble1,temp_array_dble_sls1)
+ endif
+ endif
select case( idomain )
case( IREGION_CRUST_MANTLE )
- ! region
- nspec = NSPEC_CRUST_MANTLE
+ ! checks number of elements
+ if( nspec /= NSPEC_CRUST_MANTLE ) &
+ call exit_MPI(myrank,'error in permutation nspec should be NSPEC_CRUST_MANTLE')
allocate(temp_array_real(NGLLX,NGLLY,NGLLZ,nspec))
@@ -934,8 +935,9 @@
endif
case( IREGION_OUTER_CORE )
- ! region
- nspec = NSPEC_OUTER_CORE
+ ! checks number of elements
+ if( nspec /= NSPEC_OUTER_CORE ) &
+ call exit_MPI(myrank,'error in permutation nspec should be NSPEC_OUTER_CORE')
if(ABSORBING_CONDITIONS .and. NCHUNKS /= 6 ) then
allocate(temp_array_real(NGLLX,NGLLY,NGLLZ,nspec))
@@ -946,8 +948,9 @@
endif
case( IREGION_INNER_CORE )
- ! region
- nspec = NSPEC_INNER_CORE
+ ! checks number of elements
+ if( nspec /= NSPEC_INNER_CORE ) &
+ call exit_MPI(myrank,'error in permutation nspec should be NSPEC_INNER_CORE')
allocate(temp_array_real(NGLLX,NGLLY,NGLLZ,nspec))
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-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_model.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -183,4 +183,7 @@
write(IMAIN,*)
write(IMAIN,*) 'Central cube is at a radius of ',R_CENTRAL_CUBE/1000.d0,' km'
+ ! flushes I/O buffer
+ call flush_IMAIN()
+
end subroutine sm_output_info
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/exit_mpi.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/exit_mpi.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/exit_mpi.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -105,7 +105,32 @@
end subroutine exit_MPI_without_rank
+!-------------------------------------------------------------------------------------------------
+!
+! I/O wrapper function
+!
+!-------------------------------------------------------------------------------------------------
+ subroutine flush_IMAIN()
+
+ implicit none
+
+ include "constants.h"
+
+ ! only master process writes out to main output file
+ ! file I/O in fortran is buffered by default
+ !
+ ! note: Fortran2003 includes a FLUSH statement
+ ! which is implemented by most compilers by now
+ !
+ ! otherwise:
+ ! a) comment out the line below
+ ! b) try to use instead: call flush(IMAIN)
+
+ flush(IMAIN)
+
+ end subroutine flush_IMAIN
+
!-------------------------------------------------------------------------------------------------
!
! MPI wrapper functions
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_simulation_stability.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_simulation_stability.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_simulation_stability.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -294,7 +294,10 @@
endif
write(IMAIN,*)
-
+
+ ! flushes file buffer for main output file (IMAIN)
+ call flush_IMAIN()
+
! write time stamp file to give information about progression of simulation
write(outputname,"('/timestamp',i6.6)") it
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_coupling.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -65,18 +65,18 @@
! for surface elements exactly on the CMB
do ispec2D = 1,nspec_top !NSPEC2D_TOP(IREGION_OUTER_CORE)
+
ispec = ibelm_top_outer_core(ispec2D)
+ ispec_selected = ibelm_bottom_crust_mantle(ispec2D)
! only for DOFs exactly on the CMB (top of these elements)
- k = NGLLZ
+ k = NGLLZ
+ ! get displacement on the solid side using pointwise matching
+ k_corresp = 1
+
do j = 1,NGLLY
do i = 1,NGLLX
-
- ! get displacement on the solid side using pointwise matching
- ispec_selected = ibelm_bottom_crust_mantle(ispec2D)
-
! corresponding points are located at the bottom of the mantle
- k_corresp = 1
iglob_cm = ibool_crust_mantle(i,j,k_corresp,ispec_selected)
displ_x = displ_crust_mantle(1,iglob_cm)
@@ -164,18 +164,19 @@
! for surface elements exactly on the ICB
do ispec2D = 1, nspec_bottom ! NSPEC2D_BOTTOM(IREGION_OUTER_CORE)
+
ispec = ibelm_bottom_outer_core(ispec2D)
+ ispec_selected = ibelm_top_inner_core(ispec2D)
! only for DOFs exactly on the ICB (bottom of these elements)
k = 1
+ ! get displacement on the solid side using pointwise matching
+ k_corresp = NGLLZ
+
do j = 1,NGLLY
do i = 1,NGLLX
- ! get displacement on the solid side using pointwise matching
- ispec_selected = ibelm_top_inner_core(ispec2D)
-
! corresponding points are located at the bottom of the mantle
- k_corresp = NGLLZ
iglob_ic = ibool_inner_core(i,j,k_corresp,ispec_selected)
displ_x = displ_inner_core(1,iglob_ic)
@@ -272,16 +273,15 @@
do ispec2D = 1,nspec_bottom ! NSPEC2D_BOTTOM(IREGION_CRUST_MANTLE)
ispec = ibelm_bottom_crust_mantle(ispec2D)
+ ispec_selected = ibelm_top_outer_core(ispec2D)
! only for DOFs exactly on the CMB (bottom of these elements)
k = 1
+ ! get velocity potential on the fluid side using pointwise matching
+ k_corresp = NGLLZ
+
do j = 1,NGLLY
do i = 1,NGLLX
-
- ! get velocity potential on the fluid side using pointwise matching
- ispec_selected = ibelm_top_outer_core(ispec2D)
- k_corresp = NGLLZ
-
! get normal at the CMB
nx = normal_top_outer_core(1,i,j,ispec2D)
ny = normal_top_outer_core(2,i,j,ispec2D)
@@ -378,16 +378,15 @@
do ispec2D = 1,nspec_top ! NSPEC2D_TOP(IREGION_INNER_CORE)
ispec = ibelm_top_inner_core(ispec2D)
+ ispec_selected = ibelm_bottom_outer_core(ispec2D)
! only for DOFs exactly on the ICB (top of these elements)
k = NGLLZ
+ ! get velocity potential on the fluid side using pointwise matching
+ k_corresp = 1
+
do j = 1,NGLLY
do i = 1,NGLLX
-
- ! get velocity potential on the fluid side using pointwise matching
- ispec_selected = ibelm_bottom_outer_core(ispec2D)
- k_corresp = 1
-
! get normal at the ICB
nx = normal_bottom_outer_core(1,i,j,ispec2D)
ny = normal_bottom_outer_core(2,i,j,ispec2D)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_stacey_outer_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_stacey_outer_core.f90 2012-08-13 00:58:02 UTC (rev 20566)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_stacey_outer_core.f90 2012-08-14 03:23:41 UTC (rev 20567)
@@ -225,6 +225,7 @@
endif
! ymax
+
if (SIMULATION_TYPE == 3 .and. nspec2D_ymax_outer_core > 0) then
call read_abs(7,absorb_ymax_outer_core,reclen_ymax_outer_core,NSTEP-it+1)
endif
@@ -271,6 +272,8 @@
call write_abs(7,absorb_ymax_outer_core,reclen_ymax_outer_core,it)
endif
+ ! zmin
+
! for surface elements exactly on the ICB
if (SIMULATION_TYPE == 3 .and. nspec2D_zmin_outer_core > 0) then
call read_abs(8,absorb_zmin_outer_core,reclen_zmin,NSTEP-it+1)
More information about the CIG-COMMITS
mailing list