[cig-commits] r22474 - in seismo/3D/SPECFEM3D_GLOBE: branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D trunk/src/auxiliaries trunk/src/create_header_file trunk/src/meshfem3D trunk/src/shared trunk/src/specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Jun 30 21:15:01 PDT 2013
Author: dkomati1
Date: 2013-06-30 21:15:01 -0700 (Sun, 30 Jun 2013)
New Revision: 22474
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_1D.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_2D.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_corners_chunks.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_faces_chunks.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_AVS_DX.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_GMT_global.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_1D_buffers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_cutplanes_eta.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_global.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_models.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/write_AVS_DX_global_chunks_data.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/create_name_database.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_stacey_crust_mantle.f90
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/initialize_simulation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_surf_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_1D_buffers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_eta.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_global.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_1066a.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/sort_array_coordinates.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/stretching_function.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/write_AVS_DX_global_chunks_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/create_name_database.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_crust_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_outer_core.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/initialize_simulation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/multiply_arrays_source.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90
Log:
more obvious merges
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_1D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_1D.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_1D.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -94,7 +94,6 @@
NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
nglob
-! computed in read_compute_parameters
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: r_bottom,r_top
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_2D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_2D.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_2D.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -84,7 +84,6 @@
ROTATE_SEISMOGRAMS_RT,HONOR_1D_SPHERICAL_MOHO,WRITE_SEISMOGRAMS_BY_MASTER,&
SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE
-! computed in read_compute_parameters
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: r_bottom,r_top
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_corners_chunks.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_corners_chunks.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_corners_chunks.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -93,7 +93,6 @@
character(len=150) filename,prname
-! computed in read_compute_parameters
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: r_bottom,r_top
@@ -142,7 +141,7 @@
ratio_sampling_array, ner, doubling_index,r_bottom,r_top,this_region_has_a_doubling,rmins,rmaxs,CASE_3D, &
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
- DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
+ DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.,NOISE_TOMOGRAPHY)
print *
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_faces_chunks.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_faces_chunks.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_faces_chunks.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -97,7 +97,6 @@
character(len=150) filename,prname
-! computed in read_compute_parameters
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: r_bottom,r_top
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_AVS_DX.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_AVS_DX.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -170,14 +170,12 @@
integer proc_p1,proc_p2
-! computed in read_compute_parameters
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: r_bottom,r_top
logical, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: this_region_has_a_doubling
double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: rmins,rmaxs
-
logical :: CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA
integer, dimension(NB_SQUARE_CORNERS,NB_CUT_CASE) :: DIFF_NSPEC1D_RADIAL
integer, dimension(NB_SQUARE_EDGES_ONEDIR,NB_CUT_CASE) :: DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_GMT_global.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_GMT_global.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_GMT_global.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -621,7 +621,6 @@
endif
-
enddo !i
enddo !j
@@ -708,7 +707,6 @@
if( thetaval <= 0.01 ) istamp2 = ieoff
endif
-
enddo !i
enddo !j
@@ -941,7 +939,6 @@
close(11)
if(iframe == 1) close(12)
-
! end of loop and test on all the time steps for all the movie images
endif
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_1D_buffers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_1D_buffers.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_1D_buffers.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -230,6 +230,7 @@
! nb of elements in this 1D buffer
ispeccount=0
+
do ispec=1,nspec
! remove central cube for chunk buffers
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_cutplanes_eta.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_cutplanes_eta.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_MPI_cutplanes_eta.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -103,6 +103,7 @@
! nb of elements in this cut-plane
ispecc1=0
+
do ispec=1,nspec
if(iMPIcut_eta(1,ispec)) then
ispecc1=ispecc1+1
@@ -169,6 +170,7 @@
! nb of elements in this cut-plane
ispecc2=0
+
do ispec=1,nspec
if(iMPIcut_eta(2,ispec)) then
ispecc2=ispecc2+1
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_global.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_global.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_global.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -25,7 +25,6 @@
!
!=====================================================================
-
subroutine get_global(nspec,xp,yp,zp,iglob,loc,ifseg,nglob,npointot)
! this routine MUST be in double precision to avoid sensitivity
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 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/meshfem3D_models.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -25,7 +25,6 @@
!
!=====================================================================
-
subroutine meshfem3D_models_broadcast(myrank,NSPEC, &
MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,&
R80,R220,R670,RCMB,RICB, &
@@ -661,7 +660,6 @@
lon = phi * RADIANS_TO_DEGREES
if( lon > 180.0d0 ) lon = lon - 360.0d0
-
!---
!
! ADD YOUR MODEL HERE
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 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -238,7 +238,7 @@
! mass matrices
!
! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
- ! on Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
!
! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/write_AVS_DX_global_chunks_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/write_AVS_DX_global_chunks_data.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/write_AVS_DX_global_chunks_data.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -86,7 +86,6 @@
double precision x,y,z,theta,phi_dummy,cost,p20,ell,factor
real(kind=CUSTOM_REAL) dvp,dvs
-
! writing points
open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXpointschunks.txt',status='unknown')
open(unit=11,file=prname(1:len_trim(prname))//'AVS_DXpointschunks_stability.txt',status='unknown')
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/create_name_database.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/create_name_database.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/create_name_database.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -25,11 +25,12 @@
!
!=====================================================================
-subroutine create_name_database(prname,iproc,iregion_code,LOCAL_PATH)
+ subroutine create_name_database(prname,iproc,iregion_code,LOCAL_PATH)
! create the name of the database for the mesher and the solver
implicit none
+
integer iproc,iregion_code
! name of the database file
@@ -37,10 +38,11 @@
! create the name for the database of the current slide and region
write(procname,"('/proc',i6.6,'_reg',i1,'_')") iproc,iregion_code
+
! create full name with path
prname = trim(LOCAL_PATH) // procname
-end subroutine create_name_database
+ end subroutine create_name_database
subroutine create_name_database_adios(prname,iregion_code,LOCAL_PATH)
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 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -25,15 +25,6 @@
!
!=====================================================================
-! preprocessing definition: #define _HANDOPT : turns hand-optimized code on
-! #undef _HANDOPT : turns hand-optimized code off
-! or compile with: -D_HANDOPT
-!#define _HANDOPT
-
-! note: these hand optimizations 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%.
-! the original routines are commented with "! way 1", the hand-optimized routines with "! way 2"
-
subroutine compute_element_iso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
xstore,ystore,zstore, &
@@ -402,17 +393,17 @@
endif ! end of section with gravity terms
! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl)
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl)
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl)
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl)
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl)
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl)
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
enddo ! NGLLX
enddo ! NGLLY
@@ -426,9 +417,6 @@
!--------------------------------------------------------------------------------------------------
!
-
-
-
subroutine compute_element_tiso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
xstore,ystore,zstore, &
@@ -1008,17 +996,18 @@
endif ! end of section with gravity terms
! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl)
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl)
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl)
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl)
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl)
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl)
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
+
enddo ! NGLLX
enddo ! NGLLY
enddo ! NGLLZ
@@ -1030,7 +1019,6 @@
!--------------------------------------------------------------------------------------------
!
-
subroutine compute_element_aniso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
xstore,ystore,zstore, &
@@ -1440,17 +1428,18 @@
endif ! end of section with gravity terms
! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl)
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl)
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl)
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl)
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl)
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl)
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
+
enddo ! NGLLX
enddo ! NGLLY
enddo ! NGLLZ
@@ -1566,7 +1555,6 @@
end subroutine compute_element_att_stress
-
!
!--------------------------------------------------------------------------------------------
!
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_stacey_crust_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_stacey_crust_mantle.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_stacey_crust_mantle.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -73,6 +73,7 @@
! crust & mantle
+
! xmin
! if two chunks exclude this face for one of them
if(NCHUNKS_VAL == 1 .or. ichunk == CHUNK_AC) then
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 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_stacey_outer_core.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -25,7 +25,6 @@
!
!=====================================================================
-
subroutine compute_stacey_outer_core()
use constants_solver
@@ -170,8 +169,6 @@
5) ! <= xmax
endif
-
-
if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. nspec2D_xmax_outer_core > 0 ) then
call write_abs(5,absorb_xmax_outer_core,reclen_xmax_outer_core,it)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/initialize_simulation.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/initialize_simulation.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -276,6 +276,7 @@
else
write(IMAIN,*) ' no general mantle anisotropy'
endif
+
write(IMAIN,*)
write(IMAIN,*)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -50,6 +50,7 @@
! standard include of the MPI library
include 'mpif.h'
+
include "precision.h"
integer nspec,nglob
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -483,7 +483,6 @@
OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY,MODEL, &
NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current, &
iorientation,phi,chn,sisname)
-
endif ! OUTPUT_SEISMOS_SAC_ALPHANUM .or. OUTPUT_SEISMOS_SAC_BINARY
! ASCII output format
@@ -493,7 +492,6 @@
NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current, &
SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,myrank, &
iorientation,sisname,sisname_big_file)
-
endif ! OUTPUT_SEISMOS_ASCII_TEXT
enddo ! do iorientation
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_surf_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_surf_data.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_surf_data.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -192,7 +192,7 @@
write(em_mesh_file,'(a,i1,a)') trim(outdir)//'/'//trim(filename)//'_element.surf'
command_name='rm -f '//trim(pt_mesh_file)//' '//trim(em_mesh_file)//' '//trim(mesh_file)
- !call system(trim(command_name))
+ call system(trim(command_name))
call open_file_fd(trim(pt_mesh_file)//char(0),pfd)
call open_file_fd(trim(em_mesh_file)//char(0),efd)
@@ -340,7 +340,7 @@
print *, ' '
print *, 'cat mesh files ...'
print *, trim(command_name)
- !call system(trim(command_name))
+ call system(trim(command_name))
print *, 'Done writing '//trim(mesh_file)
print *, ' '
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -465,7 +465,7 @@
print *, ' '
print *, 'cat mesh files: '
print *, trim(command_name)
- !call system(trim(command_name))
+ call system(trim(command_name))
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -397,15 +397,15 @@
! converts values into radians
! colatitude [0, PI]
- LAT_SOURCE = (90. - LAT_SOURCE)*PI/180.0
+ LAT_SOURCE = (90.0 - LAT_SOURCE)*DEGREES_TO_RADIANS
! longitude [-PI, PI]
if( LON_SOURCE < -180.0 ) LON_SOURCE = LON_SOURCE + 360.0
if( LON_SOURCE > 180.0 ) LON_SOURCE = LON_SOURCE - 360.0
- LON_SOURCE = LON_SOURCE *PI/180.0
+ LON_SOURCE = LON_SOURCE * DEGREES_TO_RADIANS
! mute radius in rad
- RADIUS_TO_MUTE = RADIUS_TO_MUTE*PI/180.0
+ RADIUS_TO_MUTE = RADIUS_TO_MUTE * DEGREES_TO_RADIANS
endif
print *,'--------'
@@ -458,7 +458,7 @@
! approximate wavefront travel distance in degrees
! (~3.5 km/s wave speed for surface waves)
- distance = 3.5 * ((it-1)*DT-t0) / 6371.0 * 180./PI
+ distance = 3.5 * ((it-1)*DT-t0) / 6371.0 * RADIANS_TO_DEGREES
! approximate distance to source (in degrees)
! (shrinks if waves travel back from antipode)
@@ -479,7 +479,7 @@
print*,'muting radius: ',0.7 * distance,'(degrees)'
! new radius of mute area (in rad)
- RADIUS_TO_MUTE = 0.7 * distance * PI/180.
+ RADIUS_TO_MUTE = 0.7 * distance * DEGREES_TO_RADIANS
else
! mute_factor used at the beginning for scaling displacement values
if( STARTTIME_TO_MUTE > TINYVAL ) then
@@ -588,17 +588,17 @@
! checks source longitude range
if( LON_SOURCE - RADIUS_TO_MUTE < -PI .or. LON_SOURCE + RADIUS_TO_MUTE > PI ) then
! source close to 180. longitudes, shifts range to [0, 2PI]
- if( phival < 0.0 ) phival = phival + 2.0*PI
+ if( phival < 0.0 ) phival = phival + TWO_PI
if( LON_SOURCE < 0.0 ) then
- dist_lon = phival - (LON_SOURCE + 2.0*PI)
+ dist_lon = phival - (LON_SOURCE + TWO_PI)
else
dist_lon = phival - LON_SOURCE
endif
else
! source well between range to [-PI, PI]
! shifts phival to be like LON_SOURCE between [-PI,PI]
- if( phival > PI ) phival = phival - 2.0*PI
- if( phival < -PI ) phival = phival + 2.0*PI
+ if( phival > PI ) phival = phival - TWO_PI
+ if( phival < -PI ) phival = phival + TWO_PI
dist_lon = phival - LON_SOURCE
endif
@@ -623,7 +623,6 @@
endif
-
enddo !i
enddo !j
@@ -710,7 +709,6 @@
if( thetaval <= 0.01 ) istamp2 = ieoff
endif
-
enddo !i
enddo !j
@@ -767,7 +765,7 @@
if( max_absol < max_average ) then
! distance (in degree) of surface waves travelled
- distance = 3.5 * ((it-1)*DT-t0) / 6371.0 * 180./PI
+ distance = 3.5 * ((it-1)*DT-t0) / 6371.0 * RADIANS_TO_DEGREES
if( distance > 10.0 .and. distance <= 20.0 ) then
! smooth transition between 10 and 20 degrees
! sets positive and negative maximum
@@ -914,13 +912,13 @@
! converts the geocentric colatitude to a geographic colatitude
if(.not. ASSUME_PERFECT_SPHERE) then
- thetaval = PI/2.0d0 - &
+ thetaval = PI_OVER_TWO - &
datan(1.006760466d0*dcos(dble(thetaval))/dmax1(TINYVAL,dble(sin(thetaval))))
endif
! gets geographic latitude and longitude in degrees
- lat = sngl(90.d0 - thetaval*180.0/PI)
- long = sngl(phival*180.0/PI)
+ lat = sngl(90.d0 - thetaval*RADIANS_TO_DEGREES)
+ long = sngl(phival*RADIANS_TO_DEGREES)
if(long > 180.0) long = long-360.0
endif
@@ -943,7 +941,6 @@
close(11)
if(iframe == 1) close(12)
-
! end of loop and test on all the time steps for all the movie images
endif
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -184,7 +184,7 @@
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,&
- static_memory_size,NGLOB1D_RADIAL_TEMP,&
+ static_memory_size,NGLOB1D_RADIAL_TEMP, &
NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NSPEC2D_TOP,NSPEC2D_BOTTOM, &
NSPEC2DMAX_YMIN_YMAX,NSPEC2DMAX_XMIN_XMAX, &
NPROC_XI,NPROC_ETA, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_1D_buffers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_1D_buffers.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_1D_buffers.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -89,7 +89,8 @@
idoubling(ispec) == IFLAG_BOTTOM_CENTRAL_CUBE .or. &
idoubling(ispec) == IFLAG_TOP_CENTRAL_CUBE .or. &
idoubling(ispec) == IFLAG_IN_FICTITIOUS_CUBE) cycle
- ! corner detection here
+
+ ! corner detection here
if(iMPIcut_xi(1,ispec) .and. iMPIcut_eta(1,ispec)) then
ispeccount=ispeccount+1
! loop on all the points
@@ -198,13 +199,12 @@
! corner detection here
if(iMPIcut_xi(1,ispec) .and. iMPIcut_eta(2,ispec)) then
- ispeccount=ispeccount+1
+ ispeccount=ispeccount+1
-! loop on all the points
- ix = 1
- iy = NGLLY
- do iz=1,NGLLZ
-
+ ! loop on all the points
+ ix = 1
+ iy = NGLLY
+ do iz=1,NGLLZ
! select point, if not already selected
if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_eta.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_eta.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_eta.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -84,7 +84,7 @@
! nb of global points shared with the other slice
npoin2D_eta = 0
-! nb of elements in this cut-plane
+ ! nb of elements in this cut-plane
ispecc1=0
do ispec=1,nspec
@@ -131,7 +131,7 @@
! nb of global points shared with the other slice
npoin2D_eta = 0
-! nb of elements in this cut-plane
+ ! nb of elements in this cut-plane
ispecc2=0
do ispec=1,nspec
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_global.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_global.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_global.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -38,7 +38,7 @@
include "constants.h"
-! parameters
+ ! input parameters
integer, intent(in) :: npointot,nspec
double precision, intent(in) :: xp(npointot),yp(npointot),zp(npointot)
@@ -67,57 +67,58 @@
enddo
enddo
- ifseg(:)=.false.
+ ifseg(:) = .false.
- nseg=1
- ifseg(1)=.true.
- ninseg(1)=npointot
+ nseg = 1
+ ifseg(1) = .true.
+ ninseg(1) = npointot
-do j=1,NDIM
-
+ do j=1,NDIM
! sort within each segment
ioff=1
do iseg=1,nseg
- if(j == 1) then
- call rank(xp(ioff),ind,ninseg(iseg))
- else if(j == 2) then
- call rank(yp(ioff),ind,ninseg(iseg))
- else
- call rank(zp(ioff),ind,ninseg(iseg))
- endif
- call swap_all(loc(ioff),xp(ioff),yp(ioff),zp(ioff),iwork,work,ind,ninseg(iseg))
- ioff=ioff+ninseg(iseg)
+ if(j == 1) then
+ call rank(xp(ioff),ind,ninseg(iseg))
+ else if(j == 2) then
+ call rank(yp(ioff),ind,ninseg(iseg))
+ else
+ call rank(zp(ioff),ind,ninseg(iseg))
+ endif
+
+ call swap_all(loc(ioff),xp(ioff),yp(ioff),zp(ioff),iwork,work,ind,ninseg(iseg))
+
+ ioff=ioff+ninseg(iseg)
enddo
-! check for jumps in current coordinate
-! compare the coordinates of the points within a small tolerance
+ ! check for jumps in current coordinate
+ ! compare the coordinates of the points within a small tolerance
if(j == 1) then
- do i=2,npointot
- if(dabs(xp(i)-xp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
- enddo
+ do i=2,npointot
+ if(dabs(xp(i)-xp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
+ enddo
else if(j == 2) then
- do i=2,npointot
- if(dabs(yp(i)-yp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
- enddo
+ do i=2,npointot
+ if(dabs(yp(i)-yp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
+ enddo
else
- do i=2,npointot
- if(dabs(zp(i)-zp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
- enddo
+ do i=2,npointot
+ if(dabs(zp(i)-zp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
+ enddo
endif
-! count up number of different segments
+ ! count up number of different segments
nseg=0
do i=1,npointot
- if(ifseg(i)) then
+ if(ifseg(i)) then
nseg=nseg+1
ninseg(nseg)=1
- else
+ else
ninseg(nseg)=ninseg(nseg)+1
- endif
+ endif
enddo
-enddo
+ enddo
-! assign global node numbers (now sorted lexicographically)
+ ! assign global node numbers (now sorted lexicographically)
ig=0
do i=1,npointot
if(ifseg(i)) ig=ig+1
@@ -126,11 +127,8 @@
nglob=ig
-! deallocate arrays
- deallocate(ind)
- deallocate(ninseg)
- deallocate(iwork)
- deallocate(work)
+ ! deallocate arrays
+ deallocate(ind,ninseg,iwork,work)
end subroutine get_global
@@ -149,9 +147,10 @@
include "constants.h"
- integer :: nspec,nglob
+ integer,intent(in) :: nspec,nglob
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ ! local parameters
! mask to sort ibool
integer, dimension(:), allocatable :: mask_ibool
integer, dimension(:,:,:,:), allocatable :: copy_ibool_ori
@@ -275,20 +274,20 @@
W(:) = A(:)
do i=1,n
- IA(i)=IW(ind(i))
- A(i)=W(ind(i))
+ IA(i) = IW(ind(i))
+ A(i) = W(ind(i))
enddo
W(:) = B(:)
do i=1,n
- B(i)=W(ind(i))
+ B(i) = W(ind(i))
enddo
W(:) = C(:)
do i=1,n
- C(i)=W(ind(i))
+ C(i) = W(ind(i))
enddo
end subroutine swap_all
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -25,7 +25,6 @@
!
!=====================================================================
-
module meshfem3D_models_par
!---
@@ -632,6 +631,7 @@
end subroutine meshfem3D_crust_broadcast
+
!
!-------------------------------------------------------------------------------------------------
!
@@ -773,7 +773,7 @@
double precision xmesh,ymesh,zmesh,r
! the 21 coefficients for an anisotropic medium in reduced notation
- double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33, &
+ double precision :: c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33, &
c34,c35,c36,c44,c45,c46,c55,c56,c66
! local parameters
@@ -1049,17 +1049,17 @@
implicit none
- integer iregion_code
+ integer :: iregion_code
! note: r is the exact radius (and not r_prem with tolerance)
- double precision xmesh,ymesh,zmesh,r
- double precision vpv,vph,vsv,vsh,rho,eta_aniso,dvp
+ double precision :: xmesh,ymesh,zmesh,r
+ double precision :: vpv,vph,vsv,vsh,rho,eta_aniso,dvp
! the 21 coefficients for an anisotropic medium in reduced notation
- double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33, &
+ double precision :: c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33, &
c34,c35,c36,c44,c45,c46,c55,c56,c66
- logical elem_in_crust
- double precision moho
+ logical :: elem_in_crust
+ double precision :: moho
! local parameters
double precision :: r_dummy,theta,phi
@@ -1075,10 +1075,11 @@
! gets point's position theta/phi, lat/lon
call xyz_2_rthetaphi_dble(xmesh,ymesh,zmesh,r_dummy,theta,phi)
call reduce(theta,phi)
- lat = (PI/2.0d0-theta)*180.0d0/PI
- lon = phi*180.0d0/PI
- if(lon>180.0d0) lon = lon-360.0d0
+ lat = (PI_OVER_TWO - theta) * RADIANS_TO_DEGREES
+ lon = phi * RADIANS_TO_DEGREES
+ if( lon > 180.0d0 ) lon = lon - 360.0d0
+
!---
!
! ADD YOUR MODEL HERE
@@ -1247,16 +1248,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
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_1066a.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_1066a.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_1066a.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -96,15 +96,14 @@
! compressional wave speed vp: km/s
! shear wave speed vs: km/s
- integer iregion_code
+ double precision :: x,rho,vp,vs,Qmu,Qkappa
+ integer :: iregion_code
- double precision x,rho,vp,vs,Qmu,Qkappa
+ ! local parameters
+ double precision :: r,frac,scaleval
+ integer :: i
- integer i
-
- double precision r,frac,scaleval
-
-! compute real physical radius in meters
+ ! compute real physical radius in meters
r = x * R_EARTH
i = 1
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -187,9 +187,12 @@
! other terms needed in the solid regions only
if(iregion_code /= IREGION_OUTER_CORE) then
- if(.not. (ANISOTROPIC_3D_MANTLE .and. iregion_code == IREGION_CRUST_MANTLE)) write(27) muvstore
+ ! note: muvstore needed for Q_mu shear attenuation in inner core
+ if(.not. (ANISOTROPIC_3D_MANTLE .and. iregion_code == IREGION_CRUST_MANTLE)) then
+ write(27) muvstore
+ endif
-! save anisotropy in the mantle only
+ ! save anisotropy in the mantle only
if(TRANSVERSE_ISOTROPY) then
if(iregion_code == IREGION_CRUST_MANTLE .and. .not. ANISOTROPIC_3D_MANTLE) then
write(27) kappahstore
@@ -198,7 +201,7 @@
endif
endif
-! save anisotropy in the inner core only
+ ! save anisotropy in the inner core only
if(ANISOTROPIC_INNER_CORE .and. iregion_code == IREGION_INNER_CORE) then
write(27) c11store
write(27) c33store
@@ -207,8 +210,6 @@
write(27) c44store
endif
-
-
if(ANISOTROPIC_3D_MANTLE .and. iregion_code == IREGION_CRUST_MANTLE) then
write(27) c11store
write(27) c12store
@@ -235,7 +236,7 @@
endif
-! Stacey
+ ! Stacey
if(ABSORBING_CONDITIONS) then
if(iregion_code == IREGION_CRUST_MANTLE) then
@@ -247,14 +248,14 @@
endif
-! mass matrices
-!
-! in the case of stacey boundary conditions, add C*delta/2 contribution to the mass matrix
-! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
-! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
-!
-! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
-! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+ ! mass matrices
+ !
+ ! in the case of stacey boundary conditions, add C*deltat/2 contribution to the mass matrix
+ ! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
+ ! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
+ ! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS .and. iregion_code == IREGION_CRUST_MANTLE) then
write(27) rmassx
write(27) rmassy
@@ -262,10 +263,10 @@
write(27) rmassz
-! additional ocean load mass matrix if oceans and if we are in the crust
+ ! additional ocean load mass matrix if oceans and if we are in the crust
if(OCEANS .and. iregion_code == IREGION_CRUST_MANTLE) write(27) rmass_ocean_load
- close(27)
+ close(27) ! solver_data.bin
open(unit=27,file=prname(1:len_trim(prname))//'solver_data_2.bin',status='unknown',form='unformatted',action='write')
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/sort_array_coordinates.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/sort_array_coordinates.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/sort_array_coordinates.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -117,7 +117,7 @@
enddo
enddo
-! assign global node numbers (now sorted lexicographically)
+ ! assign global node numbers (now sorted lexicographically)
ig=0
do i=1,npointot
if(ifseg(i)) ig=ig+1
@@ -138,12 +138,13 @@
!
implicit none
- integer n
- double precision A(n)
- integer IND(n)
+ integer :: n
+ double precision,dimension(n) :: A
+ integer,dimension(n) :: IND
- integer i,j,l,ir,indx
- double precision q
+ ! local parameters
+ integer :: i,j,l,ir,indx
+ double precision :: q
do j=1,n
IND(j)=j
@@ -164,7 +165,7 @@
ind(ir)=ind(1)
ir=ir-1
if (ir == 1) then
- ind(1)=indx
+ ind(1) = indx
return
endif
ENDIF
@@ -196,14 +197,14 @@
!
implicit none
- integer n
+ integer :: n
+ integer,dimension(n) :: IND
+ integer,dimension(n) :: IA,IB,IW
+ double precision,dimension(n) :: A,B,C,W
- integer IND(n)
- integer IA(n),IB(n),IW(n)
- double precision A(n),B(n),C(n),W(n)
+ ! local parameter
+ integer :: i
- integer i
-
do i=1,n
W(i)=A(i)
IW(i)=IA(i)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/stretching_function.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/stretching_function.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/stretching_function.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -25,7 +25,7 @@
!
!=====================================================================
-subroutine stretching_function(r_top,r_bottom,ner,stretch_tab)
+ subroutine stretching_function(r_top,r_bottom,ner,stretch_tab)
! define stretch_tab which contains r_top and r_bottom for each element layer in the crust for 3D models.
!
@@ -76,14 +76,14 @@
stretch_tab(2,i) = stretch_tab(1,i+1)
enddo
-end subroutine stretching_function
+ end subroutine stretching_function
!
!-------------------------------------------------------------------------------------------------
!
-subroutine stretching_function_regional(r_top,r_bottom,ner,stretch_tab)
+ subroutine stretching_function_regional(r_top,r_bottom,ner,stretch_tab)
! define stretch_tab which contains r_top and r_bottom for each element layer in the crust for 3D models.
!
@@ -144,6 +144,6 @@
stretch_tab(2,2) = 6336000.d0 ! bottom second layer
stretch_tab(2,3) = r_bottom ! bottom third layer
-end subroutine stretching_function_regional
+ end subroutine stretching_function_regional
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/write_AVS_DX_global_chunks_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/write_AVS_DX_global_chunks_data.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/write_AVS_DX_global_chunks_data.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -27,42 +27,53 @@
! create AVS or DX 2D data for the faces of the global chunks,
! to be recombined in postprocessing
- subroutine write_AVS_DX_global_chunks_data(myrank,prname,nspec,iboun, &
- ibool,idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool, &
- npointot,rhostore,kappavstore,muvstore,nspl,rspl,espl,espl2, &
- ELLIPTICITY,ISOTROPIC_3D_MANTLE, &
- RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R120,R80,RMOHO, &
- RMIDDLE_CRUST,ROCEAN,iregion_code)
+ subroutine write_AVS_DX_global_chunks_data(myrank,prname,nspec,iboun,ibool, &
+ idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool, &
+ npointot,rhostore,kappavstore,muvstore,nspl,rspl,espl,espl2, &
+ ELLIPTICITY,ISOTROPIC_3D_MANTLE, &
+ RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R120,R80,RMOHO, &
+ RMIDDLE_CRUST,ROCEAN,iregion_code)
implicit none
include "constants.h"
- integer nspec,myrank
- integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
+ integer :: myrank
+ ! processor identification
+ character(len=150) :: prname
+
+ integer :: nspec
+
+ logical iboun(6,nspec)
+
+ integer,dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
integer idoubling(nspec)
- logical iboun(6,nspec),ELLIPTICITY,ISOTROPIC_3D_MANTLE
+ double precision,dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
- double precision RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771, &
- R400,R120,R80,RMOHO,RMIDDLE_CRUST,ROCEAN
+ integer :: npointot
+ ! numbering of global AVS or DX points
+ integer num_ibool_AVS_DX(npointot)
+ ! logical mask used to output global points only once
+ logical mask_ibool(npointot)
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
-
real(kind=CUSTOM_REAL) kappavstore(NGLLX,NGLLY,NGLLZ,nspec)
real(kind=CUSTOM_REAL) muvstore(NGLLX,NGLLY,NGLLZ,nspec)
real(kind=CUSTOM_REAL) rhostore(NGLLX,NGLLY,NGLLZ,nspec)
-! logical mask used to output global points only once
- integer npointot
- logical mask_ibool(npointot)
+ ! for ellipticity
+ integer nspl
+ double precision rspl(NR),espl(NR),espl2(NR)
-! numbering of global AVS or DX points
- integer num_ibool_AVS_DX(npointot)
+ logical ELLIPTICITY,ISOTROPIC_3D_MANTLE
+ double precision RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771, &
+ R400,R120,R80,RMOHO,RMIDDLE_CRUST,ROCEAN
+
+ integer iregion_code
+
+ ! local parameters
integer ispec
integer i,j,k,np
integer, dimension(8) :: iglobval
@@ -75,16 +86,6 @@
double precision x,y,z,theta,phi_dummy,cost,p20,ell,factor
real(kind=CUSTOM_REAL) dvp,dvs
-! for ellipticity
- integer nspl
- double precision rspl(NR),espl(NR),espl2(NR)
-
-! processor identification
- character(len=150) prname
-
- integer iregion_code
-
-
! writing points
open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXpointschunks.txt',status='unknown')
open(unit=11,file=prname(1:len_trim(prname))//'AVS_DXpointschunks_stability.txt',status='unknown')
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/create_name_database.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/create_name_database.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/create_name_database.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -44,3 +44,21 @@
end subroutine create_name_database
+subroutine create_name_database_adios(prname,iregion_code,LOCAL_PATH)
+
+ ! create the name of the database for the mesher and the solver
+
+ implicit none
+
+ integer iregion_code
+
+! name of the database file
+ character(len=150) prname,procname,LOCAL_PATH
+
+! create the name for the database of the current slide and region
+ write(procname,"('/reg',i1,'_')") iregion_code
+
+! create full name with path
+ prname = trim(LOCAL_PATH) // procname
+
+end subroutine create_name_database_adios
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_add_sources.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -98,11 +98,11 @@
stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
endif
- ! we use a force in a single direction along one of the components:
- ! x/y/z or E/N/Z-direction would correspond to 1/2/3 = COMPONENT_FORCE_SOURCE
- ! e.g. nu_source(3,:) here would be a source normal to the surface (z-direction).
- accel_crust_mantle(:,iglob) = accel_crust_mantle(:,iglob) &
- + sngl( nu_source(COMPONENT_FORCE_SOURCE,:,isource) ) * stf_used
+ ! we use a force in a single direction along one of the components:
+ ! x/y/z or E/N/Z-direction would correspond to 1/2/3 = COMPONENT_FORCE_SOURCE
+ ! e.g. nu_source(3,:) here would be a source normal to the surface (z-direction).
+ accel_crust_mantle(:,iglob) = accel_crust_mantle(:,iglob) &
+ + sngl( nu_source(COMPONENT_FORCE_SOURCE,:,isource) ) * stf_used
else
if(USE_LDDRK)then
@@ -112,28 +112,28 @@
stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
endif
- ! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- stf_used = sngl(stf)
- else
- stf_used = stf
- endif
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
- ! add source array
- ispec = ispec_selected_source(isource)
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool_crust_mantle(i,j,k,ispec)
+ ! add source array
+ ispec = ispec_selected_source(isource)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool_crust_mantle(i,j,k,ispec)
- accel_crust_mantle(:,iglob) = accel_crust_mantle(:,iglob) &
- + sourcearrays(:,i,j,k,isource)*stf_used
+ accel_crust_mantle(:,iglob) = accel_crust_mantle(:,iglob) &
+ + sourcearrays(:,i,j,k,isource)*stf_used
enddo
enddo
enddo
- endif ! USE_FORCE_POINT_SOURCE
+ endif ! USE_FORCE_POINT_SOURCE
endif
@@ -243,15 +243,15 @@
irec_local = 0
do irec = 1,nrec
- ! adds source (only if this proc carries the source)
- if(myrank == islice_selected_rec(irec)) then
- irec_local = irec_local + 1
+ ! adds source (only if this proc carries the source)
+ if(myrank == islice_selected_rec(irec)) then
+ irec_local = irec_local + 1
- ! adds source contributions
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool_crust_mantle(i,j,k,ispec_selected_rec(irec))
+ ! adds source contributions
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool_crust_mantle(i,j,k,ispec_selected_rec(irec))
! adds adjoint source acting at this time step (it):
@@ -319,7 +319,7 @@
enddo
endif
- enddo
+ enddo
end subroutine compute_add_sources_adjoint
@@ -434,11 +434,11 @@
enddo
enddo
- endif ! USE_FORCE_POINT_SOURCE
+ endif ! USE_FORCE_POINT_SOURCE
- endif
+ endif
- enddo
+ enddo
end subroutine compute_add_sources_backward
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -131,9 +131,9 @@
gammazl = gammaz(i,j,k,ispec)
! compute the jacobian
- jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
- - xiyl*(etaxl*gammazl-etazl*gammaxl) &
- + xizl*(etaxl*gammayl-etayl*gammaxl))
+ jacobianl = 1.0_CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
duxdxl = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
duxdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
@@ -465,7 +465,7 @@
gammazl = gammaz(i,j,k,ispec)
! compute the jacobian
- jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ jacobianl = 1.0_CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
- xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ xizl*(etaxl*gammayl-etayl*gammaxl))
@@ -569,31 +569,31 @@
cosphifour = cosphisq * cosphisq
sinphifour = sinphisq * sinphisq
- costwotheta = cos(2.0*theta)
- sintwotheta = sin(2.0*theta)
- costwophi = cos(2.0*phi)
- sintwophi = sin(2.0*phi)
+ costwotheta = cos(2.0_CUSTOM_REAL*theta)
+ sintwotheta = sin(2.0_CUSTOM_REAL*theta)
+ costwophi = cos(2.0_CUSTOM_REAL*phi)
+ sintwophi = sin(2.0_CUSTOM_REAL*phi)
- cosfourtheta = cos(4.0*theta)
- cosfourphi = cos(4.0*phi)
+ cosfourtheta = cos(4.0_CUSTOM_REAL*theta)
+ cosfourphi = cos(4.0_CUSTOM_REAL*phi)
costwothetasq = costwotheta * costwotheta
costwophisq = costwophi * costwophi
sintwophisq = sintwophi * sintwophi
- etaminone = eta_aniso - 1.0
- twoetaminone = 2.0 * eta_aniso - 1.0
+ etaminone = eta_aniso - 1.0_CUSTOM_REAL
+ twoetaminone = 2.0_CUSTOM_REAL * eta_aniso - 1.0_CUSTOM_REAL
! precompute some products to reduce the CPU time
- two_eta_aniso = 2.0*eta_aniso
- four_eta_aniso = 4.0*eta_aniso
- six_eta_aniso = 6.0*eta_aniso
+ two_eta_aniso = 2.0_CUSTOM_REAL*eta_aniso
+ four_eta_aniso = 4.0_CUSTOM_REAL*eta_aniso
+ six_eta_aniso = 6.0_CUSTOM_REAL*eta_aniso
- two_rhovsvsq = 2.0*rhovsvsq
- two_rhovshsq = 2.0*rhovshsq
- four_rhovsvsq = 4.0*rhovsvsq
- four_rhovshsq = 4.0*rhovshsq
+ two_rhovsvsq = 2.0_CUSTOM_REAL*rhovsvsq
+ two_rhovshsq = 2.0_CUSTOM_REAL*rhovshsq
+ four_rhovsvsq = 4.0_CUSTOM_REAL*rhovsvsq
+ four_rhovshsq = 4.0_CUSTOM_REAL*rhovshsq
! way 2: pre-compute temporary values
@@ -983,7 +983,7 @@
gammazl = gammaz(i,j,k,ispec)
! compute the jacobian
- jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ jacobianl = 1.0_CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
- xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ xizl*(etaxl*gammayl-etayl*gammaxl))
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_crust_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_crust_mantle.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_crust_mantle.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -218,7 +218,7 @@
if (SAVE_FORWARD .and. nspec2D_xmax_crust_mantle > 0 ) &
call write_abs(1,absorb_xmax_crust_mantle,reclen_xmax_crust_mantle,it)
- endif
+ endif ! NCHUNKS_VAL == 1 .or. ichunk == CHUNK_AB
! ymin
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_outer_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_outer_core.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_outer_core.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -112,6 +112,7 @@
! file access (by process rank modulo 8) showed that the following,
! simple approach is still fastest. (assuming that files are accessed on a local scratch disk)
+ ! outer core
! xmin
! if two chunks exclude this face for one of them
if(NCHUNKS_VAL == 1 .or. ichunk == CHUNK_AC) then
@@ -188,21 +189,21 @@
do ispec2D=1,nspec2D_ymin_outer_core
- ispec=ibelm_ymin_outer_core(ispec2D)
+ ispec=ibelm_ymin_outer_core(ispec2D)
- ! exclude elements that are not on absorbing edges
- if(nkmin_eta_outer_core(1,ispec2D) == 0 .or. nimin_outer_core(1,ispec2D) == 0) cycle
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta_outer_core(1,ispec2D) == 0 .or. nimin_outer_core(1,ispec2D) == 0) cycle
- j=1
- do k=nkmin_eta_outer_core(1,ispec2D),NGLLZ
- do i=nimin_outer_core(1,ispec2D),nimax_outer_core(1,ispec2D)
- iglob=ibool_outer_core(i,j,k,ispec)
+ j=1
+ do k=nkmin_eta_outer_core(1,ispec2D),NGLLZ
+ do i=nimin_outer_core(1,ispec2D),nimax_outer_core(1,ispec2D)
+ iglob=ibool_outer_core(i,j,k,ispec)
- sn = veloc_outer_core(iglob)/vp_outer_core(i,j,k,ispec)
+ sn = veloc_outer_core(iglob)/vp_outer_core(i,j,k,ispec)
- weight=jacobian2D_ymin_outer_core(i,k,ispec2D)*wgllwgll_xz(i,k)
+ weight=jacobian2D_ymin_outer_core(i,k,ispec2D)*wgllwgll_xz(i,k)
- accel_outer_core(iglob) = accel_outer_core(iglob) - weight*sn
+ accel_outer_core(iglob) = accel_outer_core(iglob) - weight*sn
if (SAVE_FORWARD) then
absorb_ymin_outer_core(i,k,ispec2D) = weight*sn
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -33,9 +33,6 @@
! Also, t_shift is added as a new parameter to be written on sac headers!
! by Ebru Bozdag
- !subroutine get_event_info_parallel(myrank,yr,jda,ho,mi,sec,tshift_cmt, &
- ! elat,elon,depth,mb,ename,cmt_lat,cmt_lon,cmt_depth,cmt_hdur,NSOURCES)
-
subroutine get_event_info_parallel(myrank,yr,jda,ho,mi,sec,&
event_name,tshift_cmt,t_shift, &
elat,elon,depth,mb,cmt_lat, &
@@ -72,16 +69,16 @@
!integer, parameter :: LENGTH_REGION_NAME = 150
!character(len=LENGTH_REGION_NAME) region
-! get event information for SAC header on the master
+ ! get event information for SAC header on the master
if(myrank == 0) then
+ ! note: mb as (body wave) moment magnitude is not used any further,
+ ! see comment in write_output_SAC() routine
+
call get_event_info_serial(yr,jda,ho,mi,sec,event_name,tshift_cmt,t_shift, &
elat,elon,depth,mb, &
cmt_lat,cmt_lon,cmt_depth,cmt_hdur,NSOURCES)
- !call get_event_info_serial(yr,jda,ho,mi,sec,tshift_cmt,elat,elon,depth,mb,region, &
- ! cmt_lat,cmt_lon,cmt_depth,cmt_hdur,NSOURCES,LENGTH_REGION_NAME)
-
! create the event name
!write(ename(1:12),'(a12)') region(1:12)
@@ -132,10 +129,6 @@
elat_pde,elon_pde,depth_pde,mb,&
cmt_lat,cmt_lon,cmt_depth,cmt_hdur,NSOURCES)
-
- !subroutine get_event_info_serial(yr,jda,ho,mi,sec,tshift_cmt,elat,elon,depth,mb,region,&
- ! cmt_lat,cmt_lon,cmt_depth,cmt_hdur,NSOURCES,LENGTH_REGION_NAME)
-
implicit none
include "constants.h"
@@ -176,21 +169,9 @@
!
call get_value_string(CMTSOLUTION, 'solver.CMTSOLUTION','DATA/CMTSOLUTION')
- open(unit=821,file=CMTSOLUTION,iostat=ios,status='old',action='read')
+ open(unit=IIN,file=trim(CMTSOLUTION),status='old',action='read',iostat=ios)
if(ios /= 0) stop 'error opening CMTSOLUTION file (in get_event_info_serial)'
- !icounter = 0
- !do while(ios == 0)
- ! read(821,"(a)",iostat=ios) dummystring
- ! if(ios == 0) icounter = icounter + 1
- !enddo
- !close(821)
- !if(mod(icounter,NLINES_PER_CMTSOLUTION_SOURCE) /= 0) &
- ! stop 'total number of lines in CMTSOLUTION file should be a multiple of NLINES_PER_CMTSOLUTION_SOURCE'
- !NSOURCES = icounter / NLINES_PER_CMTSOLUTION_SOURCE
- !if(NSOURCES < 1) stop 'need at least one source in CMTSOLUTION file'
- !open(unit=821,file=CMTSOLUTION,status='old',action='read')
-
! example header line of CMTSOLUTION file
!PDE 2003 09 25 19 50 08.93 41.78 144.08 18.0 7.9 8.0 Hokkaido, Japan
! which is: event_id, date,origin time,latitude,longitude,depth, mb, MS, region
@@ -199,40 +180,40 @@
do isource=1,NSOURCES
! read header with event information
- read(821,*) datasource,yr,mo,da,ho,mi,sec,elat_pde,elon_pde,depth_pde,mb,ms
+ read(IIN,*) datasource,yr,mo,da,ho,mi,sec,elat_pde,elon_pde,depth_pde,mb,ms
jda=julian_day(yr,mo,da)
! ignore line with event name
- read(821,"(a)") string
+ read(IIN,"(a)") string
read(string(12:len_trim(string)),*) e_n(isource)
! read time shift
- read(821,"(a)") string
+ read(IIN,"(a)") string
read(string(12:len_trim(string)),*) t_s(isource)
! read half duration
- read(821,"(a)") string
+ read(IIN,"(a)") string
read(string(15:len_trim(string)),*) hdur(isource)
! read latitude
- read(821,"(a)") string
+ read(IIN,"(a)") string
read(string(10:len_trim(string)),*) lat(isource)
! read longitude
- read(821,"(a)") string
+ read(IIN,"(a)") string
read(string(11:len_trim(string)),*) lon(isource)
! read depth
- read(821,"(a)") string
+ read(IIN,"(a)") string
read(string(7:len_trim(string)),*) depth(isource)
! ignore the last 6 lines with moment tensor info
- read(821,"(a)") string
- read(821,"(a)") string
- read(821,"(a)") string
- read(821,"(a)") string
- read(821,"(a)") string
- read(821,"(a)") string
+ read(IIN,"(a)") string
+ read(IIN,"(a)") string
+ read(IIN,"(a)") string
+ read(IIN,"(a)") string
+ read(IIN,"(a)") string
+ read(IIN,"(a)") string
enddo
! sets tshift_cmt to zero
@@ -257,50 +238,7 @@
t_shift = minval(t_s(1:NSOURCES))
endif
- close(821)
+ close(IIN)
-
-
-! ! read header with event information
-! read(821,*) datasource,yr,mo,da,ho,mi,sec,elat,elon,depth,mb,ms,region
-!
-! jda=julian_day(yr,mo,da)
-!
-! ! ignore line with event name
-! read(821,"(a)") string
-!
-! ! read time shift
-! read(821,"(a)") string
-! read(string(12:len_trim(string)),*) tshift_cmt
-!
-! if (NSOURCES == 1) then
-!
-! ! read half duration
-! read(821,"(a)") string
-! read(string(15:len_trim(string)),*) cmt_hdur
-!
-! ! read latitude
-! read(821,"(a)") string
-! read(string(10:len_trim(string)),*) cmt_lat
-!
-! ! read longitude
-! read(821,"(a)") string
-! read(string(11:len_trim(string)),*) cmt_lon
-!
-! ! read depth
-! read(821,"(a)") string
-! read(string(7:len_trim(string)),*) cmt_depth
-!
-! else
-!
-! cmt_hdur=-1e8
-! cmt_lat=-1e8
-! cmt_lon=-1e8
-! cmt_depth=-1e8
-!
-! endif
-!
-! close(821)
-
end subroutine get_event_info_serial
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/initialize_simulation.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/initialize_simulation.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -291,38 +291,39 @@
write(IMAIN,*) 'model:'
+ ! model mesh parameters
if(ISOTROPIC_3D_MANTLE) then
- write(IMAIN,*) ' incorporates 3-D lateral variations'
+ write(IMAIN,*) ' incorporating 3-D lateral variations'
else
write(IMAIN,*) ' no 3-D lateral variations'
endif
if(HETEROGEN_3D_MANTLE) then
- write(IMAIN,*) ' incorporates heterogeneities in the mantle'
+ write(IMAIN,*) ' incorporating heterogeneities in the mantle'
else
write(IMAIN,*) ' no heterogeneities in the mantle'
endif
if(CRUSTAL) then
- write(IMAIN,*) ' incorporates crustal variations'
+ write(IMAIN,*) ' incorporating crustal variations'
else
write(IMAIN,*) ' no crustal variations'
endif
if(ONE_CRUST) then
- write(IMAIN,*) ' uses one layer only in PREM crust'
+ write(IMAIN,*) ' using one layer only in PREM crust'
else
- write(IMAIN,*) ' uses unmodified 1D crustal model with two layers'
+ write(IMAIN,*) ' using unmodified 1D crustal model with two layers'
endif
if(TRANSVERSE_ISOTROPY) then
- write(IMAIN,*) ' incorporates transverse isotropy'
+ write(IMAIN,*) ' incorporating transverse isotropy'
else
write(IMAIN,*) ' no transverse isotropy'
endif
if(ANISOTROPIC_INNER_CORE_VAL) then
- write(IMAIN,*) ' incorporates anisotropic inner core'
+ write(IMAIN,*) ' incorporating anisotropic inner core'
else
write(IMAIN,*) ' no inner-core anisotropy'
endif
if(ANISOTROPIC_3D_MANTLE_VAL) then
- write(IMAIN,*) ' incorporates anisotropic mantle'
+ write(IMAIN,*) ' incorporating anisotropic mantle'
else
write(IMAIN,*) ' no general mantle anisotropy'
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -116,12 +116,13 @@
double precision etax,etay,etaz
double precision gammax,gammay,gammaz
-! timer MPI
+ ! timer MPI
double precision time_start,tCPU
-! use dynamic allocation
+ ! use dynamic allocation
double precision, dimension(:), allocatable :: final_distance
double precision, dimension(:,:), allocatable :: final_distance_all
+
double precision distmin,final_distance_max
! receiver information
@@ -171,7 +172,7 @@
write(IMAIN,*)
endif
-! define topology of the control element
+ ! define topology of the control element
call hex_nodes(iaddx,iaddy,iaddr)
if(myrank == 0) then
@@ -182,7 +183,7 @@
write(IMAIN,*)
endif
-! allocate memory for arrays using number of stations
+ ! allocate memory for arrays using number of stations
allocate(epidist(nrec), &
ix_initial_guess(nrec), &
iy_initial_guess(nrec), &
@@ -242,9 +243,7 @@
endif
-
-
-! broadcast the information read on the master to the nodes
+ ! broadcast the information read on the master to the nodes
call MPI_BCAST(station_name,nrec*MAX_LENGTH_STATION_NAME,MPI_CHARACTER,0,MPI_COMM_WORLD,ier)
call MPI_BCAST(network_name,nrec*MAX_LENGTH_NETWORK_NAME,MPI_CHARACTER,0,MPI_COMM_WORLD,ier)
@@ -253,44 +252,44 @@
call MPI_BCAST(stele,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
call MPI_BCAST(stbur,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! loop on all the stations to locate them in the mesh
+ ! loop on all the stations to locate them in the mesh
do irec=1,nrec
-! set distance to huge initial value
+ ! set distance to huge initial value
distmin = HUGEVAL
-! convert geographic latitude stlat (degrees) to geocentric colatitude theta (radians)
+ ! convert geographic latitude stlat (degrees) to geocentric colatitude theta (radians)
if(ASSUME_PERFECT_SPHERE) then
- theta = PI/2.0d0 - stlat(irec)*PI/180.0d0
+ theta = PI_OVER_TWO - stlat(irec)*DEGREES_TO_RADIANS
else
- theta = PI/2.0d0 - atan(0.99329534d0*dtan(stlat(irec)*PI/180.0d0))
+ theta = PI_OVER_TWO - atan(0.99329534d0*dtan(stlat(irec)*DEGREES_TO_RADIANS))
endif
- phi = stlon(irec)*PI/180.0d0
+ phi = stlon(irec)*DEGREES_TO_RADIANS
call reduce(theta,phi)
-! compute epicentral distance
+ ! compute epicentral distance
epidist(irec) = acos(cos(theta)*cos(theta_source) + &
- sin(theta)*sin(theta_source)*cos(phi-phi_source))*180.0d0/PI
+ sin(theta)*sin(theta_source)*cos(phi-phi_source))*RADIANS_TO_DEGREES
-! print some information about stations
+ ! print some information about stations
if(myrank == 0) &
write(IMAIN,*) 'Station #',irec,': ',station_name(irec)(1:len_trim(station_name(irec))), &
'.',network_name(irec)(1:len_trim(network_name(irec))), &
' epicentral distance: ',sngl(epidist(irec)),' degrees'
-! record three components for each station
+ ! record three components for each station
do iorientation = 1,3
-! North
+ ! North
if(iorientation == 1) then
stazi = 0.d0
stdip = 0.d0
-! East
+ ! East
else if(iorientation == 2) then
stazi = 90.d0
stdip = 0.d0
-! Vertical
+ ! Vertical
else if(iorientation == 3) then
stazi = 0.d0
stdip = - 90.d0
@@ -298,20 +297,20 @@
call exit_MPI(myrank,'incorrect orientation')
endif
-! get the orientation of the seismometer
- thetan=(90.0d0+stdip)*PI/180.0d0
- phin=stazi*PI/180.0d0
+ ! get the orientation of the seismometer
+ thetan=(90.0d0+stdip)*DEGREES_TO_RADIANS
+ phin=stazi*DEGREES_TO_RADIANS
-! we use the same convention as in Harvard normal modes for the orientation
+ ! we use the same convention as in Harvard normal modes for the orientation
-! vertical component
+ ! vertical component
n(1) = cos(thetan)
-! N-S component
+ ! N-S component
n(2) = - sin(thetan)*cos(phin)
-! E-W component
+ ! E-W component
n(3) = sin(thetan)*sin(phin)
-! get the Cartesian components of n in the model: nu
+ ! get the Cartesian components of n in the model: nu
sint = sin(theta)
cost = cos(theta)
sinp = sin(phi)
@@ -362,7 +361,7 @@
+(y_target(irec)-dble(ystore(iglob)))**2 &
+(z_target(irec)-dble(zstore(iglob)))**2)
-! keep this point if it is closer to the receiver
+ ! keep this point if it is closer to the receiver
if(dist < distmin) then
distmin = dist
ispec_selected_rec(irec) = ispec
@@ -375,13 +374,13 @@
enddo
enddo
-! end of loop on all the spectral elements in current slice
+ ! end of loop on all the spectral elements in current slice
enddo
! end of loop on all the stations
enddo
-! create RECORDHEADER file with usual format for normal-mode codes
+ ! create RECORDHEADER file with usual format for normal-mode codes
if(myrank == 0) then
! get the base pathname for output files
@@ -552,13 +551,13 @@
dy = - (y - y_target(irec))
dz = - (z - z_target(irec))
-! compute increments
-! gamma does not change since we know the receiver is exactly on the surface
+ ! compute increments
+ ! gamma does not change since we know the receiver is exactly on the surface
dxi = xix*dx + xiy*dy + xiz*dz
deta = etax*dx + etay*dy + etaz*dz
if(RECEIVERS_CAN_BE_BURIED) dgamma = gammax*dx + gammay*dy + gammaz*dz
-! update values
+ ! update values
xi = xi + dxi
eta = eta + deta
if(RECEIVERS_CAN_BE_BURIED) gamma = gamma + dgamma
@@ -575,15 +574,15 @@
if (gamma > 1.10d0) gamma = 1.10d0
if (gamma < -1.10d0) gamma = -1.10d0
-! end of non linear iterations
+ ! end of non linear iterations
enddo
-! impose receiver exactly at the surface after final iteration
+ ! impose receiver exactly at the surface after final iteration
if(.not. RECEIVERS_CAN_BE_BURIED) gamma = 1.d0
-! compute final coordinates of point found
+ ! compute final coordinates of point found
call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
! store xi,eta and x,y,z of point found
xi_receiver_subset(irec_in_this_subset) = xi
@@ -625,7 +624,7 @@
call MPI_GATHER(z_found_subset,nrec_SUBSET_current_size,MPI_DOUBLE_PRECISION,z_found_all,nrec_SUBSET_current_size, &
MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! this is executed by main process only
+ ! this is executed by main process only
if(myrank == 0) then
! check that the gather operation went well
@@ -727,7 +726,6 @@
! writes out actual receiver location to vtk file
write(IOVTK,*) sngl(x_found(irec)), sngl(y_found(irec)), sngl(z_found(irec))
endif
-
enddo
! compute maximal distance for all the receivers
@@ -764,7 +762,9 @@
epidist(1:nrec) = epidist_found(1:nrec)
! write the list of stations and associated epicentral distance
- open(unit=27,file=trim(OUTPUT_FILES)//'/output_list_stations.txt',status='unknown')
+ open(unit=27,file=trim(OUTPUT_FILES)//'/output_list_stations.txt', &
+ status='unknown',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening file output_list_stations.txt')
write(27,*)
write(27,*) 'total number of stations: ',nrec
write(27,*)
@@ -777,7 +777,9 @@
! write out a filtered station list
if( NCHUNKS /= 6 ) then
- open(unit=27,file=trim(OUTPUT_FILES)//'/STATIONS_FILTERED',status='unknown')
+ open(unit=27,file=trim(OUTPUT_FILES)//'/STATIONS_FILTERED', &
+ status='unknown',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening file STATIONS_FILTERED')
! loop on all the stations to read station information
do irec = 1,nrec
write(27,'(a8,1x,a3,6x,f8.4,1x,f9.4,1x,f6.1,1x,f6.1)') trim(station_name(irec)),&
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -41,7 +41,7 @@
implicit none
-! standard include of the MPI library
+ ! standard include of the MPI library
include 'mpif.h'
include "constants.h"
@@ -104,11 +104,11 @@
double precision dist,typical_size
double precision xi,eta,gamma,dx,dy,dz,dxi,deta
-! topology of the control points of the surface element
+ ! topology of the control points of the surface element
integer iax,iay,iaz
integer iaddx(NGNOD),iaddy(NGNOD),iaddr(NGNOD)
-! coordinates of the control points of the surface element
+ ! coordinates of the control points of the surface element
double precision xelm(NGNOD),yelm(NGNOD),zelm(NGNOD)
integer iter_loop
@@ -171,7 +171,7 @@
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)
+ ! mask source region (mask values are between 0 and 1, with 0 around sources)
real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable :: mask_source
! **************
@@ -203,7 +203,7 @@
call MPI_BCAST(moment_tensor,6*NSOURCES,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
call MPI_BCAST(min_tshift_cmt_original,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! define topology of the control element
+ ! define topology of the control element
call hex_nodes(iaddx,iaddy,iaddr)
! initializes source mask
@@ -215,15 +215,15 @@
! get MPI starting time for all sources
time_start = MPI_WTIME()
-! loop on all the sources
-! gather source information in subsets to reduce memory requirements
+ ! loop on all the sources
+ ! gather source information in subsets to reduce memory requirements
-! loop over subsets of sources
+ ! loop over subsets of sources
do isources_already_done = 0, NSOURCES, NSOURCES_SUBSET_MAX
-! the size of the subset can be the maximum size, or less (if we are in the last subset,
-! or if there are fewer sources than the maximum size of a subset)
- NSOURCES_SUBSET_current_size = min(NSOURCES_SUBSET_MAX, NSOURCES - isources_already_done)
+ ! the size of the subset can be the maximum size, or less (if we are in the last subset,
+ ! or if there are fewer sources than the maximum size of a subset)
+ NSOURCES_SUBSET_current_size = min(NSOURCES_SUBSET_MAX, NSOURCES - isources_already_done)
! allocate arrays specific to each subset
allocate(final_distance_source_subset(NSOURCES_SUBSET_current_size), &
@@ -244,92 +244,92 @@
z_found_source(NSOURCES_SUBSET_current_size),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating temporary source arrays')
-! make sure we clean the subset array before the gather
- ispec_selected_source_subset(:) = 0
+ ! make sure we clean the subset array before the gather
+ ispec_selected_source_subset(:) = 0
-! loop over sources within this subset
- do isource_in_this_subset = 1,NSOURCES_SUBSET_current_size
+ ! loop over sources within this subset
+ do isource_in_this_subset = 1,NSOURCES_SUBSET_current_size
-! mapping from source number in current subset to real source number in all the subsets
- isource = isource_in_this_subset + isources_already_done
+ ! mapping from source number in current subset to real source number in all the subsets
+ isource = isource_in_this_subset + isources_already_done
-! convert geographic latitude lat (degrees) to geocentric colatitude theta (radians)
- if(ASSUME_PERFECT_SPHERE) then
- theta = PI/2.0d0 - lat(isource)*PI/180.0d0
- else
- theta = PI/2.0d0 - atan(0.99329534d0*dtan(lat(isource)*PI/180.0d0))
- endif
+ ! convert geographic latitude lat (degrees) to geocentric colatitude theta (radians)
+ if(ASSUME_PERFECT_SPHERE) then
+ theta = PI_OVER_TWO - lat(isource)*DEGREES_TO_RADIANS
+ else
+ theta = PI_OVER_TWO - atan(0.99329534d0*dtan(lat(isource)*DEGREES_TO_RADIANS))
+ endif
- phi = long(isource)*PI/180.0d0
- call reduce(theta,phi)
+ phi = long(isource)*DEGREES_TO_RADIANS
+ call reduce(theta,phi)
-! get the moment tensor
- Mrr = moment_tensor(1,isource)
- Mtt = moment_tensor(2,isource)
- Mpp = moment_tensor(3,isource)
- Mrt = moment_tensor(4,isource)
- Mrp = moment_tensor(5,isource)
- Mtp = moment_tensor(6,isource)
+ ! get the moment tensor
+ Mrr = moment_tensor(1,isource)
+ Mtt = moment_tensor(2,isource)
+ Mpp = moment_tensor(3,isource)
+ Mrt = moment_tensor(4,isource)
+ Mrp = moment_tensor(5,isource)
+ Mtp = moment_tensor(6,isource)
-! convert from a spherical to a Cartesian representation of the moment tensor
- st=dsin(theta)
- ct=dcos(theta)
- sp=dsin(phi)
- cp=dcos(phi)
+ ! convert from a spherical to a Cartesian representation of the moment tensor
+ st=dsin(theta)
+ ct=dcos(theta)
+ sp=dsin(phi)
+ cp=dcos(phi)
- Mxx(isource)=st*st*cp*cp*Mrr+ct*ct*cp*cp*Mtt+sp*sp*Mpp &
- +2.0d0*st*ct*cp*cp*Mrt-2.0d0*st*sp*cp*Mrp-2.0d0*ct*sp*cp*Mtp
- Myy(isource)=st*st*sp*sp*Mrr+ct*ct*sp*sp*Mtt+cp*cp*Mpp &
- +2.0d0*st*ct*sp*sp*Mrt+2.0d0*st*sp*cp*Mrp+2.0d0*ct*sp*cp*Mtp
- Mzz(isource)=ct*ct*Mrr+st*st*Mtt-2.0d0*st*ct*Mrt
- Mxy(isource)=st*st*sp*cp*Mrr+ct*ct*sp*cp*Mtt-sp*cp*Mpp &
- +2.0d0*st*ct*sp*cp*Mrt+st*(cp*cp-sp*sp)*Mrp+ct*(cp*cp-sp*sp)*Mtp
- Mxz(isource)=st*ct*cp*Mrr-st*ct*cp*Mtt &
- +(ct*ct-st*st)*cp*Mrt-ct*sp*Mrp+st*sp*Mtp
- Myz(isource)=st*ct*sp*Mrr-st*ct*sp*Mtt &
- +(ct*ct-st*st)*sp*Mrt+ct*cp*Mrp-st*cp*Mtp
+ Mxx(isource)=st*st*cp*cp*Mrr+ct*ct*cp*cp*Mtt+sp*sp*Mpp &
+ +2.0d0*st*ct*cp*cp*Mrt-2.0d0*st*sp*cp*Mrp-2.0d0*ct*sp*cp*Mtp
+ Myy(isource)=st*st*sp*sp*Mrr+ct*ct*sp*sp*Mtt+cp*cp*Mpp &
+ +2.0d0*st*ct*sp*sp*Mrt+2.0d0*st*sp*cp*Mrp+2.0d0*ct*sp*cp*Mtp
+ Mzz(isource)=ct*ct*Mrr+st*st*Mtt-2.0d0*st*ct*Mrt
+ Mxy(isource)=st*st*sp*cp*Mrr+ct*ct*sp*cp*Mtt-sp*cp*Mpp &
+ +2.0d0*st*ct*sp*cp*Mrt+st*(cp*cp-sp*sp)*Mrp+ct*(cp*cp-sp*sp)*Mtp
+ Mxz(isource)=st*ct*cp*Mrr-st*ct*cp*Mtt &
+ +(ct*ct-st*st)*cp*Mrt-ct*sp*Mrp+st*sp*Mtp
+ Myz(isource)=st*ct*sp*Mrr-st*ct*sp*Mtt &
+ +(ct*ct-st*st)*sp*Mrt+ct*cp*Mrp-st*cp*Mtp
-! record three components for each station
- do iorientation = 1,3
+ ! record three components for each station
+ do iorientation = 1,3
-! North
- if(iorientation == 1) then
- stazi = 0.d0
- stdip = 0.d0
-! East
- else if(iorientation == 2) then
- stazi = 90.d0
- stdip = 0.d0
-! Vertical
- else if(iorientation == 3) then
- stazi = 0.d0
- stdip = - 90.d0
- else
- call exit_MPI(myrank,'incorrect orientation')
- endif
+ ! North
+ if(iorientation == 1) then
+ stazi = 0.d0
+ stdip = 0.d0
+ ! East
+ else if(iorientation == 2) then
+ stazi = 90.d0
+ stdip = 0.d0
+ ! Vertical
+ else if(iorientation == 3) then
+ stazi = 0.d0
+ stdip = - 90.d0
+ else
+ call exit_MPI(myrank,'incorrect orientation')
+ endif
-! get the orientation of the seismometer
- thetan=(90.0d0+stdip)*PI/180.0d0
- phin=stazi*PI/180.0d0
+ ! get the orientation of the seismometer
+ thetan=(90.0d0+stdip)*DEGREES_TO_RADIANS
+ phin=stazi*DEGREES_TO_RADIANS
-! we use the same convention as in Harvard normal modes for the orientation
+ ! we use the same convention as in Harvard normal modes for the orientation
-! vertical component
- n(1) = dcos(thetan)
-! N-S component
- n(2) = - dsin(thetan)*dcos(phin)
-! E-W component
- n(3) = dsin(thetan)*dsin(phin)
+ ! vertical component
+ n(1) = dcos(thetan)
+ ! N-S component
+ n(2) = - dsin(thetan)*dcos(phin)
+ ! E-W component
+ n(3) = dsin(thetan)*dsin(phin)
-! get the Cartesian components of n in the model: nu
- nu_source(iorientation,1,isource) = n(1)*st*cp+n(2)*ct*cp-n(3)*sp
- nu_source(iorientation,2,isource) = n(1)*st*sp+n(2)*ct*sp+n(3)*cp
- nu_source(iorientation,3,isource) = n(1)*ct-n(2)*st
+ ! get the Cartesian components of n in the model: nu
+ nu_source(iorientation,1,isource) = n(1)*st*cp+n(2)*ct*cp-n(3)*sp
+ nu_source(iorientation,2,isource) = n(1)*st*sp+n(2)*ct*sp+n(3)*cp
+ nu_source(iorientation,3,isource) = n(1)*ct-n(2)*st
- enddo
+ enddo
-! normalized source radius
- r0 = R_UNIT_SPHERE
+ ! normalized source radius
+ r0 = R_UNIT_SPHERE
if(ELLIPTICITY) then
if(TOPOGRAPHY) then
@@ -421,19 +421,19 @@
enddo
enddo
-! calculates a gaussian mask around source point
- if( SAVE_SOURCE_MASK .and. SIMULATION_TYPE == 3 ) then
- call calc_mask_source(mask_source,ispec,NSPEC,typical_size, &
- x_target_source,y_target_source,z_target_source, &
- ibool,xstore,ystore,zstore,NGLOB)
- endif
+ ! calculates a gaussian mask around source point
+ if( SAVE_SOURCE_MASK .and. SIMULATION_TYPE == 3 ) then
+ call calc_mask_source(mask_source,ispec,NSPEC,typical_size, &
+ x_target_source,y_target_source,z_target_source, &
+ ibool,xstore,ystore,zstore,NGLOB)
+ endif
-! end of loop on all the elements in current slice
- enddo
+ ! end of loop on all the elements in current slice
+ enddo
-! *******************************************
-! find the best (xi,eta,gamma) for the source
-! *******************************************
+ ! *******************************************
+ ! find the best (xi,eta,gamma) for the source
+ ! *******************************************
! if we have not located a target element, the source is not in this slice
! therefore use first element only for fictitious iterative search
@@ -444,29 +444,29 @@
iz_initial_guess_source = 2
endif
- ! for point sources, the location will be exactly at a GLL point
- ! otherwise this tries to find best location
- if( USE_FORCE_POINT_SOURCE ) then
- ! store xi,eta,gamma and x,y,z of point found
- ! note: they have range [1.0d0,NGLLX/Y/Z], used for point sources
- ! see e.g. in compute_add_sources.f90
- xi_source_subset(isource_in_this_subset) = dble(ix_initial_guess_source)
- eta_source_subset(isource_in_this_subset) = dble(iy_initial_guess_source)
- gamma_source_subset(isource_in_this_subset) = dble(iz_initial_guess_source)
+ ! for point sources, the location will be exactly at a GLL point
+ ! otherwise this tries to find best location
+ if( USE_FORCE_POINT_SOURCE ) then
+ ! store xi,eta,gamma and x,y,z of point found
+ ! note: they have range [1.0d0,NGLLX/Y/Z], used for point sources
+ ! see e.g. in compute_add_sources.f90
+ xi_source_subset(isource_in_this_subset) = dble(ix_initial_guess_source)
+ eta_source_subset(isource_in_this_subset) = dble(iy_initial_guess_source)
+ gamma_source_subset(isource_in_this_subset) = dble(iz_initial_guess_source)
- iglob = ibool(ix_initial_guess_source,iy_initial_guess_source, &
- iz_initial_guess_source,ispec_selected_source_subset(isource_in_this_subset))
- x_found_source(isource_in_this_subset) = xstore(iglob)
- y_found_source(isource_in_this_subset) = ystore(iglob)
- z_found_source(isource_in_this_subset) = zstore(iglob)
+ iglob = ibool(ix_initial_guess_source,iy_initial_guess_source, &
+ iz_initial_guess_source,ispec_selected_source_subset(isource_in_this_subset))
+ x_found_source(isource_in_this_subset) = xstore(iglob)
+ y_found_source(isource_in_this_subset) = ystore(iglob)
+ z_found_source(isource_in_this_subset) = zstore(iglob)
- ! compute final distance between asked and found (converted to km)
- final_distance_source_subset(isource_in_this_subset) = &
- dsqrt((x_target_source-x_found_source(isource_in_this_subset))**2 + &
- (y_target_source-y_found_source(isource_in_this_subset))**2 + &
- (z_target_source-z_found_source(isource_in_this_subset))**2)*R_EARTH/1000.d0
+ ! compute final distance between asked and found (converted to km)
+ final_distance_source_subset(isource_in_this_subset) = &
+ dsqrt((x_target_source-x_found_source(isource_in_this_subset))**2 + &
+ (y_target_source-y_found_source(isource_in_this_subset))**2 + &
+ (z_target_source-z_found_source(isource_in_this_subset))**2)*R_EARTH/1000.d0
- else
+ else
! use initial guess in xi, eta and gamma
xi = xigll(ix_initial_guess_source)
@@ -513,21 +513,22 @@
enddo
- ! iterate to solve the non linear system
- do iter_loop = 1,NUM_ITER
+ ! iterate to solve the non linear system
+ do iter_loop = 1,NUM_ITER
- ! recompute jacobian for the new point
- call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+ ! recompute jacobian for the new point
+ call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
- ! compute distance to target location
- dx = - (x - x_target_source)
- dy = - (y - y_target_source)
- dz = - (z - z_target_source)
+ ! compute distance to target location
+ 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
- dgamma = gammax*dx + gammay*dy + gammaz*dz
+ ! compute increments
+ dxi = xix*dx + xiy*dy + xiz*dz
+ deta = etax*dx + etay*dy + etaz*dz
+ dgamma = gammax*dx + gammay*dy + gammaz*dz
! update values
xi = xi + dxi
@@ -663,20 +664,22 @@
endif
write(IMAIN,*) ' time shift: ',tshift_cmt(isource),' seconds'
- ! writes out actual source position to vtk file
- write(IOVTK,*) sngl(x_found_source(isource_in_this_subset)), &
- sngl(y_found_source(isource_in_this_subset)), &
- sngl(z_found_source(isource_in_this_subset))
+ ! writes out actual source position to vtk file
+ write(IOVTK,*) sngl(x_found_source(isource_in_this_subset)), &
+ sngl(y_found_source(isource_in_this_subset)), &
+ sngl(z_found_source(isource_in_this_subset))
- ! get latitude, longitude and depth of the source that will be used
- call xyz_2_rthetaphi_dble(x_found_source(isource_in_this_subset),y_found_source(isource_in_this_subset), &
- z_found_source(isource_in_this_subset),r_found_source,theta_source(isource),phi_source(isource))
- call reduce(theta_source(isource),phi_source(isource))
+ ! get latitude, longitude and depth of the source that will be used
+ call xyz_2_rthetaphi_dble(x_found_source(isource_in_this_subset), &
+ y_found_source(isource_in_this_subset), &
+ z_found_source(isource_in_this_subset), &
+ r_found_source,theta_source(isource),phi_source(isource))
+ call reduce(theta_source(isource),phi_source(isource))
- ! convert geocentric to geographic colatitude
- colat_source = PI/2.0d0 &
- - 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
+ ! convert geocentric to geographic colatitude
+ colat_source = PI_OVER_TWO &
+ - 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,*)
write(IMAIN,*) 'original (requested) position of the source:'
@@ -808,17 +811,17 @@
enddo ! end of loop over all source subsets
-! display maximum error in location estimate
+ ! display maximum error in location estimate
if(myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) 'maximum error in location of the sources: ',sngl(maxval(final_distance_source)),' km'
write(IMAIN,*)
endif
-
-! main process broadcasts the results to all the slices
+ ! 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)
+
call MPI_BCAST(xi_source,NSOURCES,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
call MPI_BCAST(eta_source,NSOURCES,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
call MPI_BCAST(gamma_source,NSOURCES,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/multiply_arrays_source.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/multiply_arrays_source.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/multiply_arrays_source.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -38,17 +38,19 @@
include "constants.h"
-! source arrays
+ ! source arrays
double precision, dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrayd
double precision, dimension(NGLLX,NGLLY,NGLLZ) :: G11,G12,G13,G21,G22,G23,G31,G32,G33
double precision, dimension(NGLLX) :: hxis,hpxis
double precision, dimension(NGLLY) :: hetas,hpetas
double precision, dimension(NGLLZ) :: hgammas,hpgammas
- integer k,l,m
+ integer :: k,l,m
- integer ir,it,iv
+ ! local parameters
+ integer :: ir,it,iv
+ ! initializes
sourcearrayd(:,k,l,m) = ZERO
do iv=1,NGLLZ
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90 2013-07-01 03:05:54 UTC (rev 22473)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90 2013-07-01 04:15:01 UTC (rev 22474)
@@ -277,7 +277,9 @@
end subroutine write_seismograms
-!=====================================================================
+!
+!-------------------------------------------------------------------------------------------------
+!
subroutine write_one_seismogram(one_seismogram,irec, &
station_name,network_name,stlat,stlon,stele,stbur,nrec, &
@@ -350,6 +352,7 @@
!----------------------------------------------------------------
call band_instrument_code(DT,bic)
+
if (ROTATE_SEISMOGRAMS_RT) then ! iorientation 1=N,2=E,3=Z,4=R,5=T
ior_start=3 ! starting from Z
ior_end =5 ! ending with T => ZRT
@@ -358,8 +361,6 @@
ior_end =3 ! ending with Z => NEZ
endif
- !do iorientation = 1,NDIM
- !do iorientation = 1,5 ! BS BS changed from 3 (NEZ) to 5 (NEZRT) components
do iorientation = ior_start,ior_end ! BS BS changed according to ROTATE_SEISMOGRAMS_RT
if(iorientation == 1) then
@@ -397,8 +398,8 @@
phi = backaz
endif
- cphi=cos(phi*pi/180)
- sphi=sin(phi*pi/180)
+ cphi=cos(phi*DEGREES_TO_RADIANS)
+ sphi=sin(phi*DEGREES_TO_RADIANS)
! BS BS do the rotation of the components and put result in
! new variable seismogram_tmp
@@ -443,7 +444,6 @@
! SAC output format
if (OUTPUT_SEISMOS_SAC_ALPHANUM .or. OUTPUT_SEISMOS_SAC_BINARY) then
-
call write_output_SAC(seismogram_tmp,irec, &
station_name,network_name,stlat,stlon,stele,stbur,nrec, &
ANGULAR_WIDTH_XI_IN_DEGREES,NEX_XI,DT,hdur,it_end, &
@@ -453,126 +453,124 @@
OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY,MODEL, &
NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current, &
iorientation,phi,chn,sisname)
-
endif ! OUTPUT_SEISMOS_SAC_ALPHANUM .or. OUTPUT_SEISMOS_SAC_BINARY
! ASCII output format
if(OUTPUT_SEISMOS_ASCII_TEXT) then
-
call write_output_ASCII(seismogram_tmp, &
DT,hdur,OUTPUT_FILES, &
NTSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current, &
SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,myrank, &
iorientation,sisname,sisname_big_file)
-
endif ! OUTPUT_SEISMOS_ASCII_TEXT
enddo ! do iorientation
end subroutine write_one_seismogram
-!=====================================================================
+!
+!-------------------------------------------------------------------------------------------------
+!
! write adjoint seismograms to text files
- subroutine write_adj_seismograms(seismograms,number_receiver_global, &
+ subroutine write_adj_seismograms(seismograms,number_receiver_global, &
nrec_local,it,nit_written,DT,NSTEP, &
NTSTEP_BETWEEN_OUTPUT_SEISMOS,hdur,LOCAL_PATH)
- implicit none
+ implicit none
- include "constants.h"
+ include "constants.h"
- integer nrec_local,NSTEP,NTSTEP_BETWEEN_OUTPUT_SEISMOS,it,nit_written
- integer, dimension(nrec_local) :: number_receiver_global
- real(kind=CUSTOM_REAL), dimension(9,nrec_local,NSTEP) :: seismograms
- double precision hdur,DT
- character(len=150) LOCAL_PATH
+ integer :: nrec_local,NSTEP,NTSTEP_BETWEEN_OUTPUT_SEISMOS,it,nit_written
+ integer, dimension(nrec_local) :: number_receiver_global
+ real(kind=CUSTOM_REAL), dimension(9,nrec_local,NSTEP) :: seismograms
+ double precision :: hdur,DT
+ character(len=150) :: LOCAL_PATH
- integer irec,irec_local
- integer iorientation,isample
+ integer :: irec,irec_local
+ integer :: iorientation,isample
- character(len=4) chn
- character(len=150) clean_LOCAL_PATH,final_LOCAL_PATH
- character(len=256) sisname
- character(len=2) bic
+ character(len=4) :: chn
+ character(len=150) :: clean_LOCAL_PATH,final_LOCAL_PATH
+ character(len=256) :: sisname
+ character(len=2) :: bic
- call band_instrument_code(DT,bic)
+ call band_instrument_code(DT,bic)
- do irec_local = 1,nrec_local
+ do irec_local = 1,nrec_local
-! get global number of that receiver
- irec = number_receiver_global(irec_local)
+ ! get global number of that receiver
+ irec = number_receiver_global(irec_local)
- do iorientation = 1,9
-
- if(iorientation == 1) then
+ do iorientation = 1,9
+ if(iorientation == 1) then
chn = 'SNN'
- else if(iorientation == 2) then
+ else if(iorientation == 2) then
chn = 'SEE'
- else if(iorientation == 3) then
+ else if(iorientation == 3) then
chn = 'SZZ'
- else if(iorientation == 4) then
+ else if(iorientation == 4) then
chn = 'SNE'
- else if(iorientation == 5) then
+ else if(iorientation == 5) then
chn = 'SNZ'
- else if(iorientation == 6) then
+ else if(iorientation == 6) then
chn = 'SEZ'
- else if(iorientation == 7) then
+ else if(iorientation == 7) then
!chn = 'LHN'
chn = bic(1:2)//'N'
- else if(iorientation == 8) then
+ else if(iorientation == 8) then
chn = bic(1:2)//'E'
- else if(iorientation == 9) then
+ else if(iorientation == 9) then
chn = bic(1:2)//'Z'
- endif
+ endif
-! create the name of the seismogram file for each slice
-! file name includes the name of the station, the network and the component
- write(sisname,"(a,i6.6,'.',a,'.',a3,'.sem')") 'S',irec,'NT',chn
+ ! create the name of the seismogram file for each slice
+ ! file name includes the name of the station, the network and the component
+ write(sisname,"(a,i6.6,'.',a,'.',a3,'.sem')") 'S',irec,'NT',chn
-! suppress white spaces if any
- clean_LOCAL_PATH = adjustl(LOCAL_PATH)
+ ! suppress white spaces if any
+ clean_LOCAL_PATH = adjustl(LOCAL_PATH)
-! create full final local path
- final_LOCAL_PATH = clean_LOCAL_PATH(1:len_trim(clean_LOCAL_PATH)) // '/'
+ ! create full final local path
+ final_LOCAL_PATH = clean_LOCAL_PATH(1:len_trim(clean_LOCAL_PATH)) // '/'
-! save seismograms in text format with no subsampling.
-! Because we do not subsample the output, this can result in large files
-! if the simulation uses many time steps. However, subsampling the output
-! here would result in a loss of accuracy when one later convolves
-! the results with the source time function
- if(it <= NTSTEP_BETWEEN_OUTPUT_SEISMOS) then
- !open new file
- open(unit=IOUT,file=final_LOCAL_PATH(1:len_trim(final_LOCAL_PATH))//sisname(1:len_trim(sisname)),&
- status='unknown',action='write')
- else if(it > NTSTEP_BETWEEN_OUTPUT_SEISMOS) then
- !append to existing file
- open(unit=IOUT,file=final_LOCAL_PATH(1:len_trim(final_LOCAL_PATH))//sisname(1:len_trim(sisname)),&
- status='old',position='append',action='write')
- endif
-! make sure we never write more than the maximum number of time steps
-! subtract half duration of the source to make sure travel time is correct
- do isample = nit_written+1,min(it,NSTEP)
-! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- write(IOUT,*) sngl(dble(isample-1)*DT - hdur),' ',seismograms(iorientation,irec_local,isample-nit_written)
- else
- write(IOUT,*) dble(isample-1)*DT - hdur,' ',seismograms(iorientation,irec_local,isample-nit_written)
- endif
- enddo
+ ! save seismograms in text format with no subsampling.
+ ! Because we do not subsample the output, this can result in large files
+ ! if the simulation uses many time steps. However, subsampling the output
+ ! here would result in a loss of accuracy when one later convolves
+ ! the results with the source time function
+ if(it <= NTSTEP_BETWEEN_OUTPUT_SEISMOS) then
+ !open new file
+ open(unit=IOUT,file=final_LOCAL_PATH(1:len_trim(final_LOCAL_PATH))//sisname(1:len_trim(sisname)),&
+ status='unknown',action='write')
+ else if(it > NTSTEP_BETWEEN_OUTPUT_SEISMOS) then
+ !append to existing file
+ open(unit=IOUT,file=final_LOCAL_PATH(1:len_trim(final_LOCAL_PATH))//sisname(1:len_trim(sisname)),&
+ status='old',position='append',action='write')
+ endif
+ ! make sure we never write more than the maximum number of time steps
+ ! subtract half duration of the source to make sure travel time is correct
+ do isample = nit_written+1,min(it,NSTEP)
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ write(IOUT,*) sngl(dble(isample-1)*DT - hdur),' ',seismograms(iorientation,irec_local,isample-nit_written)
+ else
+ write(IOUT,*) dble(isample-1)*DT - hdur,' ',seismograms(iorientation,irec_local,isample-nit_written)
+ endif
+ enddo
- close(IOUT)
+ close(IOUT)
+ enddo
+ enddo
- enddo
+ end subroutine write_adj_seismograms
- enddo
+!
+!-------------------------------------------------------------------------------------------------
+!
- end subroutine write_adj_seismograms
-
-!=====================================================================
-
subroutine band_instrument_code(DT,bic)
! This subroutine is to choose the appropriate band and instrument codes for channel names of seismograms
! based on the IRIS convention (first two letters of channel codes which were LH(Z/E/N) previously).
More information about the CIG-COMMITS
mailing list