[cig-commits] r17926 - in seismo/2D/SPECFEM2D/trunk: . DATA EXAMPLES/DATA_to_sort_older_examples EXAMPLES/DATA_to_sort_older_examples/Paul_external_9node_mesh EXAMPLES/M2_UPPA EXAMPLES/Tape2007 EXAMPLES/Tromp2005 EXAMPLES/fluid_solid
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Sun Feb 20 19:36:42 PST 2011
Author: danielpeter
Date: 2011-02-20 19:36:42 -0800 (Sun, 20 Feb 2011)
New Revision: 17926
Added:
seismo/2D/SPECFEM2D/trunk/check_stability.F90
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/DATA/Par_file.in
seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/Par_file_attenuation_2D
seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/Paul_external_9node_mesh/Par_file_CylMeshEnfouiCMor_9n
seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/SOURCE.adj
seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/SOURCE_attenuation_2D
seismo/2D/SPECFEM2D/trunk/EXAMPLES/M2_UPPA/Par_file_M2_UPPA
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_132rec
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005_s100
seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/Par_file_fluid_solid
seismo/2D/SPECFEM2D/trunk/Makefile.in
seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt
seismo/2D/SPECFEM2D/trunk/checkgrid.F90
seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_poro_fluid.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_poro_solid.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_viscoelastic.f90
seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90
seismo/2D/SPECFEM2D/trunk/locate_source_force.F90
seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/plotpost.F90
seismo/2D/SPECFEM2D/trunk/prepare_initialfield.F90
seismo/2D/SPECFEM2D/trunk/prepare_source_time_function.f90
seismo/2D/SPECFEM2D/trunk/read_databases.f90
seismo/2D/SPECFEM2D/trunk/read_parameter_file.F90
seismo/2D/SPECFEM2D/trunk/read_source_file.f90
seismo/2D/SPECFEM2D/trunk/save_databases.f90
seismo/2D/SPECFEM2D/trunk/set_sources.f90
seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
updates source parameters in SPECFEM2D/; adds new file check_stability.F90
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2011-02-21 03:36:42 UTC (rev 17926)
@@ -47,7 +47,7 @@
deltat = 1.d-3 # duration of a time step
# source parameters
-NSOURCE = 1 # number of sources [source info read in CMTSOLUTION file]
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# constants for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file.in 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file.in 2011-02-21 03:36:42 UTC (rev 17926)
@@ -47,7 +47,7 @@
deltat = 1.d-3 # duration of a time step
# source parameters
-NSOURCE = 1 # number of sources [source info read in CMTSOLUTION file]
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# constants for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/Par_file_attenuation_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/Par_file_attenuation_2D 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/Par_file_attenuation_2D 2011-02-21 03:36:42 UTC (rev 17926)
@@ -47,7 +47,7 @@
deltat = 7.5e-4 # duration of a time step
# source parameters
-NSOURCE = 1 # number of sources [source info read in CMTSOLUTION file]
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# constants for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/Paul_external_9node_mesh/Par_file_CylMeshEnfouiCMor_9n
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/Paul_external_9node_mesh/Par_file_CylMeshEnfouiCMor_9n 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/Paul_external_9node_mesh/Par_file_CylMeshEnfouiCMor_9n 2011-02-21 03:36:42 UTC (rev 17926)
@@ -46,7 +46,7 @@
deltat = 1.d-4 # duration of a time step
# source parameters
-NSOURCE = 1 # number of sources [source info read in CMTSOLUTION file]
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# constants for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/SOURCE.adj
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/SOURCE.adj 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/SOURCE.adj 2011-02-21 03:36:42 UTC (rev 17926)
@@ -5,7 +5,7 @@
source_type = 1 # elastic force or acoustic pressure = 1 or moment tensor = 2
time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10 # dominant source frequency (Hz) if not Dirac or Heaviside
-tshift = 0. # offset of the source, irrelevant if NSOURCE=1
+tshift = 0. # offset of the source, irrelevant if NSOURCES=1
angleforce = -90. # angle of the source (for a force only)
Mxx = 1. # Mxx component (for a moment tensor source only)
Mzz = 1. # Mzz component (for a moment tensor source only)
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/SOURCE_attenuation_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/SOURCE_attenuation_2D 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/DATA_to_sort_older_examples/SOURCE_attenuation_2D 2011-02-21 03:36:42 UTC (rev 17926)
@@ -5,7 +5,7 @@
source_type = 1 # elastic force or acoustic pressure = 1 or moment tensor = 2
time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 18.0 # dominant source frequency (Hz) if not Dirac or Heaviside
-tshift_src = 0.0 # time shift when multi sources (if one source, must be zero)
+tshift = 0.0 # time shift when multi sources (if one source, must be zero)
angleforce = 0. # angle of the source (for a force only)
Mxx = 1. # Mxx component (for a moment tensor source only)
Mzz = 1. # Mzz component (for a moment tensor source only)
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/M2_UPPA/Par_file_M2_UPPA
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/M2_UPPA/Par_file_M2_UPPA 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/M2_UPPA/Par_file_M2_UPPA 2011-02-21 03:36:42 UTC (rev 17926)
@@ -47,7 +47,7 @@
deltat = 1.d-3 # duration of a time step
# source parameters
-NSOURCE = 1 # number of sources [source info read in CMTSOLUTION file]
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# constants for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_132rec
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_132rec 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_132rec 2011-02-21 03:36:42 UTC (rev 17926)
@@ -47,7 +47,7 @@
deltat = 6.0d-2 # duration of a time step
# source parameters
-NSOURCE = 1 # number of sources [source info read in CMTSOLUTION file]
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# constants for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007/Par_file_Tape2007_onerec 2011-02-21 03:36:42 UTC (rev 17926)
@@ -47,7 +47,7 @@
deltat = 6.0d-2 # duration of a time step
# source parameters
-NSOURCE = 1 # number of sources [source info read in CMTSOLUTION file]
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# constants for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005 2011-02-21 03:36:42 UTC (rev 17926)
@@ -47,7 +47,7 @@
deltat = 2.0d-2 # duration of a time step
# source parameters
-NSOURCE = 1 # number of sources [source info read in CMTSOLUTION file]
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# constants for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005_s100
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005_s100 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tromp2005/Par_file_Tromp2005_s100 2011-02-21 03:36:42 UTC (rev 17926)
@@ -47,7 +47,7 @@
deltat = 2.0d-4 # duration of a time step
# source parameters
-NSOURCE = 1 # number of sources [source info read in CMTSOLUTION file]
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# constants for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/Par_file_fluid_solid
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/Par_file_fluid_solid 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/fluid_solid/Par_file_fluid_solid 2011-02-21 03:36:42 UTC (rev 17926)
@@ -47,7 +47,7 @@
deltat = 0.25d-3 # duration of a time step
# source parameters
-NSOURCE = 1 # number of sources [source info read in CMTSOLUTION file]
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# constants for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/Makefile.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile.in 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/Makefile.in 2011-02-21 03:36:42 UTC (rev 17926)
@@ -125,6 +125,7 @@
$O/attenuation_compute_param.o \
$O/calendar.o \
$O/checkgrid.o \
+ $O/check_stability.o \
$O/compute_arrays_source.o \
$O/compute_Bielak_conditions.o \
$O/compute_curl_one_element.o \
@@ -360,6 +361,9 @@
$O/compute_gradient_attenuation.o: compute_gradient_attenuation.f90 constants.h
${F90} $(FLAGS_NO_CHECK) -c -o $O/compute_gradient_attenuation.o compute_gradient_attenuation.f90
+$O/check_stability.o: check_stability.F90 constants.h
+ ${F90} $(FLAGS_NO_CHECK) -c -o $O/check_stability.o check_stability.F90
+
$O/calendar.o: calendar.f90
${F90} $(FLAGS_CHECK) -c -o $O/calendar.o calendar.f90
Modified: seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt
===================================================================
--- seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt 2011-02-21 03:36:42 UTC (rev 17926)
@@ -379,7 +379,7 @@
In section "# source parameters":
The code now support multi sources.
-NSOURCE is the number of source.
+NSOURCES is the number of source.
Parameters of the sources are displayed in the file SOURCE, which must be
in the directory DATA. The components of a moment tensor source must be given in N.m,
not in dyne.cm as in the DATA/CMTSOLUTION source file of the 3D version of the code.
Added: seismo/2D/SPECFEM2D/trunk/check_stability.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/check_stability.F90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/check_stability.F90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -0,0 +1,302 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 6.1
+! ------------------------------
+!
+! Copyright Universite de Pau, CNRS and INRIA, France,
+! and Princeton University / California Institute of Technology, USA.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
+! Christina Morency, cmorency aT princeton DOT edu
+!
+! This software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
+! using a spectral-element method (SEM).
+!
+! This software is governed by the CeCILL license under French law and
+! abiding by the rules of distribution of free software. You can use,
+! modify and/or redistribute the software under the terms of the CeCILL
+! license as circulated by CEA, CNRS and INRIA at the following URL
+! "http://www.cecill.info".
+!
+! As a counterpart to the access to the source code and rights to copy,
+! modify and redistribute granted by the license, users are provided only
+! with a limited warranty and the software's author, the holder of the
+! economic rights, and the successive licensors have only limited
+! liability.
+!
+! In this respect, the user's attention is drawn to the risks associated
+! with loading, using, modifying and/or developing or reproducing the
+! software by the user in light of its specific status of free software,
+! that may mean that it is complicated to manipulate, and that also
+! therefore means that it is reserved for developers and experienced
+! professionals having in-depth computer knowledge. Users are therefore
+! encouraged to load and test the software's suitability as regards their
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and, more generally, to use and operate it in the
+! same conditions as regards security.
+!
+! The full text of the license is available in file "LICENSE".
+!
+!========================================================================
+
+
+ subroutine check_stability(myrank,time,it,NSTEP,npoin, &
+ any_elastic_glob,any_elastic,displ_elastic, &
+ any_poroelastic_glob,any_poroelastic, &
+ displs_poroelastic,displw_poroelastic, &
+ any_acoustic_glob,any_acoustic,potential_acoustic, &
+ year_start,month_start,time_start)
+
+! checks simulation stability and outputs timerun infos
+
+ implicit none
+ include "constants.h"
+#ifdef USE_MPI
+ include "mpif.h"
+#endif
+
+ integer :: myrank,it,NSTEP
+ integer :: npoin
+
+ double precision :: time
+
+ logical :: any_elastic_glob,any_elastic
+ real(kind=CUSTOM_REAL), dimension(3,npoin) :: displ_elastic
+
+ logical :: any_poroelastic_glob,any_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displs_poroelastic,displw_poroelastic
+
+ logical :: any_acoustic_glob,any_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin) :: potential_acoustic
+
+ double precision :: time_start
+ integer :: year_start,month_start
+
+ ! local parameters
+ double precision displnorm_all,displnorm_all_glob
+ ! timer to count elapsed time
+ double precision :: time_end
+ integer :: year_end,month_end
+ double precision :: tCPU,t_remain,t_total
+ integer :: ihours,iminutes,iseconds,int_tCPU, &
+ ihours_remain,iminutes_remain,iseconds_remain,int_t_remain, &
+ ihours_total,iminutes_total,iseconds_total,int_t_total
+ ! to determine date and time at which the run will finish
+ character(len=8) datein
+ character(len=10) timein
+ character(len=5) :: zone
+ integer, dimension(8) :: time_values
+ character(len=3), dimension(12) :: month_name
+ character(len=3), dimension(0:6) :: weekday_name
+ data month_name /'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/
+ data weekday_name /'Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat'/
+ integer :: year,mon,day,hr,minutes,timestamp,julian_day_number,day_of_week
+ integer, external :: idaywk
+#ifdef USE_MPI
+ integer :: ier
+#endif
+
+ ! user output
+ if (myrank == 0) then
+ write(IOUT,*)
+ if(time >= 1.d-3 .and. time < 1000.d0) then
+ write(IOUT,"('Time step number ',i7,' t = ',f9.4,' s out of ',i7)") it,time,NSTEP
+ else
+ write(IOUT,"('Time step number ',i7,' t = ',1pe12.6,' s out of ',i7)") it,time,NSTEP
+ endif
+ write(IOUT,*) 'We have done ',sngl(100.d0*dble(it-1)/dble(NSTEP-1)),'% of the total'
+ endif
+
+
+ ! elastic wavefield
+ if(any_elastic_glob) then
+ if(any_elastic) then
+ displnorm_all = maxval(sqrt(displ_elastic(1,:)**2 &
+ + displ_elastic(2,:)**2 &
+ + displ_elastic(3,:)**2))
+ else
+ displnorm_all = 0.d0
+ endif
+
+ displnorm_all_glob = displnorm_all
+#ifdef USE_MPI
+ call MPI_ALLREDUCE (displnorm_all, displnorm_all_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+#endif
+
+ if (myrank == 0) &
+ write(IOUT,*) 'Max norm of vector field in solid (elastic) = ',displnorm_all_glob
+
+ ! check stability of the code in solid, exit if unstable
+ ! negative values can occur with some compilers when the unstable value is greater
+ ! than the greatest possible floating-point number of the machine
+ if(displnorm_all_glob > STABILITY_THRESHOLD .or. displnorm_all_glob < 0) &
+ call exit_MPI('code became unstable and blew up in solid (elastic)')
+
+ endif
+
+ ! poroelastic wavefield
+ if(any_poroelastic_glob) then
+ if(any_poroelastic) then
+ displnorm_all = maxval(sqrt(displs_poroelastic(1,:)**2 &
+ + displs_poroelastic(2,:)**2))
+ else
+ displnorm_all = 0.d0
+ endif
+
+ displnorm_all_glob = displnorm_all
+#ifdef USE_MPI
+ call MPI_ALLREDUCE (displnorm_all, displnorm_all_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+#endif
+
+ if (myrank == 0) &
+ write(IOUT,*) 'Max norm of vector field in solid (poroelastic) = ',displnorm_all_glob
+
+ ! check stability of the code in solid, exit if unstable
+ ! negative values can occur with some compilers when the unstable value is greater
+ ! than the greatest possible floating-point number of the machine
+ if(displnorm_all_glob > STABILITY_THRESHOLD .or. displnorm_all_glob < 0) &
+ call exit_MPI('code became unstable and blew up in solid (poroelastic)')
+
+ if(any_poroelastic) then
+ displnorm_all = maxval(sqrt(displw_poroelastic(1,:)**2 &
+ + displw_poroelastic(2,:)**2))
+ else
+ displnorm_all = 0.d0
+ endif
+
+ displnorm_all_glob = displnorm_all
+#ifdef USE_MPI
+ call MPI_ALLREDUCE (displnorm_all, displnorm_all_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+#endif
+
+ if (myrank == 0) &
+ write(IOUT,*) 'Max norm of vector field in fluid (poroelastic) = ',displnorm_all_glob
+
+ ! check stability of the code in solid, exit if unstable
+ ! negative values can occur with some compilers when the unstable value is greater
+ ! than the greatest possible floating-point number of the machine
+ if(displnorm_all_glob > STABILITY_THRESHOLD .or. displnorm_all_glob < 0) &
+ call exit_MPI('code became unstable and blew up in fluid (poroelastic)')
+
+ endif
+
+
+ ! acoustic wavefield
+ if(any_acoustic_glob) then
+ if(any_acoustic) then
+ displnorm_all = maxval(abs(potential_acoustic(:)))
+ else
+ displnorm_all = 0.d0
+ endif
+
+ displnorm_all_glob = displnorm_all
+#ifdef USE_MPI
+ call MPI_ALLREDUCE (displnorm_all, displnorm_all_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+#endif
+
+ if (myrank == 0) &
+ write(IOUT,*) 'Max absolute value of scalar field in fluid (acoustic) = ',displnorm_all_glob
+
+ ! check stability of the code in fluid, exit if unstable
+ ! negative values can occur with some compilers when the unstable value is greater
+ ! than the greatest possible floating-point number of the machine
+ if(displnorm_all_glob > STABILITY_THRESHOLD .or. displnorm_all_glob < 0) &
+ call exit_MPI('code became unstable and blew up in fluid (acoustic)')
+
+ endif
+
+ ! count elapsed wall-clock time
+ call date_and_time(datein,timein,zone,time_values)
+ ! time_values(1): year
+ ! time_values(2): month of the year
+ ! time_values(3): day of the month
+ ! time_values(5): hour of the day
+ ! time_values(6): minutes of the hour
+ ! time_values(7): seconds of the minute
+ ! time_values(8): milliseconds of the second
+ ! this fails if we cross the end of the month
+ time_end = 86400.d0*time_values(3) + 3600.d0*time_values(5) + &
+ 60.d0*time_values(6) + time_values(7) + time_values(8) / 1000.d0
+ month_end = time_values(2)
+ year_end = time_values(1)
+
+ ! elapsed time since beginning of the simulation
+ if (myrank == 0) then
+ if(month_end == month_start .and. year_end == year_start) then
+ tCPU = time_end - time_start
+ int_tCPU = int(tCPU)
+ ihours = int_tCPU / 3600
+ iminutes = (int_tCPU - 3600*ihours) / 60
+ iseconds = int_tCPU - 3600*ihours - 60*iminutes
+ write(IOUT,*) 'Elapsed time in seconds = ',tCPU
+ write(IOUT,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
+ write(IOUT,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
+
+ ! compute estimated remaining simulation time
+ t_remain = (NSTEP - it) * (tCPU/dble(it))
+ int_t_remain = int(t_remain)
+ ihours_remain = int_t_remain / 3600
+ iminutes_remain = (int_t_remain - 3600*ihours_remain) / 60
+ iseconds_remain = int_t_remain - 3600*ihours_remain - 60*iminutes_remain
+ write(IOUT,*) 'Time steps remaining = ',NSTEP - it
+ write(IOUT,*) 'Estimated remaining time in seconds = ',t_remain
+ write(IOUT,"(' Estimated remaining time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
+ ihours_remain,iminutes_remain,iseconds_remain
+
+ ! compute estimated total simulation time
+ t_total = t_remain + tCPU
+ int_t_total = int(t_total)
+ ihours_total = int_t_total / 3600
+ iminutes_total = (int_t_total - 3600*ihours_total) / 60
+ iseconds_total = int_t_total - 3600*ihours_total - 60*iminutes_total
+ write(IOUT,*) 'Estimated total run time in seconds = ',t_total
+ write(IOUT,"(' Estimated total run time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
+ ihours_total,iminutes_total,iseconds_total
+
+ if(it < NSTEP) then
+ ! compute date and time at which the run should finish
+ ! (useful for long runs); for simplicity only minutes
+ ! are considered, seconds are ignored; in any case the prediction is not
+ ! accurate down to seconds because of system and network fluctuations
+ year = time_values(1)
+ mon = time_values(2)
+ day = time_values(3)
+ hr = time_values(5)
+ minutes = time_values(6)
+
+ ! get timestamp in minutes of current date and time
+ call convtime(timestamp,year,mon,day,hr,minutes)
+
+ ! add remaining minutes
+ timestamp = timestamp + nint(t_remain / 60.d0)
+
+ ! get date and time of that future timestamp in minutes
+ call invtime(timestamp,year,mon,day,hr,minutes)
+
+ ! convert to Julian day to get day of the week
+ call calndr(day,mon,year,julian_day_number)
+ day_of_week = idaywk(julian_day_number)
+
+ write(IOUT,"(' The run will finish approximately on: ',a3,' ',a3,' ',i2.2,', ',i4.4,' ',i2.2,':',i2.2)") &
+ weekday_name(day_of_week),month_name(mon),day,year,hr,minutes
+
+ endif
+ write(IOUT,*)
+ else
+ write(IOUT,*) 'The calendar has crossed the end of the month during the simulation,'
+ write(IOUT,*) 'cannot produce accurate CPU time estimates any more.'
+ write(IOUT,*)
+ endif
+ endif
+
+ if (myrank == 0) write(IOUT,*)
+
+ end subroutine check_stability
+
Modified: seismo/2D/SPECFEM2D/trunk/checkgrid.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -42,14 +42,15 @@
!
!========================================================================
- subroutine checkgrid(vpext,vsext,rhoext,density,poroelastcoef,porosity,tortuosity,permeability,ibool,kmato, &
- coord,npoin,vpImin,vpImax,vpIImin,vpIImax, &
- assign_external_model,nspec,UPPER_LIMIT_DISPLAY,numat,deltat, &
- f0,tshift_src,initialfield,time_function_type, &
- coorg,xinterp,zinterp,shapeint,knods,simulation_title, &
- npgeo,pointsdisp,ngnod,any_elastic,any_poroelastic,all_anisotropic, &
- myrank,nproc,NSOURCE,poroelastic, &
- freq0,Q0,TURN_VISCATTENUATION_ON)
+ subroutine checkgrid(vpext,vsext,rhoext,density,poroelastcoef, &
+ porosity,tortuosity,permeability,ibool,kmato, &
+ coord,npoin,vpImin,vpImax,vpIImin,vpIImax, &
+ assign_external_model,nspec,UPPER_LIMIT_DISPLAY,numat,deltat, &
+ f0,initialfield,time_function_type, &
+ coorg,xinterp,zinterp,shapeint,knods,simulation_title, &
+ npgeo,pointsdisp,ngnod,any_elastic,any_poroelastic,all_anisotropic, &
+ myrank,nproc,NSOURCES,poroelastic, &
+ freq0,Q0,TURN_VISCATTENUATION_ON)
! check the mesh, stability and number of points per wavelength
@@ -64,18 +65,7 @@
! for instance to analyze Cuthill-McKee mesh partitioning etc.
integer :: UPPER_LIMIT_DISPLAY
-! color palette
- integer, parameter :: NUM_COLORS = 236
- double precision, dimension(NUM_COLORS) :: red,green,blue
-
-#ifdef USE_MPI
- integer :: icol
-#endif
-
- integer i,j,ispec,material,npoin,nspec,numat,NSOURCE
-
- integer, dimension(NSOURCE) :: time_function_type
-
+ integer :: npoin,nspec,numat
integer, dimension(nspec) :: kmato
logical, dimension(nspec) :: poroelastic
integer, dimension(NGLLX,NGLLX,nspec) :: ibool
@@ -88,8 +78,34 @@
double precision coord(NDIM,npoin)
- double precision vpImin,vpImax,vsmin,vsmax,densmin,densmax,vpImax_local,vpImin_local,vsmin_local
- double precision vpIImin,vpIImax,vpIImax_local,vpIImin_local
+ integer :: NSOURCES
+ integer, dimension(NSOURCES) :: time_function_type
+ double precision, dimension(NSOURCES) :: f0
+
+ integer :: pointsdisp,npgeo,ngnod
+
+ integer :: knods(ngnod,nspec)
+
+ double precision :: xinterp(pointsdisp,pointsdisp),zinterp(pointsdisp,pointsdisp)
+ double precision :: shapeint(ngnod,pointsdisp,pointsdisp)
+
+ double precision :: coorg(NDIM,npgeo)
+
+! title of the plot
+ character(len=60) :: simulation_title
+
+ double precision :: vpImin,vpImax
+ double precision :: vpIImin,vpIImax
+ double precision :: deltat
+
+ logical :: assign_external_model,initialfield,any_elastic,any_poroelastic,all_anisotropic, &
+ TURN_VISCATTENUATION_ON
+
+ integer :: myrank,nproc
+
+ ! local parameters
+ double precision vpIImax_local,vpIImin_local
+ double precision vsmin,vsmax,densmin,densmax,vpImax_local,vpImin_local,vsmin_local
double precision kappa_s,kappa_f,kappa_fr,mu_s,mu_fr,denst_s,denst_f,denst,phi,tort,cpIloc,cpIIloc,csloc
double precision D_biot,H_biot,C_biot,M_biot,cpIsquare,cpIIsquare,cssquare
double precision f0min,f0max,freq0,Q0,w_c,eta_f,perm
@@ -97,22 +113,22 @@
double precision distance_min,distance_max,distance_min_local,distance_max_local
double precision courant_stability_number_max,lambdaPImin,lambdaPImax,lambdaPIImin,lambdaPIImax, &
lambdaSmin,lambdaSmax
- double precision deltat,distance_1,distance_2,distance_3,distance_4
- double precision, dimension(NSOURCE) :: f0,tshift_src
- logical assign_external_model,initialfield,any_elastic,any_poroelastic,all_anisotropic, &
- TURN_VISCATTENUATION_ON
+ double precision distance_1,distance_2,distance_3,distance_4
! for the stability condition
! maximum polynomial degree for which we can compute the stability condition
integer, parameter :: NGLLX_MAX_STABILITY = 15
double precision :: percent_GLL(NGLLX_MAX_STABILITY)
- integer pointsdisp,npgeo,ngnod,is,ir,in,nnum
+! color palette
+ integer, parameter :: NUM_COLORS = 236
+ double precision, dimension(NUM_COLORS) :: red,green,blue
double precision :: xmax,zmax,height,usoffset,sizex,sizez,courant_stability_number
double precision :: x1,z1,x2,z2,ratio_page,xmin,zmin,lambdaS_local,lambdaPI_local
#ifdef USE_MPI
+ integer :: icol
double precision :: vpImin_glob,vpImax_glob,vsmin_glob,vsmax_glob,densmin_glob,densmax_glob
double precision :: vpIImin_glob,vpIImax_glob
double precision :: distance_min_glob,distance_max_glob
@@ -128,28 +144,21 @@
integer, dimension(:), allocatable :: RGB_recv
real, dimension(nspec) :: greyscale_send
real, dimension(:), allocatable :: greyscale_recv
- integer :: nspec_recv
- integer :: num_ispec
- integer :: myrank, iproc, nproc
- integer :: ier
+ integer :: nspec_recv
+ integer :: num_ispec
+ integer :: iproc
+ integer :: ier
+ integer :: i,j,ispec,material
+ integer :: is,ir,in,nnum
#ifdef USE_MPI
integer, dimension(MPI_STATUS_SIZE) :: request_mpi_status
#endif
- integer knods(ngnod,nspec)
+ ! check
+ if(UPPER_LIMIT_DISPLAY > nspec) &
+ call exit_MPI('cannot have UPPER_LIMIT_DISPLAY > nspec in checkgrid.F90')
- double precision xinterp(pointsdisp,pointsdisp),zinterp(pointsdisp,pointsdisp)
- double precision shapeint(ngnod,pointsdisp,pointsdisp)
-
- double precision coorg(NDIM,npgeo)
-
-! title of the plot
- character(len=60) simulation_title
-
-
- if(UPPER_LIMIT_DISPLAY > nspec) stop 'cannot have UPPER_LIMIT_DISPLAY > nspec in checkgrid.F90'
-
#ifndef USE_MPI
allocate(coorg_recv(1,1))
allocate(RGB_recv(1))
@@ -347,26 +356,44 @@
any_elastic_glob = any_elastic
any_poroelastic_glob = any_poroelastic
#ifdef USE_MPI
- call MPI_ALLREDUCE (vpImin, vpImin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (vpImax, vpImax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (vpIImin, vpIImin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (vpIImax, vpIImax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (vsmin, vsmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (vsmax, vsmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (densmin, densmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (densmax, densmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (distance_min, distance_min_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (distance_max, distance_max_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (vpImin, vpImin_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (vpImax, vpImax_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (vpIImin, vpIImin_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (vpIImax, vpIImax_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (vsmin, vsmin_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (vsmax, vsmax_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (densmin, densmin_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (densmax, densmax_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (distance_min, distance_min_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (distance_max, distance_max_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
call MPI_ALLREDUCE (courant_stability_number_max, courant_stability_max_glob, 1, MPI_DOUBLE_PRECISION, &
- MPI_MAX, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (lambdaPImin, lambdaPImin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (lambdaPImax, lambdaPImax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (lambdaPIImin, lambdaPIImin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (lambdaPIImax, lambdaPIImax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (lambdaSmin, lambdaSmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (lambdaSmax, lambdaSmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (any_elastic, any_elastic_glob, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ier)
- call MPI_ALLREDUCE (any_poroelastic, any_poroelastic_glob, 1, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ier)
+ MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (lambdaPImin, lambdaPImin_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (lambdaPImax, lambdaPImax_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (lambdaPIImin, lambdaPIImin_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (lambdaPIImax, lambdaPIImax_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (lambdaSmin, lambdaSmin_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MIN, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (lambdaSmax, lambdaSmax_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MAX, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (any_elastic, any_elastic_glob, 1, MPI_LOGICAL, &
+ MPI_LOR, MPI_COMM_WORLD, ier)
+ call MPI_ALLREDUCE (any_poroelastic, any_poroelastic_glob, 1, MPI_LOGICAL, &
+ MPI_LOR, MPI_COMM_WORLD, ier)
vpImin = vpImin_glob
vpImax = vpImax_glob
vpIImin = vpIImin_glob
@@ -416,23 +443,27 @@
if(.not. initialfield) then
f0max = -HUGEVAL
f0min = HUGEVAL
- write(IOUT,*) ' USER_T0 = ',USER_T0
- do i = 1,NSOURCE
+! write(IOUT,*) ' USER_T0 = ',USER_T0
+
+ do i = 1,NSOURCES
+
+ ! excludes Dirac and Heaviside sources
if(time_function_type(i) /= 4 .and. time_function_type(i) /= 5) then
- write(IOUT,*) ' Onset time = ',tshift_src(i)
- write(IOUT,*) ' Fundamental period = ',1.d0/f0(i)
- write(IOUT,*) ' Fundamental frequency = ',f0(i)
+! write(IOUT,*) ' Onset time = ',t0+tshift_src(i)
+! write(IOUT,*) ' Fundamental period = ',1.d0/f0(i)
+! write(IOUT,*) ' Fundamental frequency = ',f0(i)
+! ! checks source onset time
+! if( t0+tshift_src(i) <= 1.d0/f0(i)) then
+! call exit_MPI('Onset time too small')
+! else
+! write(IOUT,*) ' --> onset time ok'
+! endif
+
! sets min/max frequency
if(f0(i) > f0max) f0max = f0(i)
if(f0(i) < f0min) f0min = f0(i)
- ! checks source onset time
- if(tshift_src(i) <= 1.d0/f0(i)) then
- call exit_MPI('Onset time too small')
- else
- write(IOUT,*) ' --> onset time ok'
- endif
- if( i == NSOURCE ) then
+ if( i == NSOURCES ) then
write(IOUT,*) '----'
write(IOUT,*) ' Nb pts / lambdaPImin_fmax max = ',lambdaPImax/(2.5d0*f0min)
write(IOUT,*) ' Nb pts / lambdaPImin_fmax min = ',lambdaPImin/(2.5d0*f0max)
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -46,7 +46,7 @@
anyabs,assign_external_model,ibool,kmato,numabs, &
elastic,poroelastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,b_potential_dot_dot_acoustic,b_potential_acoustic, &
- density,elastcoef,xix,xiz,gammax,gammaz,jacobian, &
+ density,poroelastcoef,xix,xiz,gammax,gammaz,jacobian, &
vpext,rhoext,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
@@ -64,10 +64,10 @@
integer :: npoin,nspec,nelemabs,numat,it,NSTEP,SIMULATION_TYPE
integer :: nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax
- integer, dimension(nspec_xmin) :: ib_left
- integer, dimension(nspec_xmax) :: ib_right
- integer, dimension(nspec_zmin) :: ib_bottom
- integer, dimension(nspec_zmax) :: ib_top
+ integer, dimension(nelemabs) :: ib_left
+ integer, dimension(nelemabs) :: ib_right
+ integer, dimension(nelemabs) :: ib_bottom
+ integer, dimension(nelemabs) :: ib_top
logical :: anyabs,assign_external_model
logical :: SAVE_FORWARD
@@ -80,10 +80,12 @@
logical, dimension(nspec) :: elastic,poroelastic
logical, dimension(4,nelemabs) :: codeabs
- real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
- real(kind=CUSTOM_REAL), dimension(npoin) :: b_potential_dot_dot_acoustic,b_potential_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin) :: &
+ potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin) :: &
+ b_potential_dot_dot_acoustic,b_potential_acoustic
double precision, dimension(2,numat) :: density
- double precision, dimension(4,3,numat) :: elastcoef
+ double precision, dimension(4,3,numat) :: poroelastcoef
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,rhoext
@@ -229,8 +231,8 @@
ispec = numabs(ispecabs)
! get elastic parameters of current spectral element
- lambdal_relaxed = elastcoef(1,1,kmato(ispec))
- mul_relaxed = elastcoef(2,1,kmato(ispec))
+ lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
+ mul_relaxed = poroelastcoef(2,1,kmato(ispec))
kappal = lambdal_relaxed + TWO*mul_relaxed/3._CUSTOM_REAL
rhol = density(1,kmato(ispec))
@@ -411,3 +413,338 @@
end subroutine compute_forces_acoustic
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine compute_forces_acoustic_2(npoin,nspec,nelemabs,numat,it,NSTEP, &
+ anyabs,assign_external_model,ibool,kmato,numabs, &
+ elastic,poroelastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
+ potential_acoustic, &
+ density,poroelastcoef,xix,xiz,gammax,gammaz,jacobian, &
+ vpext,rhoext,hprime_xx,hprimewgll_xx, &
+ hprime_zz,hprimewgll_zz,wxgll,wzgll, &
+ ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ jbegin_left,jend_left,jbegin_right,jend_right, &
+ SIMULATION_TYPE,SAVE_FORWARD,nspec_xmin,nspec_xmax,&
+ nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top, &
+ b_absorb_acoustic_left,b_absorb_acoustic_right, &
+ b_absorb_acoustic_bottom,b_absorb_acoustic_top)
+
+! compute forces for the acoustic elements
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: npoin,nspec,nelemabs,numat,it,NSTEP,SIMULATION_TYPE
+
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+ integer, dimension(nspec) :: kmato
+ integer, dimension(nelemabs) :: numabs,ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ jbegin_left,jend_left,jbegin_right,jend_right
+
+ logical, dimension(nspec) :: elastic,poroelastic
+ logical, dimension(4,nelemabs) :: codeabs
+
+ real(kind=CUSTOM_REAL), dimension(npoin) :: &
+ potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
+
+ double precision, dimension(2,numat) :: density
+ double precision, dimension(4,3,numat) :: poroelastcoef
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,rhoext
+
+ logical :: anyabs,assign_external_model
+ logical :: SAVE_FORWARD
+
+! derivatives of Lagrange polynomials
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
+
+! Gauss-Lobatto-Legendre weights
+ real(kind=CUSTOM_REAL), dimension(NGLLX) :: wxgll
+ real(kind=CUSTOM_REAL), dimension(NGLLZ) :: wzgll
+
+ integer :: nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax
+ integer, dimension(nelemabs) :: ib_left
+ integer, dimension(nelemabs) :: ib_right
+ integer, dimension(nelemabs) :: ib_bottom
+ integer, dimension(nelemabs) :: ib_top
+
+ double precision, dimension(NGLLZ,nspec_xmin,NSTEP) :: b_absorb_acoustic_left
+ double precision, dimension(NGLLZ,nspec_xmax,NSTEP) :: b_absorb_acoustic_right
+ double precision, dimension(NGLLX,nspec_zmax,NSTEP) :: b_absorb_acoustic_top
+ double precision, dimension(NGLLX,nspec_zmin,NSTEP) :: b_absorb_acoustic_bottom
+
+!---
+!--- local variables
+!---
+
+ integer :: ispec,i,j,k,iglob,ispecabs,ibegin,iend,jbegin,jend
+
+! spatial derivatives
+ real(kind=CUSTOM_REAL) :: dux_dxi,dux_dgamma,dux_dxl,dux_dzl
+ real(kind=CUSTOM_REAL) :: weight,xxi,zxi,xgamma,zgamma,jacobian1D
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: tempx1,tempx2
+
+! Jacobian matrix and determinant
+ real(kind=CUSTOM_REAL) :: xixl,xizl,gammaxl,gammazl,jacobianl
+
+! material properties of the elastic medium
+ real(kind=CUSTOM_REAL) :: mul_relaxed,lambdal_relaxed,kappal,cpl,rhol
+
+ integer :: ifirstelem,ilastelem
+
+ ifirstelem = 1
+ ilastelem = nspec
+
+! loop over spectral elements
+ do ispec = ifirstelem,ilastelem
+
+!---
+!--- acoustic spectral element
+!---
+ if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
+
+ rhol = density(1,kmato(ispec))
+
+ ! first double loop over GLL points to compute and store gradients
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+
+ ! derivative along x and along z
+ dux_dxi = ZERO
+ dux_dgamma = ZERO
+
+ ! first double loop over GLL points to compute and store gradients
+ ! we can merge the two loops because NGLLX == NGLLZ
+ do k = 1,NGLLX
+ dux_dxi = dux_dxi + potential_acoustic(ibool(k,j,ispec))*hprime_xx(i,k)
+ dux_dgamma = dux_dgamma + potential_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
+ enddo
+
+ xixl = xix(i,j,ispec)
+ xizl = xiz(i,j,ispec)
+ gammaxl = gammax(i,j,ispec)
+ gammazl = gammaz(i,j,ispec)
+
+ ! derivatives of potential
+ dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
+ dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
+ jacobianl = jacobian(i,j,ispec)
+
+ ! if external density model
+ if(assign_external_model) rhol = rhoext(i,j,ispec)
+
+ ! for acoustic medium
+ ! also add GLL integration weights
+ tempx1(i,j) = wzgll(j)*jacobianl*(xixl*dux_dxl + xizl*dux_dzl) / rhol
+ tempx2(i,j) = wxgll(i)*jacobianl*(gammaxl*dux_dxl + gammazl*dux_dzl) / rhol
+ enddo
+ enddo
+
+!
+! second double-loop over GLL to compute all the terms
+!
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+
+ iglob = ibool(i,j,ispec)
+
+ ! along x direction and z direction
+ ! and assemble the contributions
+ do k = 1,NGLLX
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - &
+ (tempx1(k,j)*hprimewgll_xx(k,i) + tempx2(i,k)*hprimewgll_zz(k,j))
+ enddo
+
+ enddo ! second loop over the GLL points
+ enddo
+
+ endif ! end of test if acoustic element
+
+ enddo ! end of loop over all spectral elements
+
+!
+!--- absorbing boundaries
+!
+ if(anyabs) then
+
+ do ispecabs=1,nelemabs
+
+ ispec = numabs(ispecabs)
+
+ ! Sommerfeld condition if acoustic
+ if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
+
+ ! get elastic parameters of current spectral element
+ lambdal_relaxed = poroelastcoef(1,1,kmato(ispec))
+ mul_relaxed = poroelastcoef(2,1,kmato(ispec))
+ kappal = lambdal_relaxed + TWO*mul_relaxed/3._CUSTOM_REAL
+ rhol = density(1,kmato(ispec))
+
+ cpl = sqrt(kappal/rhol)
+
+ !--- left absorbing boundary
+ if(codeabs(ILEFT,ispecabs)) then
+ i = 1
+ jbegin = jbegin_left(ispecabs)
+ jend = jend_left(ispecabs)
+ do j = jbegin,jend
+ iglob = ibool(i,j,ispec)
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
+ xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
+ zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xgamma**2 + zgamma**2)
+ weight = jacobian1D * wzgll(j)
+
+ if( SIMULATION_TYPE == 1 ) then
+ ! adds absorbing boundary contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ elseif(SIMULATION_TYPE == 2) then
+ ! adds (previously) stored contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - b_absorb_acoustic_left(j,ib_left(ispecabs),NSTEP-it+1)
+ endif
+
+ if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
+ ! saves contribution
+ b_absorb_acoustic_left(j,ib_left(ispecabs),it) = &
+ potential_dot_acoustic(iglob)*weight/cpl/rhol
+ endif
+
+ enddo
+
+ endif ! end of left absorbing boundary
+
+ !--- right absorbing boundary
+ if(codeabs(IRIGHT,ispecabs)) then
+ i = NGLLX
+ jbegin = jbegin_right(ispecabs)
+ jend = jend_right(ispecabs)
+ do j = jbegin,jend
+ iglob = ibool(i,j,ispec)
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
+ xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
+ zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xgamma**2 + zgamma**2)
+ weight = jacobian1D * wzgll(j)
+
+ if( SIMULATION_TYPE == 1 ) then
+ ! adds absorbing boundary contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ elseif(SIMULATION_TYPE == 2) then
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - b_absorb_acoustic_right(j,ib_right(ispecabs),NSTEP-it+1)
+ endif
+
+ if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
+ ! saves contribution
+ b_absorb_acoustic_right(j,ib_right(ispecabs),it) = &
+ potential_dot_acoustic(iglob)*weight/cpl/rhol
+ endif
+ enddo
+ endif ! end of right absorbing boundary
+
+ !--- bottom absorbing boundary
+ if(codeabs(IBOTTOM,ispecabs)) then
+ j = 1
+ ibegin = ibegin_bottom(ispecabs)
+ iend = iend_bottom(ispecabs)
+ ! exclude corners to make sure there is no contradiction on the normal
+ if(codeabs(ILEFT,ispecabs)) ibegin = 2
+ if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
+ do i = ibegin,iend
+ iglob = ibool(i,j,ispec)
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
+ xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
+ zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xxi**2 + zxi**2)
+ weight = jacobian1D * wxgll(i)
+
+ if( SIMULATION_TYPE == 1 ) then
+ ! adds absorbing boundary contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ elseif(SIMULATION_TYPE == 2) then
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - b_absorb_acoustic_bottom(i,ib_bottom(ispecabs),NSTEP-it+1)
+ endif
+
+ if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
+ ! saves contribution
+ b_absorb_acoustic_bottom(i,ib_bottom(ispecabs),it) = &
+ potential_dot_acoustic(iglob)*weight/cpl/rhol
+ endif
+ enddo
+ endif ! end of bottom absorbing boundary
+
+ !--- top absorbing boundary
+ if(codeabs(ITOP,ispecabs)) then
+ j = NGLLZ
+ ibegin = ibegin_top(ispecabs)
+ iend = iend_top(ispecabs)
+ ! exclude corners to make sure there is no contradiction on the normal
+ if(codeabs(ILEFT,ispecabs)) ibegin = 2
+ if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
+ do i = ibegin,iend
+ iglob = ibool(i,j,ispec)
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
+ xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
+ zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xxi**2 + zxi**2)
+ weight = jacobian1D * wxgll(i)
+
+ if( SIMULATION_TYPE == 1 ) then
+ ! adds absorbing boundary contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ elseif(SIMULATION_TYPE == 2) then
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - b_absorb_acoustic_top(i,ib_top(ispecabs),NSTEP-it+1)
+ endif
+
+ if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
+ ! saves contribution
+ b_absorb_acoustic_top(i,ib_top(ispecabs),it) = &
+ potential_dot_acoustic(iglob)*weight/cpl/rhol
+ endif
+ enddo
+ endif ! end of top absorbing boundary
+
+ endif ! acoustic ispec
+ enddo
+ endif ! end of absorbing boundaries
+
+ end subroutine compute_forces_acoustic_2
+
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_poro_fluid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_poro_fluid.f90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_poro_fluid.f90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -58,7 +58,7 @@
b_viscodampx,b_viscodampz,&
ibegin_bottom_poro,iend_bottom_poro,ibegin_top_poro,iend_top_poro, &
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
- C_k,M_k,NSOURCE,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
+ C_k,M_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
b_absorb_poro_w_left,b_absorb_poro_w_right,b_absorb_poro_w_bottom,b_absorb_poro_w_top,&
nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top,f0,freq0,Q0)
@@ -67,16 +67,16 @@
implicit none
include "constants.h"
- integer :: NSOURCE, i_source
- integer, dimension(NSOURCE) ::ispec_selected_source,source_type,is_proc_source
+ integer :: NSOURCES, i_source
+ integer, dimension(NSOURCES) ::ispec_selected_source,source_type,is_proc_source
integer :: npoin,nspec,nelemabs,numat,it,NSTEP
integer :: nrec,SIMULATION_TYPE,myrank
integer, dimension(nrec) :: ispec_selected_rec,which_proc_receiver
integer :: nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax
- integer, dimension(nspec_xmin) :: ib_left
- integer, dimension(nspec_xmax) :: ib_right
- integer, dimension(nspec_zmin) :: ib_bottom
- integer, dimension(nspec_zmax) :: ib_top
+ integer, dimension(nelemabs) :: ib_left
+ integer, dimension(nelemabs) :: ib_right
+ integer, dimension(nelemabs) :: ib_bottom
+ integer, dimension(nelemabs) :: ib_top
logical :: anyabs,initialfield,TURN_ATTENUATION_ON
logical :: SAVE_FORWARD
@@ -99,8 +99,8 @@
double precision, dimension(numat) :: porosity,tortuosity
double precision, dimension(4,3,numat) :: poroelastcoef
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
- real(kind=CUSTOM_REAL), dimension(NSOURCE,NSTEP) :: source_time_function
- real(kind=CUSTOM_REAL), dimension(NSOURCE,NDIM,NGLLX,NGLLZ) :: sourcearray
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NSTEP) :: source_time_function
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLZ) :: sourcearray
real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays
real(kind=CUSTOM_REAL), dimension(npoin) :: C_k,M_k
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLZ,nspec_xmin,NSTEP) :: b_absorb_poro_w_left
@@ -797,7 +797,7 @@
! --- add the source
if(.not. initialfield) then
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
! if this processor carries the source and the source element is poroelastic
if (is_proc_source(i_source) == 1 .and. poroelastic(ispec_selected_source(i_source))) then
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_poro_solid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_poro_solid.f90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_poro_solid.f90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -58,7 +58,7 @@
b_viscodampx,b_viscodampz,&
ibegin_bottom_poro,iend_bottom_poro,ibegin_top_poro,iend_top_poro, &
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
- mufr_k,B_k,NSOURCE,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
+ mufr_k,B_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
b_absorb_poro_s_left,b_absorb_poro_s_right,b_absorb_poro_s_bottom,b_absorb_poro_s_top,&
nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top,f0,freq0,Q0)
@@ -67,16 +67,16 @@
implicit none
include "constants.h"
- integer :: NSOURCE, i_source
- integer, dimension(NSOURCE) :: ispec_selected_source,source_type,is_proc_source
+ integer :: NSOURCES, i_source
+ integer, dimension(NSOURCES) :: ispec_selected_source,source_type,is_proc_source
integer :: npoin,nspec,nelemabs,numat,it,NSTEP
integer :: nrec,SIMULATION_TYPE,myrank
integer, dimension(nrec) :: ispec_selected_rec,which_proc_receiver
integer :: nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax
- integer, dimension(nspec_xmin) :: ib_left
- integer, dimension(nspec_xmax) :: ib_right
- integer, dimension(nspec_zmin) :: ib_bottom
- integer, dimension(nspec_zmax) :: ib_top
+ integer, dimension(nelemabs) :: ib_left
+ integer, dimension(nelemabs) :: ib_right
+ integer, dimension(nelemabs) :: ib_bottom
+ integer, dimension(nelemabs) :: ib_top
logical :: anyabs,initialfield,TURN_ATTENUATION_ON
logical :: SAVE_FORWARD
@@ -99,8 +99,8 @@
double precision, dimension(numat) :: porosity,tortuosity
double precision, dimension(4,3,numat) :: poroelastcoef
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
- real(kind=CUSTOM_REAL), dimension(NSOURCE,NSTEP) :: source_time_function
- real(kind=CUSTOM_REAL), dimension(NSOURCE,NDIM,NGLLX,NGLLZ) :: sourcearray
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NSTEP) :: source_time_function
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLZ) :: sourcearray
real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays
real(kind=CUSTOM_REAL), dimension(npoin) :: mufr_k,B_k
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLZ,nspec_xmin,NSTEP) :: b_absorb_poro_s_left
@@ -820,7 +820,7 @@
! --- add the source
if(.not. initialfield) then
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
! if this processor carries the source and the source element is poroelastic
if (is_proc_source(i_source) == 1 .and. poroelastic(ispec_selected_source(i_source))) then
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_viscoelastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_viscoelastic.f90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_viscoelastic.f90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -57,7 +57,7 @@
deltat,coord,add_Bielak_conditions, &
x0_source, z0_source, A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset,f0, &
v0x_left,v0z_left,v0x_right,v0z_right,v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right,t0x_bot,t0z_bot,&
- nleft,nright,nbot,over_critical_angle,NSOURCE,nrec,SIMULATION_TYPE,SAVE_FORWARD,b_absorb_elastic_left,&
+ nleft,nright,nbot,over_critical_angle,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,b_absorb_elastic_left,&
b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top,nspec_xmin,nspec_xmax,&
nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top,mu_k,kappa_k)
@@ -68,24 +68,24 @@
include "constants.h"
logical :: p_sv
- integer :: NSOURCE, i_source
+ integer :: NSOURCES, i_source
integer :: npoin,nspec,myrank,nelemabs,numat,it,NSTEP
- integer, dimension(NSOURCE) :: ispec_selected_source,is_proc_source,source_type
+ integer, dimension(NSOURCES) :: ispec_selected_source,is_proc_source,source_type
integer :: nrec,SIMULATION_TYPE
integer, dimension(nrec) :: ispec_selected_rec,which_proc_receiver
integer :: nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax
- integer, dimension(nspec_xmin) :: ib_left
- integer, dimension(nspec_xmax) :: ib_right
- integer, dimension(nspec_zmin) :: ib_bottom
- integer, dimension(nspec_zmax) :: ib_top
+ integer, dimension(nelemabs) :: ib_left
+ integer, dimension(nelemabs) :: ib_right
+ integer, dimension(nelemabs) :: ib_bottom
+ integer, dimension(nelemabs) :: ib_top
logical :: anyabs,assign_external_model,initialfield,TURN_ATTENUATION_ON,add_Bielak_conditions
logical :: SAVE_FORWARD
double precision :: deltatcube,deltatfourth,twelvedeltat,fourdeltatsquare
- double precision, dimension(NSOURCE) :: angleforce
+ double precision, dimension(NSOURCES) :: angleforce
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
integer, dimension(nspec) :: kmato
@@ -102,8 +102,8 @@
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,vsext,rhoext
double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c15ext,c13ext,c33ext,c35ext,c55ext
- real(kind=CUSTOM_REAL), dimension(NSOURCE,NSTEP) :: source_time_function
- real(kind=CUSTOM_REAL), dimension(NSOURCE,NDIM,NGLLX,NGLLZ) :: sourcearray
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NSTEP) :: source_time_function
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLZ) :: sourcearray
real(kind=CUSTOM_REAL), dimension(3,npoin) :: b_accel_elastic,b_displ_elastic
real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays
@@ -847,7 +847,7 @@
! --- add the source if it is a moment tensor
if(.not. initialfield) then
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
! if this processor carries the source and the source element is elastic
if (is_proc_source(i_source) == 1 .and. elastic(ispec_selected_source(i_source))) then
@@ -882,7 +882,7 @@
endif !if(source_type(i_source) == 2)
endif ! if this processor carries the source and the source element is elastic
- enddo ! do i_source=1,NSOURCE
+ enddo ! do i_source=1,NSOURCES
if(SIMULATION_TYPE == 2) then ! adjoint wavefield
Modified: seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -43,8 +43,8 @@
!========================================================================
subroutine enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
- potential_acoustic,acoustic_surface, &
- ibool,nelem_acoustic_surface,npoin,nspec)
+ potential_acoustic,acoustic_surface, &
+ ibool,nelem_acoustic_surface,npoin,nspec)
! free surface for an acoustic medium
! if acoustic, the free surface condition is a Dirichlet condition for the potential,
@@ -56,11 +56,12 @@
integer :: nelem_acoustic_surface,npoin,nspec
- integer, dimension(5,max(1,nelem_acoustic_surface)) :: acoustic_surface
+ integer, dimension(5,nelem_acoustic_surface) :: acoustic_surface
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
- real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin) :: &
+ potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
!---
!--- local variables
@@ -70,17 +71,17 @@
do ispec_acoustic_surface = 1, nelem_acoustic_surface
- ispec = acoustic_surface(1,ispec_acoustic_surface)
- do j = acoustic_surface(4,ispec_acoustic_surface), acoustic_surface(5,ispec_acoustic_surface)
- do i = acoustic_surface(2,ispec_acoustic_surface), acoustic_surface(3,ispec_acoustic_surface)
- iglob = ibool(i,j,ispec)
- potential_acoustic(iglob) = ZERO
- potential_dot_acoustic(iglob) = ZERO
- potential_dot_dot_acoustic(iglob) = ZERO
+ ispec = acoustic_surface(1,ispec_acoustic_surface)
+
+ do j = acoustic_surface(4,ispec_acoustic_surface), acoustic_surface(5,ispec_acoustic_surface)
+ do i = acoustic_surface(2,ispec_acoustic_surface), acoustic_surface(3,ispec_acoustic_surface)
+ iglob = ibool(i,j,ispec)
+ potential_acoustic(iglob) = ZERO
+ potential_dot_acoustic(iglob) = ZERO
+ potential_dot_dot_acoustic(iglob) = ZERO
+ enddo
+ enddo
- enddo
- enddo
-
enddo
end subroutine enforce_acoustic_free_surface
Modified: seismo/2D/SPECFEM2D/trunk/locate_source_force.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -114,7 +114,8 @@
do i = 2,NGLLX-1
iglob = ibool(i,j,ispec)
- dist = sqrt((x_source-dble(coord(1,iglob)))**2 + (z_source-dble(coord(2,iglob)))**2)
+ dist = sqrt((x_source-dble(coord(1,iglob)))**2 &
+ + (z_source-dble(coord(2,iglob)))**2)
! keep this point if it is closer to the source
if(dist < distmin) then
@@ -133,7 +134,8 @@
#ifdef USE_MPI
! global minimum distance computed over all processes
- call MPI_ALLREDUCE (distmin, dist_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ierror)
+ call MPI_ALLREDUCE (distmin, dist_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MIN, MPI_COMM_WORLD, ierror)
#else
dist_glob = distmin
@@ -144,10 +146,13 @@
if ( abs(dist_glob - distmin) < TINYVAL ) is_proc_source = 1
#ifdef USE_MPI
- ! determining the number of processes that contain the source (useful when the source is located on an interface)
- call MPI_ALLREDUCE (is_proc_source, nb_proc_source, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierror)
+ ! determining the number of processes that contain the source
+ ! (useful when the source is located on an interface)
+ call MPI_ALLREDUCE (is_proc_source, nb_proc_source, 1, MPI_INTEGER, &
+ MPI_SUM, MPI_COMM_WORLD, ierror)
#else
+
nb_proc_source = is_proc_source
#endif
@@ -157,7 +162,8 @@
! when several processes contain the source, we elect one of them (minimum rank).
if ( nb_proc_source > 1 ) then
- call MPI_ALLGATHER(is_proc_source, 1, MPI_INTEGER, allgather_is_proc_source(1), 1, MPI_INTEGER, MPI_COMM_WORLD, ierror)
+ call MPI_ALLGATHER(is_proc_source, 1, MPI_INTEGER, allgather_is_proc_source(1), &
+ 1, MPI_INTEGER, MPI_COMM_WORLD, ierror)
locate_is_proc_source = maxloc(allgather_is_proc_source) - 1
if ( myrank /= locate_is_proc_source(1) ) then
@@ -181,37 +187,39 @@
do iter_loop = 1,NUM_ITER
! recompute jacobian for the new point
- call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian,coorg,knods,ispec_selected_source,ngnod,nspec,npgeo, &
- .true.)
+ call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian, &
+ coorg,knods,ispec_selected_source,ngnod,nspec,npgeo, &
+ .true.)
! compute distance to target location
- dx = - (x - x_source)
- dz = - (z - z_source)
+ dx = - (x - x_source)
+ dz = - (z - z_source)
! compute increments
- dxi = xix*dx + xiz*dz
- dgamma = gammax*dx + gammaz*dz
+ dxi = xix*dx + xiz*dz
+ dgamma = gammax*dx + gammaz*dz
! update values
- xi = xi + dxi
- gamma = gamma + dgamma
+ xi = xi + dxi
+ gamma = gamma + dgamma
! impose that we stay in that element
! (useful if user gives a source outside the mesh for instance)
! we can go slightly outside the [1,1] segment since with finite elements
! the polynomial solution is defined everywhere
! this can be useful for convergence of itertive scheme with distorted elements
- if (xi > 1.10d0) xi = 1.10d0
- if (xi < -1.10d0) xi = -1.10d0
- if (gamma > 1.10d0) gamma = 1.10d0
- if (gamma < -1.10d0) gamma = -1.10d0
+ if (xi > 1.10d0) xi = 1.10d0
+ if (xi < -1.10d0) xi = -1.10d0
+ if (gamma > 1.10d0) gamma = 1.10d0
+ if (gamma < -1.10d0) gamma = -1.10d0
! end of non linear iterations
enddo
! compute final coordinates of point found
- call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian,coorg,knods,ispec_selected_source,ngnod,nspec,npgeo, &
- .true.)
+ call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian, &
+ coorg,knods,ispec_selected_source,ngnod,nspec,npgeo, &
+ .true.)
! store xi,gamma of point found
xi_source = xi
@@ -229,6 +237,9 @@
write(IOUT,*) ' original x: ',sngl(x_source)
write(IOUT,*) ' original z: ',sngl(z_source)
write(IOUT,*) 'closest estimate found: ',sngl(final_distance),' m away'
+#ifdef USE_MPI
+ write(IOUT,*) ' in rank ',myrank
+#endif
write(IOUT,*) ' in element ',ispec_selected_source
write(IOUT,*) ' at xi,gamma coordinates = ',xi_source,gamma_source
write(IOUT,*)
Modified: seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -46,7 +46,8 @@
!---- locate_source_moment_tensor finds the correct position of the moment-tensor source
!----
- subroutine locate_source_moment_tensor(ibool,coord,nspec,npoin,xigll,zigll,x_source,z_source, &
+ subroutine locate_source_moment_tensor(ibool,coord,nspec,npoin, &
+ xigll,zigll,x_source,z_source, &
ispec_selected_source,is_proc_source,nb_proc_source,nproc,myrank, &
xi_source,gamma_source,coorg,knods,ngnod,npgeo,ipass)
@@ -114,7 +115,8 @@
do i = 2,NGLLX-1
iglob = ibool(i,j,ispec)
- dist = sqrt((x_source-dble(coord(1,iglob)))**2 + (z_source-dble(coord(2,iglob)))**2)
+ dist = sqrt((x_source-dble(coord(1,iglob)))**2 &
+ + (z_source-dble(coord(2,iglob)))**2)
! keep this point if it is closer to the source
if(dist < distmin) then
@@ -132,7 +134,8 @@
#ifdef USE_MPI
! global minimum distance computed over all processes
- call MPI_ALLREDUCE (distmin, dist_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ierror)
+ call MPI_ALLREDUCE (distmin, dist_glob, 1, MPI_DOUBLE_PRECISION, &
+ MPI_MIN, MPI_COMM_WORLD, ierror)
#else
dist_glob = distmin
@@ -143,8 +146,10 @@
if ( dist_glob == distmin ) is_proc_source = 1
#ifdef USE_MPI
- ! determining the number of processes that contain the source (useful when the source is located on an interface)
- call MPI_ALLREDUCE (is_proc_source, nb_proc_source, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierror)
+ ! determining the number of processes that contain the source
+ ! (useful when the source is located on an interface)
+ call MPI_ALLREDUCE (is_proc_source, nb_proc_source, 1, MPI_INTEGER, &
+ MPI_SUM, MPI_COMM_WORLD, ierror)
#else
nb_proc_source = is_proc_source
@@ -156,7 +161,8 @@
! when several processes contain the source, we elect one of them (minimum rank).
if ( nb_proc_source > 1 ) then
- call MPI_ALLGATHER(is_proc_source, 1, MPI_INTEGER, allgather_is_proc_source(1), 1, MPI_INTEGER, MPI_COMM_WORLD, ierror)
+ call MPI_ALLGATHER(is_proc_source, 1, MPI_INTEGER, allgather_is_proc_source(1), &
+ 1, MPI_INTEGER, MPI_COMM_WORLD, ierror)
locate_is_proc_source = maxloc(allgather_is_proc_source) - 1
if ( myrank /= locate_is_proc_source(1) ) then
@@ -180,8 +186,9 @@
do iter_loop = 1,NUM_ITER
! recompute jacobian for the new point
- call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian,coorg,knods,ispec_selected_source,ngnod,nspec,npgeo, &
- .true.)
+ call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian, &
+ coorg,knods,ispec_selected_source,ngnod,nspec,npgeo, &
+ .true.)
! compute distance to target location
dx = - (x - x_source)
@@ -209,8 +216,9 @@
enddo
! compute final coordinates of point found
- call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian,coorg,knods,ispec_selected_source,ngnod,nspec,npgeo, &
- .true.)
+ call recompute_jacobian(xi,gamma,x,z,xix,xiz,gammax,gammaz,jacobian, &
+ coorg,knods,ispec_selected_source,ngnod,nspec,npgeo, &
+ .true.)
! store xi,gamma of point found
xi_source = xi
@@ -228,6 +236,9 @@
write(IOUT,*) ' original x: ',sngl(x_source)
write(IOUT,*) ' original z: ',sngl(z_source)
write(IOUT,*) 'closest estimate found: ',sngl(final_distance),' m away'
+#ifdef USE_MPI
+ write(IOUT,*) ' in rank ',myrank
+#endif
write(IOUT,*) ' in element ',ispec_selected_source
write(IOUT,*) ' at xi,gamma coordinates = ',xi_source,gamma_source
write(IOUT,*)
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -418,7 +418,7 @@
print *,'Parameter file successfully read... '
! reads in source descriptions
- call read_source_file(NSOURCE,deltat,f0_attenuation)
+ call read_source_file(NSOURCES)
! reads in tangential detection
if (force_normal_to_surface .or. rec_normal_to_surface) then
@@ -505,7 +505,7 @@
! check if we are in the last layer, which contains topography,
! and modify the position of the source accordingly if it is located exactly at the surface
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
if(source_surf(i_source) .and. ilayer == number_of_layers) &
zs(i_source) = value_spline(xs(i_source),xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
enddo
@@ -911,7 +911,7 @@
nnodes_tangential_curve,nodes_tangential_curve)
! print position of the source
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
print *
print *,'Position (x,z) of the source = ',xs(i_source),zs(i_source)
print *
Modified: seismo/2D/SPECFEM2D/trunk/plotpost.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/plotpost.F90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/plotpost.F90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -46,7 +46,7 @@
xinterp,zinterp,shapeint,Uxinterp,Uzinterp,flagrange,density,porosity,tortuosity,&
poroelastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,nelem_acoustic_surface, acoustic_edges, &
- simulation_title,npoin,npgeo,vpmin,vpmax,nrec,NSOURCE, &
+ simulation_title,npoin,npgeo,vpmin,vpmax,nrec,NSOURCES, &
colors,numbers,subsamp,imagetype,interpol,meshvect,modelvect, &
boundvect,assign_external_model,cutsnaps,sizemax_arrows,nelemabs,numat,pointsdisp, &
nspec,ngnod,coupled_acoustic_elastic,coupled_acoustic_poro,coupled_elastic_poro, &
@@ -92,7 +92,7 @@
double precision, dimension(NUM_COLORS) :: red,green,blue
integer it,nrec,nelemabs,numat,pointsdisp,pointsdisp_loop,nspec
- integer i,npoin,npgeo,ngnod,NSOURCE
+ integer i,npoin,npgeo,ngnod,NSOURCES
integer kmato(nspec),knods(ngnod,nspec)
integer ibool(NGLLX,NGLLZ,nspec)
@@ -105,7 +105,7 @@
double precision density(2,numat),poroelastcoef(4,3,numat),porosity(numat),tortuosity(numat)
double precision dt,timeval
- double precision, dimension(NSOURCE) :: x_source,z_source
+ double precision, dimension(NSOURCES) :: x_source,z_source
double precision displ(3,npoin),coord(NDIM,npoin)
double precision vpext(NGLLX,NGLLZ,nspec)
@@ -3011,9 +3011,9 @@
!
!---- write position of the source
!
- do i=1,NSOURCE
+ do i=1,NSOURCES
if(i == 1) write(24,*) '% beginning of source line'
- if(i == NSOURCE) write(24,*) '% end of source line'
+ if(i == NSOURCES) write(24,*) '% end of source line'
xw = x_source(i)
zw = z_source(i)
xw = (xw-xmin)*ratio_page + orig_x
Modified: seismo/2D/SPECFEM2D/trunk/prepare_initialfield.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/prepare_initialfield.F90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/prepare_initialfield.F90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -45,7 +45,7 @@
subroutine prepare_initialfield(myrank,any_acoustic,any_poroelastic,over_critical_angle, &
- NSOURCE,source_type,angleforce,x_source,z_source,f0, &
+ NSOURCES,source_type,angleforce,x_source,z_source,f0, &
npoin,numat,poroelastcoef,density,coord, &
angleforce_refl,c_inc,c_refl,cploc,csloc,time_offset, &
A_plane, B_plane, C_plane, &
@@ -60,9 +60,9 @@
integer :: myrank
logical :: any_acoustic,any_poroelastic
- integer :: NSOURCE
- integer, dimension(NSOURCE) :: source_type
- double precision, dimension(NSOURCE) :: angleforce,x_source,z_source,f0
+ integer :: NSOURCES
+ integer, dimension(NSOURCES) :: source_type
+ double precision, dimension(NSOURCES) :: angleforce,x_source,z_source,f0
integer :: npoin,numat
double precision, dimension(4,3,numat) :: poroelastcoef
@@ -117,7 +117,7 @@
write(IOUT,*)
! only implemented for one source
- if(NSOURCE > 1) call exit_MPI('calculation of the initial wave is only implemented for one source')
+ if(NSOURCES > 1) call exit_MPI('calculation of the initial wave is only implemented for one source')
if (source_type(1) == 1) then
write(IOUT,*) 'initial P wave of', angleforce(1)*180.d0/pi, 'degrees introduced.'
else if (source_type(1) == 2) then
@@ -317,7 +317,7 @@
subroutine prepare_initialfield_paco(myrank,nelemabs,left_bound,right_bound,bot_bound, &
numabs,codeabs,ibool,nspec, &
- source_type,NSOURCE,c_inc,c_refl, &
+ source_type,NSOURCES,c_inc,c_refl, &
count_bottom,count_left,count_right)
implicit none
@@ -335,8 +335,8 @@
integer :: nspec
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
- integer :: NSOURCE
- integer :: source_type(NSOURCE)
+ integer :: NSOURCES
+ integer :: source_type(NSOURCES)
double precision :: c_inc,c_refl
Modified: seismo/2D/SPECFEM2D/trunk/prepare_source_time_function.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/prepare_source_time_function.f90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/prepare_source_time_function.f90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -44,7 +44,7 @@
!========================================================================
- subroutine prepare_source_time_function(myrank,NSTEP,NSOURCE,source_time_function, &
+ subroutine prepare_source_time_function(myrank,NSTEP,NSOURCES,source_time_function, &
time_function_type,f0,tshift_src,factor,aval, &
t0,nb_proc_source,deltat)
@@ -55,19 +55,19 @@
integer :: myrank,NSTEP
- integer :: NSOURCE
- integer, dimension(NSOURCE) :: time_function_type
- double precision, dimension(NSOURCE) :: f0,tshift_src,factor
- double precision, dimension(NSOURCE) :: aval
+ integer :: NSOURCES
+ integer, dimension(NSOURCES) :: time_function_type
+ double precision, dimension(NSOURCES) :: f0,tshift_src,factor
+ double precision, dimension(NSOURCES) :: aval
double precision :: t0
- integer,dimension(NSOURCE) :: nb_proc_source
+ integer,dimension(NSOURCES) :: nb_proc_source
double precision :: deltat
- real(kind=CUSTOM_REAL),dimension(NSOURCE,NSTEP) :: source_time_function
+ real(kind=CUSTOM_REAL),dimension(NSOURCES,NSTEP) :: source_time_function
! local parameters
double precision :: stf_used,time
- double precision, dimension(NSOURCE) :: hdur,hdur_gauss
+ double precision, dimension(NSOURCES) :: hdur,hdur_gauss
double precision, external :: netlib_specfun_erf
integer :: it,i_source
@@ -81,42 +81,45 @@
endif
! ! loop on all the sources
- ! do i_source=1,NSOURCE
+ ! do i_source=1,NSOURCES
! loop on all the time steps
do it = 1,NSTEP
+ ! note: t0 is the simulation start time, tshift_src is the time shift of the source
+ ! relative to this start time
+
! compute current time
time = (it-1)*deltat
stf_used = 0.d0
! loop on all the sources
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
if( time_function_type(i_source) == 1 ) then
! Ricker (second derivative of a Gaussian) source time function
source_time_function(i_source,it) = - factor(i_source) * &
- (ONE-TWO*aval(i_source)*(time-tshift_src(i_source))**2) * &
- exp(-aval(i_source)*(time-tshift_src(i_source))**2)
+ (ONE-TWO*aval(i_source)*(time-t0-tshift_src(i_source))**2) * &
+ exp(-aval(i_source)*(time-t0-tshift_src(i_source))**2)
! source_time_function(i_source,it) = - factor(i_source) * &
! TWO*aval(i_source)*sqrt(aval(i_source))*&
- ! (time-tshift_src(i_source))/pi * exp(-aval(i_source)*(time-tshift_src(i_source))**2)
+ ! (time-t0-tshift_src(i_source))/pi * exp(-aval(i_source)*(time-t0-tshift_src(i_source))**2)
else if( time_function_type(i_source) == 2 ) then
! first derivative of a Gaussian source time function
source_time_function(i_source,it) = - factor(i_source) * &
- TWO*aval(i_source)*(time-tshift_src(i_source)) * &
- exp(-aval(i_source)*(time-tshift_src(i_source))**2)
+ TWO*aval(i_source)*(time-t0-tshift_src(i_source)) * &
+ exp(-aval(i_source)*(time-t0-tshift_src(i_source))**2)
else if(time_function_type(i_source) == 3 .or. time_function_type(i_source) == 4) then
! Gaussian or Dirac (we use a very thin Gaussian instead) source time function
source_time_function(i_source,it) = factor(i_source) * &
- exp(-aval(i_source)*(time-tshift_src(i_source))**2)
+ exp(-aval(i_source)*(time-t0-tshift_src(i_source))**2)
else if(time_function_type(i_source) == 5) then
@@ -124,7 +127,7 @@
hdur(i_source) = 1.d0 / f0(i_source)
hdur_gauss(i_source) = hdur(i_source) * 5.d0 / 3.d0
source_time_function(i_source,it) = factor(i_source) * 0.5d0*(1.0d0 + &
- netlib_specfun_erf(SOURCE_DECAY_MIMIC_TRIANGLE*(time-tshift_src(i_source))/hdur_gauss(i_source)))
+ netlib_specfun_erf(SOURCE_DECAY_MIMIC_TRIANGLE*(time-t0-tshift_src(i_source))/hdur_gauss(i_source)))
else
call exit_MPI('unknown source time function')
@@ -135,7 +138,7 @@
enddo
! output relative time in third column, in case user wants to check it as well
- ! if (myrank == 0 .and. i_source==1 ) write(55,*) sngl(time-tshift_src(1)),real(source_time_function(1,it),4),sngl(time)
+ ! if (myrank == 0 .and. i_source==1 ) write(55,*) sngl(time-t0-tshift_src(1)),real(source_time_function(1,it),4),sngl(time)
if (myrank == 0) then
! note: earliest start time of the simulation is: (it-1)*deltat - t0
write(55,*) sngl(time-t0),sngl(stf_used),sngl(time)
@@ -150,10 +153,8 @@
! than one if the nearest point is on the interface between several partitions with an explosive source.
! since source contribution is linear, the source_time_function is cut down by that number (it would have been similar
! if we just had elected one of those processes).
- do i_source=1,NSOURCE
-
+ do i_source=1,NSOURCES
source_time_function(i_source,:) = source_time_function(i_source,:) / nb_proc_source(i_source)
-
enddo
end subroutine prepare_source_time_function
Modified: seismo/2D/SPECFEM2D/trunk/read_databases.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/read_databases.f90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/read_databases.f90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -53,7 +53,7 @@
seismotype,imagetype,assign_external_model,READ_EXTERNAL_SEP_FILE, &
outputgrid,OUTPUT_ENERGY,TURN_ATTENUATION_ON, &
TURN_VISCATTENUATION_ON,Q0,freq0,p_sv, &
- NSTEP,deltat,NTSTEP_BETWEEN_OUTPUT_SEISMO,NSOURCE)
+ NSTEP,deltat,NTSTEP_BETWEEN_OUTPUT_SEISMO,NSOURCES)
! starts reading in parameters from input Database file
@@ -75,7 +75,7 @@
double precision :: Q0,freq0
double precision :: deltat
- integer :: NSTEP,NSOURCE
+ integer :: NSTEP,NSOURCES
integer :: NTSTEP_BETWEEN_OUTPUT_INFO,NTSTEP_BETWEEN_OUTPUT_SEISMO
! local parameters
@@ -189,7 +189,7 @@
!---- read source information
read(IIN,"(a80)") datlin
- read(IIN,*) NSOURCE
+ read(IIN,*) NSOURCES
! output formats
230 format('./OUTPUT_FILES/Database',i5.5)
@@ -236,7 +236,7 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine read_databases_sources(NSOURCE,source_type,time_function_type, &
+ subroutine read_databases_sources(NSOURCES,source_type,time_function_type, &
x_source,z_source,Mxx,Mzz,Mxz,f0,tshift_src,factor,angleforce)
! reads source parameters
@@ -244,17 +244,30 @@
implicit none
include "constants.h"
- integer :: NSOURCE
- integer, dimension(NSOURCE) :: source_type,time_function_type
- double precision, dimension(NSOURCE) :: x_source,z_source, &
+ integer :: NSOURCES
+ integer, dimension(NSOURCES) :: source_type,time_function_type
+ double precision, dimension(NSOURCES) :: x_source,z_source, &
Mxx,Mzz,Mxz,f0,tshift_src,factor,angleforce
! local parameters
integer :: i_source
character(len=80) :: datlin
+ ! initializes
+ source_type(:) = 0
+ time_function_type(:) = 0
+ x_source(:) = 0.d0
+ z_source(:) = 0.d0
+ Mxx(:) = 0.d0
+ Mzz(:) = 0.d0
+ Mxz(:) = 0.d0
+ f0(:) = 0.d0
+ tshift_src(:) = 0.d0
+ factor(:) = 0.d0
+ angleforce(:) = 0.d0
+
! reads in source info from Database file
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
read(IIN,"(a80)") datlin
read(IIN,*) source_type(i_source),time_function_type(i_source), &
x_source(i_source),z_source(i_source),f0(i_source),tshift_src(i_source), &
@@ -315,10 +328,14 @@
double precision, dimension(:), allocatable :: coorgread
character(len=80) :: datlin
+ ! initializes
+ coorg(:,:) = 0.d0
+
! reads the spectral macrobloc nodal coordinates
- ipoin = 0
read(IIN,"(a80)") datlin
+ ! reads in values
+ ipoin = 0
allocate(coorgread(NDIM))
do ip = 1,npgeo
! reads coordinates
@@ -385,10 +402,16 @@
integer, dimension(:), allocatable :: knods_read
character(len=80) :: datlin
+ ! initializes
+ kmato(:) = 0
+ knods(:,:) = 0
+
! reads spectral macrobloc data
- n = 0
read(IIN,"(a80)") datlin
+
+ ! reads in values
allocate(knods_read(ngnod))
+ n = 0
do ispec = 1,nspec
read(IIN,*) n,kmato_read,(knods_read(k), k=1,ngnod)
if(ipass == 1) then
@@ -450,6 +473,11 @@
! local parameters
integer :: num_interface,ie,my_interfaces_read
+ ! initializes
+ my_neighbours(:) = 0
+ my_nelmnts_neighbours(:) = 0
+ my_interfaces(:,:,:) = 0
+
! reads in interfaces
do num_interface = 1, ninterface
read(IIN,*) my_neighbours(num_interface), my_nelmnts_neighbours(num_interface)
@@ -498,9 +526,21 @@
logical :: codeabsread(4)
character(len=80) :: datlin
+ ! initializes
+ codeabs(:,:) = .false.
+ ibegin_bottom(:) = 0
+ iend_bottom(:) = 0
+ ibegin_top(:) = 0
+ iend_top(:) = 0
+ jbegin_left(:) = 0
+ jend_left(:) = 0
+ jbegin_right(:) = 0
+ jend_right(:) = 0
+
! reads in absorbing edges
read(IIN,"(a80)") datlin
-
+
+ ! reads in values
if( anyabs ) then
do inum = 1,nelemabs
read(IIN,*) numabsread,codeabsread(1),codeabsread(2),codeabsread(3),&
@@ -557,6 +597,9 @@
integer :: inum,acoustic_edges_read
character(len=80) :: datlin
+ ! initializes
+ acoustic_edges(:,:) = 0
+
! reads in any possible free surface edges
read(IIN,"(a80)") datlin
@@ -621,6 +664,14 @@
solid_poro_poro_ispec_read,solid_poro_elastic_ispec_read
character(len=80) :: datlin
+ ! initializes
+ fluid_solid_acoustic_ispec(:) = 0
+ fluid_solid_elastic_ispec(:) = 0
+ fluid_poro_acoustic_ispec(:) = 0
+ fluid_poro_poroelastic_ispec(:) = 0
+ solid_poro_elastic_ispec(:) = 0
+ solid_poro_poroelastic_ispec(:) = 0
+
! reads acoustic elastic coupled edges
read(IIN,"(a80)") datlin
@@ -704,6 +755,9 @@
integer :: i
character(len=80) :: datlin
+ ! initializes
+ nodes_tangential_curve(:,:) = 0.d0
+
! reads tangential detection curve
read(IIN,"(a80)") datlin
read(IIN,*) force_normal_to_surface,rec_normal_to_surface
Modified: seismo/2D/SPECFEM2D/trunk/read_parameter_file.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/read_parameter_file.F90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/read_parameter_file.F90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -78,7 +78,7 @@
integer :: nt
double precision :: deltat
- integer :: NSOURCE
+ integer :: NSOURCES
logical :: force_normal_to_surface
! variables used for attenuation
@@ -180,7 +180,7 @@
call read_value_double_precision(IIN,IGNORE_JUNK,deltat)
! read source infos
- call read_value_integer(IIN,IGNORE_JUNK,NSOURCE)
+ call read_value_integer(IIN,IGNORE_JUNK,NSOURCES)
call read_value_logical(IIN,IGNORE_JUNK,force_normal_to_surface)
! read constants for attenuation
Modified: seismo/2D/SPECFEM2D/trunk/read_source_file.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/read_source_file.f90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/read_source_file.f90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -54,34 +54,33 @@
contains
- subroutine read_source_file(NSOURCE,deltat,f0_attenuation)
+ subroutine read_source_file(NSOURCES)
! reads in source file DATA/SOURCE
implicit none
include "constants.h"
- integer :: NSOURCE
- double precision :: deltat,f0_attenuation
+ integer :: NSOURCES
! local parameters
- integer :: ios,icounter,i_source,nsources
+ integer :: ios,icounter,i_source,num_sources
character(len=150) dummystring
integer, parameter :: IIN_SOURCE = 22
! allocates memory arrays
- allocate(source_surf(NSOURCE))
- allocate(xs(NSOURCE))
- allocate(zs(NSOURCE))
- allocate(source_type(NSOURCE))
- allocate(time_function_type(NSOURCE))
- allocate(f0(NSOURCE))
- allocate(tshift_src(NSOURCE))
- allocate(angleforce(NSOURCE))
- allocate(Mxx(NSOURCE))
- allocate(Mxz(NSOURCE))
- allocate(Mzz(NSOURCE))
- allocate(factor(NSOURCE))
+ allocate(source_surf(NSOURCES))
+ allocate(xs(NSOURCES))
+ allocate(zs(NSOURCES))
+ allocate(source_type(NSOURCES))
+ allocate(time_function_type(NSOURCES))
+ allocate(f0(NSOURCES))
+ allocate(tshift_src(NSOURCES))
+ allocate(angleforce(NSOURCES))
+ allocate(Mxx(NSOURCES))
+ allocate(Mxz(NSOURCES))
+ allocate(Mzz(NSOURCES))
+ allocate(factor(NSOURCES))
! counts lines
open(unit=IIN_SOURCE,file='DATA/SOURCE',iostat=ios,status='old',action='read')
@@ -94,18 +93,20 @@
enddo
close(IIN_SOURCE)
+ ! checks counter
if(mod(icounter,NLINES_PER_SOURCE) /= 0) &
stop 'total number of lines in SOURCE file should be a multiple of NLINES_PER_SOURCE'
- nsources = icounter / NLINES_PER_SOURCE
+ ! total number of sources
+ num_sources = icounter / NLINES_PER_SOURCE
- if(nsources < 1) stop 'need at least one source in SOURCE file'
- if(nsources /= NSOURCE) &
+ if(num_sources < 1) stop 'need at least one source in SOURCE file'
+ if(num_sources /= NSOURCES) &
stop 'total number of sources read is different than declared in Par_file'
! reads in source parameters
open(unit=IIN_SOURCE,file='DATA/SOURCE',status='old',action='read')
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
call read_value_logical(IIN_SOURCE,IGNORE_JUNK,source_surf(i_source))
call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,xs(i_source))
call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,zs(i_source))
@@ -119,23 +120,9 @@
call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,Mxz(i_source))
call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,factor(i_source))
- ! note: this is slightly different than in specfem2D.f90,
- ! tshift_src will be set outside of this next if statement, i.e. it will be set for all sources
- ! regardless of their type (it just makes a distinction between type 5 sources and the rest)
+ ! note: we will further process source info in solver,
+ ! here we just read in the given specifics and show them
- ! if Dirac source time function, use a very thin Gaussian instead
- ! if Heaviside source time function, use a very thin error function instead
- if(time_function_type(i_source) == 4 .or. time_function_type(i_source) == 5) then
- f0(i_source) = 1.d0 / (10.d0 * deltat)
- endif
-
- ! time delay of the source in seconds, use a 20 % security margin (use 2 / f0 if error function)
- if(time_function_type(i_source)== 5) then
- tshift_src(i_source) = 2.0d0 / f0(i_source) + tshift_src(i_source)
- else
- tshift_src(i_source) = 1.20d0 / f0(i_source) + tshift_src(i_source)
- endif
-
print *
print *,'Source', i_source
print *,'Position xs, zs = ',xs(i_source),zs(i_source)
@@ -148,15 +135,9 @@
print *,'Mxz of the source if moment tensor = ',Mxz(i_source)
print *,'Multiplying factor = ',factor(i_source)
print *
- enddo ! do i_source=1,NSOURCE
+ enddo ! do i_source=1,NSOURCES
close(IIN_SOURCE)
-
- ! if source is not a Dirac or Heavyside then f0_attenuation is f0 of the first source
- if(.not. (time_function_type(1) == 4 .or. time_function_type(1) == 5)) then
- f0_attenuation = f0(1)
- endif
-
end subroutine read_source_file
end module source_file
Modified: seismo/2D/SPECFEM2D/trunk/save_databases.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/save_databases.f90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/save_databases.f90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -138,10 +138,10 @@
write(15,*) 'nt deltat'
write(15,*) nt,deltat
- write(15,*) 'NSOURCE'
- write(15,*) NSOURCE
+ write(15,*) 'NSOURCES'
+ write(15,*) NSOURCES
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
write(15,*) 'source', i_source
write(15,*) source_type(i_source),time_function_type(i_source), &
xs(i_source),zs(i_source),f0(i_source),tshift_src(i_source), &
Modified: seismo/2D/SPECFEM2D/trunk/set_sources.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/set_sources.f90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/set_sources.f90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -44,7 +44,7 @@
!========================================================================
- subroutine set_sources(myrank,NSOURCE,source_type,time_function_type, &
+ subroutine set_sources(myrank,NSOURCES,source_type,time_function_type, &
x_source,z_source,Mxx,Mzz,Mxz,f0,tshift_src,factor,angleforce,aval, &
t0,initialfield,ipass,deltat)
@@ -54,11 +54,11 @@
include "constants.h"
integer :: myrank
- integer :: NSOURCE
- integer, dimension(NSOURCE) :: source_type,time_function_type
- double precision, dimension(NSOURCE) :: x_source,z_source, &
+ integer :: NSOURCES
+ integer, dimension(NSOURCES) :: source_type,time_function_type
+ double precision, dimension(NSOURCES) :: x_source,z_source, &
Mxx,Mzz,Mxz,f0,tshift_src,factor,angleforce
- double precision, dimension(NSOURCE) :: aval
+ double precision, dimension(NSOURCES) :: aval
double precision :: t0
double precision :: deltat
integer :: ipass
@@ -66,59 +66,74 @@
! local parameters
integer :: i_source
-
+ double precision, dimension(NSOURCES) :: t0_source,hdur
+ double precision :: min_tshift_src_original
+
! checks the input
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
! checks source type
if(.not. initialfield) then
if (source_type(i_source) == 1) then
if ( myrank == 0 .and. ipass == 1 ) then
-
+ ! user output
write(IOUT,212) x_source(i_source),z_source(i_source),f0(i_source),tshift_src(i_source), &
factor(i_source),angleforce(i_source)
-
endif
else if(source_type(i_source) == 2) then
if ( myrank == 0 .and. ipass == 1 ) then
-
+ ! user output
write(IOUT,222) x_source(i_source),z_source(i_source),f0(i_source),tshift_src(i_source), &
factor(i_source),Mxx(i_source),Mzz(i_source),Mxz(i_source)
-
endif
else
call exit_MPI('Unknown source type number !')
endif
endif
- ! note: this is slightly different than in meshfem2D.f90,
- ! t0 will only be set within this if statement, i.e. only for type 4 or 5 sources
- ! (since we set f0 to a new values for these two types of sources)
-
! if Dirac source time function, use a very thin Gaussian instead
! if Heaviside source time function, use a very thin error function instead
- if(time_function_type(i_source) == 4 .or. time_function_type(i_source) == 5) then
+ if(time_function_type(i_source) == 4 .or. time_function_type(i_source) == 5) &
f0(i_source) = 1.d0 / (10.d0 * deltat)
- ! time delay of the source in seconds, use a 20 % security margin (use 2 / f0 if error function)
- if(time_function_type(i_source)== 5) then
- tshift_src(i_source) = 2.0d0 / f0(i_source) + tshift_src(i_source)
- else
- tshift_src(i_source) = 1.20d0 / f0(i_source) + tshift_src(i_source)
- endif
+ ! checks source frequency
+ if( abs(f0(i_source)) < TINYVAL ) then
+ call exit_MPI('Error source frequency is zero')
endif
+
+ ! half-duration of source
+ hdur(i_source) = 1.d0 / f0(i_source)
+ ! sets source start times, shifted by the given (non-zero) time-shift
+ if(time_function_type(i_source)== 5) then
+ t0_source(i_source) = 2.0d0 * hdur(i_source) + tshift_src(i_source)
+ else
+ t0_source(i_source) = 1.20d0 * hdur(i_source) + tshift_src(i_source)
+ endif
+
! for the source time function
aval(i_source) = PI*PI*f0(i_source)*f0(i_source)
! convert angle from degrees to radians
angleforce(i_source) = angleforce(i_source) * PI / 180.d0
+
+ enddo ! do i_source=1,NSOURCES
- enddo ! do i_source=1,NSOURCE
-
! initializes simulation start time
- t0 = tshift_src(1)
-
+ if( NSOURCES == 1 ) then
+ ! simulation start time
+ t0 = t0_source(1)
+ ! sets source time shift relative to simulation start time
+ min_tshift_src_original = tshift_src(1)
+ tshift_src(1) = 0.d0
+ else
+ ! starts with earliest start time
+ t0 = minval( t0_source(:) )
+ ! sets source time shifts relative to simulation start time
+ min_tshift_src_original = minval( tshift_src(:) )
+ tshift_src(:) = t0_source(:) - t0
+ endif
+
! checks if user set USER_T0 to fix simulation start time
! note: USER_T0 has to be positive
if( USER_T0 > 0.d0 ) then
@@ -129,41 +144,42 @@
! notifies user
if( myrank == 0 .and. ipass == 1) then
write(IOUT,*)
- write(IOUT,*) ' USER_T0: ',USER_T0,'initial t0_start: ',t0
+ write(IOUT,*) ' using USER_T0 . . . . . . . . . = ',USER_T0
+ write(IOUT,*) ' original t0 . . . . . . . . . = ',t0
+ write(IOUT,*) ' min_tshift_src_original . . . = ',min_tshift_src_original
+ write(IOUT,*)
endif
! checks if automatically set t0 is too small
! note: times in seismograms are shifted by t0(1)
- if( t0 <= USER_T0 ) then
+ if( t0 <= USER_T0 + min_tshift_src_original ) then
+
! sets new simulation start time such that
! simulation starts at t = - t0 = - USER_T0
t0 = USER_T0
! notifies user
if( myrank == 0 .and. ipass == 1) then
- write(IOUT,*) ' fix new simulation start time. . . . . = ', - t0
- write(IOUT,*)
+ write(IOUT,*) ' fix new simulation start time . = ', - t0
endif
! loops over all sources
- do i_source=1,NSOURCE
- ! gets the given, initial time shifts
+ do i_source=1,NSOURCES
+ ! sets the given, initial time shifts
if( time_function_type(i_source) == 5 ) then
- tshift_src(i_source) = tshift_src(i_source) - 2.0d0 / f0(i_source)
+ tshift_src(i_source) = t0_source(i_source) - 2.0d0 * hdur(i_source)
else
- tshift_src(i_source) = tshift_src(i_source) - 1.20d0 / f0(i_source)
+ tshift_src(i_source) = t0_source(i_source) - 1.20d0 * hdur(i_source)
endif
-
- ! sets new t0 according to simulation start time,
- ! using absolute time shifts for each source such that
- ! a zero time shift would have the maximum gaussian at time t = (it-1)*DT - t0_start = 0
- tshift_src(i_source) = USER_T0 + tshift_src(i_source)
-
+ ! user output
if( myrank == 0 .and. ipass == 1) then
- write(IOUT,*) ' source ',i_source,'uses tshift . . . . . = ',tshift_src(i_source)
+ write(IOUT,*) ' source ',i_source,'uses tshift = ',tshift_src(i_source)
endif
-
enddo
+ ! user output
+ if( myrank == 0 .and. ipass == 1) then
+ write(IOUT,*)
+ endif
else
! start time needs to be at least t0 for numerical stability
@@ -184,6 +200,36 @@
call exit_MPI('error negative USER_T0 parameter in constants.h')
endif
+ ! checks onset times
+ if(.not. initialfield) then
+
+ ! loops over sources
+ do i_source = 1,NSOURCES
+
+ ! excludes Dirac and Heaviside sources
+ if(time_function_type(i_source) /= 4 .and. time_function_type(i_source) /= 5) then
+
+ ! user output
+ if( myrank == 0 .and. ipass == 1 ) then
+ write(IOUT,*) ' Onset time. . . . . . = ',t0+tshift_src(i_source)
+ write(IOUT,*) ' Fundamental period. . = ',1.d0/f0(i_source)
+ write(IOUT,*) ' Fundamental frequency = ',f0(i_source)
+ endif
+
+ ! checks source onset time
+ if( t0+tshift_src(i_source) <= 1.d0/f0(i_source)) then
+ call exit_MPI('Onset time too small')
+ else
+ if( myrank == 0 .and. ipass == 1 ) then
+ write(IOUT,*) ' --> onset time ok'
+ endif
+ endif
+ endif
+ enddo
+
+ endif
+
+
! output formats
212 format(//,5x,'Source Type. . . . . . . . . . . . . . = Collocated Force',/5x, &
'X-position (meters). . . . . . . . . . =',1pe20.10,/5x, &
Modified: seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -43,7 +43,7 @@
!
!========================================================================
-subroutine setup_sources_receivers(NSOURCE,initialfield,source_type,&
+ subroutine setup_sources_receivers(NSOURCES,initialfield,source_type,&
coord,ibool,npoin,nspec,nelem_acoustic_surface,acoustic_surface,elastic,poroelastic, &
x_source,z_source,ispec_selected_source,ispec_selected_rec, &
is_proc_source,nb_proc_source,ipass,&
@@ -57,7 +57,7 @@
include "constants.h"
logical :: initialfield
- integer :: NSOURCE
+ integer :: NSOURCES
integer :: npgeo,ngnod,myrank,ipass,nproc
integer :: npoin,nspec,nelem_acoustic_surface
@@ -77,10 +77,10 @@
character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec) :: network_name
! for sources
- integer, dimension(NSOURCE) :: source_type
- integer, dimension(NSOURCE) :: ispec_selected_source,is_proc_source,nb_proc_source,iglob_source
- real(kind=CUSTOM_REAL), dimension(NSOURCE,NDIM,NGLLX,NGLLZ) :: sourcearray
- double precision, dimension(NSOURCE) :: x_source,z_source,xi_source,gamma_source,Mxx,Mzz,Mxz
+ integer, dimension(NSOURCES) :: source_type
+ integer, dimension(NSOURCES) :: ispec_selected_source,is_proc_source,nb_proc_source,iglob_source
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLZ) :: sourcearray
+ double precision, dimension(NSOURCES) :: x_source,z_source,xi_source,gamma_source,Mxx,Mzz,Mxz
logical, dimension(nspec) :: elastic,poroelastic
integer, dimension(ngnod,nspec) :: knods
@@ -95,26 +95,27 @@
! Local variables
integer i_source,ispec,ispec_acoustic_surface
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
- if(source_type(i_source) == 1) then
+ if(source_type(i_source) == 1) then
- ! collocated force source
- call locate_source_force(ibool,coord,nspec,npoin,xigll,zigll,x_source(i_source),z_source(i_source), &
+ ! collocated force source
+ call locate_source_force(ibool,coord,nspec,npoin,xigll,zigll,x_source(i_source),z_source(i_source), &
ispec_selected_source(i_source),is_proc_source(i_source),nb_proc_source(i_source),&
nproc,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,ngnod,npgeo,ipass,&
iglob_source(i_source))
-! check that acoustic source is not exactly on the free surface because pressure is zero there
- if(is_proc_source(i_source) == 1) then
- do ispec_acoustic_surface = 1,nelem_acoustic_surface
- ispec = acoustic_surface(1,ispec_acoustic_surface)
- ixmin = acoustic_surface(2,ispec_acoustic_surface)
- ixmax = acoustic_surface(3,ispec_acoustic_surface)
- izmin = acoustic_surface(4,ispec_acoustic_surface)
- izmax = acoustic_surface(5,ispec_acoustic_surface)
- if( .not. elastic(ispec) .and. .not. poroelastic(ispec) .and. ispec == ispec_selected_source(i_source) ) then
- if ( (izmin==1 .and. izmax==1 .and. ixmin==1 .and. ixmax==NGLLX .and. &
+ ! check that acoustic source is not exactly on the free surface because pressure is zero there
+ if(is_proc_source(i_source) == 1) then
+ do ispec_acoustic_surface = 1,nelem_acoustic_surface
+ ispec = acoustic_surface(1,ispec_acoustic_surface)
+ ixmin = acoustic_surface(2,ispec_acoustic_surface)
+ ixmax = acoustic_surface(3,ispec_acoustic_surface)
+ izmin = acoustic_surface(4,ispec_acoustic_surface)
+ izmax = acoustic_surface(5,ispec_acoustic_surface)
+ if( .not. elastic(ispec) .and. .not. poroelastic(ispec) .and. &
+ ispec == ispec_selected_source(i_source) ) then
+ if ( (izmin==1 .and. izmax==1 .and. ixmin==1 .and. ixmax==NGLLX .and. &
gamma_source(i_source) < -0.99d0) .or.&
(izmin==NGLLZ .and. izmax==NGLLZ .and. ixmin==1 .and. ixmax==NGLLX .and. &
gamma_source(i_source) > 0.99d0) .or.&
@@ -130,37 +131,40 @@
gamma_source(i_source) > 0.99d0 .and. xi_source(i_source) < -0.99d0) .or.&
(izmin==NGLLZ .and. izmax==NGLLZ .and. ixmin==NGLLX .and. ixmax==NGLLX .and. &
gamma_source(i_source) > 0.99d0 .and. xi_source(i_source) > 0.99d0) ) then
- call exit_MPI('an acoustic source cannot be located exactly on the free surface because pressure is zero there')
+ call exit_MPI('an acoustic source cannot be located exactly '// &
+ 'on the free surface because pressure is zero there')
+ endif
endif
- endif
- enddo
- endif
+ enddo
+ endif
- else if(source_type(i_source) == 2) then
- ! moment-tensor source
- call locate_source_moment_tensor(ibool,coord,nspec,npoin,xigll,zigll,x_source(i_source),z_source(i_source), &
+ else if(source_type(i_source) == 2) then
+ ! moment-tensor source
+ call locate_source_moment_tensor(ibool,coord,nspec,npoin,xigll,zigll,x_source(i_source),z_source(i_source), &
ispec_selected_source(i_source),is_proc_source(i_source),nb_proc_source(i_source),&
nproc,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,ngnod,npgeo,ipass)
- ! compute source array for moment-tensor source
- call compute_arrays_source(ispec_selected_source(i_source),xi_source(i_source),gamma_source(i_source),&
+ ! compute source array for moment-tensor source
+ call compute_arrays_source(ispec_selected_source(i_source),xi_source(i_source),gamma_source(i_source),&
sourcearray(i_source,:,:,:), &
Mxx(i_source),Mzz(i_source),Mxz(i_source),xix,xiz,gammax,gammaz,xigll,zigll,nspec)
- else if(.not.initialfield) then
- call exit_MPI('incorrect source type')
- endif
+ else if(.not.initialfield) then
+
+ call exit_MPI('incorrect source type')
+
+ endif
- ! locate receivers in the mesh
- call locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank,&
+ ! locate receivers in the mesh
+ call locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank,&
st_xval,st_zval,ispec_selected_rec, &
xi_receiver,gamma_receiver,station_name,network_name,x_source(i_source),z_source(i_source),&
coorg,knods,ngnod,npgeo,ipass, &
x_final_receiver, z_final_receiver)
- enddo ! do i_source=1,NSOURCE
+ enddo ! do i_source=1,NSOURCES
-end subroutine setup_sources_receivers
+ end subroutine setup_sources_receivers
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2011-02-21 01:39:04 UTC (rev 17925)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2011-02-21 03:36:42 UTC (rev 17926)
@@ -328,7 +328,7 @@
! character(len=80) datlin
- integer NSOURCE,i_source
+ integer NSOURCES,i_source
integer, dimension(:), allocatable :: source_type,time_function_type
double precision, dimension(:), allocatable :: x_source,z_source,xi_source,gamma_source,&
Mxx,Mzz,Mxz,f0,tshift_src,factor,angleforce
@@ -385,11 +385,14 @@
double precision :: mul_relaxed,lambdal_relaxed,lambdalplus2mul_relaxed,kappal
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_elastic,veloc_elastic,displ_elastic
- double precision, dimension(:,:), allocatable :: coord, flagrange,xinterp,zinterp,Uxinterp,Uzinterp,vector_field_display
+ double precision, dimension(:,:), allocatable :: &
+ coord, flagrange,xinterp,zinterp,Uxinterp,Uzinterp,vector_field_display
! material properties of the poroelastic medium (solid phase:s and fluid phase [defined as w=phi(u_f-u_s)]: w)
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accels_poroelastic,velocs_poroelastic,displs_poroelastic
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accelw_poroelastic,velocw_poroelastic,displw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
+ accels_poroelastic,velocs_poroelastic,displs_poroelastic
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
+ accelw_poroelastic,velocw_poroelastic,displw_poroelastic
double precision, dimension(:), allocatable :: porosity,tortuosity
double precision, dimension(:,:), allocatable :: density,permeability
@@ -403,20 +406,21 @@
double precision, dimension(:,:), allocatable :: anisotropy
! for acoustic medium
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
+ potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
! inverse mass matrices
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_elastic,rmass_inverse_acoustic
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
+ rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic
! to evaluate cpI, cpII, and cs, and rI (poroelastic medium)
real(kind=CUSTOM_REAL) :: rhol_s,rhol_f,rhol_bar,phil,tortl
real(kind=CUSTOM_REAL) :: mul_s,kappal_s
real(kind=CUSTOM_REAL) :: kappal_f
-! double precision :: etal_f
real(kind=CUSTOM_REAL) :: mul_fr,kappal_fr
real(kind=CUSTOM_REAL) :: D_biot,H_biot,C_biot,M_biot,B_biot,cpIsquare,cpIIsquare,cssquare
- real(kind=CUSTOM_REAL) :: ratio,dd1 !gamma1,gamma2,gamma3,gamma4,afactor,bfactor,cfactor
+ real(kind=CUSTOM_REAL) :: ratio,dd1
double precision, dimension(:,:,:), allocatable :: vpext,vsext,rhoext
double precision, dimension(:,:,:), allocatable :: Qp_attenuationext,Qs_attenuationext
@@ -434,7 +438,6 @@
integer, dimension(:), allocatable :: ispec_selected_source,iglob_source,&
is_proc_source,nb_proc_source
- double precision displnorm_all,displnorm_all_glob
double precision, dimension(:), allocatable :: aval
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: source_time_function
double precision, external :: netlib_specfun_erf
@@ -608,24 +611,14 @@
double precision, external :: hgll
! timer to count elapsed time
- double precision :: time_start,time_end
- integer :: year_start,year_end,month_start,month_end
+ double precision :: time_start
+ integer :: year_start,month_start
- double precision :: tCPU,t_remain,t_total
- integer :: ihours,iminutes,iseconds,int_tCPU, &
- ihours_remain,iminutes_remain,iseconds_remain,int_t_remain, &
- ihours_total,iminutes_total,iseconds_total,int_t_total
! to determine date and time at which the run will finish
character(len=8) datein
character(len=10) timein
character(len=5) :: zone
integer, dimension(8) :: time_values
- character(len=3), dimension(12) :: month_name
- character(len=3), dimension(0:6) :: weekday_name
- data month_name /'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/
- data weekday_name /'Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat'/
- integer :: year,mon,day,hr,minutes,timestamp,julian_day_number,day_of_week
- integer, external :: idaywk
! for MPI and partitioning
integer :: ier
@@ -803,40 +796,40 @@
seismotype,imagetype,assign_external_model,READ_EXTERNAL_SEP_FILE, &
outputgrid,OUTPUT_ENERGY,TURN_ATTENUATION_ON, &
TURN_VISCATTENUATION_ON,Q0,freq0,p_sv, &
- NSTEP,deltat,NTSTEP_BETWEEN_OUTPUT_SEISMO,NSOURCE)
+ NSTEP,deltat,NTSTEP_BETWEEN_OUTPUT_SEISMO,NSOURCES)
!
!--- source information
!
if(ipass == 1) then
- allocate( source_type(NSOURCE) )
- allocate( time_function_type(NSOURCE) )
- allocate( x_source(NSOURCE) )
- allocate( z_source(NSOURCE) )
- allocate( f0(NSOURCE) )
- allocate( tshift_src(NSOURCE) )
- allocate( factor(NSOURCE) )
- allocate( angleforce(NSOURCE) )
- allocate( Mxx(NSOURCE) )
- allocate( Mxz(NSOURCE) )
- allocate( Mzz(NSOURCE) )
- allocate( aval(NSOURCE) )
- allocate( ispec_selected_source(NSOURCE) )
- allocate( iglob_source(NSOURCE) )
- allocate( source_courbe_eros(NSOURCE) )
- allocate( xi_source(NSOURCE) )
- allocate( gamma_source(NSOURCE) )
- allocate( is_proc_source(NSOURCE) )
- allocate( nb_proc_source(NSOURCE) )
- allocate( sourcearray(NSOURCE,NDIM,NGLLX,NGLLZ) )
+ allocate( source_type(NSOURCES) )
+ allocate( time_function_type(NSOURCES) )
+ allocate( x_source(NSOURCES) )
+ allocate( z_source(NSOURCES) )
+ allocate( f0(NSOURCES) )
+ allocate( tshift_src(NSOURCES) )
+ allocate( factor(NSOURCES) )
+ allocate( angleforce(NSOURCES) )
+ allocate( Mxx(NSOURCES) )
+ allocate( Mxz(NSOURCES) )
+ allocate( Mzz(NSOURCES) )
+ allocate( aval(NSOURCES) )
+ allocate( ispec_selected_source(NSOURCES) )
+ allocate( iglob_source(NSOURCES) )
+ allocate( source_courbe_eros(NSOURCES) )
+ allocate( xi_source(NSOURCES) )
+ allocate( gamma_source(NSOURCES) )
+ allocate( is_proc_source(NSOURCES) )
+ allocate( nb_proc_source(NSOURCES) )
+ allocate( sourcearray(NSOURCES,NDIM,NGLLX,NGLLZ) )
endif
! reads in source infos
- call read_databases_sources(NSOURCE,source_type,time_function_type, &
+ call read_databases_sources(NSOURCES,source_type,time_function_type, &
x_source,z_source,Mxx,Mzz,Mxz,f0,tshift_src,factor,angleforce)
! sets source parameters
- call set_sources(myrank,NSOURCE,source_type,time_function_type, &
+ call set_sources(myrank,NSOURCES,source_type,time_function_type, &
x_source,z_source,Mxx,Mzz,Mxz,f0,tshift_src,factor,angleforce,aval, &
t0,initialfield,ipass,deltat)
@@ -844,6 +837,11 @@
!---- read attenuation information
!
call read_databases_atten(N_SLS,f0_attenuation)
+
+ ! if source is not a Dirac or Heavyside then f0_attenuation is f0 of the first source
+ if(.not. (time_function_type(1) == 4 .or. time_function_type(1) == 5)) then
+ f0_attenuation = f0(1)
+ endif
!
@@ -1081,6 +1079,12 @@
allocate(jend_left_poro(nelemabs))
allocate(jbegin_right_poro(nelemabs))
allocate(jend_right_poro(nelemabs))
+
+ allocate(ib_left(nelemabs))
+ allocate(ib_right(nelemabs))
+ allocate(ib_bottom(nelemabs))
+ allocate(ib_top(nelemabs))
+
endif
!
@@ -1093,16 +1097,16 @@
if( anyabs ) then
- nspec_xmin = ZERO
- nspec_xmax = ZERO
- nspec_zmin = ZERO
- nspec_zmax = ZERO
- if(ipass == 1) then
- allocate(ib_left(nelemabs))
- allocate(ib_right(nelemabs))
- allocate(ib_bottom(nelemabs))
- allocate(ib_top(nelemabs))
- endif
+ nspec_xmin = 0
+ nspec_xmax = 0
+ nspec_zmin = 0
+ nspec_zmax = 0
+! if(ipass == 1) then
+! allocate(ib_left(nelemabs))
+! allocate(ib_right(nelemabs))
+! allocate(ib_bottom(nelemabs))
+! allocate(ib_top(nelemabs))
+! endif
do inum = 1,nelemabs
if (codeabs(IBOTTOM,inum)) then
nspec_zmin = nspec_zmin + 1
@@ -1177,12 +1181,12 @@
else
- if(.not. allocated(ib_left)) then
- allocate(ib_left(1))
- allocate(ib_right(1))
- allocate(ib_bottom(1))
- allocate(ib_top(1))
- endif
+! if(.not. allocated(ib_left)) then
+! allocate(ib_left(1))
+! allocate(ib_right(1))
+! allocate(ib_bottom(1))
+! allocate(ib_top(1))
+! endif
if(.not. allocated(b_absorb_elastic_left)) then
allocate(b_absorb_elastic_left(1,1,1,1))
@@ -1622,8 +1626,8 @@
allocate(hgammar_store(nrec,NGLLZ))
! allocate Lagrange interpolators for sources
- allocate(hxis_store(NSOURCE,NGLLX))
- allocate(hgammas_store(NSOURCE,NGLLZ))
+ allocate(hxis_store(NSOURCES,NGLLX))
+ allocate(hgammas_store(NSOURCES,NGLLZ))
! allocate other global arrays
allocate(coord(NDIM,npoin))
@@ -1843,7 +1847,7 @@
!---- define actual location of source and receivers
- call setup_sources_receivers(NSOURCE,initialfield,source_type,&
+ call setup_sources_receivers(NSOURCES,initialfield,source_type,&
coord,ibool,npoin,nspec,nelem_acoustic_surface,acoustic_surface,elastic,poroelastic, &
x_source,z_source,ispec_selected_source,ispec_selected_rec, &
is_proc_source,nb_proc_source,ipass,&
@@ -1951,7 +1955,7 @@
! for the source
if (force_normal_to_surface) then
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
if (is_proc_source(i_source) == 1) then
distmin = HUGEVAL
do i = 1, nnodes_tangential_curve
@@ -2004,7 +2008,7 @@
angleforce(i_source) = angleforce_recv
#endif
endif ! if (is_proc_source(i_source) == 1)
- enddo ! do i_source=1,NSOURCE
+ enddo ! do i_source=1,NSOURCES
endif ! if (force_normal_to_surface)
! CHRIS --- how to deal with multiple source. Use first source now. ---
@@ -2162,7 +2166,7 @@
enddo
! define and store Lagrange interpolators at all the sources
- do i = 1,NSOURCE
+ do i = 1,NSOURCES
call lagrange_any(xi_source(i),NGLLX,xigll,hxis,hpxis)
call lagrange_any(gamma_source(i),NGLLZ,zigll,hgammas,hpgammas)
hxis_store(i,:) = hxis(:)
@@ -2689,13 +2693,14 @@
else
stop 'incorrect value of DISPLAY_SUBSET_OPTION'
endif
- call checkgrid(vpext,vsext,rhoext,density,poroelastcoef,porosity,tortuosity,permeability,ibool,kmato, &
+ call checkgrid(vpext,vsext,rhoext,density,poroelastcoef, &
+ porosity,tortuosity,permeability,ibool,kmato, &
coord,npoin,vpImin,vpImax,vpIImin,vpIImax, &
assign_external_model,nspec,UPPER_LIMIT_DISPLAY,numat,deltat, &
- f0,tshift_src,initialfield,time_function_type, &
+ f0,initialfield,time_function_type, &
coorg,xinterp,zinterp,shape2D_display,knods,simulation_title, &
npgeo,pointsdisp,ngnod,any_elastic,any_poroelastic,all_anisotropic, &
- myrank,nproc,NSOURCE,poroelastic, &
+ myrank,nproc,NSOURCES,poroelastic, &
freq0,Q0,TURN_VISCATTENUATION_ON)
! convert receiver angle to radians
@@ -2988,7 +2993,7 @@
! Calculation of the initial field for a plane wave
call prepare_initialfield(myrank,any_acoustic,any_poroelastic,over_critical_angle, &
- NSOURCE,source_type,angleforce,x_source,z_source,f0, &
+ NSOURCES,source_type,angleforce,x_source,z_source,f0, &
npoin,numat,poroelastcoef,density,coord, &
angleforce_refl,c_inc,c_refl,cploc,csloc,time_offset, &
A_plane, B_plane, C_plane, &
@@ -3002,7 +3007,7 @@
call prepare_initialfield_paco(myrank,nelemabs,left_bound,right_bound,bot_bound, &
numabs,codeabs,ibool,nspec, &
- source_type,NSOURCE,c_inc,c_refl, &
+ source_type,NSOURCES,c_inc,c_refl, &
count_bottom,count_left,count_right)
allocate(v0x_left(count_left,NSTEP))
@@ -3069,11 +3074,11 @@
! compute the source time function and store it in a text file
if(.not. initialfield) then
- allocate(source_time_function(NSOURCE,NSTEP))
+ allocate(source_time_function(NSOURCES,NSTEP))
source_time_function(:,:) = 0.0_CUSTOM_REAL
! computes source time function array
- call prepare_source_time_function(myrank,NSTEP,NSOURCE,source_time_function, &
+ call prepare_source_time_function(myrank,NSTEP,NSOURCES,source_time_function, &
time_function_type,f0,tshift_src,factor,aval, &
t0,nb_proc_source,deltat)
else
@@ -3666,16 +3671,16 @@
!
if (myrank == 0) write(IOUT,400)
-! count elapsed wall-clock time
+ ! count elapsed wall-clock time
call date_and_time(datein,timein,zone,time_values)
-! time_values(1): year
-! time_values(2): month of the year
-! time_values(3): day of the month
-! time_values(5): hour of the day
-! time_values(6): minutes of the hour
-! time_values(7): seconds of the minute
-! time_values(8): milliseconds of the second
-! this fails if we cross the end of the month
+ ! time_values(1): year
+ ! time_values(2): month of the year
+ ! time_values(3): day of the month
+ ! time_values(5): hour of the day
+ ! time_values(6): minutes of the hour
+ ! time_values(7): seconds of the minute
+ ! time_values(8): milliseconds of the second
+ ! this fails if we cross the end of the month
time_start = 86400.d0*time_values(3) + 3600.d0*time_values(5) + &
60.d0*time_values(6) + time_values(7) + time_values(8) / 1000.d0
month_start = time_values(2)
@@ -3986,55 +3991,95 @@
!-----------------------------------------
if(any_acoustic) then
- potential_acoustic = potential_acoustic + deltat*potential_dot_acoustic + deltatsquareover2*potential_dot_dot_acoustic
- potential_dot_acoustic = potential_dot_acoustic + deltatover2*potential_dot_dot_acoustic
+ potential_acoustic = potential_acoustic &
+ + deltat*potential_dot_acoustic &
+ + deltatsquareover2*potential_dot_dot_acoustic
+ potential_dot_acoustic = potential_dot_acoustic &
+ + deltatover2*potential_dot_dot_acoustic
potential_dot_dot_acoustic = ZERO
- if(SIMULATION_TYPE == 2) then ! Adjoint calculation
- b_potential_acoustic = b_potential_acoustic + b_deltat*b_potential_dot_acoustic + &
- b_deltatsquareover2*b_potential_dot_dot_acoustic
- b_potential_dot_acoustic = b_potential_dot_acoustic + b_deltatover2*b_potential_dot_dot_acoustic
- b_potential_dot_dot_acoustic = ZERO
- endif
+ if(SIMULATION_TYPE == 2) then ! Adjoint calculation
+ b_potential_acoustic = b_potential_acoustic &
+ + b_deltat*b_potential_dot_acoustic &
+ + b_deltatsquareover2*b_potential_dot_dot_acoustic
+ b_potential_dot_acoustic = b_potential_dot_acoustic &
+ + b_deltatover2*b_potential_dot_dot_acoustic
+ b_potential_dot_dot_acoustic = ZERO
+ endif
! free surface for an acoustic medium
if ( nelem_acoustic_surface > 0 ) then
call enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
- potential_acoustic,acoustic_surface,ibool,nelem_acoustic_surface,npoin,nspec)
+ potential_acoustic,acoustic_surface, &
+ ibool,nelem_acoustic_surface,npoin,nspec)
- if(SIMULATION_TYPE == 2) then ! Adjoint calculation
- call enforce_acoustic_free_surface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
- b_potential_acoustic,acoustic_surface,ibool,nelem_acoustic_surface,npoin,nspec)
- endif
+ if(SIMULATION_TYPE == 2) then ! Adjoint calculation
+ call enforce_acoustic_free_surface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
+ b_potential_acoustic,acoustic_surface, &
+ ibool,nelem_acoustic_surface,npoin,nspec)
+ endif
endif
! *********************************************************
! ************* compute forces for the acoustic elements
! *********************************************************
- call compute_forces_acoustic(npoin,nspec,nelemabs,numat,it,NSTEP, &
+! call compute_forces_acoustic(npoin,nspec,nelemabs,numat,it,NSTEP, &
+! anyabs,assign_external_model,ibool,kmato,numabs, &
+! elastic,poroelastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
+! potential_acoustic,b_potential_dot_dot_acoustic,b_potential_acoustic, &
+! density,poroelastcoef,xix,xiz,gammax,gammaz,jacobian, &
+! vpext,rhoext,hprime_xx,hprimewgll_xx, &
+! hprime_zz,hprimewgll_zz,wxgll,wzgll, &
+! ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+! jbegin_left,jend_left,jbegin_right,jend_right, &
+! SIMULATION_TYPE,SAVE_FORWARD,b_absorb_acoustic_left,&
+! b_absorb_acoustic_right,b_absorb_acoustic_bottom,&
+! b_absorb_acoustic_top,nspec_xmin,nspec_xmax,&
+! nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top)
+
+
+ call compute_forces_acoustic_2(npoin,nspec,nelemabs,numat,it,NSTEP, &
anyabs,assign_external_model,ibool,kmato,numabs, &
elastic,poroelastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
- potential_acoustic,b_potential_dot_dot_acoustic,b_potential_acoustic, &
+ potential_acoustic, &
density,poroelastcoef,xix,xiz,gammax,gammaz,jacobian, &
vpext,rhoext,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
- jbegin_left,jend_left,jbegin_right,jend_right,SIMULATION_TYPE,SAVE_FORWARD,b_absorb_acoustic_left,&
- b_absorb_acoustic_right,b_absorb_acoustic_bottom,&
- b_absorb_acoustic_top,nspec_xmin,nspec_xmax,&
- nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top)
+ jbegin_left,jend_left,jbegin_right,jend_right, &
+ SIMULATION_TYPE,SAVE_FORWARD,nspec_xmin,nspec_xmax,&
+ nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top, &
+ b_absorb_acoustic_left,b_absorb_acoustic_right, &
+ b_absorb_acoustic_bottom,b_absorb_acoustic_top)
+ if( SIMULATION_TYPE == 2 ) then
+ call compute_forces_acoustic_2(npoin,nspec,nelemabs,numat,it,NSTEP, &
+ anyabs,assign_external_model,ibool,kmato,numabs, &
+ elastic,poroelastic,codeabs,b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
+ b_potential_acoustic, &
+ density,poroelastcoef,xix,xiz,gammax,gammaz,jacobian, &
+ vpext,rhoext,hprime_xx,hprimewgll_xx, &
+ hprime_zz,hprimewgll_zz,wxgll,wzgll, &
+ ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
+ jbegin_left,jend_left,jbegin_right,jend_right, &
+ SIMULATION_TYPE,SAVE_FORWARD,nspec_xmin,nspec_xmax,&
+ nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top, &
+ b_absorb_acoustic_left,b_absorb_acoustic_right, &
+ b_absorb_acoustic_bottom,b_absorb_acoustic_top)
+ endif
- if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1) then
+
+ if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1) then
+
!--- left absorbing boundary
- if(nspec_xmin >0) then
- do ispec = 1,nspec_xmin
- do i=1,NGLLZ
- write(65) b_absorb_acoustic_left(i,ispec,it)
- enddo
- enddo
- endif
+ if(nspec_xmin >0) then
+ do ispec = 1,nspec_xmin
+ do i=1,NGLLZ
+ write(65) b_absorb_acoustic_left(i,ispec,it)
+ enddo
+ enddo
+ endif
!--- right absorbing boundary
if(nspec_xmax >0) then
@@ -4263,7 +4308,7 @@
! --- add the source
if(.not. initialfield) then
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
! if this processor carries the source and the source element is acoustic
if (is_proc_source(i_source) == 1 .and. .not. elastic(ispec_selected_source(i_source)) .and. &
.not. poroelastic(ispec_selected_source(i_source))) then
@@ -4299,7 +4344,7 @@
endif
endif ! if this processor carries the source and the source element is acoustic
- enddo ! do i_source=1,NSOURCE
+ enddo ! do i_source=1,NSOURCES
if(SIMULATION_TYPE == 2) then ! adjoint wavefield
irec_local = 0
@@ -4366,13 +4411,15 @@
! free surface for an acoustic medium
if ( nelem_acoustic_surface > 0 ) then
- call enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
- potential_acoustic,acoustic_surface,ibool,nelem_acoustic_surface,npoin,nspec)
+ call enforce_acoustic_free_surface(potential_dot_dot_acoustic,potential_dot_acoustic, &
+ potential_acoustic,acoustic_surface, &
+ ibool,nelem_acoustic_surface,npoin,nspec)
- if(SIMULATION_TYPE == 2) then
- call enforce_acoustic_free_surface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
- b_potential_acoustic,acoustic_surface,ibool,nelem_acoustic_surface,npoin,nspec)
- endif
+ if(SIMULATION_TYPE == 2) then
+ call enforce_acoustic_free_surface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
+ b_potential_acoustic,acoustic_surface, &
+ ibool,nelem_acoustic_surface,npoin,nspec)
+ endif
endif
@@ -4400,7 +4447,7 @@
A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset, f0(1),&
v0x_left(1,it),v0z_left(1,it),v0x_right(1,it),v0z_right(1,it),v0x_bot(1,it),v0z_bot(1,it), &
t0x_left(1,it),t0z_left(1,it),t0x_right(1,it),t0z_right(1,it),t0x_bot(1,it),t0z_bot(1,it), &
- count_left,count_right,count_bottom,over_critical_angle,NSOURCE,nrec,SIMULATION_TYPE,SAVE_FORWARD, &
+ count_left,count_right,count_bottom,over_critical_angle,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD, &
b_absorb_elastic_left,b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top, &
nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top,mu_k,kappa_k)
@@ -4878,7 +4925,7 @@
! --- add the source if it is a collocated force
if(.not. initialfield) then
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
! if this processor carries the source and the source element is elastic
if (is_proc_source(i_source) == 1 .and. elastic(ispec_selected_source(i_source))) then
@@ -4939,7 +4986,7 @@
endif
endif ! if this processor carries the source and the source element is elastic
- enddo ! do i_source=1,NSOURCE
+ enddo ! do i_source=1,NSOURCES
endif ! if not using an initial field
endif !if(any_elastic)
@@ -5019,7 +5066,7 @@
b_viscodampx,b_viscodampz,&
ibegin_bottom_poro,iend_bottom_poro,ibegin_top_poro,iend_top_poro, &
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
- mufr_k,B_k,NSOURCE,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
+ mufr_k,B_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
b_absorb_poro_s_left,b_absorb_poro_s_right,b_absorb_poro_s_bottom,b_absorb_poro_s_top,&
nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top,f0(1),freq0,Q0)
@@ -5042,7 +5089,7 @@
b_viscodampx,b_viscodampz,&
ibegin_bottom_poro,iend_bottom_poro,ibegin_top_poro,iend_top_poro, &
jbegin_left_poro,jend_left_poro,jbegin_right_poro,jend_right_poro,&
- C_k,M_k,NSOURCE,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
+ C_k,M_k,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,&
b_absorb_poro_w_left,b_absorb_poro_w_right,b_absorb_poro_w_bottom,b_absorb_poro_w_top,&
nspec_xmin,nspec_xmax,nspec_zmin,nspec_zmax,ib_left,ib_right,ib_bottom,ib_top,f0(1),freq0,Q0)
@@ -5526,7 +5573,7 @@
! --- add the source if it is a collocated force
if(.not. initialfield) then
- do i_source=1,NSOURCE
+ do i_source=1,NSOURCES
! if this processor carries the source and the source element is elastic
if (is_proc_source(i_source) == 1 .and. poroelastic(ispec_selected_source(i_source))) then
@@ -5576,7 +5623,7 @@
endif
endif ! if this processor carries the source and the source element is elastic
- enddo ! do i_source=1,NSOURCE
+ enddo ! do i_source=1,NSOURCES
endif ! if not using an initial field
endif !if(any_poroelastic)
@@ -5723,21 +5770,22 @@
! acoustic medium
if(any_acoustic) then
- write(outputname,'(a,i6.6,a)') 'lastframe_acoustic',myrank,'.bin'
- open(unit=55,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
- do j=1,npoin
- read(55) b_potential_acoustic(j),&
- b_potential_dot_acoustic(j),&
- b_potential_dot_dot_acoustic(j)
- enddo
- close(55)
+ write(outputname,'(a,i6.6,a)') 'lastframe_acoustic',myrank,'.bin'
+ open(unit=55,file='OUTPUT_FILES/'//outputname,status='old',action='read',form='unformatted')
+ do j=1,npoin
+ read(55) b_potential_acoustic(j),&
+ b_potential_dot_acoustic(j),&
+ b_potential_dot_dot_acoustic(j)
+ enddo
+ close(55)
! free surface for an acoustic medium
- if ( nelem_acoustic_surface > 0 ) then
- call enforce_acoustic_free_surface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
- b_potential_acoustic,acoustic_surface,ibool,nelem_acoustic_surface,npoin,nspec)
+ if ( nelem_acoustic_surface > 0 ) then
+ call enforce_acoustic_free_surface(b_potential_dot_dot_acoustic,b_potential_dot_acoustic, &
+ b_potential_acoustic,acoustic_surface, &
+ ibool,nelem_acoustic_surface,npoin,nspec)
+ endif
endif
- endif
! elastic medium
if(any_elastic) then
@@ -5838,175 +5886,14 @@
!---- display time step and max of norm of displacement
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5 .or. it == NSTEP) then
-
- if (myrank == 0) then
- write(IOUT,*)
- if(time >= 1.d-3 .and. time < 1000.d0) then
- write(IOUT,"('Time step number ',i7,' t = ',f9.4,' s out of ',i7)") it,time,NSTEP
- else
- write(IOUT,"('Time step number ',i7,' t = ',1pe12.6,' s out of ',i7)") it,time,NSTEP
- endif
- write(IOUT,*) 'We have done ',sngl(100.d0*dble(it-1)/dble(NSTEP-1)),'% of the total'
- endif
-
- if(any_elastic_glob) then
- if(any_elastic) then
- displnorm_all = maxval(sqrt(displ_elastic(1,:)**2 + displ_elastic(2,:)**2 + displ_elastic(3,:)**2))
- else
- displnorm_all = 0.d0
- endif
- displnorm_all_glob = displnorm_all
-#ifdef USE_MPI
- call MPI_ALLREDUCE (displnorm_all, displnorm_all_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
-#endif
- if (myrank == 0) write(IOUT,*) 'Max norm of vector field in solid (elastic) = ',displnorm_all_glob
-! check stability of the code in solid, exit if unstable
-! negative values can occur with some compilers when the unstable value is greater
-! than the greatest possible floating-point number of the machine
- if(displnorm_all_glob > STABILITY_THRESHOLD .or. displnorm_all_glob < 0) &
- call exit_MPI('code became unstable and blew up in solid (elastic)')
- endif
-
- if(any_poroelastic_glob) then
- if(any_poroelastic) then
- displnorm_all = maxval(sqrt(displs_poroelastic(1,:)**2 + displs_poroelastic(2,:)**2))
- else
- displnorm_all = 0.d0
- endif
- displnorm_all_glob = displnorm_all
-#ifdef USE_MPI
- call MPI_ALLREDUCE (displnorm_all, displnorm_all_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
-#endif
- if (myrank == 0) write(IOUT,*) 'Max norm of vector field in solid (poroelastic) = ',displnorm_all_glob
-! check stability of the code in solid, exit if unstable
-! negative values can occur with some compilers when the unstable value is greater
-! than the greatest possible floating-point number of the machine
- if(displnorm_all_glob > STABILITY_THRESHOLD .or. displnorm_all_glob < 0) &
- call exit_MPI('code became unstable and blew up in solid (poroelastic)')
-
- if(any_poroelastic) then
- displnorm_all = maxval(sqrt(displw_poroelastic(1,:)**2 + displw_poroelastic(2,:)**2))
- else
- displnorm_all = 0.d0
- endif
- displnorm_all_glob = displnorm_all
-#ifdef USE_MPI
- call MPI_ALLREDUCE (displnorm_all, displnorm_all_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
-#endif
- if (myrank == 0) write(IOUT,*) 'Max norm of vector field in fluid (poroelastic) = ',displnorm_all_glob
-! check stability of the code in solid, exit if unstable
-! negative values can occur with some compilers when the unstable value is greater
-! than the greatest possible floating-point number of the machine
- if(displnorm_all_glob > STABILITY_THRESHOLD .or. displnorm_all_glob < 0) &
- call exit_MPI('code became unstable and blew up in fluid (poroelastic)')
- endif
-
- if(any_acoustic_glob) then
- if(any_acoustic) then
- displnorm_all = maxval(abs(potential_acoustic(:)))
- else
- displnorm_all = 0.d0
- endif
- displnorm_all_glob = displnorm_all
-#ifdef USE_MPI
- call MPI_ALLREDUCE (displnorm_all, displnorm_all_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
-#endif
- if (myrank == 0) write(IOUT,*) 'Max absolute value of scalar field in fluid (acoustic) = ',displnorm_all_glob
-! check stability of the code in fluid, exit if unstable
-! negative values can occur with some compilers when the unstable value is greater
-! than the greatest possible floating-point number of the machine
- if(displnorm_all_glob > STABILITY_THRESHOLD .or. displnorm_all_glob < 0) &
- call exit_MPI('code became unstable and blew up in fluid (acoustic)')
- endif
-
-! count elapsed wall-clock time
- call date_and_time(datein,timein,zone,time_values)
-! time_values(1): year
-! time_values(2): month of the year
-! time_values(3): day of the month
-! time_values(5): hour of the day
-! time_values(6): minutes of the hour
-! time_values(7): seconds of the minute
-! time_values(8): milliseconds of the second
-! this fails if we cross the end of the month
- time_end = 86400.d0*time_values(3) + 3600.d0*time_values(5) + &
- 60.d0*time_values(6) + time_values(7) + time_values(8) / 1000.d0
- month_end = time_values(2)
- year_end = time_values(1)
-
-! elapsed time since beginning of the simulation
- if (myrank == 0) then
- if(month_end == month_start .and. year_end == year_start) then
- tCPU = time_end - time_start
- int_tCPU = int(tCPU)
- ihours = int_tCPU / 3600
- iminutes = (int_tCPU - 3600*ihours) / 60
- iseconds = int_tCPU - 3600*ihours - 60*iminutes
- write(IOUT,*) 'Elapsed time in seconds = ',tCPU
- write(IOUT,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
- write(IOUT,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
-
- ! compute estimated remaining simulation time
- t_remain = (NSTEP - it) * (tCPU/dble(it))
- int_t_remain = int(t_remain)
- ihours_remain = int_t_remain / 3600
- iminutes_remain = (int_t_remain - 3600*ihours_remain) / 60
- iseconds_remain = int_t_remain - 3600*ihours_remain - 60*iminutes_remain
- write(IOUT,*) 'Time steps remaining = ',NSTEP - it
- write(IOUT,*) 'Estimated remaining time in seconds = ',t_remain
- write(IOUT,"(' Estimated remaining time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
- ihours_remain,iminutes_remain,iseconds_remain
-
- ! compute estimated total simulation time
- t_total = t_remain + tCPU
- int_t_total = int(t_total)
- ihours_total = int_t_total / 3600
- iminutes_total = (int_t_total - 3600*ihours_total) / 60
- iseconds_total = int_t_total - 3600*ihours_total - 60*iminutes_total
- write(IOUT,*) 'Estimated total run time in seconds = ',t_total
- write(IOUT,"(' Estimated total run time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") &
- ihours_total,iminutes_total,iseconds_total
-
- if(it < NSTEP) then
-
- ! compute date and time at which the run should finish (useful for long runs); for simplicity only minutes
- ! are considered, seconds are ignored; in any case the prediction is not
- ! accurate down to seconds because of system and network fluctuations
- year = time_values(1)
- mon = time_values(2)
- day = time_values(3)
- hr = time_values(5)
- minutes = time_values(6)
-
- ! get timestamp in minutes of current date and time
- call convtime(timestamp,year,mon,day,hr,minutes)
-
- ! add remaining minutes
- timestamp = timestamp + nint(t_remain / 60.d0)
-
- ! get date and time of that future timestamp in minutes
- call invtime(timestamp,year,mon,day,hr,minutes)
-
- ! convert to Julian day to get day of the week
- call calndr(day,mon,year,julian_day_number)
- day_of_week = idaywk(julian_day_number)
-
- write(IOUT,"(' The run will finish approximately on: ',a3,' ',a3,' ',i2.2,', ',i4.4,' ',i2.2,':',i2.2)") &
- weekday_name(day_of_week),month_name(mon),day,year,hr,minutes
-
- endif
-
- write(IOUT,*)
- else
- write(IOUT,*) 'The calendar has crossed the end of the month during the simulation,'
- write(IOUT,*) 'cannot produce accurate CPU time estimates any more.'
- write(IOUT,*)
- endif
+ call check_stability(myrank,time,it,NSTEP,npoin, &
+ any_elastic_glob,any_elastic,displ_elastic, &
+ any_poroelastic_glob,any_poroelastic, &
+ displs_poroelastic,displw_poroelastic, &
+ any_acoustic_glob,any_acoustic,potential_acoustic, &
+ year_start,month_start,time_start)
endif
- if (myrank == 0) write(IOUT,*)
- endif
-
! loop on all the receivers to compute and store the seismograms
do irecloc = 1,nrecloc
@@ -6503,7 +6390,7 @@
Uxinterp,Uzinterp,flagrange,density,porosity,tortuosity,&
poroelastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,nelem_acoustic_surface,acoustic_edges, &
- simulation_title,npoin,npgeo,vpImin,vpImax,nrec,NSOURCE, &
+ simulation_title,npoin,npgeo,vpImin,vpImax,nrec,NSOURCES, &
colors,numbers,subsamp,imagetype,interpol,meshvect,modelvect, &
boundvect,assign_external_model,cutsnaps,sizemax_arrows,nelemabs,numat,pointsdisp, &
nspec,ngnod,coupled_acoustic_elastic,coupled_acoustic_poro,coupled_elastic_poro, &
@@ -6539,7 +6426,7 @@
Uxinterp,Uzinterp,flagrange,density,porosity,tortuosity,&
poroelastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,nelem_acoustic_surface,acoustic_edges, &
- simulation_title,npoin,npgeo,vpImin,vpImax,nrec,NSOURCE, &
+ simulation_title,npoin,npgeo,vpImin,vpImax,nrec,NSOURCES, &
colors,numbers,subsamp,imagetype,interpol,meshvect,modelvect, &
boundvect,assign_external_model,cutsnaps,sizemax_arrows,nelemabs,numat,pointsdisp, &
nspec,ngnod,coupled_acoustic_elastic,coupled_acoustic_poro,coupled_elastic_poro, &
@@ -6575,7 +6462,7 @@
Uxinterp,Uzinterp,flagrange,density,porosity,tortuosity,&
poroelastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,nelem_acoustic_surface,acoustic_edges, &
- simulation_title,npoin,npgeo,vpImin,vpImax,nrec,NSOURCE, &
+ simulation_title,npoin,npgeo,vpImin,vpImax,nrec,NSOURCES, &
colors,numbers,subsamp,imagetype,interpol,meshvect,modelvect, &
boundvect,assign_external_model,cutsnaps,sizemax_arrows,nelemabs,numat,pointsdisp, &
nspec,ngnod,coupled_acoustic_elastic,coupled_acoustic_poro,coupled_elastic_poro, &
More information about the CIG-COMMITS
mailing list