[cig-commits] r17241 - in seismo/3D/SPECFEM3D/trunk: . USER_MANUAL UTILS/oldstuff
pieyre at geodynamics.org
pieyre at geodynamics.org
Thu Oct 7 06:50:43 PDT 2010
Author: pieyre
Date: 2010-10-07 06:50:42 -0700 (Thu, 07 Oct 2010)
New Revision: 17241
Added:
seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/check_buffers_2D.f90
seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/get_MPI_cutplanes_eta.f90
seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/get_MPI_cutplanes_xi.f90
seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/model_interface_bedrock.f90
Removed:
seismo/3D/SPECFEM3D/trunk/check_buffers_2D.f90
seismo/3D/SPECFEM3D/trunk/get_MPI_cutplanes_eta.f90
seismo/3D/SPECFEM3D/trunk/get_MPI_cutplanes_xi.f90
seismo/3D/SPECFEM3D/trunk/model_interface_bedrock.f90
Modified:
seismo/3D/SPECFEM3D/trunk/Makefile.in
seismo/3D/SPECFEM3D/trunk/USER_MANUAL/manual_SPECFEM3D.tex
seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90
Log:
modified some files which didn't compile with gfortran ! updated the manual, moved old source code files in UTILS/oldstuff and modified the Makefile accordingly
Modified: seismo/3D/SPECFEM3D/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/Makefile.in 2010-10-07 10:09:09 UTC (rev 17240)
+++ seismo/3D/SPECFEM3D/trunk/Makefile.in 2010-10-07 13:50:42 UTC (rev 17241)
@@ -75,8 +75,6 @@
$O/get_coupling_surfaces.o \
$O/get_model.o \
$O/get_MPI.o \
- $O/get_MPI_cutplanes_eta.o \
- $O/get_MPI_cutplanes_xi.o \
$O/get_attenuation_model.o \
$O/get_cmt.o \
$O/get_element_face.o \
@@ -195,7 +193,6 @@
@COND_PYRE_FALSE at DEFAULT = \
@COND_PYRE_FALSE@ generate_databases \
@COND_PYRE_FALSE@ specfem3D \
- at COND_PYRE_FALSE@ check_buffers_2D \
@COND_PYRE_FALSE@ combine_AVS_DX \
@COND_PYRE_FALSE@ combine_vol_data \
@COND_PYRE_FALSE@ combine_surf_data \
@@ -229,7 +226,6 @@
@COND_PYRE_FALSE@ ${FCLINK} -o xspecfem3D $(XSPECFEM_OBJECTS) $(COND_MPI_OBJECTS) $(MPILIBS)
@COND_PYRE_FALSE@
-check_buffers_2D: xcheck_buffers_2D
combine_AVS_DX: xcombine_AVS_DX
convolve_source_timefunction: xconvolve_source_timefunction
create_header_file: xcreate_header_file
@@ -250,9 +246,6 @@
xcombine_AVS_DX: $O/combine_AVS_DX.o $(LIBSPECFEM)
${FCCOMPILE_CHECK} -o xcombine_AVS_DX $O/combine_AVS_DX.o $(LIBSPECFEM)
-xcheck_buffers_2D: $O/check_buffers_2D.o $(LIBSPECFEM)
- ${FCCOMPILE_CHECK} -o xcheck_buffers_2D $O/check_buffers_2D.o $(LIBSPECFEM)
-
xcombine_vol_data: $O/combine_vol_data.o $O/write_c_binary.o $O/read_parameter_file.o $O/read_value_parameters.o $O/get_value_parameters.o $O/param_reader.o
${FCCOMPILE_CHECK} -o xcombine_vol_data $O/combine_vol_data.o $O/write_c_binary.o $O/read_parameter_file.o $O/read_value_parameters.o $O/get_value_parameters.o $O/param_reader.o
@@ -262,7 +255,7 @@
clean:
rm -f $O/* *.o *.gnu *.mod OUTPUT_FILES/timestamp* OUTPUT_FILES/starttime*txt work.pc* \
- xgenerate_databases xspecfem3D xcombine_AVS_DX xcheck_buffers_2D \
+ xgenerate_databases xspecfem3D xcombine_AVS_DX \
xconvolve_source_timefunction xcreate_header_file \
xcreate_movie_shakemap_AVS_DX_GMT xcombine_vol_data xcombine_surf_data
@@ -439,9 +432,6 @@
$O/save_header_file.o: constants.h save_header_file.f90
${FCCOMPILE_CHECK} -c -o $O/save_header_file.o save_header_file.f90
-$O/check_buffers_2D.o: constants.h check_buffers_2D.f90
- ${FCCOMPILE_CHECK} -c -o $O/check_buffers_2D.o check_buffers_2D.f90
-
$O/read_parameter_file.o: constants.h read_parameter_file.f90
${FCCOMPILE_CHECK} -c -o $O/read_parameter_file.o read_parameter_file.f90
@@ -475,12 +465,6 @@
$O/get_flags_boundaries.o: constants.h get_flags_boundaries.f90
${FCCOMPILE_CHECK} -c -o $O/get_flags_boundaries.o get_flags_boundaries.f90
-$O/get_MPI_cutplanes_xi.o: constants.h get_MPI_cutplanes_xi.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_MPI_cutplanes_xi.o get_MPI_cutplanes_xi.f90
-
-$O/get_MPI_cutplanes_eta.o: constants.h get_MPI_cutplanes_eta.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_MPI_cutplanes_eta.o get_MPI_cutplanes_eta.f90
-
$O/get_cmt.o: constants.h get_cmt.f90
${FCCOMPILE_CHECK} -c -o $O/get_cmt.o get_cmt.f90
Modified: seismo/3D/SPECFEM3D/trunk/USER_MANUAL/manual_SPECFEM3D.tex
===================================================================
--- seismo/3D/SPECFEM3D/trunk/USER_MANUAL/manual_SPECFEM3D.tex 2010-10-07 10:09:09 UTC (rev 17240)
+++ seismo/3D/SPECFEM3D/trunk/USER_MANUAL/manual_SPECFEM3D.tex 2010-10-07 13:50:42 UTC (rev 17241)
@@ -319,7 +319,7 @@
\par\end{flushleft}
You are now ready to compile the internal mesher. This is an alternative to CUBIT
-for mesh generation of relatively simple geological models. The mesher is no longer
+for the mesh generation of relatively simple geological models. The mesher is no longer
dedicated to Southern California and more flexiblity is provided in this version of the package.
In the directory \texttt{meshfem3D} type `\texttt{make meshfem3D}'. If all paths and flags
@@ -346,14 +346,6 @@
the \texttt{Par\_file}:
\begin{description}
-%% \item [{\texttt{SIMULATION\_TYPE}}] is set to 1 for forward simulations,
-%% 2 for adjoint simulations (see Section \ref{sec:Adjoint-simulation-finite})
-%% and 3 for kernel simulations (see Section \ref{sec:Finite-Frequency-Kernels}).
-%% \item [{\texttt{SAVE\_FORWARD}}] is only set to \texttt{.true.} for a forward
-%% simulation with the last frame of the simulation saved, as part of
-%% the finite-frequency kernel calculations (see Section \ref{sec:Finite-Frequency-Kernels}).
-%% For a regular forward simulation, leave \texttt{SIMULATION\_TYPE}
-%% and \texttt{SAVE\_FORWARD} at their default values.
\item [{\texttt{LATITUDE\_MIN}}] Minimum latitude in the block (negative
for South).
\item [{\texttt{LATITUDE\_MAX}}] Maximum latitude in the block.
@@ -388,14 +380,6 @@
by \texttt{NGLLX} in the \texttt{constants.h} file. We generally use
$\mbox{\texttt{NGLLX\/}}=5$, for a total of $5^{3}=125$ points per
elements. We suggest not to change this value.
-%% \item [{\texttt{\noun{UTM\_PROJECTION\_ZONE}}}] UTM projection zone in
-%% which your model resides, only valid when \texttt{SUPPRESS\_UTM\_}~\\
-%% \texttt{PROJECTION} is \texttt{.false.}.
-%% \item [{\texttt{SUPPRESS\_UTM\_PROJECTION}}] set to be \texttt{.false.}
-%% when your model range is specified in the geographical coordinates,
-%% and needs to be \texttt{.true.} when your model is specified in a
-%% cartesian coordinates. \noun{UTM projection zone in which your simulation
-%% region resides.}
\item [{$\nexeta$}] The number of spectral elements along the other side
of the block. This number \textit{must} be 8~$\times$~a multiple
of $\nproceta$ defined below.
@@ -415,79 +399,11 @@
if \texttt{USE\_REGULAR\_MESH} is set to \texttt{.true.}).
\item [{\texttt{NZ\_DOUBLING\_2}}] The position of the second doubling layer (only interpreted
if \texttt{USE\_REGULAR\_MESH} is set to \texttt{.true.} and if \texttt{NDOUBLINGS} is set to \texttt{2}).
-
-%%\item [{\texttt{MODEL}}] Must be set to one of the following:
-%%\begin{description}
-%% \item [{\texttt{SoCal}}] Isotropic, southern California layercake model
-%% developed by \citet{DrHe90}.
-%% \item [{\texttt{Harvard\_LA}}] 3D model based upon the high-resolution
-%% Los Angeles basin model developed by \citet{SuSh03} the Salton Trough
-%% model developed by \citet{lovelyetal06}, the regional tomographic
-%% model of \citet{hauksson2000}, and the Moho map determined by \citet{zhu&kanamori2000}.
-%% \end{description}
-%% \item [{\texttt{OCEANS}}] Set to \texttt{.true.} if the effect of the oceans
-%% on seismic wave propagation should be incorporated based upon the
-%% approximate treatment discussed in \citet{KoTr02b}. This feature
-%% is inexpensive from a numerical perspective, both in terms of memory
-%% requirements and CPU time. This approximation is accurate at periods
-%% of roughly 20~s and longer. At shorter periods the effect of water
-%% phases/reverberations is not taken into account, even when the flag
-%% is on.
-%% \item [{\texttt{TOPOGRAPHY}}] Set to \texttt{.true.} if topography and
-%% bathymetry should be incorporated based upon model ETOPO5 \citep{Etopo5}.
-%% This feature adds no cost to the simulation.
-%% \item [{\texttt{ATTENUATION}}] Set to \texttt{.true.} if attenuation should
-%% be incorporated. Turning this feature on increases the memory requirements
-%% significantly (roughly by a factor of~1.5), and is numerically fairly
-%% expensive. See \citet{KoTr99,KoTr02a} for a discussion on the implementation
-%% of attenuation based upon standard linear solids.
-%% \item [{\texttt{USE\_OLSEN\_ATTENUATION}}] Set to \texttt{.true.} if you
-%% want to use the attenuation model that scaled from the velocity model
-%% using Olsen's empirical relation (reference).
-%% \item [{\texttt{ABSORBING\_CONDITIONS}}] Set to \texttt{.true.} to turn
-%% on Clayton-Enquist absorbing boundary conditions (see \citet{KoTr99}).
-%% \item [{\texttt{RECORD\_LENGTH\_IN\_MINUTES}}] Choose the desired record
-%% length of the synthetic seismograms (in minutes). This controls the
-%% length of the numerical simulation, i.e., twice the record length
-%% requires twice as much CPU time. This feature is not used at the time
-%% of meshing but is required for the solver, i.e., you may change this
-%% parameter after running the mesher.
-%% \item [{\texttt{MOVIE\_SURFACE}}] Set to \texttt{.false.}, unless you want
-%% to create a movie of seismic wave propagation on the Earth's surface.
-%% Turning this option on generates large output files. See Section \ref{sec:Movies}
-%% for a discussion on the generation of movies. This feature is not
-%% used at the time of meshing but is relevant for the solver.
-%% \item [{\texttt{MOVIE\_VOLUME}}] Set to \texttt{.false.}, unless you want
-%% to create a movie of seismic wave propagation in the Earth's interior.
-%% Turning this option on generates huge output files. See Section \ref{sec:Movies}
-%% for a discussion on the generation of movies. This feature is not
-%% used at the time of meshing but is relevant for the solver.
-%% \item [{\texttt{NTSTEP\_BETWEEN\_FRAMES}}] Determines the number of timesteps
-%% between movie frames. Typically you want to save a snapshot every
-%% 100 timesteps. The smaller you make this number the more output will
-%% be generated! See Section \ref{sec:Movies} for a discussion on the
-%% generation of movies. This feature is not used at the time of meshing
-%% but is relevant for the solver.
-%% \item [{\texttt{CREATE\_SHAKEMAP}}] Set this flag to \texttt{.true.} to
-%% create a ShakeMap\textregistered, i.e., a peak ground velocity map
-%% of the maximum absolute value of the two horizontal components of the velocity vector.
-%% \item [{\texttt{SAVE\_DISPLACEMENT}}] Set this flag to \texttt{.true.}
-%% if you want to save the displacement instead of velocity for the movie
-%% frames.
-%% \item [{\texttt{USE\_HIGHRES\_FOR\_MOVIES}}] Set this flag to \texttt{.true.}
-%% if you want to save the values at all the NGLL grid points for the
-%% movie frames.
\item [{\texttt{CREATE\_ABAQUS\_FILES}}] Set this flag to \texttt{.true.} to save Abaqus FEA \url{www.simulia.com}
mesh files for subsequent viewing. Turning the flag on generates files in the \texttt{LOCAL\_PATH} directory.
See Section~\ref{sec:Mesh-graphics} for a discussion of mesh viewing features.
\item [{\texttt{CREATE\_DX\_FILES}}] Set this flag to \texttt{.true.} to save OpenDX \url{www.opendx.org}
mesh files for subsequent viewing.
-
-%% \item [{\texttt{SAVE\_MESH\_FILES}}] Set this flag to \texttt{.true.} to
-%% save AVS \url{www.avs.com}, OpenDX \url{www.opendx.org}, or ParaView \url{www.paraview.org}
-%% mesh files for subsequent viewing. Turning the flag on generates large
-%% (distributed) files in the \texttt{LOCAL\_PATH} directory. See Section~\ref{sec:Mesh-graphics}
-%% for a discussion of mesh viewing features.
\item [{\texttt{LOCAL\_PATH}}] Directory in which the databases generated
by the mesher will be written. Generally one uses a directory on the
local disk of the compute nodes, although on some machines these databases
@@ -514,26 +430,6 @@
\begin{lyxcode}
\texttt{XI\_begin} \texttt{XI\_end} \texttt{ETA\_begin} \texttt{ETA\_end} \texttt{material\_ID}
\end{lyxcode}
-
-%% \item [{\texttt{NTSTEP\_BETWEEN\_OUTPUT\_INFO}}] This parameter specifies
-%% the interval at which basic information about a run is written to
-%% the file system (\texttt{timestamp{*}} files in the \texttt{OUTPUT\_FILES}
-%% directory). If you have access to a fast machine, set \texttt{NTSTEP\_BETWEEN\_OUTPUT\_INFO}
-%% to a relatively high value (e.g., at least 100, or even 1000 or more)
-%% to avoid writing output text files too often. This feature is not
-%% used at the time of meshing. One can set this parameter to a larger
-%% value than the number of time steps to avoid writing output during
-%% the run.
-%% \item [{\texttt{NTSTEP\_BETWEEN\_OUTPUT\_SEISMOS}}] This parameter specifies
-%% the interval at which synthetic seismograms are written in the \texttt{LOCAL\_PATH}
-%% directory. If a run crashes, you may still find usable (but shorter
-%% than requested) seismograms in this directory. On a fast machine set
-%% \texttt{NTSTEP\_BETWEEN\_OUTPUT\_SEISMOS} to a relatively high value
-%% to avoid writing to the seismograms too often. This feature is not
-%% used at the time of meshing.
-%% \item [{\texttt{PRINT\_SOURCE\_TIME\_FUNCTION}}] Turn this flag on to print
-%% information about the source time function in the file \texttt{OUTPUT\_FILES/plot\_source\_time\_function.txt}.
-%% This feature is not used at the time of meshing.
\end{description}
The topography of the model is defined as a set of elevation
@@ -577,18 +473,7 @@
!~integer,~parameter~::~IMAIN~=~ISTANDARD\_OUTPUT
\end{lyxcode}
-%% Another file generated by the mesher is the header file \texttt{OUTPUT\_FILES/values\_from\_mesher.h}.
-%% This file specifies a number of constants and flags needed by the
-%% solver. These values are passed statically to the solver for reasons
-%% of speed. Some useful statistics about the mesh are also provided
-%% in this file.
-%% For a given model, set of nodes, and set of parameters in \texttt{Par\_file},
-%% one only needs to run the mesher once and for all, even if one wants
-%% to run several simulations with different sources and/or receivers
-%% (the source and receiver information is used in the solver only).
-
-
%% \section{Checking the MPI Buffers (Optional)}
%% The mesher writes MPI communication tables in the \texttt{OUTPUT\_FILES}
@@ -753,7 +638,17 @@
This feature is only relevant for the solver.
\end{description}
+%% Another file generated by \texttt{xgenerate\_databases} is the header file \texttt{OUTPUT\_FILES/values\_from\_mesher.h}.
+%% This file specifies a number of constants and flags needed by the
+%% solver. These values are passed statically to the solver for reasons
+%% of speed. Some useful statistics about the mesh are also provided
+%% in this file.
+%% For a given model, set of nodes, and set of parameters in \texttt{Par\_file},
+%% one only needs to run the mesher once and for all, even if one wants
+%% to run several simulations with different sources and/or receivers
+%% (the source and receiver information is used in the solver only).
+
\chapter{\label{cha:Running-the-Solver}Running the Solver \texttt{xspecfem3D}}
Now that you have successfully generated the databases, you are ready to compile
Copied: seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/check_buffers_2D.f90 (from rev 17240, seismo/3D/SPECFEM3D/trunk/check_buffers_2D.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/check_buffers_2D.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/check_buffers_2D.f90 2010-10-07 13:50:42 UTC (rev 17241)
@@ -0,0 +1,326 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! code to check that all the internal MPI buffers are okay along xi and eta
+! we compare the coordinates of the points in the buffers
+
+ program check_buffers_2D
+
+ implicit none
+
+ include "constants.h"
+
+ integer ithisproc,iotherproc
+
+ integer ipoin
+
+ integer npoin2d_xi_save,npoin2d_xi_mesher,npoin2d_xi
+ integer npoin2d_eta_save,npoin2d_eta_mesher,npoin2d_eta
+
+! for addressing of the slices
+ integer iproc_xi,iproc_eta,iproc
+ integer iproc_read
+ integer, dimension(:,:), allocatable :: addressing
+
+ double precision diff
+
+! 2-D addressing and buffers for summation between slices
+ integer, dimension(:), allocatable :: iboolleft_xi,iboolright_xi, &
+ iboolleft_eta,iboolright_eta
+
+! coordinates of the points to compare
+ double precision, dimension(:), allocatable :: xleft_xi,yleft_xi,zleft_xi, &
+ xright_xi,yright_xi,zright_xi,xleft_eta,yleft_eta,zleft_eta, &
+ xright_eta,yright_eta,zright_eta
+
+! parameters read from parameter file
+ integer NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT, &
+ NER_MOHO_16,NER_BOTTOM_MOHO,NEX_XI,NEX_ETA, &
+ NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,SIMULATION_TYPE
+ integer NSOURCES, NTSTEP_BETWEEN_READ_ADJSRC
+
+ double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
+ double precision DT,LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX,HDUR_MOVIE
+ double precision THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM
+
+ logical HARVARD_3D_GOCAD_MODEL,ATTENUATION,USE_OLSEN_ATTENUATION, &
+ OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL, &
+ BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS,SAVE_FORWARD
+ logical ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
+ logical MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
+ USE_HIGHRES_FOR_MOVIES,SUPPRESS_UTM_PROJECTION,USE_REGULAR_MESH
+ integer NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
+
+ character(len=256) OUTPUT_FILES,LOCAL_PATH,MODEL
+
+! parameters deduced from parameters read from file
+ integer NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
+ integer NER
+
+! now this is for all the regions
+ integer NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB
+
+! processor identification
+ character(len=256) prname,prname_other
+
+! ************** PROGRAM STARTS HERE **************
+
+ print *
+ print *,'Check all MPI buffers along xi and eta'
+ print *
+
+! read the parameter file
+ call read_parameter_file(LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX, &
+ UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
+ NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT,NER_MOHO_16,NER_BOTTOM_MOHO, &
+ NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,DT, &
+ ATTENUATION,USE_OLSEN_ATTENUATION,HARVARD_3D_GOCAD_MODEL,LOCAL_PATH,NSOURCES, &
+ THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
+ OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL,ANISOTROPY, &
+ BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS, &
+ MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
+ NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
+ SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION, &
+ NTSTEP_BETWEEN_OUTPUT_INFO,SUPPRESS_UTM_PROJECTION,MODEL,USE_REGULAR_MESH,SIMULATION_TYPE,SAVE_FORWARD, &
+ NTSTEP_BETWEEN_READ_ADJSRC)
+
+! compute other parameters based upon values read
+ call compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
+ NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM, &
+ NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,USE_REGULAR_MESH)
+
+! get the base pathname for output files
+ call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
+
+ print *
+ print *,'There are ',NPROC,' slices numbered from 0 to ',NPROC-1
+ print *,'There are ',NPROC_XI,' slices along xi'
+ print *,'There are ',NPROC_ETA,' slices along eta'
+ print *
+
+! dynamic memory allocation for arrays
+ allocate(addressing(0:NPROC_XI-1,0:NPROC_ETA-1))
+
+! open file with global slice number addressing
+ print *,'reading slice addressing'
+ open(unit=34,file=trim(OUTPUT_FILES)//'/addressing.txt',status='old',action='read')
+ do iproc = 0,NPROC-1
+ read(34,*) iproc_read,iproc_xi,iproc_eta
+ if(iproc_read /= iproc) stop 'incorrect slice number read'
+ addressing(iproc_xi,iproc_eta) = iproc
+ enddo
+ close(34)
+
+! dynamic memory allocation for arrays
+ allocate(iboolleft_xi(NPOIN2DMAX_XMIN_XMAX))
+ allocate(iboolright_xi(NPOIN2DMAX_XMIN_XMAX))
+ allocate(iboolleft_eta(NPOIN2DMAX_YMIN_YMAX))
+ allocate(iboolright_eta(NPOIN2DMAX_YMIN_YMAX))
+ allocate(xleft_xi(NPOIN2DMAX_XMIN_XMAX))
+ allocate(yleft_xi(NPOIN2DMAX_XMIN_XMAX))
+ allocate(zleft_xi(NPOIN2DMAX_XMIN_XMAX))
+ allocate(xright_xi(NPOIN2DMAX_XMIN_XMAX))
+ allocate(yright_xi(NPOIN2DMAX_XMIN_XMAX))
+ allocate(zright_xi(NPOIN2DMAX_XMIN_XMAX))
+ allocate(xleft_eta(NPOIN2DMAX_YMIN_YMAX))
+ allocate(yleft_eta(NPOIN2DMAX_YMIN_YMAX))
+ allocate(zleft_eta(NPOIN2DMAX_YMIN_YMAX))
+ allocate(xright_eta(NPOIN2DMAX_YMIN_YMAX))
+ allocate(yright_eta(NPOIN2DMAX_YMIN_YMAX))
+ allocate(zright_eta(NPOIN2DMAX_YMIN_YMAX))
+
+! double loop on NPROC_XI and NPROC_ETA
+ do iproc_eta=0,NPROC_ETA-1
+
+ print *,'checking row ',iproc_eta
+
+ do iproc_xi=0,NPROC_XI-2
+
+ print *,'checking slice ixi = ',iproc_xi,' in that row'
+
+ ithisproc = addressing(iproc_xi,iproc_eta)
+ iotherproc = addressing(iproc_xi+1,iproc_eta)
+
+! create the name for the database of the current slide
+ call create_serial_name_database(prname,ithisproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
+ call create_serial_name_database(prname_other,iotherproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
+
+! read 2-D addressing for summation between slices along xi with MPI
+
+! read iboolright_xi of this slice
+ write(*,*) 'reading MPI buffer iboolright_xi slice ',ithisproc
+ open(unit=34,file=prname(1:len_trim(prname))//'iboolright_xi.txt',status='old',action='read')
+ npoin2D_xi = 1
+ 360 continue
+ read(34,*) iboolright_xi(npoin2D_xi), &
+ xright_xi(npoin2D_xi),yright_xi(npoin2D_xi),zright_xi(npoin2D_xi)
+ if(iboolright_xi(npoin2D_xi) > 0) then
+ npoin2D_xi = npoin2D_xi + 1
+ goto 360
+ endif
+ npoin2D_xi = npoin2D_xi - 1
+ write(*,*) 'found ',npoin2D_xi,' points in iboolright_xi slice ',ithisproc
+ read(34,*) npoin2D_xi_mesher
+ if(npoin2D_xi > NPOIN2DMAX_XMIN_XMAX .or. npoin2D_xi /= npoin2D_xi_mesher) then
+ stop 'incorrect iboolright_xi read'
+ endif
+ close(34)
+
+! save to compare to other side
+ npoin2D_xi_save = npoin2D_xi
+
+! read iboolleft_xi of other slice
+ write(*,*) 'reading MPI buffer iboolleft_xi slice ',iotherproc
+ open(unit=34,file=prname_other(1:len_trim(prname_other))//'iboolleft_xi.txt',status='old',action='read')
+ npoin2D_xi = 1
+ 350 continue
+ read(34,*) iboolleft_xi(npoin2D_xi), &
+ xleft_xi(npoin2D_xi),yleft_xi(npoin2D_xi),zleft_xi(npoin2D_xi)
+ if(iboolleft_xi(npoin2D_xi) > 0) then
+ npoin2D_xi = npoin2D_xi + 1
+ goto 350
+ endif
+ npoin2D_xi = npoin2D_xi - 1
+ write(*,*) 'found ',npoin2D_xi,' points in iboolleft_xi slice ',iotherproc
+ read(34,*) npoin2D_xi_mesher
+ if(npoin2D_xi > NPOIN2DMAX_XMIN_XMAX .or. npoin2D_xi /= npoin2D_xi_mesher) then
+ stop 'incorrect iboolleft_xi read'
+ endif
+ close(34)
+
+ if(npoin2D_xi_save == npoin2D_xi) then
+ print *,'okay, same size for both buffers'
+ else
+ stop 'wrong buffer size'
+ endif
+
+! check the coordinates of all the points in the buffer
+! to see if it is correctly sorted
+ do ipoin = 1,npoin2D_xi
+ diff = dmax1(dabs(xleft_xi(ipoin)-xright_xi(ipoin)), &
+ dabs(yleft_xi(ipoin)-yright_xi(ipoin)),dabs(zleft_xi(ipoin)-zright_xi(ipoin)))
+ if(diff > 0.0000001d0) then
+ print *,'different: ',ipoin,iboolleft_xi(ipoin),iboolright_xi(ipoin),diff
+ stop 'error: different'
+ endif
+ enddo
+
+ enddo
+ enddo
+
+
+! double loop on NPROC_XI and NPROC_ETA
+ do iproc_xi=0,NPROC_XI-1
+
+ print *,'checking row ',iproc_xi
+
+ do iproc_eta=0,NPROC_ETA-2
+
+ print *,'checking slice ieta = ',iproc_eta,' in that row'
+
+ ithisproc = addressing(iproc_xi,iproc_eta)
+ iotherproc = addressing(iproc_xi,iproc_eta+1)
+
+! create the name for the database of the current slide
+ call create_serial_name_database(prname,ithisproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
+ call create_serial_name_database(prname_other,iotherproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
+
+! read 2-D addressing for summation between slices along xi with MPI
+
+! read iboolright_eta of this slice
+ write(*,*) 'reading MPI buffer iboolright_eta slice ',ithisproc
+ open(unit=34,file=prname(1:len_trim(prname))//'iboolright_eta.txt',status='old',action='read')
+ npoin2D_eta = 1
+ 460 continue
+ read(34,*) iboolright_eta(npoin2D_eta), &
+ xright_eta(npoin2D_eta),yright_eta(npoin2D_eta),zright_eta(npoin2D_eta)
+ if(iboolright_eta(npoin2D_eta) > 0) then
+ npoin2D_eta = npoin2D_eta + 1
+ goto 460
+ endif
+ npoin2D_eta = npoin2D_eta - 1
+ write(*,*) 'found ',npoin2D_eta,' points in iboolright_eta slice ',ithisproc
+ read(34,*) npoin2D_eta_mesher
+ if(npoin2D_eta > NPOIN2DMAX_YMIN_YMAX .or. npoin2D_eta /= npoin2D_eta_mesher) then
+ stop 'incorrect iboolright_eta read'
+ endif
+ close(34)
+
+! save to compare to other side
+ npoin2D_eta_save = npoin2D_eta
+
+! read iboolleft_eta of other slice
+ write(*,*) 'reading MPI buffer iboolleft_eta slice ',iotherproc
+ open(unit=34,file=prname_other(1:len_trim(prname_other))//'iboolleft_eta.txt',status='old',action='read')
+ npoin2D_eta = 1
+ 450 continue
+ read(34,*) iboolleft_eta(npoin2D_eta), &
+ xleft_eta(npoin2D_eta),yleft_eta(npoin2D_eta),zleft_eta(npoin2D_eta)
+ if(iboolleft_eta(npoin2D_eta) > 0) then
+ npoin2D_eta = npoin2D_eta + 1
+ goto 450
+ endif
+ npoin2D_eta = npoin2D_eta - 1
+ write(*,*) 'found ',npoin2D_eta,' points in iboolleft_eta slice ',iotherproc
+ read(34,*) npoin2D_eta_mesher
+ if(npoin2D_eta > NPOIN2DMAX_YMIN_YMAX .or. npoin2D_eta /= npoin2D_eta_mesher) then
+ stop 'incorrect iboolleft_eta read'
+ endif
+ close(34)
+
+ if(npoin2D_eta_save == npoin2D_eta) then
+ print *,'okay, same size for both buffers'
+ else
+ stop 'wrong buffer size'
+ endif
+
+! check the coordinates of all the points in the buffer
+! to see if it is correctly sorted
+ do ipoin = 1,npoin2D_eta
+ diff = dmax1(dabs(xleft_eta(ipoin)-xright_eta(ipoin)), &
+ dabs(yleft_eta(ipoin)-yright_eta(ipoin)),dabs(zleft_eta(ipoin)-zright_eta(ipoin)))
+ if(diff > 0.0000001d0) then
+ print *,'different: ',ipoin,iboolleft_eta(ipoin),iboolright_eta(ipoin),diff
+ stop 'error: different'
+ endif
+ enddo
+
+ enddo
+ enddo
+
+ print *
+ print *,'done'
+ print *
+
+ end program check_buffers_2D
+
Copied: seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/get_MPI_cutplanes_eta.f90 (from rev 17240, seismo/3D/SPECFEM3D/trunk/get_MPI_cutplanes_eta.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/get_MPI_cutplanes_eta.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/get_MPI_cutplanes_eta.f90 2010-10-07 13:50:42 UTC (rev 17241)
@@ -0,0 +1,179 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine get_MPI_cutplanes_eta(myrank,prname,nspec,iMPIcut_eta,ibool, &
+ xstore,ystore,zstore,mask_ibool,npointot, &
+ NSPEC2D_A_XI,NSPEC2D_B_XI)
+
+! this routine detects cut planes along eta
+! In principle the left cut plane of the first slice
+! and the right cut plane of the last slice are not used
+! in the solver except if we want to have periodic conditions
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec,myrank
+ integer NSPEC2D_A_XI,NSPEC2D_B_XI
+
+ logical iMPIcut_eta(2,nspec)
+
+ integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
+
+ double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+
+! logical mask used to create arrays iboolleft_eta and iboolright_eta
+ integer npointot
+ logical mask_ibool(npointot)
+
+! global element numbering
+ integer ispec
+
+! MPI cut-plane element numbering
+ integer ispecc1,ispecc2,npoin2D_eta,ix,iy,iz
+ integer nspec2Dtheor1,nspec2Dtheor2
+
+! processor identification
+ character(len=256) prname
+
+! theoretical number of surface elements in the buffers
+! cut planes along eta=constant correspond to XI faces
+ nspec2Dtheor1 = NSPEC2D_A_XI
+ nspec2Dtheor2 = NSPEC2D_B_XI
+
+! write the MPI buffers for the left and right edges of the slice
+! and the position of the points to check that the buffers are fine
+
+!
+! determine if the element falls on the left MPI cut plane
+!
+
+! global point number and coordinates left MPI cut-plane
+ open(unit=10,file=prname(1:len_trim(prname))//'iboolleft_eta.txt',status='unknown')
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+! nb of global points shared with the other slice
+ npoin2D_eta = 0
+
+! nb of elements in this cut-plane
+ ispecc1=0
+
+ do ispec=1,nspec
+ if(iMPIcut_eta(1,ispec)) then
+
+ ispecc1=ispecc1+1
+
+! loop on all the points in that 2-D element, including edges
+ iy = 1
+ do ix=1,NGLLX
+ 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.
+ npoin2D_eta = npoin2D_eta + 1
+
+ write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
+ endif
+
+ enddo
+ enddo
+
+ endif
+ enddo
+
+! put flag to indicate end of the list of points
+ write(10,*) '0 0 0. 0. 0.'
+
+! write total number of points
+ write(10,*) npoin2D_eta
+
+ close(10)
+
+! compare number of surface elements detected to analytical value
+ if(ispecc1 /= nspec2Dtheor1 .and. ispecc1 /= nspec2Dtheor2) &
+ call exit_MPI(myrank,'error MPI cut-planes detection in eta=left')
+
+!
+! determine if the element falls on the right MPI cut plane
+!
+
+! global point number and coordinates right MPI cut-plane
+ open(unit=10,file=prname(1:len_trim(prname))//'iboolright_eta.txt',status='unknown')
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+! nb of global points shared with the other slice
+ npoin2D_eta = 0
+
+! nb of elements in this cut-plane
+ ispecc2=0
+
+ do ispec=1,nspec
+ if(iMPIcut_eta(2,ispec)) then
+
+ ispecc2=ispecc2+1
+
+! loop on all the points in that 2-D element, including edges
+ iy = NGLLY
+ do ix=1,NGLLX
+ 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.
+ npoin2D_eta = npoin2D_eta + 1
+
+ write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
+ endif
+
+ enddo
+ enddo
+
+ endif
+ enddo
+
+! put flag to indicate end of the list of points
+ write(10,*) '0 0 0. 0. 0.'
+
+! write total number of points
+ write(10,*) npoin2D_eta
+
+ close(10)
+
+! compare number of surface elements detected to analytical value
+ if(ispecc2 /= nspec2Dtheor1 .and. ispecc2 /= nspec2Dtheor2) &
+ call exit_MPI(myrank,'error MPI cut-planes detection in eta=right')
+
+ end subroutine get_MPI_cutplanes_eta
+
Copied: seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/get_MPI_cutplanes_xi.f90 (from rev 17240, seismo/3D/SPECFEM3D/trunk/get_MPI_cutplanes_xi.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/get_MPI_cutplanes_xi.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/get_MPI_cutplanes_xi.f90 2010-10-07 13:50:42 UTC (rev 17241)
@@ -0,0 +1,178 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine get_MPI_cutplanes_xi(myrank,prname,nspec,iMPIcut_xi,ibool, &
+ xstore,ystore,zstore,mask_ibool,npointot, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA)
+
+! this routine detects cut planes along xi
+! In principle the left cut plane of the first slice
+! and the right cut plane of the last slice are not used
+! in the solver except if we want to have periodic conditions
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec,myrank
+ integer NSPEC2D_A_ETA,NSPEC2D_B_ETA
+
+ logical iMPIcut_xi(2,nspec)
+
+ integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
+
+ double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+
+! logical mask used to create arrays iboolleft_xi and iboolright_xi
+ integer npointot
+ logical mask_ibool(npointot)
+
+! global element numbering
+ integer ispec
+
+! MPI cut-plane element numbering
+ integer ispecc1,ispecc2,npoin2D_xi,ix,iy,iz
+ integer nspec2Dtheor1,nspec2Dtheor2
+
+! processor identification
+ character(len=256) prname
+
+! theoretical number of surface elements in the buffers
+! cut planes along xi=constant correspond to ETA faces
+ nspec2Dtheor1 = NSPEC2D_A_ETA
+ nspec2Dtheor2 = NSPEC2D_B_ETA
+
+! write the MPI buffers for the left and right edges of the slice
+! and the position of the points to check that the buffers are fine
+
+!
+! determine if the element falls on the left MPI cut plane
+!
+
+! global point number and coordinates left MPI cut-plane
+ open(unit=10,file=prname(1:len_trim(prname))//'iboolleft_xi.txt',status='unknown')
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+! nb of global points shared with the other slice
+ npoin2D_xi = 0
+
+! nb of elements in this cut-plane
+ ispecc1=0
+
+ do ispec=1,nspec
+ if(iMPIcut_xi(1,ispec)) then
+
+ ispecc1=ispecc1+1
+
+! loop on all the points in that 2-D element, including edges
+ ix = 1
+ do iy=1,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.
+ npoin2D_xi = npoin2D_xi + 1
+
+ write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
+ endif
+
+ enddo
+ enddo
+
+ endif
+ enddo
+
+! put flag to indicate end of the list of points
+ write(10,*) '0 0 0. 0. 0.'
+
+! write total number of points
+ write(10,*) npoin2D_xi
+
+ close(10)
+
+! compare number of surface elements detected to analytical value
+ if(ispecc1 /= nspec2Dtheor1 .and. ispecc1 /= nspec2Dtheor2) &
+ call exit_MPI(myrank,'error MPI cut-planes detection in xi=left')
+
+!
+! determine if the element falls on the right MPI cut plane
+!
+
+! global point number and coordinates right MPI cut-plane
+ open(unit=10,file=prname(1:len_trim(prname))//'iboolright_xi.txt',status='unknown')
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+! nb of global points shared with the other slice
+ npoin2D_xi = 0
+
+! nb of elements in this cut-plane
+ ispecc2=0
+
+ do ispec=1,nspec
+ if(iMPIcut_xi(2,ispec)) then
+
+ ispecc2=ispecc2+1
+
+! loop on all the points in that 2-D element, including edges
+ ix = NGLLX
+ do iy=1,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.
+ npoin2D_xi = npoin2D_xi + 1
+
+ write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
+ endif
+
+ enddo
+ enddo
+
+ endif
+ enddo
+
+! put flag to indicate end of the list of points
+ write(10,*) '0 0 0. 0. 0.'
+
+! write total number of points
+ write(10,*) npoin2D_xi
+
+ close(10)
+
+! compare number of surface elements detected to analytical value
+ if(ispecc2 /= nspec2Dtheor1 .and. ispecc2 /= nspec2Dtheor2) &
+ call exit_MPI(myrank,'error MPI cut-planes detection in xi=right')
+
+ end subroutine get_MPI_cutplanes_xi
Copied: seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/model_interface_bedrock.f90 (from rev 17240, seismo/3D/SPECFEM3D/trunk/model_interface_bedrock.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/model_interface_bedrock.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/oldstuff/model_interface_bedrock.f90 2010-10-07 13:50:42 UTC (rev 17241)
@@ -0,0 +1,390 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! interface model file
+! example file only, unused so far
+
+! ! Piero
+! module bedrock
+!
+! real,dimension(:,:),allocatable :: ibedrock
+!
+! end module bedrock
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+! subroutine model_bedrock_broadcast(myrank)
+!
+!! standard routine to setup model
+!
+! use bedrock
+!
+! implicit none
+!
+! include "constants.h"
+! ! standard include of the MPI library
+! include 'mpif.h'
+!
+! integer :: myrank
+!
+! ! local parameters
+! integer :: idummy
+!
+! ! dummy to ignore compiler warnings
+! idummy = myrank
+!
+! allocate(ibedrock(NX_TOPO_ANT,NY_TOPO_ANT))
+
+! if(myrank == 0) then
+! call read_bedrock_file(ibedrock)
+! ! write(IMAIN,*)
+! ! write(IMAIN,*) 'regional bedrock file read ranges in m from ',minval(ibedrock),' to ',maxval(ibedrock)
+! ! write(IMAIN,*)
+! endif
+
+! ! broadcast the information read on the master to the nodes
+! ! call MPI_BCAST(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT,MPI_REAL,0,MPI_COMM_WORLD,ier)
+! call bcast_all_cr(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT)
+
+! end subroutine model_bedrock_broadcast
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+!
+! subroutine read_bedrock_file()
+!
+! use bedrock
+!
+! implicit none
+!
+! include "constants.h"
+!---
+!
+! ADD YOUR MODEL HERE
+!
+!---
+!
+! end subroutine read_external_model
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+! subroutine model_bedrock_store()
+!
+! use bedrock
+!
+! implicit none
+!
+! !! DK DK store the position of the six stations to be able to
+! !! DK DK exclude circles around each station to make sure they are on the bedrock
+! !! DK DK and not in the ice
+! utm_x_station(1) = 783500.6250000d0
+! utm_y_station(1) = -11828.7519531d0
+
+! utm_x_station(2) = 853644.5000000d0
+! utm_y_station(2) = -114.0138092d0
+
+! utm_x_station(3) = 863406.0000000d0
+! utm_y_station(3) = -53736.1640625d0
+
+! utm_x_station(4) = 823398.8125000d0
+! utm_y_station(4) = 29847.4511719d0
+
+! utm_x_station(5) = 863545.3750000d0
+! utm_y_station(5) = 19669.6621094d0
+
+! utm_x_station(6) = 817099.3750000d0
+! utm_y_station(6) = -24430.2871094d0
+
+! print*,myrank,'après store the position of the six stations'
+! call flush(6)
+
+! print*, myrank,minval(nodes_coords_ext_mesh(1,:))
+! call flush(6)
+
+
+! print*, myrank,maxval(nodes_coords_ext_mesh(1,:))
+! call flush(6)
+
+
+! do ispec = 1, nspec
+
+! zmesh = zstore(2,2,2,ispec)
+
+! ! if(doubling_index == IFLAG_ONE_LAYER_TOPOGRAPHY) then
+! if(any(ibelm_top == ispec)) then
+! doubling_value_found_for_Piero = IFLAG_ONE_LAYER_TOPOGRAPHY
+
+! else if(zmesh < Z_23p4km) then
+! doubling_value_found_for_Piero = IFLAG_MANTLE_BELOW_23p4km
+
+! else if(zmesh < Z_14km) then
+! doubling_value_found_for_Piero = IFLAG_14km_to_23p4km
+
+! else
+! doubling_value_found_for_Piero = IFLAG_BEDROCK_down_to_14km
+! endif
+! idoubling(ispec) = doubling_value_found_for_Piero
+
+! do k = 1, NGLLZ
+! do j = 1, NGLLY
+! do i = 1, NGLLX
+
+
+! if(idoubling(ispec) == IFLAG_ONE_LAYER_TOPOGRAPHY .or. &
+! idoubling(ispec) == IFLAG_BEDROCK_down_to_14km) then
+
+! ! since we have suppressed UTM projection for Piero Basini, UTMx is the same as long
+! ! and UTMy is the same as lat
+! long = xstore(i,j,k,ispec)
+! lat = ystore(i,j,k,ispec)
+
+! ! get coordinate of corner in model
+! icornerlong = int((long - ORIG_LONG_TOPO) / DEGREES_PER_CELL_TOPO) + 1
+! icornerlat = int((lat - ORIG_LAT_TOPO) / DEGREES_PER_CELL_TOPO) + 1
+
+! ! avoid edge effects and extend with identical point if outside model
+! if(icornerlong < 1) icornerlong = 1
+! if(icornerlong > NX_TOPO-1) icornerlong = NX_TOPO-1
+! if(icornerlat < 1) icornerlat = 1
+! if(icornerlat > NY_TOPO-1) icornerlat = NY_TOPO-1
+
+! ! compute coordinates of corner
+! long_corner = ORIG_LONG_TOPO + (icornerlong-1)*DEGREES_PER_CELL_TOPO
+! lat_corner = ORIG_LAT_TOPO + (icornerlat-1)*DEGREES_PER_CELL_TOPO
+
+! ! compute ratio for interpolation
+! ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO
+! ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO
+
+! ! avoid edge effects
+! if(ratio_xi < 0.) ratio_xi = 0.
+! if(ratio_xi > 1.) ratio_xi = 1.
+! if(ratio_eta < 0.) ratio_eta = 0.
+! if(ratio_eta > 1.) ratio_eta = 1.
+
+! ! interpolate elevation at current point
+! elevation_bedrock = &
+! ibedrock(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
+! ibedrock(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
+! ibedrock(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
+! ibedrock(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
+
+! !! DK DK exclude circles around each station to make sure they are on the bedrock
+! !! DK DK and not in the ice
+! is_around_a_station = .false.
+! do istation = 1,NUMBER_OF_STATIONS
+! if(sqrt((long - utm_x_station(istation))**2 + (lat - utm_y_station(istation))**2) < RADIUS_TO_EXCLUDE) then
+! is_around_a_station = .true.
+! exit
+! endif
+! enddo
+
+! ! define elastic parameters in the model
+
+! ! we are above the bedrock interface i.e. in the ice, and not too close to a station
+! if(zmesh >= elevation_bedrock .and. .not. is_around_a_station) then
+! vp = 3800.d0
+! vs = 1900.d0
+! rho = 900.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_ICE
+
+! ! we are below the bedrock interface i.e. in the bedrock, or close to a station
+! else
+! vp = 5800.d0
+! vs = 3200.d0
+! rho = 2600.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
+! endif
+
+! else if(idoubling(ispec) == IFLAG_14km_to_23p4km) then
+! vp = 6800.d0
+! vs = 3900.d0
+! rho = 2900.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
+
+! else if(idoubling(ispec) == IFLAG_MANTLE_BELOW_23p4km) then
+! vp = 8100.d0
+! vs = 4480.d0
+! rho = 3380.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
+
+! endif
+
+! !pll 8/06
+! if(CUSTOM_REAL == SIZE_REAL) then
+! rhostore(i,j,k,ispec) = sngl(rho)
+! vpstore(i,j,k,ispec) = sngl(vp)
+! vsstore(i,j,k,ispec) = sngl(vs)
+! else
+! rhostore(i,j,k,ispec) = rho
+! vpstore(i,j,k,ispec) = vp
+! vsstore(i,j,k,ispec) = vs
+! end if
+
+! kappastore(i,j,k,ispec) = rhostore(i,j,k,ispec)*(vpstore(i,j,k,ispec)*vpstore(i,j,k,ispec) - &
+! 4.d0*vsstore(i,j,k,ispec)*vsstore(i,j,k,ispec)/3.d0)
+! mustore(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)*&
+! vsstore(i,j,k,ispec)
+
+! ! Stacey, a completer par la suite
+! rho_vp(i,j,k,ispec) = rhostore(i,j,k,ispec)*vpstore(i,j,k,ispec)
+! rho_vs(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)
+! !end pll
+
+! ! kappastore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
+! ! (materials_ext_mesh(2,mat_ext_mesh(ispec))*materials_ext_mesh(2,mat_ext_mesh(ispec)) - &
+! ! 4.d0*materials_ext_mesh(3,mat_ext_mesh(ispec))*materials_ext_mesh(3,mat_ext_mesh(ispec))/3.d0)
+! ! mustore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
+! materials_ext_mesh(3,mat_ext_mesh(ispec))*&
+! ! x materials_ext_mesh(3,mat_ext_mesh(ispec))
+! enddo
+! enddo
+! enddo
+! enddo
+!
+! end subroutine
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+!pll
+! subroutine interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
+
+! implicit none
+
+! include "constants.h"
+
+! integer :: iflag,flag_below,flag_above
+! integer :: ispec,nspec
+! integer :: i,j,k
+! double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
+! real(kind=CUSTOM_REAL), dimension(NX_TOPO_ANT,NY_TOPO_ANT) :: ibedrock
+! integer, parameter :: NUMBER_OF_STATIONS = 1
+! double precision, parameter :: RADIUS_TO_EXCLUDE = 250.d0
+! double precision, dimension(NUMBER_OF_STATIONS) :: utm_x_station,utm_y_station
+
+! !-------------------
+
+! !for Piero
+! logical :: is_around_a_station
+! integer :: istation
+
+! ! store bedrock values
+! integer :: icornerlat,icornerlong
+! double precision :: lat,long,elevation_bedrock
+! double precision :: lat_corner,long_corner,ratio_xi,ratio_eta
+
+
+! !! DK DK store the position of the six stations to be able to
+! !! DK DK exclude circles around each station to make sure they are on the bedrock
+! !! DK DK and not in the ice
+! utm_x_station(1) = 783500.6250000d0
+! utm_y_station(1) = -11828.7519531d0
+
+! utm_x_station(2) = 853644.5000000d0
+! utm_y_station(2) = -114.0138092d0
+
+! utm_x_station(3) = 863406.0000000d0
+! utm_y_station(3) = -53736.1640625d0
+
+! utm_x_station(4) = 823398.8125000d0
+! utm_y_station(4) = 29847.4511719d0
+
+! utm_x_station(5) = 863545.3750000d0
+! utm_y_station(5) = 19669.6621094d0
+
+! utm_x_station(6) = 817099.3750000d0
+! utm_y_station(6) = -24430.2871094d0
+
+! ! since we have suppressed UTM projection for Piero Basini, UTMx is the same as long
+! ! and UTMy is the same as lat
+! long = xstore(i,j,k,ispec)
+! lat = ystore(i,j,k,ispec)
+
+! ! get coordinate of corner in model
+! icornerlong = int((long - ORIG_LONG_TOPO_ANT) / DEGREES_PER_CELL_TOPO_ANT) + 1
+! icornerlat = int((lat - ORIG_LAT_TOPO_ANT) / DEGREES_PER_CELL_TOPO_ANT) + 1
+
+! ! avoid edge effects and extend with identical point if outside model
+! if(icornerlong < 1) icornerlong = 1
+! if(icornerlong > NX_TOPO_ANT-1) icornerlong = NX_TOPO_ANT-1
+! if(icornerlat < 1) icornerlat = 1
+! if(icornerlat > NY_TOPO_ANT-1) icornerlat = NY_TOPO_ANT-1
+
+! ! compute coordinates of corner
+! long_corner = ORIG_LONG_TOPO_ANT + (icornerlong-1)*DEGREES_PER_CELL_TOPO_ANT
+! lat_corner = ORIG_LAT_TOPO_ANT + (icornerlat-1)*DEGREES_PER_CELL_TOPO_ANT
+
+! ! compute ratio for interpolation
+! ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO_ANT
+! ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO_ANT
+
+! ! avoid edge effects
+! if(ratio_xi < 0.) ratio_xi = 0.
+! if(ratio_xi > 1.) ratio_xi = 1.
+! if(ratio_eta < 0.) ratio_eta = 0.
+! if(ratio_eta > 1.) ratio_eta = 1.
+
+! ! interpolate elevation at current point
+! elevation_bedrock = &
+! ibedrock(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
+! ibedrock(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
+! ibedrock(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
+! ibedrock(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
+
+! !! DK DK exclude circles around each station to make sure they are on the bedrock
+! !! DK DK and not in the ice
+! is_around_a_station = .false.
+! do istation = 1,NUMBER_OF_STATIONS
+! if(sqrt((xstore(i,j,k,ispec) - utm_x_station(istation))**2 + (ystore(i,j,k,ispec) - &
+! utm_y_station(istation))**2) < RADIUS_TO_EXCLUDE) then
+! is_around_a_station = .true.
+! exit
+! endif
+! enddo
+
+! ! we are above the bedrock interface i.e. in the ice, and not too close to a station
+! if(zstore(i,j,k,ispec) >= elevation_bedrock .and. .not. is_around_a_station) then
+! iflag = flag_above
+! !iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_ICE
+! ! we are below the bedrock interface i.e. in the bedrock, or close to a station
+! else
+! iflag = flag_below
+! !iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
+! endif
+
+
+! end subroutine interface
Deleted: seismo/3D/SPECFEM3D/trunk/check_buffers_2D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/check_buffers_2D.f90 2010-10-07 10:09:09 UTC (rev 17240)
+++ seismo/3D/SPECFEM3D/trunk/check_buffers_2D.f90 2010-10-07 13:50:42 UTC (rev 17241)
@@ -1,326 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
-! code to check that all the internal MPI buffers are okay along xi and eta
-! we compare the coordinates of the points in the buffers
-
- program check_buffers_2D
-
- implicit none
-
- include "constants.h"
-
- integer ithisproc,iotherproc
-
- integer ipoin
-
- integer npoin2d_xi_save,npoin2d_xi_mesher,npoin2d_xi
- integer npoin2d_eta_save,npoin2d_eta_mesher,npoin2d_eta
-
-! for addressing of the slices
- integer iproc_xi,iproc_eta,iproc
- integer iproc_read
- integer, dimension(:,:), allocatable :: addressing
-
- double precision diff
-
-! 2-D addressing and buffers for summation between slices
- integer, dimension(:), allocatable :: iboolleft_xi,iboolright_xi, &
- iboolleft_eta,iboolright_eta
-
-! coordinates of the points to compare
- double precision, dimension(:), allocatable :: xleft_xi,yleft_xi,zleft_xi, &
- xright_xi,yright_xi,zright_xi,xleft_eta,yleft_eta,zleft_eta, &
- xright_eta,yright_eta,zright_eta
-
-! parameters read from parameter file
- integer NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT, &
- NER_MOHO_16,NER_BOTTOM_MOHO,NEX_XI,NEX_ETA, &
- NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,SIMULATION_TYPE
- integer NSOURCES, NTSTEP_BETWEEN_READ_ADJSRC
-
- double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
- double precision DT,LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX,HDUR_MOVIE
- double precision THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM
-
- logical HARVARD_3D_GOCAD_MODEL,ATTENUATION,USE_OLSEN_ATTENUATION, &
- OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL, &
- BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS,SAVE_FORWARD
- logical ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
- logical MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
- USE_HIGHRES_FOR_MOVIES,SUPPRESS_UTM_PROJECTION,USE_REGULAR_MESH
- integer NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
-
- character(len=256) OUTPUT_FILES,LOCAL_PATH,MODEL
-
-! parameters deduced from parameters read from file
- integer NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
- integer NER
-
-! now this is for all the regions
- integer NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
- NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB
-
-! processor identification
- character(len=256) prname,prname_other
-
-! ************** PROGRAM STARTS HERE **************
-
- print *
- print *,'Check all MPI buffers along xi and eta'
- print *
-
-! read the parameter file
- call read_parameter_file(LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX, &
- UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
- NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT,NER_MOHO_16,NER_BOTTOM_MOHO, &
- NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,DT, &
- ATTENUATION,USE_OLSEN_ATTENUATION,HARVARD_3D_GOCAD_MODEL,LOCAL_PATH,NSOURCES, &
- THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
- OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL,ANISOTROPY, &
- BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS, &
- MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
- NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
- SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION, &
- NTSTEP_BETWEEN_OUTPUT_INFO,SUPPRESS_UTM_PROJECTION,MODEL,USE_REGULAR_MESH,SIMULATION_TYPE,SAVE_FORWARD, &
- NTSTEP_BETWEEN_READ_ADJSRC)
-
-! compute other parameters based upon values read
- call compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
- NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
- NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM, &
- NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,USE_REGULAR_MESH)
-
-! get the base pathname for output files
- call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
-
- print *
- print *,'There are ',NPROC,' slices numbered from 0 to ',NPROC-1
- print *,'There are ',NPROC_XI,' slices along xi'
- print *,'There are ',NPROC_ETA,' slices along eta'
- print *
-
-! dynamic memory allocation for arrays
- allocate(addressing(0:NPROC_XI-1,0:NPROC_ETA-1))
-
-! open file with global slice number addressing
- print *,'reading slice addressing'
- open(unit=34,file=trim(OUTPUT_FILES)//'/addressing.txt',status='old',action='read')
- do iproc = 0,NPROC-1
- read(34,*) iproc_read,iproc_xi,iproc_eta
- if(iproc_read /= iproc) stop 'incorrect slice number read'
- addressing(iproc_xi,iproc_eta) = iproc
- enddo
- close(34)
-
-! dynamic memory allocation for arrays
- allocate(iboolleft_xi(NPOIN2DMAX_XMIN_XMAX))
- allocate(iboolright_xi(NPOIN2DMAX_XMIN_XMAX))
- allocate(iboolleft_eta(NPOIN2DMAX_YMIN_YMAX))
- allocate(iboolright_eta(NPOIN2DMAX_YMIN_YMAX))
- allocate(xleft_xi(NPOIN2DMAX_XMIN_XMAX))
- allocate(yleft_xi(NPOIN2DMAX_XMIN_XMAX))
- allocate(zleft_xi(NPOIN2DMAX_XMIN_XMAX))
- allocate(xright_xi(NPOIN2DMAX_XMIN_XMAX))
- allocate(yright_xi(NPOIN2DMAX_XMIN_XMAX))
- allocate(zright_xi(NPOIN2DMAX_XMIN_XMAX))
- allocate(xleft_eta(NPOIN2DMAX_YMIN_YMAX))
- allocate(yleft_eta(NPOIN2DMAX_YMIN_YMAX))
- allocate(zleft_eta(NPOIN2DMAX_YMIN_YMAX))
- allocate(xright_eta(NPOIN2DMAX_YMIN_YMAX))
- allocate(yright_eta(NPOIN2DMAX_YMIN_YMAX))
- allocate(zright_eta(NPOIN2DMAX_YMIN_YMAX))
-
-! double loop on NPROC_XI and NPROC_ETA
- do iproc_eta=0,NPROC_ETA-1
-
- print *,'checking row ',iproc_eta
-
- do iproc_xi=0,NPROC_XI-2
-
- print *,'checking slice ixi = ',iproc_xi,' in that row'
-
- ithisproc = addressing(iproc_xi,iproc_eta)
- iotherproc = addressing(iproc_xi+1,iproc_eta)
-
-! create the name for the database of the current slide
- call create_serial_name_database(prname,ithisproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
- call create_serial_name_database(prname_other,iotherproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
-
-! read 2-D addressing for summation between slices along xi with MPI
-
-! read iboolright_xi of this slice
- write(*,*) 'reading MPI buffer iboolright_xi slice ',ithisproc
- open(unit=34,file=prname(1:len_trim(prname))//'iboolright_xi.txt',status='old',action='read')
- npoin2D_xi = 1
- 360 continue
- read(34,*) iboolright_xi(npoin2D_xi), &
- xright_xi(npoin2D_xi),yright_xi(npoin2D_xi),zright_xi(npoin2D_xi)
- if(iboolright_xi(npoin2D_xi) > 0) then
- npoin2D_xi = npoin2D_xi + 1
- goto 360
- endif
- npoin2D_xi = npoin2D_xi - 1
- write(*,*) 'found ',npoin2D_xi,' points in iboolright_xi slice ',ithisproc
- read(34,*) npoin2D_xi_mesher
- if(npoin2D_xi > NPOIN2DMAX_XMIN_XMAX .or. npoin2D_xi /= npoin2D_xi_mesher) then
- stop 'incorrect iboolright_xi read'
- endif
- close(34)
-
-! save to compare to other side
- npoin2D_xi_save = npoin2D_xi
-
-! read iboolleft_xi of other slice
- write(*,*) 'reading MPI buffer iboolleft_xi slice ',iotherproc
- open(unit=34,file=prname_other(1:len_trim(prname_other))//'iboolleft_xi.txt',status='old',action='read')
- npoin2D_xi = 1
- 350 continue
- read(34,*) iboolleft_xi(npoin2D_xi), &
- xleft_xi(npoin2D_xi),yleft_xi(npoin2D_xi),zleft_xi(npoin2D_xi)
- if(iboolleft_xi(npoin2D_xi) > 0) then
- npoin2D_xi = npoin2D_xi + 1
- goto 350
- endif
- npoin2D_xi = npoin2D_xi - 1
- write(*,*) 'found ',npoin2D_xi,' points in iboolleft_xi slice ',iotherproc
- read(34,*) npoin2D_xi_mesher
- if(npoin2D_xi > NPOIN2DMAX_XMIN_XMAX .or. npoin2D_xi /= npoin2D_xi_mesher) then
- stop 'incorrect iboolleft_xi read'
- endif
- close(34)
-
- if(npoin2D_xi_save == npoin2D_xi) then
- print *,'okay, same size for both buffers'
- else
- stop 'wrong buffer size'
- endif
-
-! check the coordinates of all the points in the buffer
-! to see if it is correctly sorted
- do ipoin = 1,npoin2D_xi
- diff = dmax1(dabs(xleft_xi(ipoin)-xright_xi(ipoin)), &
- dabs(yleft_xi(ipoin)-yright_xi(ipoin)),dabs(zleft_xi(ipoin)-zright_xi(ipoin)))
- if(diff > 0.0000001d0) then
- print *,'different: ',ipoin,iboolleft_xi(ipoin),iboolright_xi(ipoin),diff
- stop 'error: different'
- endif
- enddo
-
- enddo
- enddo
-
-
-! double loop on NPROC_XI and NPROC_ETA
- do iproc_xi=0,NPROC_XI-1
-
- print *,'checking row ',iproc_xi
-
- do iproc_eta=0,NPROC_ETA-2
-
- print *,'checking slice ieta = ',iproc_eta,' in that row'
-
- ithisproc = addressing(iproc_xi,iproc_eta)
- iotherproc = addressing(iproc_xi,iproc_eta+1)
-
-! create the name for the database of the current slide
- call create_serial_name_database(prname,ithisproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
- call create_serial_name_database(prname_other,iotherproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
-
-! read 2-D addressing for summation between slices along xi with MPI
-
-! read iboolright_eta of this slice
- write(*,*) 'reading MPI buffer iboolright_eta slice ',ithisproc
- open(unit=34,file=prname(1:len_trim(prname))//'iboolright_eta.txt',status='old',action='read')
- npoin2D_eta = 1
- 460 continue
- read(34,*) iboolright_eta(npoin2D_eta), &
- xright_eta(npoin2D_eta),yright_eta(npoin2D_eta),zright_eta(npoin2D_eta)
- if(iboolright_eta(npoin2D_eta) > 0) then
- npoin2D_eta = npoin2D_eta + 1
- goto 460
- endif
- npoin2D_eta = npoin2D_eta - 1
- write(*,*) 'found ',npoin2D_eta,' points in iboolright_eta slice ',ithisproc
- read(34,*) npoin2D_eta_mesher
- if(npoin2D_eta > NPOIN2DMAX_YMIN_YMAX .or. npoin2D_eta /= npoin2D_eta_mesher) then
- stop 'incorrect iboolright_eta read'
- endif
- close(34)
-
-! save to compare to other side
- npoin2D_eta_save = npoin2D_eta
-
-! read iboolleft_eta of other slice
- write(*,*) 'reading MPI buffer iboolleft_eta slice ',iotherproc
- open(unit=34,file=prname_other(1:len_trim(prname_other))//'iboolleft_eta.txt',status='old',action='read')
- npoin2D_eta = 1
- 450 continue
- read(34,*) iboolleft_eta(npoin2D_eta), &
- xleft_eta(npoin2D_eta),yleft_eta(npoin2D_eta),zleft_eta(npoin2D_eta)
- if(iboolleft_eta(npoin2D_eta) > 0) then
- npoin2D_eta = npoin2D_eta + 1
- goto 450
- endif
- npoin2D_eta = npoin2D_eta - 1
- write(*,*) 'found ',npoin2D_eta,' points in iboolleft_eta slice ',iotherproc
- read(34,*) npoin2D_eta_mesher
- if(npoin2D_eta > NPOIN2DMAX_YMIN_YMAX .or. npoin2D_eta /= npoin2D_eta_mesher) then
- stop 'incorrect iboolleft_eta read'
- endif
- close(34)
-
- if(npoin2D_eta_save == npoin2D_eta) then
- print *,'okay, same size for both buffers'
- else
- stop 'wrong buffer size'
- endif
-
-! check the coordinates of all the points in the buffer
-! to see if it is correctly sorted
- do ipoin = 1,npoin2D_eta
- diff = dmax1(dabs(xleft_eta(ipoin)-xright_eta(ipoin)), &
- dabs(yleft_eta(ipoin)-yright_eta(ipoin)),dabs(zleft_eta(ipoin)-zright_eta(ipoin)))
- if(diff > 0.0000001d0) then
- print *,'different: ',ipoin,iboolleft_eta(ipoin),iboolright_eta(ipoin),diff
- stop 'error: different'
- endif
- enddo
-
- enddo
- enddo
-
- print *
- print *,'done'
- print *
-
- end program check_buffers_2D
-
Modified: seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90 2010-10-07 10:09:09 UTC (rev 17240)
+++ seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90 2010-10-07 13:50:42 UTC (rev 17241)
@@ -77,11 +77,11 @@
integer:: nadj_rec_local
real(kind=CUSTOM_REAL),dimension(NGLOB_ADJOINT):: b_potential_dot_dot_acoustic
!<YANGL
+ logical :: ibool_read_adj_arrays
+ integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC
! real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
- integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC
real(kind=CUSTOM_REAL),dimension(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearray
- logical :: ibool_read_adj_arrays
!>YANGL
! local parameters
@@ -272,7 +272,9 @@
! potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
! - adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j,k) / kappastore(i,j,k,ispec)
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - adj_sourcearrays(irec_local,NTSTEP_BETWEEN_READ_ADJSRC-mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC),1,i,j,k) / kappastore(i,j,k,ispec)
+ - adj_sourcearrays(irec_local,NTSTEP_BETWEEN_READ_ADJSRC &
+ -mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC),1,i,j,k) &
+ / kappastore(i,j,k,ispec)
!>YANGL
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90 2010-10-07 10:09:09 UTC (rev 17240)
+++ seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90 2010-10-07 13:50:42 UTC (rev 17241)
@@ -75,10 +75,10 @@
real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_ADJOINT):: b_accel
!<YANGL
! real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
+ logical :: ibool_read_adj_arrays
+ integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC
real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
- integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC
real(kind=CUSTOM_REAL),dimension(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearray
- logical :: ibool_read_adj_arrays
!>YANGL
! local parameters
@@ -243,7 +243,8 @@
iglob = ibool(i,j,k,ispec_selected_rec(irec))
!<YANGL
! accel(:,iglob) = accel(:,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,:,i,j,k)
- accel(:,iglob) = accel(:,iglob) + adj_sourcearrays(irec_local,NTSTEP_BETWEEN_READ_ADJSRC-mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC),:,i,j,k)
+ accel(:,iglob) = accel(:,iglob) + &
+ adj_sourcearrays(irec_local,NTSTEP_BETWEEN_READ_ADJSRC-mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC),:,i,j,k)
!>YANGL
enddo
enddo
Deleted: seismo/3D/SPECFEM3D/trunk/get_MPI_cutplanes_eta.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/get_MPI_cutplanes_eta.f90 2010-10-07 10:09:09 UTC (rev 17240)
+++ seismo/3D/SPECFEM3D/trunk/get_MPI_cutplanes_eta.f90 2010-10-07 13:50:42 UTC (rev 17241)
@@ -1,179 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine get_MPI_cutplanes_eta(myrank,prname,nspec,iMPIcut_eta,ibool, &
- xstore,ystore,zstore,mask_ibool,npointot, &
- NSPEC2D_A_XI,NSPEC2D_B_XI)
-
-! this routine detects cut planes along eta
-! In principle the left cut plane of the first slice
-! and the right cut plane of the last slice are not used
-! in the solver except if we want to have periodic conditions
-
- implicit none
-
- include "constants.h"
-
- integer nspec,myrank
- integer NSPEC2D_A_XI,NSPEC2D_B_XI
-
- logical iMPIcut_eta(2,nspec)
-
- integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
-
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
-
-! logical mask used to create arrays iboolleft_eta and iboolright_eta
- integer npointot
- logical mask_ibool(npointot)
-
-! global element numbering
- integer ispec
-
-! MPI cut-plane element numbering
- integer ispecc1,ispecc2,npoin2D_eta,ix,iy,iz
- integer nspec2Dtheor1,nspec2Dtheor2
-
-! processor identification
- character(len=256) prname
-
-! theoretical number of surface elements in the buffers
-! cut planes along eta=constant correspond to XI faces
- nspec2Dtheor1 = NSPEC2D_A_XI
- nspec2Dtheor2 = NSPEC2D_B_XI
-
-! write the MPI buffers for the left and right edges of the slice
-! and the position of the points to check that the buffers are fine
-
-!
-! determine if the element falls on the left MPI cut plane
-!
-
-! global point number and coordinates left MPI cut-plane
- open(unit=10,file=prname(1:len_trim(prname))//'iboolleft_eta.txt',status='unknown')
-
-! erase the logical mask used to mark points already found
- mask_ibool(:) = .false.
-
-! nb of global points shared with the other slice
- npoin2D_eta = 0
-
-! nb of elements in this cut-plane
- ispecc1=0
-
- do ispec=1,nspec
- if(iMPIcut_eta(1,ispec)) then
-
- ispecc1=ispecc1+1
-
-! loop on all the points in that 2-D element, including edges
- iy = 1
- do ix=1,NGLLX
- 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.
- npoin2D_eta = npoin2D_eta + 1
-
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
-
- enddo
- enddo
-
- endif
- enddo
-
-! put flag to indicate end of the list of points
- write(10,*) '0 0 0. 0. 0.'
-
-! write total number of points
- write(10,*) npoin2D_eta
-
- close(10)
-
-! compare number of surface elements detected to analytical value
- if(ispecc1 /= nspec2Dtheor1 .and. ispecc1 /= nspec2Dtheor2) &
- call exit_MPI(myrank,'error MPI cut-planes detection in eta=left')
-
-!
-! determine if the element falls on the right MPI cut plane
-!
-
-! global point number and coordinates right MPI cut-plane
- open(unit=10,file=prname(1:len_trim(prname))//'iboolright_eta.txt',status='unknown')
-
-! erase the logical mask used to mark points already found
- mask_ibool(:) = .false.
-
-! nb of global points shared with the other slice
- npoin2D_eta = 0
-
-! nb of elements in this cut-plane
- ispecc2=0
-
- do ispec=1,nspec
- if(iMPIcut_eta(2,ispec)) then
-
- ispecc2=ispecc2+1
-
-! loop on all the points in that 2-D element, including edges
- iy = NGLLY
- do ix=1,NGLLX
- 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.
- npoin2D_eta = npoin2D_eta + 1
-
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
-
- enddo
- enddo
-
- endif
- enddo
-
-! put flag to indicate end of the list of points
- write(10,*) '0 0 0. 0. 0.'
-
-! write total number of points
- write(10,*) npoin2D_eta
-
- close(10)
-
-! compare number of surface elements detected to analytical value
- if(ispecc2 /= nspec2Dtheor1 .and. ispecc2 /= nspec2Dtheor2) &
- call exit_MPI(myrank,'error MPI cut-planes detection in eta=right')
-
- end subroutine get_MPI_cutplanes_eta
-
Deleted: seismo/3D/SPECFEM3D/trunk/get_MPI_cutplanes_xi.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/get_MPI_cutplanes_xi.f90 2010-10-07 10:09:09 UTC (rev 17240)
+++ seismo/3D/SPECFEM3D/trunk/get_MPI_cutplanes_xi.f90 2010-10-07 13:50:42 UTC (rev 17241)
@@ -1,178 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine get_MPI_cutplanes_xi(myrank,prname,nspec,iMPIcut_xi,ibool, &
- xstore,ystore,zstore,mask_ibool,npointot, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA)
-
-! this routine detects cut planes along xi
-! In principle the left cut plane of the first slice
-! and the right cut plane of the last slice are not used
-! in the solver except if we want to have periodic conditions
-
- implicit none
-
- include "constants.h"
-
- integer nspec,myrank
- integer NSPEC2D_A_ETA,NSPEC2D_B_ETA
-
- logical iMPIcut_xi(2,nspec)
-
- integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
-
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
-
-! logical mask used to create arrays iboolleft_xi and iboolright_xi
- integer npointot
- logical mask_ibool(npointot)
-
-! global element numbering
- integer ispec
-
-! MPI cut-plane element numbering
- integer ispecc1,ispecc2,npoin2D_xi,ix,iy,iz
- integer nspec2Dtheor1,nspec2Dtheor2
-
-! processor identification
- character(len=256) prname
-
-! theoretical number of surface elements in the buffers
-! cut planes along xi=constant correspond to ETA faces
- nspec2Dtheor1 = NSPEC2D_A_ETA
- nspec2Dtheor2 = NSPEC2D_B_ETA
-
-! write the MPI buffers for the left and right edges of the slice
-! and the position of the points to check that the buffers are fine
-
-!
-! determine if the element falls on the left MPI cut plane
-!
-
-! global point number and coordinates left MPI cut-plane
- open(unit=10,file=prname(1:len_trim(prname))//'iboolleft_xi.txt',status='unknown')
-
-! erase the logical mask used to mark points already found
- mask_ibool(:) = .false.
-
-! nb of global points shared with the other slice
- npoin2D_xi = 0
-
-! nb of elements in this cut-plane
- ispecc1=0
-
- do ispec=1,nspec
- if(iMPIcut_xi(1,ispec)) then
-
- ispecc1=ispecc1+1
-
-! loop on all the points in that 2-D element, including edges
- ix = 1
- do iy=1,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.
- npoin2D_xi = npoin2D_xi + 1
-
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
-
- enddo
- enddo
-
- endif
- enddo
-
-! put flag to indicate end of the list of points
- write(10,*) '0 0 0. 0. 0.'
-
-! write total number of points
- write(10,*) npoin2D_xi
-
- close(10)
-
-! compare number of surface elements detected to analytical value
- if(ispecc1 /= nspec2Dtheor1 .and. ispecc1 /= nspec2Dtheor2) &
- call exit_MPI(myrank,'error MPI cut-planes detection in xi=left')
-
-!
-! determine if the element falls on the right MPI cut plane
-!
-
-! global point number and coordinates right MPI cut-plane
- open(unit=10,file=prname(1:len_trim(prname))//'iboolright_xi.txt',status='unknown')
-
-! erase the logical mask used to mark points already found
- mask_ibool(:) = .false.
-
-! nb of global points shared with the other slice
- npoin2D_xi = 0
-
-! nb of elements in this cut-plane
- ispecc2=0
-
- do ispec=1,nspec
- if(iMPIcut_xi(2,ispec)) then
-
- ispecc2=ispecc2+1
-
-! loop on all the points in that 2-D element, including edges
- ix = NGLLX
- do iy=1,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.
- npoin2D_xi = npoin2D_xi + 1
-
- write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
-
- enddo
- enddo
-
- endif
- enddo
-
-! put flag to indicate end of the list of points
- write(10,*) '0 0 0. 0. 0.'
-
-! write total number of points
- write(10,*) npoin2D_xi
-
- close(10)
-
-! compare number of surface elements detected to analytical value
- if(ispecc2 /= nspec2Dtheor1 .and. ispecc2 /= nspec2Dtheor2) &
- call exit_MPI(myrank,'error MPI cut-planes detection in xi=right')
-
- end subroutine get_MPI_cutplanes_xi
Deleted: seismo/3D/SPECFEM3D/trunk/model_interface_bedrock.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/model_interface_bedrock.f90 2010-10-07 10:09:09 UTC (rev 17240)
+++ seismo/3D/SPECFEM3D/trunk/model_interface_bedrock.f90 2010-10-07 13:50:42 UTC (rev 17241)
@@ -1,390 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
-! interface model file
-! example file only, unused so far
-
-! ! Piero
-! module bedrock
-!
-! real,dimension(:,:),allocatable :: ibedrock
-!
-! end module bedrock
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-! subroutine model_bedrock_broadcast(myrank)
-!
-!! standard routine to setup model
-!
-! use bedrock
-!
-! implicit none
-!
-! include "constants.h"
-! ! standard include of the MPI library
-! include 'mpif.h'
-!
-! integer :: myrank
-!
-! ! local parameters
-! integer :: idummy
-!
-! ! dummy to ignore compiler warnings
-! idummy = myrank
-!
-! allocate(ibedrock(NX_TOPO_ANT,NY_TOPO_ANT))
-
-! if(myrank == 0) then
-! call read_bedrock_file(ibedrock)
-! ! write(IMAIN,*)
-! ! write(IMAIN,*) 'regional bedrock file read ranges in m from ',minval(ibedrock),' to ',maxval(ibedrock)
-! ! write(IMAIN,*)
-! endif
-
-! ! broadcast the information read on the master to the nodes
-! ! call MPI_BCAST(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT,MPI_REAL,0,MPI_COMM_WORLD,ier)
-! call bcast_all_cr(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT)
-
-! end subroutine model_bedrock_broadcast
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-!
-! subroutine read_bedrock_file()
-!
-! use bedrock
-!
-! implicit none
-!
-! include "constants.h"
-!---
-!
-! ADD YOUR MODEL HERE
-!
-!---
-!
-! end subroutine read_external_model
-
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-
-! subroutine model_bedrock_store()
-!
-! use bedrock
-!
-! implicit none
-!
-! !! DK DK store the position of the six stations to be able to
-! !! DK DK exclude circles around each station to make sure they are on the bedrock
-! !! DK DK and not in the ice
-! utm_x_station(1) = 783500.6250000d0
-! utm_y_station(1) = -11828.7519531d0
-
-! utm_x_station(2) = 853644.5000000d0
-! utm_y_station(2) = -114.0138092d0
-
-! utm_x_station(3) = 863406.0000000d0
-! utm_y_station(3) = -53736.1640625d0
-
-! utm_x_station(4) = 823398.8125000d0
-! utm_y_station(4) = 29847.4511719d0
-
-! utm_x_station(5) = 863545.3750000d0
-! utm_y_station(5) = 19669.6621094d0
-
-! utm_x_station(6) = 817099.3750000d0
-! utm_y_station(6) = -24430.2871094d0
-
-! print*,myrank,'après store the position of the six stations'
-! call flush(6)
-
-! print*, myrank,minval(nodes_coords_ext_mesh(1,:))
-! call flush(6)
-
-
-! print*, myrank,maxval(nodes_coords_ext_mesh(1,:))
-! call flush(6)
-
-
-! do ispec = 1, nspec
-
-! zmesh = zstore(2,2,2,ispec)
-
-! ! if(doubling_index == IFLAG_ONE_LAYER_TOPOGRAPHY) then
-! if(any(ibelm_top == ispec)) then
-! doubling_value_found_for_Piero = IFLAG_ONE_LAYER_TOPOGRAPHY
-
-! else if(zmesh < Z_23p4km) then
-! doubling_value_found_for_Piero = IFLAG_MANTLE_BELOW_23p4km
-
-! else if(zmesh < Z_14km) then
-! doubling_value_found_for_Piero = IFLAG_14km_to_23p4km
-
-! else
-! doubling_value_found_for_Piero = IFLAG_BEDROCK_down_to_14km
-! endif
-! idoubling(ispec) = doubling_value_found_for_Piero
-
-! do k = 1, NGLLZ
-! do j = 1, NGLLY
-! do i = 1, NGLLX
-
-
-! if(idoubling(ispec) == IFLAG_ONE_LAYER_TOPOGRAPHY .or. &
-! idoubling(ispec) == IFLAG_BEDROCK_down_to_14km) then
-
-! ! since we have suppressed UTM projection for Piero Basini, UTMx is the same as long
-! ! and UTMy is the same as lat
-! long = xstore(i,j,k,ispec)
-! lat = ystore(i,j,k,ispec)
-
-! ! get coordinate of corner in model
-! icornerlong = int((long - ORIG_LONG_TOPO) / DEGREES_PER_CELL_TOPO) + 1
-! icornerlat = int((lat - ORIG_LAT_TOPO) / DEGREES_PER_CELL_TOPO) + 1
-
-! ! avoid edge effects and extend with identical point if outside model
-! if(icornerlong < 1) icornerlong = 1
-! if(icornerlong > NX_TOPO-1) icornerlong = NX_TOPO-1
-! if(icornerlat < 1) icornerlat = 1
-! if(icornerlat > NY_TOPO-1) icornerlat = NY_TOPO-1
-
-! ! compute coordinates of corner
-! long_corner = ORIG_LONG_TOPO + (icornerlong-1)*DEGREES_PER_CELL_TOPO
-! lat_corner = ORIG_LAT_TOPO + (icornerlat-1)*DEGREES_PER_CELL_TOPO
-
-! ! compute ratio for interpolation
-! ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO
-! ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO
-
-! ! avoid edge effects
-! if(ratio_xi < 0.) ratio_xi = 0.
-! if(ratio_xi > 1.) ratio_xi = 1.
-! if(ratio_eta < 0.) ratio_eta = 0.
-! if(ratio_eta > 1.) ratio_eta = 1.
-
-! ! interpolate elevation at current point
-! elevation_bedrock = &
-! ibedrock(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
-! ibedrock(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
-! ibedrock(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
-! ibedrock(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
-
-! !! DK DK exclude circles around each station to make sure they are on the bedrock
-! !! DK DK and not in the ice
-! is_around_a_station = .false.
-! do istation = 1,NUMBER_OF_STATIONS
-! if(sqrt((long - utm_x_station(istation))**2 + (lat - utm_y_station(istation))**2) < RADIUS_TO_EXCLUDE) then
-! is_around_a_station = .true.
-! exit
-! endif
-! enddo
-
-! ! define elastic parameters in the model
-
-! ! we are above the bedrock interface i.e. in the ice, and not too close to a station
-! if(zmesh >= elevation_bedrock .and. .not. is_around_a_station) then
-! vp = 3800.d0
-! vs = 1900.d0
-! rho = 900.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_ICE
-
-! ! we are below the bedrock interface i.e. in the bedrock, or close to a station
-! else
-! vp = 5800.d0
-! vs = 3200.d0
-! rho = 2600.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
-! endif
-
-! else if(idoubling(ispec) == IFLAG_14km_to_23p4km) then
-! vp = 6800.d0
-! vs = 3900.d0
-! rho = 2900.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
-
-! else if(idoubling(ispec) == IFLAG_MANTLE_BELOW_23p4km) then
-! vp = 8100.d0
-! vs = 4480.d0
-! rho = 3380.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
-
-! endif
-
-! !pll 8/06
-! if(CUSTOM_REAL == SIZE_REAL) then
-! rhostore(i,j,k,ispec) = sngl(rho)
-! vpstore(i,j,k,ispec) = sngl(vp)
-! vsstore(i,j,k,ispec) = sngl(vs)
-! else
-! rhostore(i,j,k,ispec) = rho
-! vpstore(i,j,k,ispec) = vp
-! vsstore(i,j,k,ispec) = vs
-! end if
-
-! kappastore(i,j,k,ispec) = rhostore(i,j,k,ispec)*(vpstore(i,j,k,ispec)*vpstore(i,j,k,ispec) - &
-! 4.d0*vsstore(i,j,k,ispec)*vsstore(i,j,k,ispec)/3.d0)
-! mustore(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)*&
-! vsstore(i,j,k,ispec)
-
-! ! Stacey, a completer par la suite
-! rho_vp(i,j,k,ispec) = rhostore(i,j,k,ispec)*vpstore(i,j,k,ispec)
-! rho_vs(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)
-! !end pll
-
-! ! kappastore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
-! ! (materials_ext_mesh(2,mat_ext_mesh(ispec))*materials_ext_mesh(2,mat_ext_mesh(ispec)) - &
-! ! 4.d0*materials_ext_mesh(3,mat_ext_mesh(ispec))*materials_ext_mesh(3,mat_ext_mesh(ispec))/3.d0)
-! ! mustore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
-! materials_ext_mesh(3,mat_ext_mesh(ispec))*&
-! ! x materials_ext_mesh(3,mat_ext_mesh(ispec))
-! enddo
-! enddo
-! enddo
-! enddo
-!
-! end subroutine
-
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-
-!pll
-! subroutine interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
-
-! implicit none
-
-! include "constants.h"
-
-! integer :: iflag,flag_below,flag_above
-! integer :: ispec,nspec
-! integer :: i,j,k
-! double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
-! real(kind=CUSTOM_REAL), dimension(NX_TOPO_ANT,NY_TOPO_ANT) :: ibedrock
-! integer, parameter :: NUMBER_OF_STATIONS = 1
-! double precision, parameter :: RADIUS_TO_EXCLUDE = 250.d0
-! double precision, dimension(NUMBER_OF_STATIONS) :: utm_x_station,utm_y_station
-
-! !-------------------
-
-! !for Piero
-! logical :: is_around_a_station
-! integer :: istation
-
-! ! store bedrock values
-! integer :: icornerlat,icornerlong
-! double precision :: lat,long,elevation_bedrock
-! double precision :: lat_corner,long_corner,ratio_xi,ratio_eta
-
-
-! !! DK DK store the position of the six stations to be able to
-! !! DK DK exclude circles around each station to make sure they are on the bedrock
-! !! DK DK and not in the ice
-! utm_x_station(1) = 783500.6250000d0
-! utm_y_station(1) = -11828.7519531d0
-
-! utm_x_station(2) = 853644.5000000d0
-! utm_y_station(2) = -114.0138092d0
-
-! utm_x_station(3) = 863406.0000000d0
-! utm_y_station(3) = -53736.1640625d0
-
-! utm_x_station(4) = 823398.8125000d0
-! utm_y_station(4) = 29847.4511719d0
-
-! utm_x_station(5) = 863545.3750000d0
-! utm_y_station(5) = 19669.6621094d0
-
-! utm_x_station(6) = 817099.3750000d0
-! utm_y_station(6) = -24430.2871094d0
-
-! ! since we have suppressed UTM projection for Piero Basini, UTMx is the same as long
-! ! and UTMy is the same as lat
-! long = xstore(i,j,k,ispec)
-! lat = ystore(i,j,k,ispec)
-
-! ! get coordinate of corner in model
-! icornerlong = int((long - ORIG_LONG_TOPO_ANT) / DEGREES_PER_CELL_TOPO_ANT) + 1
-! icornerlat = int((lat - ORIG_LAT_TOPO_ANT) / DEGREES_PER_CELL_TOPO_ANT) + 1
-
-! ! avoid edge effects and extend with identical point if outside model
-! if(icornerlong < 1) icornerlong = 1
-! if(icornerlong > NX_TOPO_ANT-1) icornerlong = NX_TOPO_ANT-1
-! if(icornerlat < 1) icornerlat = 1
-! if(icornerlat > NY_TOPO_ANT-1) icornerlat = NY_TOPO_ANT-1
-
-! ! compute coordinates of corner
-! long_corner = ORIG_LONG_TOPO_ANT + (icornerlong-1)*DEGREES_PER_CELL_TOPO_ANT
-! lat_corner = ORIG_LAT_TOPO_ANT + (icornerlat-1)*DEGREES_PER_CELL_TOPO_ANT
-
-! ! compute ratio for interpolation
-! ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO_ANT
-! ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO_ANT
-
-! ! avoid edge effects
-! if(ratio_xi < 0.) ratio_xi = 0.
-! if(ratio_xi > 1.) ratio_xi = 1.
-! if(ratio_eta < 0.) ratio_eta = 0.
-! if(ratio_eta > 1.) ratio_eta = 1.
-
-! ! interpolate elevation at current point
-! elevation_bedrock = &
-! ibedrock(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
-! ibedrock(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
-! ibedrock(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
-! ibedrock(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
-
-! !! DK DK exclude circles around each station to make sure they are on the bedrock
-! !! DK DK and not in the ice
-! is_around_a_station = .false.
-! do istation = 1,NUMBER_OF_STATIONS
-! if(sqrt((xstore(i,j,k,ispec) - utm_x_station(istation))**2 + (ystore(i,j,k,ispec) - &
-! utm_y_station(istation))**2) < RADIUS_TO_EXCLUDE) then
-! is_around_a_station = .true.
-! exit
-! endif
-! enddo
-
-! ! we are above the bedrock interface i.e. in the ice, and not too close to a station
-! if(zstore(i,j,k,ispec) >= elevation_bedrock .and. .not. is_around_a_station) then
-! iflag = flag_above
-! !iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_ICE
-! ! we are below the bedrock interface i.e. in the bedrock, or close to a station
-! else
-! iflag = flag_below
-! !iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
-! endif
-
-
-! end subroutine interface
More information about the CIG-COMMITS
mailing list