diff --git a/DATA/CMTSOLUTION b/DATA/CMTSOLUTION deleted file mode 120000 index 6bbf21120..000000000 --- a/DATA/CMTSOLUTION +++ /dev/null @@ -1 +0,0 @@ -../EXAMPLES/applications/homogeneous_halfspace/DATA/CMTSOLUTION \ No newline at end of file diff --git a/DATA/CMTSOLUTION b/DATA/CMTSOLUTION new file mode 100644 index 000000000..3d78a2e2f --- /dev/null +++ b/DATA/CMTSOLUTION @@ -0,0 +1,13 @@ +PDE 1999 01 01 00 00 00.00 67000 67000 -25000 4.2 4.2 homog_test +event name: homog_test +time shift: 0.0000 +half duration: 5.0 +latorUTM: 67000.0 +longorUTM: 67000.0 +depth: 30.0 +Mrr: -7.600000e+27 +Mtt: 7.700000e+27 +Mpp: -2.000000e+26 +Mrt: -2.500000e+28 +Mrp: 4.000000e+26 +Mtp: -2.500000e+27 diff --git a/DATA/Par_file b/DATA/Par_file deleted file mode 120000 index 9b1a3472b..000000000 --- a/DATA/Par_file +++ /dev/null @@ -1 +0,0 @@ -../EXAMPLES/applications/homogeneous_halfspace/DATA/Par_file \ No newline at end of file diff --git a/DATA/Par_file b/DATA/Par_file new file mode 100644 index 000000000..01a5e8972 --- /dev/null +++ b/DATA/Par_file @@ -0,0 +1,398 @@ +#----------------------------------------------------------- +# +# Simulation input parameters +# +#----------------------------------------------------------- + +# forward or adjoint simulation +# 1 = forward, 2 = adjoint, 3 = both simultaneously +SIMULATION_TYPE = 1 +# 0 = earthquake simulation, 1/2/3 = three steps in noise simulation +NOISE_TOMOGRAPHY = 0 +SAVE_FORWARD = .false. + +# solve a full FWI inverse problem from a single calling program with no I/Os, storing everything in memory, +# or run a classical forward or adjoint problem only and save the seismograms and/or sensitivity kernels to disk (with costlier I/Os) +INVERSE_FWI_FULL_PROBLEM = .false. + +# UTM projection parameters +# Use a negative zone number for the Southern hemisphere: +# The Northern hemisphere corresponds to zones +1 to +60, +# The Southern hemisphere corresponds to zones -1 to -60. +UTM_PROJECTION_ZONE = 11 +SUPPRESS_UTM_PROJECTION = .true. + +# number of MPI processors +NPROC = 4 + +# time step parameters +NSTEP = 5000 +DT = 0.05 + +# set to true to use local-time stepping (LTS) +LTS_MODE = .false. + +# Partitioning algorithm for decompose_mesh +# choose partitioner: 1==SCOTCH (default), 2==METIS, 3==PATOH, 4==ROWS_PART +PARTITIONING_TYPE = 1 + +#----------------------------------------------------------- +# +# LDDRK time scheme +# +#----------------------------------------------------------- +USE_LDDRK = .false. +INCREASE_CFL_FOR_LDDRK = .false. +RATIO_BY_WHICH_TO_INCREASE_IT = 1.4 + +#----------------------------------------------------------- +# +# Mesh +# +#----------------------------------------------------------- + +# Number of nodes for 2D and 3D shape functions for hexahedra. +# We use either 8-node mesh elements (bricks) or 27-node elements. +# If you use our internal mesher, the only option is 8-node bricks (27-node elements are not supported). +NGNOD = 8 + +# models: +# available options are: +# default (model parameters described by mesh properties) +# 1D models available are: +# 1d_prem,1d_socal,1d_cascadia +# 3D models available are: +# aniso,external,gll,salton_trough,tomo,SEP,coupled,... +MODEL = default + +# path for external tomographic models files +TOMOGRAPHY_PATH = DATA/tomo_files/ +# if you are using a SEP model (oil-industry format) +SEP_MODEL_DIRECTORY = DATA/my_SEP_model/ + +#----------------------------------------------------------- + +# parameters describing the model +APPROXIMATE_OCEAN_LOAD = .false. +TOPOGRAPHY = .false. +ATTENUATION = .false. +ANISOTROPY = .false. +GRAVITY = .false. + +# in case of attenuation, reference frequency in Hz at which the velocity values in the velocity model are given (unused otherwise) +ATTENUATION_f0_REFERENCE = 18.d0 + +# attenuation period range over which we try to mimic a constant Q factor +MIN_ATTENUATION_PERIOD = 999999998.d0 +MAX_ATTENUATION_PERIOD = 999999999.d0 +# ignore this range and ask the code to compute it automatically instead based on the estimated resolution of the mesh (use this unless you know what you are doing) +COMPUTE_FREQ_BAND_AUTOMATIC = .true. + +# Olsen's constant for Q_mu = constant * V_s attenuation rule +USE_OLSEN_ATTENUATION = .false. +OLSEN_ATTENUATION_RATIO = 0.05 + +#----------------------------------------------------------- +# +# Absorbing boundary conditions +# +#----------------------------------------------------------- + +# C-PML boundary conditions for a regional simulation +# (if set to .false., and STACEY_ABSORBING_CONDITIONS is also set to .false., you get a free surface instead +# in the case of elastic or viscoelastic mesh elements, and a rigid surface in the case of acoustic (fluid) elements +PML_CONDITIONS = .false. + +# C-PML top surface +PML_INSTEAD_OF_FREE_SURFACE = .false. + +# C-PML dominant frequency +f0_FOR_PML = 0.05555 + +# parameters used to rotate C-PML boundary conditions by a given angle (not completed yet) +# ROTATE_PML_ACTIVATE = .false. +# ROTATE_PML_ANGLE = 0. + +# absorbing boundary conditions for a regional simulation +# (if set to .false., and PML_CONDITIONS is also set to .false., you get a free surface instead +# in the case of elastic or viscoelastic mesh elements, and a rigid surface in the case of acoustic (fluid) elements +STACEY_ABSORBING_CONDITIONS = .true. + +# absorbing top surface (defined in mesh as 'free_surface_file') +STACEY_INSTEAD_OF_FREE_SURFACE = .false. + +# When STACEY_ABSORBING_CONDITIONS is set to .true. : +# absorbing conditions are defined in xmin, xmax, ymin, ymax and zmin +# this option BOTTOM_FREE_SURFACE can be set to .true. to +# make zmin free surface instead of absorbing condition +BOTTOM_FREE_SURFACE = .false. + +#----------------------------------------------------------- +# +# undoing attenuation and/or PMLs for sensitivity kernel calculations +# +#----------------------------------------------------------- + +# to undo attenuation and/or PMLs for sensitivity kernel calculations or forward runs with SAVE_FORWARD +# use the flag below. It performs undoing of attenuation and/or of PMLs in an exact way for sensitivity kernel calculations +# but requires disk space for temporary storage, and uses a significant amount of memory used as buffers for temporary storage. +# When that option is on the second parameter indicates how often the code dumps restart files to disk (if in doubt, use something between 100 and 1000). +UNDO_ATTENUATION_AND_OR_PML = .false. +NT_DUMP_ATTENUATION = 500 + +#----------------------------------------------------------- +# +# Visualization +# +#----------------------------------------------------------- + +# save AVS or OpenDX movies +# MOVIE_TYPE = 1 to show the top surface +# MOVIE_TYPE = 2 to show all the external faces of the mesh +CREATE_SHAKEMAP = .false. +MOVIE_SURFACE = .false. +MOVIE_TYPE = 1 +MOVIE_VOLUME = .false. +SAVE_DISPLACEMENT = .false. +MOVIE_VOLUME_STRESS = .false. +USE_HIGHRES_FOR_MOVIES = .false. +NTSTEP_BETWEEN_FRAMES = 200 +HDUR_MOVIE = 0.0 + +# save AVS or OpenDX mesh files to check the mesh +SAVE_MESH_FILES = .true. + +# path to store the local database file on each node +LOCAL_PATH = OUTPUT_FILES/DATABASES_MPI + +# interval at which we output time step info and max of norm of displacement +NTSTEP_BETWEEN_OUTPUT_INFO = 500 + +#----------------------------------------------------------- +# +# Sources +# +#----------------------------------------------------------- + +# sources and receivers Z coordinates given directly (i.e. as their true position) instead of as their depth +USE_SOURCES_RECEIVERS_Z = .false. + +# use a (tilted) FORCESOLUTION force point source (or several) instead of a CMTSOLUTION moment-tensor source. +# This can be useful e.g. for oil industry foothills simulations or asteroid simulations +# in which the source is a vertical force, normal force, tilted force, impact etc. +# If this flag is turned on, the FORCESOLUTION file must be edited by giving: +# - the corresponding time-shift parameter, +# - the half duration (hdur, in s) for Gaussian/Step function, dominant frequency (f0, in Hz) for Ricker, +# - the coordinates of the source, +# - the source time function type (0=Gaussian function, 1=Ricker wavelet, 2=Step function), +# - the magnitude of the force source, +# - the components of a (non necessarily unitary) direction vector for the force source in the E/N/Z_UP basis. +# The direction vector is made unitary internally in the code and thus only its direction matters here; +# its norm is ignored and the norm of the force used is the factor force source times the source time function. +USE_FORCE_POINT_SOURCE = .false. + +# set to true to use a Ricker source time function instead of the source time functions set by default +# to represent a (tilted) FORCESOLUTION force point source or a CMTSOLUTION moment-tensor source. +USE_RICKER_TIME_FUNCTION = .false. + +# use an external source time function +# you must add a file with your source time function and the file name path +# relative to working directory at the end of CMTSOLUTION or FORCESOLUTION file +# (with multiple sources, one file per source is required). +# This file must have a single column containing the amplitudes of the source at all time steps; +# time step size used must be equal to DT as defined at the beginning of this Par_file. +# Be sure when this option is .false. to remove the name of stf file in CMTSOLUTION or FORCESOLUTION +USE_EXTERNAL_SOURCE_FILE = .false. + +# print source time function +PRINT_SOURCE_TIME_FUNCTION = .false. + +# source encoding +# (for acoustic simulations only for now) determines source encoding factor +/-1 depending on sign of moment tensor +# (see e.g. Krebs et al., 2009. Fast full-wavefield seismic inversion using encoded sources, Geophysics, 74 (6), WCC177-WCC188.) +USE_SOURCE_ENCODING = .false. + +#----------------------------------------------------------- +# +# Seismograms +# +#----------------------------------------------------------- + +# interval in time steps for writing of seismograms +NTSTEP_BETWEEN_OUTPUT_SEISMOS = 10000 + +# set to n to reduce the sampling rate of output seismograms by a factor of n +# defaults to 1, which means no down-sampling +NTSTEP_BETWEEN_OUTPUT_SAMPLE = 1 + +# decide if we save displacement, velocity, acceleration and/or pressure in forward runs (they can be set to true simultaneously) +# currently pressure seismograms are implemented in acoustic (i.e. fluid) elements only +SAVE_SEISMOGRAMS_DISPLACEMENT = .true. +SAVE_SEISMOGRAMS_VELOCITY = .false. +SAVE_SEISMOGRAMS_ACCELERATION = .false. +SAVE_SEISMOGRAMS_PRESSURE = .false. # currently implemented in acoustic (i.e. fluid) elements only + +# option to save strain seismograms +# this option is useful for strain Green's tensor +SAVE_SEISMOGRAMS_STRAIN = .false. + +# save seismograms also when running the adjoint runs for an inverse problem +# (usually they are unused and not very meaningful, leave this off in almost all cases) +SAVE_SEISMOGRAMS_IN_ADJOINT_RUN = .true. + +# save seismograms in binary or ASCII format (binary is smaller but may not be portable between machines) +USE_BINARY_FOR_SEISMOGRAMS = .false. + +# output seismograms in Seismic Unix format (binary with 240-byte-headers) +SU_FORMAT = .false. + +# output seismograms in ASDF (requires asdf-library) +ASDF_FORMAT = .false. + +# output seismograms in HDF5 (requires hdf5-library and WRITE_SEISMOGRAMS_BY_MAIN) +HDF5_FORMAT = .false. + +# decide if main process writes all the seismograms or if all processes do it in parallel +WRITE_SEISMOGRAMS_BY_MAIN = .false. + +# save all seismograms in one large combined file instead of one file per seismogram +# to avoid overloading shared non-local file systems such as LUSTRE or GPFS for instance +SAVE_ALL_SEISMOS_IN_ONE_FILE = .false. + +# use a trick to increase accuracy of pressure seismograms in fluid (acoustic) elements: +# use the second derivative of the source for the source time function instead of the source itself, +# and then record -potential_acoustic() as pressure seismograms instead of -potential_dot_dot_acoustic(); +# this is mathematically equivalent, but numerically significantly more accurate because in the explicit +# Newmark time scheme acceleration is accurate at zeroth order while displacement is accurate at second order, +# thus in fluid elements potential_dot_dot_acoustic() is accurate at zeroth order while potential_acoustic() +# is accurate at second order and thus contains significantly less numerical noise. +USE_TRICK_FOR_BETTER_PRESSURE = .false. + +#----------------------------------------------------------- +# +# Energy calculation +# +#----------------------------------------------------------- + +# to plot energy curves, for instance to monitor how CPML absorbing layers behave; +# should be turned OFF in most cases because a bit expensive +OUTPUT_ENERGY = .false. +# every how many time steps we compute energy (which is a bit expensive to compute) +NTSTEP_BETWEEN_OUTPUT_ENERGY = 10 + +#----------------------------------------------------------- +# +# Adjoint kernel outputs +# +#----------------------------------------------------------- + +# interval in time steps for reading adjoint traces +# 0 = read the whole adjoint sources at start time +NTSTEP_BETWEEN_READ_ADJSRC = 0 + +# read adjoint sources using ASDF (requires asdf-library) +READ_ADJSRC_ASDF = .false. + +# this parameter must be set to .true. to compute anisotropic kernels +# in crust and mantle (related to the 21 Cij in geographical coordinates) +# default is .false. to compute isotropic kernels (related to alpha and beta) +ANISOTROPIC_KL = .false. + +# compute transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho) +# rather than fully anisotropic kernels in case ANISOTROPIC_KL is set to .true. +SAVE_TRANSVERSE_KL = .false. + +# this parameter must be set to .true. to compute anisotropic kernels for +# cost function using velocity observable rather than displacement +ANISOTROPIC_VELOCITY_KL = .false. + +# outputs approximate Hessian for preconditioning +APPROXIMATE_HESS_KL = .false. + +# save Moho mesh and compute Moho boundary kernels +SAVE_MOHO_MESH = .false. + +#----------------------------------------------------------- +# +# Coupling with an injection technique (DSM, AxiSEM, or FK) +# +#----------------------------------------------------------- +COUPLE_WITH_INJECTION_TECHNIQUE = .false. +INJECTION_TECHNIQUE_TYPE = 3 # 1 = DSM, 2 = AxiSEM, 3 = FK +MESH_A_CHUNK_OF_THE_EARTH = .false. +TRACTION_PATH = DATA/AxiSEM_tractions/3/ +FKMODEL_FILE = FKmodel +RECIPROCITY_AND_KH_INTEGRAL = .false. # does not work yet + +#----------------------------------------------------------- +# +# Run modes +# +#----------------------------------------------------------- + +# Simultaneous runs +# added the ability to run several calculations (several earthquakes) +# in an embarrassingly-parallel fashion from within the same run; +# this can be useful when using a very large supercomputer to compute +# many earthquakes in a catalog, in which case it can be better from +# a batch job submission point of view to start fewer and much larger jobs, +# each of them computing several earthquakes in parallel. +# To turn that option on, set parameter NUMBER_OF_SIMULTANEOUS_RUNS to a value greater than 1. +# To implement that, we create NUMBER_OF_SIMULTANEOUS_RUNS MPI sub-communicators, +# each of them being labeled "my_local_mpi_comm_world", and we use them +# in all the routines in "src/shared/parallel.f90", except in MPI_ABORT() because in that case +# we need to kill the entire run. +# When that option is on, of course the number of processor cores used to start +# the code in the batch system must be a multiple of NUMBER_OF_SIMULTANEOUS_RUNS, +# all the individual runs must use the same number of processor cores, +# which as usual is NPROC in the Par_file, +# and thus the total number of processor cores to request from the batch system +# should be NUMBER_OF_SIMULTANEOUS_RUNS * NPROC. +# All the runs to perform must be placed in directories called run0001, run0002, run0003 and so on +# (with exactly four digits). +# +# Imagine you have 10 independent calculations to do, each of them on 100 cores; you have three options: +# +# 1/ submit 10 jobs to the batch system +# +# 2/ submit a single job on 1000 cores to the batch, and in that script create a sub-array of jobs to start 10 jobs, +# each running on 100 cores (see e.g. http://www.schedmd.com/slurmdocs/job_array.html ) +# +# 3/ submit a single job on 1000 cores to the batch, start SPECFEM3D on 1000 cores, create 10 sub-communicators, +# cd into one of 10 subdirectories (called e.g. run0001, run0002,... run0010) depending on the sub-communicator +# your MPI rank belongs to, and run normally on 100 cores using that sub-communicator. +# +# The option below implements 3/. +# +NUMBER_OF_SIMULTANEOUS_RUNS = 1 + +# if we perform simultaneous runs in parallel, if only the source and receivers vary between these runs +# but not the mesh nor the model (velocity and density) then we can also read the mesh and model files +# from a single run in the beginning and broadcast them to all the others; for a large number of simultaneous +# runs for instance when solving inverse problems iteratively this can DRASTICALLY reduce I/Os to disk in the solver +# (by a factor equal to NUMBER_OF_SIMULTANEOUS_RUNS), and reducing I/Os is crucial in the case of huge runs. +# Thus, always set this option to .true. if the mesh and the model are the same for all simultaneous runs. +# In that case there is no need to duplicate the mesh and model file database (the content of the DATABASES_MPI +# directories) in each of the run0001, run0002,... directories, it is sufficient to have one in run0001 +# and the code will broadcast it to the others) +BROADCAST_SAME_MESH_AND_MODEL = .true. + +#----------------------------------------------------------- + +# set to true to use GPUs +GPU_MODE = .false. + +# ADIOS Options for I/Os +ADIOS_ENABLED = .false. +ADIOS_FOR_DATABASES = .false. +ADIOS_FOR_MESH = .false. +ADIOS_FOR_FORWARD_ARRAYS = .false. +ADIOS_FOR_KERNELS = .false. +ADIOS_FOR_UNDO_ATTENUATION = .false. + +# HDF5 Database I/O +# (note the flags for HDF5 and ADIOS are mutually exclusive, only one can be used) +HDF5_ENABLED = .false. +HDF5_FOR_MOVIES = .false. +HDF5_IO_NODES = 0 # HDF5 IO server with number of IO dedicated procs + diff --git a/DATA/STATIONS b/DATA/STATIONS deleted file mode 120000 index 7c5e8b5d6..000000000 --- a/DATA/STATIONS +++ /dev/null @@ -1 +0,0 @@ -../EXAMPLES/applications/homogeneous_halfspace/DATA/STATIONS \ No newline at end of file diff --git a/DATA/STATIONS b/DATA/STATIONS new file mode 100644 index 000000000..1ee47ff70 --- /dev/null +++ b/DATA/STATIONS @@ -0,0 +1,4 @@ +X20 DB 67000.00 22732.14 0.0 0.0 +X30 DB 67000.00 34696.43 0.0 0.0 +X40 DB 67000.00 46660.71 0.0 0.0 +X50 DB 67000.00 58625.00 0.0 0.0 diff --git a/doc/USER_MANUAL/05_running_the_solver.tex b/doc/USER_MANUAL/05_running_the_solver.tex index ef3e70409..641160343 100644 --- a/doc/USER_MANUAL/05_running_the_solver.tex +++ b/doc/USER_MANUAL/05_running_the_solver.tex @@ -193,7 +193,7 @@ \chapter{Running the Solver \texttt{xspecfem3D}}\label{cha:Running-the-Solver} \item \texttt{latorUTM:} Set the latitude or UTM $x$ coordinate. \item \texttt{longorUTM:} Set the longitude or UTM $y$ coordinate. \item \texttt{depth:} Set the depth of the source (in km). -\item \texttt{source time function:} Set the type of source-time function: 0 = Gaussian function, 1 = Ricker function, 2 = Heaviside (step) function, 3 = monochromatic function, 4 = Gaussian function as defined in \citet{Meschede2011}. +\item \texttt{source time function:} Set the type of source-time function: 0 = Gaussian function, 1 = Ricker function, 2 = Heaviside (step) function, 3 = monochromatic function, 4 = Gaussian function as defined in \citet{Meschede2011}, 5 = Brune function, and 6 = Smoothed Brune function. For the Brune and Smoothed Brune source time functions \texttt{hdurorf0} is the rise time. The Brune and Smoothed Brune functions are currently implemented only for the viscoelastic simulations. When {\texttt{USE\_RICKER\_TIME\_FUNCTION}} is turned on in the main parameter file \texttt{DATA/Par\_file}, it will override this source time function type selection and always use a Ricker wavelet. Note that we use the standard definition of a Ricker, for a dominant frequency $f_0$: diff --git a/m4 b/m4 index e490e14fb..230983ede 160000 --- a/m4 +++ b/m4 @@ -1 +1 @@ -Subproject commit e490e14fb13595428d39055304bcf0ee7ab94806 +Subproject commit 230983edee3f18d9a5e975608b6374751d024f75 diff --git a/src/specfem3D/comp_source_time_function.f90 b/src/specfem3D/comp_source_time_function.f90 index 822eb3b38..4062454e2 100644 --- a/src/specfem3D/comp_source_time_function.f90 +++ b/src/specfem3D/comp_source_time_function.f90 @@ -263,6 +263,58 @@ end function comp_source_time_function_d2rck !------------------------------------------------------------------------------------------------- ! + double precision function comp_source_time_function_brune(t,f0) + + use constants, only: PI + + implicit none + + double precision, intent(in) :: t,f0 + + ! local variables + double precision :: ft + + ! Brune source-time function + if(t .lt. 0.d0)then + comp_source_time_function_brune = 0.d0 + else + ft = f0*t + comp_source_time_function_brune = 1.d0 - exp( -ft ) * (1.0d0+ft) + endif + + end function comp_source_time_function_brune + +! +!------------------------------------------------------------------------------------------------- +! + + double precision function comp_source_time_function_smooth_brune(t,f0) + + implicit none + + double precision, intent(in) :: t,f0 + + ! local variables + double precision,parameter :: tau0=2.31d0 + double precision :: f,ft + + ! Brune source-time function + ft = f0*t + if (ft .lt. 0.d0) then + comp_source_time_function_smooth_brune = 0.d0 + elseif (ft .ge. 0.d0 .and. ft .lt. tau0)then + comp_source_time_function_smooth_brune = 1.d0 - exp(-ft)*( 1.0d0 + ft + & + 0.5d0*ft**2 - (1.5d0*ft**3)/tau0 + (1.5d0*ft**4)/(tau0**2) - & + (0.5d0*ft**5)/(tau0**3) ) + else ! (ft .gt. tau0)then + comp_source_time_function_smooth_brune = 1.d0 - exp( -ft ) * (1.0d0+ft) + endif + + end function comp_source_time_function_smooth_brune + +! +!------------------------------------------------------------------------------------------------- +! double precision function comp_source_time_function_mono(t,f0) diff --git a/src/specfem3D/compute_add_sources_viscoelastic.F90 b/src/specfem3D/compute_add_sources_viscoelastic.F90 index 4311eaebf..176fa98c0 100644 --- a/src/specfem3D/compute_add_sources_viscoelastic.F90 +++ b/src/specfem3D/compute_add_sources_viscoelastic.F90 @@ -838,7 +838,7 @@ double precision function get_stf_viscoelastic(time_source_dble,isource,it_tmp_e ! returns source time function value for specified time - use constants, only: USE_MONOCHROMATIC_CMT_SOURCE + use constants, only: PI,USE_MONOCHROMATIC_CMT_SOURCE use specfem_par, only: USE_FORCE_POINT_SOURCE,USE_RICKER_TIME_FUNCTION, & hdur,hdur_Gaussian,force_stf @@ -857,6 +857,7 @@ double precision function get_stf_viscoelastic(time_source_dble,isource,it_tmp_e double precision, external :: comp_source_time_function,comp_source_time_function_rickr, & comp_source_time_function_gauss,comp_source_time_function_gauss_2, & + comp_source_time_function_brune,comp_source_time_function_smooth_brune, & comp_source_time_function_mono,comp_source_time_function_ext ! external source time function @@ -890,6 +891,18 @@ double precision function get_stf_viscoelastic(time_source_dble,isource,it_tmp_e case (4) ! Gaussian source time function by Meschede et al. (2011) stf = comp_source_time_function_gauss_2(time_source_dble,hdur(isource)) + case (5) + ! Brune source time function + ! hdur is the rise time + ! Frequency parameter: + f0=2.0d0*PI/hdur(isource) + stf = comp_source_time_function_brune(time_source_dble,f0) + case (6) + ! Smoothed Brune source time function + ! hdur is the rise time + ! Frequency parameter: + f0=2.0d0*PI/hdur(isource) + stf = comp_source_time_function_smooth_brune(time_source_dble,f0) case default stop 'unsupported force_stf value!' end select diff --git a/src/specfem3D/get_force.f90 b/src/specfem3D/get_force.f90 index 49de3eb0a..d9a1fc07e 100644 --- a/src/specfem3D/get_force.f90 +++ b/src/specfem3D/get_force.f90 @@ -250,6 +250,20 @@ subroutine get_force(FORCESOLUTION,tshift_src,hdur,lat,long,depth,NSOURCES, & ! null half-duration indicates a Dirac ! replace with very short Gaussian function if (hdur(isource) < 5. * DT ) hdur(isource) = 5. * DT + case (5) + ! Brune source time function + ! half-duration is the rise time + ! (see constants.h: TINYVAL = 1.d-9 ) + if (hdur(isource) < TINYVAL ) then + stop 'Error set force period, make sure all forces have a non-zero rise time' + endif + case (6) + ! Smoothed Brune source time function + ! half-duration is the rise time + ! (see constants.h: TINYVAL = 1.d-9 ) + if (hdur(isource) < TINYVAL ) then + stop 'Error set force period, make sure all forces have a non-zero rise time' + endif case default stop 'unsupported source time function type (force_stf) value!' end select diff --git a/src/specfem3D/locate_source.F90 b/src/specfem3D/locate_source.F90 index f664d6a66..107f9b994 100644 --- a/src/specfem3D/locate_source.F90 +++ b/src/specfem3D/locate_source.F90 @@ -505,6 +505,18 @@ subroutine locate_source() ! Gaussian by Meschede et al. (2011) write(IMAIN,*) ' using Gaussian source time function by Meschede et al. (2011), eq.(2)' write(IMAIN,*) ' tau: ',hdur(isource),' seconds' + case (5) + ! Brune + write(IMAIN,*) ' using Brune source time function' + ! prints rise time for point forces + write(IMAIN,*) + write(IMAIN,*) ' rise time: ',hdur(isource),' seconds' + case (6) + ! Smoothed Brune + write(IMAIN,*) ' using Smoothed Brune source time function' + ! prints rise time for point forces + write(IMAIN,*) + write(IMAIN,*) ' rise time: ',hdur(isource),' seconds' case default stop 'unsupported force_stf value!' end select diff --git a/src/specfem3D/setup_sources_receivers.f90 b/src/specfem3D/setup_sources_receivers.f90 index 9aa75a804..3aafddbbc 100644 --- a/src/specfem3D/setup_sources_receivers.f90 +++ b/src/specfem3D/setup_sources_receivers.f90 @@ -461,6 +461,14 @@ subroutine setup_stf_constants() case (4) ! Gaussian source time function by Meschede et al. (2011) t0 = min(t0,1.5d0 * (tshift_src(isource) - hdur(isource))) + case (5) + ! Brune + ! This needs to be CHECKED!!! + t0 = min(t0,1.5d0 * (tshift_src(isource) - hdur(isource))) + case (6) + ! Smotthed Brune + ! This needs to be CHECKED!!! + t0 = min(t0,1.5d0 * (tshift_src(isource) - hdur(isource))) case default stop 'unsupported force_stf value!' end select