[cig-commits] r20842 - in seismo/3D/SPECFEM3D/trunk: doc/USER_MANUAL examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files examples/Mount_StHelens/in_data_files examples/homogeneous_halfspace/in_data_files examples/homogeneous_poroelastic/in_data_files examples/layered_halfspace/in_data_files examples/meshfem3D_examples/many_interfaces examples/meshfem3D_examples/simple_model/in_data_files examples/meshfem3D_examples/socal1D/in_data_files examples/tomographic_model/in_data_files examples/waterlayered_halfspace/in_data_files in_data_files src/generate_databases src/shared src/specfem3D
joseph.charles at geodynamics.org
joseph.charles at geodynamics.org
Wed Oct 17 09:52:57 PDT 2012
Author: joseph.charles
Date: 2012-10-17 09:52:57 -0700 (Wed, 17 Oct 2012)
New Revision: 20842
Modified:
seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.pdf
seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex
seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/Par_file
seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/in_data_files/FORCESOLUTION
seismo/3D/SPECFEM3D/trunk/in_data_files/Par_file
seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
seismo/3D/SPECFEM3D/trunk/src/shared/get_force.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90
Log:
adds the hdur parameter in FORCESOLUTION files and the use of a quasi-Heaviside source time function for runs involving a FORCESOLUTION file
Modified: seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex
===================================================================
--- seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/doc/USER_MANUAL/manual_SPECFEM3D.tex 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1059,10 +1059,12 @@
\item [{\texttt{USE\_FORCE\_POINT\_SOURCE}}] Turn this flag on to use a \texttt{FORCESOLUTION} force source
instead of a \texttt{CMTSOLUTION} moment-tensor source. This can be useful e.g. for oil industry
foothills simulations in which the source is a vertical force, normal force,
-tilted force, or an impact etc. You will need to give the East, North and vertical components of an
+tilted force, or an impact etc. By default a quasi-Heaviside time function is used to represent the force point source but a Ricker
+source time function can be computed instead by changing the value of \texttt{USE\_RICKER\_IPATI} in the main constants file
+\texttt{constants.h} located in the \texttt{src/shared/} subdirectory.
+Note that in the \texttt{FORCESOLUTION} file, you will need to edit the East, North and vertical components of an
arbitrary (non-unitary) direction vector of the force vector; thus refer to Appendix~\ref{cha:Coordinates}
-for the orientation of the reference frame.
-This vector is made unitary internally in the code and thus only its direction matters here;
+for the orientation of the reference frame. This vector is made unitary internally in the code and thus only its direction matters here;
its norm is ignored and the norm of the force used is the factor force source times the source time function.
\item [{\texttt{PRINT\_SOURCE\_TIME\_FUNCTION}}] Turn this flag on to print
information about the source time function in the file \texttt{in\_out\_files/OUTPUT\_FILES/plot\_source\_time\_function.txt}.
@@ -1105,7 +1107,7 @@
\begin{description}
\item [{\texttt{Par\_file}}] the main parameter file which was discussed in detail in the previous
Chapter~\ref{cha:Creating-Distributed-Databases},
-\item [{\texttt{CMTSOLUTION}} or {\texttt{FORCESOLUTION}}] the earthquake source parameter file or the force source parameter file or, and
+\item [{\texttt{CMTSOLUTION}} or {\texttt{FORCESOLUTION}}] the earthquake source parameter file or the force source parameter file, and
\item [{\texttt{STATIONS}}] the stations file.
\end{description}
@@ -1149,7 +1151,7 @@
\end{figure}
{\small \par}
\end{lyxcode}
-The \texttt{CMTSOLUTION} should be edited in the following way:
+The \texttt{CMTSOLUTION} file should be edited in the following way:
\begin{itemize}
\item Set the \texttt{time shift} parameter equal to $0.0$ (the solver
@@ -1242,13 +1244,14 @@
\vspace{1cm}
-\noindent The \texttt{FORCESOLUTION} should be edited in the following way:
+\noindent The \texttt{FORCESOLUTION} file should be edited in the following way:
\begin{itemize}
\item Set the \texttt{time shift} parameter equal to $0.0$ (the solver
will not run otherwise.) The time shift parameter would simply apply
an overall time shift to the synthetics, something that can be done
in the post-processing (see Section \ref{sec:Process-data-and-syn}).
+\item Set the \texttt{half duration} parameter of the source time function.
\item Set the latitude, longitude, depth of the source in geographical
coordinates.
\item Set the magnitude of the force source.
@@ -1262,9 +1265,9 @@
$N_{\mathrm{sources}}$ entries, one for each subevent (i.e.,
concatenate $N_{\mathrm{sources}}$ \texttt{FORCESOLUTION} files to a
single \texttt{FORCESOLUTION} file). At least one entry (not
-necessarily the first) must have a zero time shift, and all the other
-entries must have non-negative time shift. Each subevent can have its
-own half latitude, longitude, depth, and force parameters.
+necessarily the first) must have a zero \texttt{time shift}, and all the other
+entries must have non-negative \texttt{time shift}. Each subevent can have its
+own half latitude, longitude, depth, \texttt{half duration} and force parameters.
{\small }%
\begin{figure}[H]
Modified: seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
time shift: 0.0000
+hdur: 0.5
latitude: 1320.0
longitude: 1320.0
depth: 1.0
Modified: seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ACOUSTIC/in_data_files/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -69,7 +69,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
time shift: 0.0000
+hdur: 0.5
latitude: 1320.0
longitude: 1320.0
depth: 1.0
Modified: seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/BENCHMARK_CLAERBOUT_ADJOINT/ELASTIC/in_data_files/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -69,7 +69,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
time shift: 0.0000
+hdur: 0.5
latitude: 46.197
longitude: -122.186
depth: 5.0
Modified: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/in_data_files/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -69,7 +69,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
time shift: 0.0000
+hdur: 0.5
latitude: 67000.0
longitude: 67000.0
depth: 25.0
Modified: seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/in_data_files/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -71,7 +71,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
time shift: 0.0000
+hdur: 0.5
latitude: 100.0
longitude: 150.0
depth: 0.150
Modified: seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/homogeneous_poroelastic/in_data_files/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -71,7 +71,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
time shift: 0.0000
+hdur: 0.5
latitude: 67000.0
longitude: 67000.0
depth: 25.05
Modified: seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/layered_halfspace/in_data_files/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -70,7 +70,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
time shift: 0.0000
+hdur: 0.5
latitude: 15150.0
longitude: 14925.0
depth: 10.0
Modified: seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/many_interfaces/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -70,7 +70,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
-time shift: 0.0000
+time shift: 0.0000
+hdur: 0.5
latitude: 33.6
longitude: -118.4
depth: 10.0
Modified: seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/in_data_files/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -70,7 +70,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
time shift: 0.0000
+hdur: 0.5
latitude: 34.0745
longitude: -118.3792
depth: 5.4000
Modified: seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/socal1D/in_data_files/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -70,7 +70,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
time shift: 0.0000
+hdur: 0.5
latitude: 67000.0
longitude: 67000.0
depth: 25.0
Modified: seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/tomographic_model/in_data_files/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -70,7 +70,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
time shift: 0.0000
+hdur: 0.5
latitude: 67000.0
longitude: 66999.9
depth: 1.0
Modified: seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/examples/waterlayered_halfspace/in_data_files/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -70,7 +70,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/in_data_files/FORCESOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/in_data_files/FORCESOLUTION 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/in_data_files/FORCESOLUTION 2012-10-17 16:52:57 UTC (rev 20842)
@@ -1,5 +1,6 @@
FORCE 001
time shift: 0.0000
+hdur: 0.5
latitude: 67000.0
longitude: 67000.0
depth: 25.0
Modified: seismo/3D/SPECFEM3D/trunk/in_data_files/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/in_data_files/Par_file 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/in_data_files/Par_file 2012-10-17 16:52:57 UTC (rev 20842)
@@ -71,7 +71,8 @@
# in which the source is a vertical force, normal force, inclined force, impact etc.
USE_FORCE_POINT_SOURCE = .false.
# If this flag is turned on, the FORCESOLUTION file must be edited by precising:
-# - the related time-shift parameter,
+# - the corresponding time-shift parameter,
+# - the half duration parameter of the source,
# - the coordinates of the source,
# - the magnitude of the force source,
# - the components of a (non-unitary) direction vector for the force source in the E/N/Z_UP basis.
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -407,7 +407,7 @@
write(IMAIN,*)
if(USE_FORCE_POINT_SOURCE) then
- write(IMAIN,*) 'using a force point source instead of a CMTSOLUTION source'
+ write(IMAIN,*) 'using a FORCESOLUTION source instead of a CMTSOLUTION source'
else
write(IMAIN,*) 'using a CMTSOLUTION source'
write(IMAIN,*)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2012-10-17 16:52:57 UTC (rev 20842)
@@ -194,7 +194,8 @@
! no lagrange interpolation on seismograms (we take the value on one NGLL point)
logical, parameter :: FASTER_RECEIVERS_POINTS_ONLY = .false.
-! set to use a Ricker source time function instead of a gaussian
+! set to use a Ricker source time function instead of a gaussian (when using CMTSOLUTION files)
+! or a quasi-Heaviside (when using FORCESOLUTION files)
logical, parameter :: USE_RICKER_IPATI = .false.
! use this t0 as earliest starting time rather than the automatically calculated one
@@ -341,7 +342,7 @@
! number of lines per source in CMTSOLUTION file
integer, parameter :: NLINES_PER_CMTSOLUTION_SOURCE = 13
- integer, parameter :: NLINES_PER_FORCESOLUTION_SOURCE = 9
+ integer, parameter :: NLINES_PER_FORCESOLUTION_SOURCE = 10
! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 4
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_force.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_force.f90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_force.f90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -83,6 +83,10 @@
read(1,"(a)") string
read(string(12:len_trim(string)),*) t_shift(isource)
+ ! read half duration
+ read(1,"(a)") string
+ read(string(6:len_trim(string)),*) hdur(isource)
+
! read latitude
read(1,"(a)") string
read(string(10:len_trim(string)),*) lat(isource)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -30,7 +30,7 @@
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
xi_source,eta_source,gamma_source, &
- hdur,hdur_gaussian,tshift_cmt,dt,t0, &
+ hdur,hdur_gaussian,tshift_src,dt,t0, &
sourcearrays,kappastore,ispec_is_acoustic,&
SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
@@ -65,7 +65,7 @@
integer :: NSOURCES,myrank,it
integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
- double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,tshift_cmt
+ double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,tshift_src
double precision :: dt,t0
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
@@ -118,127 +118,136 @@
if (SIMULATION_TYPE == 1 .and. nsources_local > 0) then
!way 2
- if(GPU_MODE) then
- if( NSOURCES > 0 ) then
- do isource = 1,NSOURCES
- if(USE_FORCE_POINT_SOURCE) then
- ! precomputes source time function factor
- stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function_rickr( &
- dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
- else
- if( USE_RICKER_IPATI ) then
- stf_pre_compute(isource) = comp_source_time_function_rickr( &
- dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
- else
- stf_pre_compute(isource) = comp_source_time_function_gauss( &
- dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
- endif
- endif
- enddo
- stf_used_total = stf_used_total + sum(stf_pre_compute(:))
- ! only implements SIMTYPE=1 and NOISE_TOM=0
- ! write(*,*) "fortran dt = ", dt
- ! change dt -> DT
- call compute_add_sources_ac_cuda(Mesh_pointer, phase_is_inner, &
- NSOURCES, stf_pre_compute, myrank)
- endif
+ if(GPU_MODE) then
+ if( NSOURCES > 0 ) then
+ do isource = 1,NSOURCES
+ ! precomputes source time function factor
+ if(USE_FORCE_POINT_SOURCE) then
+ if( USE_RICKER_IPATI ) then
+ stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function_rickr( &
+ dble(it-1)*DT-t0-tshift_src(isource),hdur(isource))
+ else
+ stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function( &
+ dble(it-1)*DT-t0-tshift_src(isource),hdur(isource))
+ endif
+ else
+ if( USE_RICKER_IPATI ) then
+ stf_pre_compute(isource) = comp_source_time_function_rickr( &
+ dble(it-1)*DT-t0-tshift_src(isource),hdur(isource))
+ else
+ stf_pre_compute(isource) = comp_source_time_function_gauss( &
+ dble(it-1)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
+ endif
+ endif
+ enddo
+ stf_used_total = stf_used_total + sum(stf_pre_compute(:))
+ ! only implements SIMTYPE=1 and NOISE_TOM=0
+ ! write(*,*) "fortran dt = ", dt
+ ! change dt -> DT
+ call compute_add_sources_ac_cuda(Mesh_pointer, phase_is_inner, &
+ NSOURCES, stf_pre_compute, myrank)
+ endif
- else ! .NOT. GPU_MODE
+ else ! .NOT. GPU_MODE
- ! adds acoustic sources
- do isource = 1,NSOURCES
+ ! adds acoustic sources
+ do isource = 1,NSOURCES
- ! add the source (only if this proc carries the source)
- if(myrank == islice_selected_source(isource)) then
+ ! add the source (only if this proc carries the source)
+ if(myrank == islice_selected_source(isource)) then
- ispec = ispec_selected_source(isource)
+ ispec = ispec_selected_source(isource)
- if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
- if( ispec_is_acoustic(ispec) ) then
+ if( ispec_is_acoustic(ispec) ) then
- if(USE_FORCE_POINT_SOURCE) then
+ if(USE_FORCE_POINT_SOURCE) then
- ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
- iglob = ibool(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec)
+ ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
+ iglob = ibool(nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)), &
+ ispec)
- f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing FORCESOLUTION file format
+ f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing FORCESOLUTION file format
- !if (it == 1 .and. myrank == 0) then
- ! write(IMAIN,*) 'using a source of dominant frequency ',f0
- ! write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
- ! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
- !endif
+ !if (it == 1 .and. myrank == 0) then
+ ! write(IMAIN,*) 'using a source of dominant frequency ',f0
+ ! write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+ ! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ !endif
- ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = factor_force_source(isource) * &
- comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+ if( USE_RICKER_IPATI ) then
+ stf_used = factor_force_source(isource) * &
+ comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_src(isource),f0)
+ else
+ stf_used = factor_force_source(isource) * &
+ comp_source_time_function(dble(it-1)*DT-t0-tshift_src(isource),f0)
+ endif
- ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
- ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
- ! to add minus the source to Chi_dot_dot to get plus the source in pressure:
+ ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
+ ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
+ ! to add minus the source to Chi_dot_dot to get plus the source in pressure:
- ! acoustic source for pressure gets divided by kappa
- ! source contribution
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - stf_used / kappastore(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)),ispec)
+ ! acoustic source for pressure gets divided by kappa
+ ! source contribution
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - stf_used / kappastore(nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)),ispec)
- else
+ else
- if( USE_RICKER_IPATI ) then
- stf = comp_source_time_function_rickr( &
- dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
- else
- ! gaussian source time
- stf = comp_source_time_function_gauss( &
- dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
- endif
+ if( USE_RICKER_IPATI ) then
+ stf = comp_source_time_function_rickr( &
+ dble(it-1)*DT-t0-tshift_src(isource),hdur(isource))
+ else
+ ! gaussian source time
+ stf = comp_source_time_function_gauss( &
+ dble(it-1)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
+ endif
- ! quasi-heaviside
- !stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ ! quasi-heaviside
+ !stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
- ! source encoding
- stf = stf * pm1_source_encoding(isource)
+ ! source encoding
+ stf = stf * pm1_source_encoding(isource)
- ! distinguishes between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- stf_used = sngl(stf)
- else
- stf_used = stf
- endif
+ ! distinguishes between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
- ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
- ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
- ! to add minus the source to Chi_dot_dot to get plus the source in pressure
+ ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
+ ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
+ ! to add minus the source to Chi_dot_dot to get plus the source in pressure
- ! add source array
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- ! adds source contribution
- ! note: acoustic source for pressure gets divided by kappa
- iglob = ibool(i,j,k,ispec)
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - sourcearrays(isource,1,i,j,k) * stf_used / kappastore(i,j,k,ispec)
- enddo
- enddo
- enddo
+ ! add source array
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! adds source contribution
+ ! note: acoustic source for pressure gets divided by kappa
+ iglob = ibool(i,j,k,ispec)
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ - sourcearrays(isource,1,i,j,k) * stf_used / kappastore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
- endif ! USE_FORCE_POINT_SOURCE
+ endif ! USE_FORCE_POINT_SOURCE
- stf_used_total = stf_used_total + stf_used
+ stf_used_total = stf_used_total + stf_used
- endif ! ispec_is_acoustic
- endif ! ispec_is_inner
- endif ! myrank
+ endif ! ispec_is_acoustic
+ endif ! ispec_is_inner
+ endif ! myrank
- enddo ! NSOURCES
- endif ! GPU_MODE
+ enddo ! NSOURCES
+ endif ! GPU_MODE
endif
! NOTE: adjoint sources and backward wavefield timing:
@@ -408,130 +417,138 @@
! adjoint simulations
if (SIMULATION_TYPE == 3 .and. nsources_local > 0) then
- ! on GPU
- if(GPU_MODE) then
- if( NSOURCES > 0 ) then
- do isource = 1,NSOURCES
- if(USE_FORCE_POINT_SOURCE) then
- ! precomputes source time function factors
- stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function_rickr( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
- else
- if( USE_RICKER_IPATI ) then
- stf_pre_compute(isource) = comp_source_time_function_rickr( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
- else
- stf_pre_compute(isource) = comp_source_time_function_gauss( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
- endif
- endif
- enddo
- stf_used_total = stf_used_total + sum(stf_pre_compute(:))
- ! only implements SIMTYPE=3
- call compute_add_sources_ac_s3_cuda(Mesh_pointer, phase_is_inner, &
- NSOURCES,stf_pre_compute, myrank)
- endif
+ ! on GPU
+ if(GPU_MODE) then
+ if( NSOURCES > 0 ) then
+ do isource = 1,NSOURCES
+ ! precomputes source time function factors
+ if(USE_FORCE_POINT_SOURCE) then
+ if( USE_RICKER_IPATI ) then
+ stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function_rickr( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur(isource))
+ else
+ stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur(isource))
+ endif
+ else
+ if( USE_RICKER_IPATI ) then
+ stf_pre_compute(isource) = comp_source_time_function_rickr( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur(isource))
+ else
+ stf_pre_compute(isource) = comp_source_time_function_gauss( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
+ endif
+ endif
+ enddo
+ stf_used_total = stf_used_total + sum(stf_pre_compute(:))
+ ! only implements SIMTYPE=3
+ call compute_add_sources_ac_s3_cuda(Mesh_pointer, phase_is_inner, &
+ NSOURCES,stf_pre_compute, myrank)
+ endif
- else ! .NOT. GPU_MODE
+ else ! .NOT. GPU_MODE
- ! adds acoustic sources
- do isource = 1,NSOURCES
+ ! adds acoustic sources
+ do isource = 1,NSOURCES
- ! add the source (only if this proc carries the source)
- if(myrank == islice_selected_source(isource)) then
+ ! add the source (only if this proc carries the source)
+ if(myrank == islice_selected_source(isource)) then
- ispec = ispec_selected_source(isource)
+ ispec = ispec_selected_source(isource)
- if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
- if( ispec_is_acoustic(ispec) ) then
+ if( ispec_is_acoustic(ispec) ) then
- if(USE_FORCE_POINT_SOURCE) then
+ if(USE_FORCE_POINT_SOURCE) then
- ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
- iglob = ibool(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec)
+ ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
+ iglob = ibool(nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)), &
+ ispec)
- f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
+ f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
- !if (it == 1 .and. myrank == 0) then
- ! write(IMAIN,*) 'using a source of dominant frequency ',f0
- ! write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
- ! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
- !endif
+ !if (it == 1 .and. myrank == 0) then
+ ! write(IMAIN,*) 'using a source of dominant frequency ',f0
+ ! write(IMAIN,*) 'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+ ! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ !endif
- ! gaussian source time function
- !stf_used = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ ! gaussian source time function
+ !stf_used = comp_source_time_function(dble(it-1)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
- ! we use nu_source(:,3) here because we want a source normal to the surface.
- ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = factor_force_source(isource) * comp_source_time_function_rickr( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
+ if( USE_RICKER_IPATI ) then
+ stf_used = factor_force_source(isource) * comp_source_time_function_rickr( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),f0)
+ else
+ stf_used = factor_force_source(isource) * comp_source_time_function( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),f0)
+ endif
- ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
- ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
- ! to add minus the source to Chi_dot_dot to get plus the source in pressure:
+ ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
+ ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
+ ! to add minus the source to Chi_dot_dot to get plus the source in pressure:
- ! acoustic source for pressure gets divided by kappa
- ! source contribution
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
+ ! acoustic source for pressure gets divided by kappa
+ ! source contribution
+ b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
- stf_used / kappastore(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)),ispec)
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)),ispec)
- else
+ else
- if( USE_RICKER_IPATI ) then
- stf = comp_source_time_function_rickr( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
- else
- ! gaussian source time
- stf = comp_source_time_function_gauss( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
- endif
+ if( USE_RICKER_IPATI ) then
+ stf = comp_source_time_function_rickr( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur(isource))
+ else
+ ! gaussian source time
+ stf = comp_source_time_function_gauss( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
+ endif
- ! quasi-heaviside
- !stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ ! quasi-heaviside
+ !stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
- ! source encoding
- stf = stf * pm1_source_encoding(isource)
+ ! source encoding
+ stf = stf * pm1_source_encoding(isource)
- ! distinguishes between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- stf_used = sngl(stf)
- else
- stf_used = stf
- endif
+ ! distinguishes between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
- ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
- ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
- ! to add minus the source to Chi_dot_dot to get plus the source in pressure
+ ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
+ ! the sign is negative because pressure p = - Chi_dot_dot therefore we need
+ ! to add minus the source to Chi_dot_dot to get plus the source in pressure
- ! add source array
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- ! adds source contribution
- ! note: acoustic source for pressure gets divided by kappa
- iglob = ibool(i,j,k,ispec)
- b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
- - sourcearrays(isource,1,i,j,k) * stf_used / kappastore(i,j,k,ispec)
- enddo
- enddo
- enddo
+ ! add source array
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! adds source contribution
+ ! note: acoustic source for pressure gets divided by kappa
+ iglob = ibool(i,j,k,ispec)
+ b_potential_dot_dot_acoustic(iglob) = b_potential_dot_dot_acoustic(iglob) &
+ - sourcearrays(isource,1,i,j,k) * stf_used / kappastore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
- endif ! USE_FORCE_POINT_SOURCE
+ endif ! USE_FORCE_POINT_SOURCE
- stf_used_total = stf_used_total + stf_used
+ stf_used_total = stf_used_total + stf_used
- endif ! ispec_is_elastic
- endif ! ispec_is_inner
- endif ! myrank
+ endif ! ispec_is_elastic
+ endif ! ispec_is_inner
+ endif ! myrank
- enddo ! NSOURCES
- endif ! GPU_MODE
+ enddo ! NSOURCES
+ endif ! GPU_MODE
endif
! master prints out source time function to file
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -30,7 +30,7 @@
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
xi_source,eta_source,gamma_source, &
- hdur,hdur_gaussian,tshift_cmt,dt,t0,sourcearrays, &
+ hdur,hdur_gaussian,tshift_src,dt,t0,sourcearrays, &
ispec_is_elastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
nadj_rec_local,adj_sourcearrays,b_accel, &
@@ -68,7 +68,7 @@
integer :: NSOURCES,myrank,it
integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
- double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,tshift_cmt
+ double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,tshift_src
double precision :: dt,t0
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
@@ -123,17 +123,22 @@
if(GPU_MODE) then
if( NSOURCES > 0 ) then
do isource = 1,NSOURCES
+ ! precomputes source time function factor
if(USE_FORCE_POINT_SOURCE) then
- ! precomputes source time function factor
- stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function_rickr( &
- dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
+ if( USE_RICKER_IPATI ) then
+ stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function_rickr( &
+ dble(it-1)*DT-t0-tshift_src(isource),hdur(isource))
+ else
+ stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function( &
+ dble(it-1)*DT-t0-tshift_src(isource),hdur(isource))
+ endif
else
if( USE_RICKER_IPATI ) then
stf_pre_compute(isource) = comp_source_time_function_rickr( &
- dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
+ dble(it-1)*DT-t0-tshift_src(isource),hdur(isource))
else
stf_pre_compute(isource) = comp_source_time_function( &
- dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ dble(it-1)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
endif
endif
enddo
@@ -174,8 +179,11 @@
! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
!endif
- ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+ if( USE_RICKER_IPATI) then
+ stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_src(isource),f0)
+ else
+ stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_src(isource),f0)
+ endif
! add the inclined force source array
! distinguish between single and double precision for reals
@@ -197,9 +205,9 @@
else
if( USE_RICKER_IPATI) then
- stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
+ stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_src(isource),hdur(isource))
else
- stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
endif
! distinguish between single and double precision for reals
@@ -410,17 +418,22 @@
if(GPU_MODE) then
if( NSOURCES > 0 ) then
do isource = 1,NSOURCES
+ ! precomputes source time function factors
if(USE_FORCE_POINT_SOURCE) then
- ! precomputes source time function factors
- stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function_rickr( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
+ if( USE_RICKER_IPATI ) then
+ stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function_rickr( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur(isource))
+ else
+ stf_pre_compute(isource) = factor_force_source(isource) * comp_source_time_function( &
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur(isource))
+ endif
else
if( USE_RICKER_IPATI ) then
stf_pre_compute(isource) = comp_source_time_function_rickr( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur(isource))
else
stf_pre_compute(isource) = comp_source_time_function( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
endif
endif
enddo
@@ -460,8 +473,11 @@
! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
!endif
- ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf = comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
+ if( USE_RICKER_IPATI ) then
+ stf = comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_src(isource),f0)
+ else
+ stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_src(isource),f0)
+ endif
! add the inclined force source array
! distinguish between single and double precision for reals
@@ -486,10 +502,10 @@
! (also compare to it-1 for forward simulation)
if( USE_RICKER_IPATI ) then
stf = comp_source_time_function_rickr( &
- dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
+ dble(it-1)*DT-t0-tshift_src(isource),hdur(isource))
else
stf = comp_source_time_function( &
- dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
endif
! distinguish between single and double precision for reals
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -33,7 +33,7 @@
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
xi_source,eta_source,gamma_source, &
- hdur,hdur_gaussian,tshift_cmt,dt,t0,sourcearrays, &
+ hdur,hdur_gaussian,tshift_src,dt,t0,sourcearrays, &
ispec_is_poroelastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
nadj_rec_local,adj_sourcearrays,b_accels,b_accelw, &
@@ -67,7 +67,7 @@
integer :: NSOURCES,myrank,it
integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
- double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,tshift_cmt
+ double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,tshift_src
double precision :: dt,t0
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
@@ -153,10 +153,15 @@
! print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
endif
- ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
- !stf_used = comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+ ! gaussian source time function
+ !stf_used = comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_src(isource),f0)
+ if( USE_RICKER_IPATI ) then
+ stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_src(isource),f0)
+ else
+ stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_src(isource),f0)
+ endif
+
! add the inclined force source array
! the source is applied to both solid and fluid phase: bulk source.
@@ -183,10 +188,10 @@
else
- !stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ !stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
!t0 = 1.2d0/hdur(isource)
- !stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
- stf = comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ !stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_src(isource),hdur(isource))
+ stf = comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
@@ -413,9 +418,11 @@
! write(IMAIN,*) 'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
!endif
- ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- ! should be the same than for the forward simulation (check above)
- stf = comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),f0)
+ if( USE_RICKER_IPATI ) then
+ stf = comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_src(isource),f0)
+ else
+ stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_src(isource),f0)
+ endif
! add the inclined force source array
! the source is applied to both solid and fluid phase: bulk source
@@ -446,10 +453,10 @@
! see note above: time step corresponds now to NSTEP-it
! (also compare to it-1 for forward simulation)
- !stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ !stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
!t0 = 1.2d0/hdur(isource)
- !stf = comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur(isource))
- stf = comp_source_time_function_gauss(dble(NSTEP-it)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ !stf = comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur(isource))
+ stf = comp_source_time_function_gauss(dble(NSTEP-it)*DT-t0-tshift_src(isource),hdur_gaussian(isource))
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic.f90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -282,7 +282,7 @@
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
xi_source,eta_source,gamma_source, &
- hdur,hdur_gaussian,tshift_cmt,dt,t0, &
+ hdur,hdur_gaussian,tshift_src,dt,t0, &
sourcearrays,kappastore,ispec_is_acoustic,&
SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_elastic.F90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -236,7 +236,7 @@
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
xi_source,eta_source,gamma_source, &
- hdur,hdur_gaussian,tshift_cmt,dt,t0,sourcearrays, &
+ hdur,hdur_gaussian,tshift_src,dt,t0,sourcearrays, &
ispec_is_elastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
nadj_rec_local,adj_sourcearrays,b_accel, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -206,7 +206,7 @@
ibool,ispec_is_inner,phase_is_inner, &
NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
xi_source,eta_source,gamma_source, &
- hdur,hdur_gaussian,tshift_cmt,dt,t0,sourcearrays, &
+ hdur,hdur_gaussian,tshift_src,dt,t0,sourcearrays, &
ispec_is_poroelastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
nadj_rec_local,adj_sourcearrays,b_accels_poroelastic,b_accelw_poroelastic, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -30,7 +30,7 @@
subroutine locate_source(ibool,NSOURCES,myrank,NSPEC_AB,NGLOB_AB,xstore,ystore,zstore, &
xigll,yigll,zigll,NPROC, &
- tshift_cmt,min_tshift_cmt_original,yr,jda,ho,mi,utm_x_source,utm_y_source, &
+ tshift_src,min_tshift_src_original,yr,jda,ho,mi,utm_x_source,utm_y_source, &
DT,hdur,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
islice_selected_source,ispec_selected_source, &
xi_source,eta_source,gamma_source, &
@@ -64,8 +64,8 @@
integer yr,jda,ho,mi
- double precision tshift_cmt(NSOURCES)
- double precision sec,min_tshift_cmt_original
+ double precision tshift_src(NSOURCES)
+ double precision sec,min_tshift_src_original
integer iprocloop
@@ -169,11 +169,11 @@
! read all the sources (note: each process reads the source file)
if (USE_FORCE_POINT_SOURCE) then
- call get_force(tshift_cmt,hdur,lat,long,depth,NSOURCES,min_tshift_cmt_original,factor_force_source, &
+ call get_force(tshift_src,hdur,lat,long,depth,NSOURCES,min_tshift_src_original,factor_force_source, &
comp_dir_vect_source_E,comp_dir_vect_source_N,comp_dir_vect_source_Z_UP)
else
- call get_cmt(yr,jda,ho,mi,sec,tshift_cmt,hdur,lat,long,depth,moment_tensor, &
- DT,NSOURCES,min_tshift_cmt_original)
+ call get_cmt(yr,jda,ho,mi,sec,tshift_src,hdur,lat,long,depth,moment_tensor, &
+ DT,NSOURCES,min_tshift_src_original)
endif
! define topology of the control element
@@ -752,11 +752,14 @@
! prints frequency content for point forces
f0 = hdur(isource)
- t0_ricker = 1.2d0/f0
write(IMAIN,*) ' using a source of dominant frequency ',f0
write(IMAIN,*) ' lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
write(IMAIN,*) ' lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
- write(IMAIN,*) ' t0_ricker = ',t0_ricker,'tshift_cmt = ',tshift_cmt(isource)
+ if( USE_RICKER_IPATI ) then
+ t0_ricker = 1.2d0/f0
+ write(IMAIN,*) ' t0_ricker = ',t0_ricker
+ endif
+ write(IMAIN,*) ' time shift = ',tshift_src(isource)
write(IMAIN,*)
write(IMAIN,*) ' half duration -> frequency: ',hdur(isource),' seconds**(-1)'
else
@@ -772,7 +775,7 @@
endif
write(IMAIN,*) ' half duration: ',hdur(isource),' seconds'
endif
- write(IMAIN,*) ' time shift: ',tshift_cmt(isource),' seconds'
+ write(IMAIN,*) ' time shift: ',tshift_src(isource),' seconds'
write(IMAIN,*)
write(IMAIN,*) 'original (requested) position of the source:'
write(IMAIN,*)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -74,7 +74,7 @@
use specfem_par_movie
implicit none
- double precision :: t0_acoustic,min_tshift_cmt_original
+ double precision :: t0_acoustic,min_tshift_src_original
integer :: yr,jda,ho,mi
integer :: isource,ispec,ier
@@ -90,7 +90,7 @@
xi_source(NSOURCES), &
eta_source(NSOURCES), &
gamma_source(NSOURCES), &
- tshift_cmt(NSOURCES), &
+ tshift_src(NSOURCES), &
hdur(NSOURCES), &
hdur_gaussian(NSOURCES), &
utm_x_source(NSOURCES), &
@@ -104,6 +104,12 @@
comp_dir_vect_source_N(NSOURCES), &
comp_dir_vect_source_Z_UP(NSOURCES),stat=ier)
if( ier /= 0 ) stop 'error allocating arrays for force point sources'
+ else
+ allocate(factor_force_source(1), &
+ comp_dir_vect_source_E(1), &
+ comp_dir_vect_source_N(1), &
+ comp_dir_vect_source_Z_UP(1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy FORCESOLUTION arrays'
endif
! for source encoding (acoustic sources so far only)
@@ -116,7 +122,7 @@
! xi_source, eta_source & gamma_source
call locate_source(ibool,NSOURCES,myrank,NSPEC_AB,NGLOB_AB, &
xstore,ystore,zstore,xigll,yigll,zigll,NPROC, &
- tshift_cmt,min_tshift_cmt_original,yr,jda,ho,mi,utm_x_source,utm_y_source, &
+ tshift_src,min_tshift_src_original,yr,jda,ho,mi,utm_x_source,utm_y_source, &
DT,hdur,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
islice_selected_source,ispec_selected_source, &
xi_source,eta_source,gamma_source, &
@@ -128,7 +134,7 @@
USE_FORCE_POINT_SOURCE,factor_force_source,comp_dir_vect_source_E, &
comp_dir_vect_source_N,comp_dir_vect_source_Z_UP)
- if(abs(minval(tshift_cmt)) > TINYVAL) call exit_MPI(myrank,'one tshift_cmt must be zero, others must be positive')
+ if(abs(minval(tshift_src)) > TINYVAL) call exit_MPI(myrank,'one tshift_src must be zero, others must be positive')
! filter source time function by Gaussian with hdur = HDUR_MOVIE when outputing movies or shakemaps
if (MOVIE_SURFACE .or. MOVIE_VOLUME .or. CREATE_SHAKEMAP) then
@@ -146,7 +152,7 @@
! define t0 as the earliest start time
! note: an earlier start time also reduces numerical noise due to a
! non-zero offset at the beginning of the source time function
- t0 = - 2.0d0 * minval(tshift_cmt(:) - hdur(:)) ! - 1.5d0 * minval(tshift_cmt-hdur)
+ t0 = - 2.0d0 * minval(tshift_src(:) - hdur(:)) ! - 1.5d0 * minval(tshift_src-hdur)
! uses an earlier start time if source is acoustic with a gaussian source time function
t0_acoustic = 0.0d0
@@ -155,7 +161,7 @@
ispec = ispec_selected_source(isource)
if( ispec_is_acoustic(ispec) ) then
! uses an earlier start time
- t0_acoustic = - 3.0d0 * ( tshift_cmt(isource) - hdur(isource) )
+ t0_acoustic = - 3.0d0 * ( tshift_src(isource) - hdur(isource) )
if( t0_acoustic > t0 ) t0 = t0_acoustic
endif
endif
@@ -168,10 +174,10 @@
! point force sources will start depending on the frequency given by hdur
if( USE_FORCE_POINT_SOURCE .or. USE_RICKER_IPATI ) then
! note: point force sources will give the dominant frequency in hdur,
- ! thus the main period is 1/hdur.
- ! also, these sources use a Ricker source time function instead of a gaussian.
- ! for a Ricker source time function, a start time ~1.2 * main_period is a good choice
- t0 = - 1.2d0 * minval(tshift_cmt(:) - 1.0d0/hdur(:))
+ ! thus the main period is 1/hdur.
+ ! also, these sources use a Ricker source time function instead of a gaussian.
+ ! for a Ricker source time function, a start time ~1.2 * main_period is a good choice
+ t0 = - 1.2d0 * minval(tshift_src(:) - 1.0d0/hdur(:))
endif
! checks if user set USER_T0 to fix simulation start time
@@ -184,17 +190,17 @@
! notifies user
if( myrank == 0 ) then
write(IMAIN,*) 'USER_T0: ',USER_T0
- write(IMAIN,*) 't0: ',t0,'min_tshift_cmt_original: ',min_tshift_cmt_original
+ write(IMAIN,*) 't0: ',t0,'min_tshift_src_original: ',min_tshift_src_original
write(IMAIN,*)
endif
! checks if automatically set t0 is too small
- ! note: min_tshift_cmt_original can be a positive or negative time shift (minimum from all tshift)
- if( t0 <= USER_T0 + min_tshift_cmt_original ) then
- ! by default, tshift_cmt(:) holds relative time shifts with a minimum time shift set to zero
+ ! note: min_tshift_src_original can be a positive or negative time shift (minimum from all tshift)
+ if( t0 <= USER_T0 + min_tshift_src_original ) then
+ ! by default, tshift_src(:) holds relative time shifts with a minimum time shift set to zero
! re-adds (minimum) original time shift such that sources will kick in
! according to their absolute time shift
- tshift_cmt(:) = tshift_cmt(:) + min_tshift_cmt_original
+ tshift_src(:) = tshift_src(:) + min_tshift_src_original
! sets new simulation start time such that
! simulation starts at t = - t0 = - USER_T0
@@ -211,7 +217,7 @@
if( myrank == 0 ) then
write(IMAIN,*) 'error: USER_T0 is too small'
write(IMAIN,*) ' must make one of three adjustements:'
- write(IMAIN,*) ' - increase USER_T0 to be at least: ',t0-min_tshift_cmt_original
+ write(IMAIN,*) ' - increase USER_T0 to be at least: ',t0-min_tshift_src_original
write(IMAIN,*) ' - decrease time shift in CMTSOLUTION file'
write(IMAIN,*) ' - decrease hdur in CMTSOLUTION file'
endif
@@ -973,4 +979,3 @@
endif
end subroutine setup_sources_receivers_VTKfile
-
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -102,7 +102,7 @@
double precision, dimension(:,:,:), allocatable :: nu_source
double precision, dimension(:), allocatable :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
double precision, dimension(:), allocatable :: xi_source,eta_source,gamma_source
- double precision, dimension(:), allocatable :: tshift_cmt,hdur,hdur_gaussian
+ double precision, dimension(:), allocatable :: tshift_src,hdur,hdur_gaussian
double precision, dimension(:), allocatable :: utm_x_source,utm_y_source
double precision, external :: comp_source_time_function
double precision :: t0
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90 2012-10-17 12:33:40 UTC (rev 20841)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90 2012-10-17 16:52:57 UTC (rev 20842)
@@ -190,7 +190,7 @@
etax(:,:,:,ispec),etay(:,:,:,ispec),etaz(:,:,:,ispec), &
gammax(:,:,:,ispec),gammay(:,:,:,ispec),gammaz(:,:,:,ispec))
- stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_cmt(irec),hdur_gaussian(irec))
+ stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-tshift_src(irec),hdur_gaussian(irec))
stf_deltat = stf * deltat
Mxx_der(irec_local) = Mxx_der(irec_local) + eps_s(1,1) * stf_deltat
Myy_der(irec_local) = Myy_der(irec_local) + eps_s(2,2) * stf_deltat
More information about the CIG-COMMITS
mailing list