[cig-commits] r8471 - in seismo/2D/SPECFEM2D/trunk: . DATA
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:48:50 PST 2007
Author: walter
Date: 2007-12-07 15:48:49 -0800 (Fri, 07 Dec 2007)
New Revision: 8471
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_CercleMod50_50_10Hz_Paul
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_acoustic
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_elastic
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Thomas_Lacq
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA
seismo/2D/SPECFEM2D/trunk/DATA/STATIONS
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/constants.h
seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90
seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.f90
Log:
added Heaviside source to SPECFEM2D
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2007-12-07 23:48:49 UTC (rev 8471)
@@ -5,73 +5,74 @@
#-------------------------------->
# title of job, and file that contains interface data
-title = Test InterfaceCercleMod50_50
-interfacesfile = interfaceCercleMod50_50_Paul.dat
+title = Lacq Thomas
+interfacesfile = interfaces_thomas_Lacq.dat
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
-xmax = 4000.d0 # abscissa of right side of the model
-nx = 300 # number of elements along X
+xmax = 25000.d0 # abscissa of right side of the model
+nx = 250 # number of elements along X
ngnod = 9 # noeuds de controle pour blocs (4 ou 9)
initialfield = .false. # use a plane wave as source or not
-readmodel = .false. # read external earth model or not
-ELASTIC = .false. # elastic or acoustic simulation
+read_external_model = .false. # read external earth model or not
+ELASTIC = .true. # elastic or acoustic simulation
TURN_ANISOTROPY_ON = .false. # turn anisotropy on or off
TURN_ATTENUATION_ON = .false. # turn attenuation on or off
# absorbing boundaries parameters
-absorbhaut = .true. # absorbing boundary active or not
+absorbhaut = .false. # absorbing boundary active or not
absorbbas = .true.
absorbgauche = .true.
absorbdroite = .true.
# time step parameters
-nt = 6600 # nb total de pas de temps
-dt = 2.d-4 # valeur du pas de temps
+nt = 10000 # nb total de pas de temps
+dt = 1.d-3 # valeur du pas de temps
# source parameters
source_surf = .false. # source dans le volume ou a la surface
-xs = 2000. # source location x in meters
-zs = 2850. # source location z in meters
+xs = 12500. # source location x in meters
+zs = 10000. # source location z in meters
source_type = 1 # force = 1 or explosion = 2
-time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
+time_function_type = 5 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
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)
Mxz = 0. # Mxz component (for a moment tensor source only)
-factor = 1.d-6 # amplification factor
+factor = 0.1 # amplification factor
# receiver line parameters
-enreg_surf = .false. # enregistrement volume ou surface
+enreg_surf = .false. # enregistrement volume ou surface
sismostype = 2 # record 1=displ 2=veloc 3=accel 4=pressure
-nrec = 101 # number of receivers
-xdeb = 500. # first receiver x in meters
-zdeb = 2750. # first receiver z in meters
-xfin = 3500. # last receiver x in meters
-zfin = 2750. # last receiver z in meters
+nrec = 4 # number of receivers
+xdeb = 5000. # first receiver x in meters
+zdeb = 15000. # first receiver z in meters
+xfin = 20000. # last receiver x in meters
+zfin = 15000. # last receiver z in meters
anglerec = 0.d0 # angle to rotate components at receivers
# display parameters
-itaff = 200 # display frequency in time steps
-output_postscript_snapshot = .false. # output Postscript image of the results
+itaff = 250 # display frequency in time steps
+output_postscript_snapshot = .true. # output Postscript image of the results
output_PNM_image = .false. # output PNM image of the results
-vecttype = 2 # display 1=displ 2=veloc 3=accel
-cutvect = 1. # amplitude min en % pour vector plots
-meshvect = .true. # display mesh on vector plots or not
-modelvect = .false. # display velocity model on vector plots
-boundvect = .false. # display boundary conditions on plots
-interpol = .false. # interpolation of the display or not
-pointsdisp = 6 # points for interpolation of display (set to 1 for lower-left corner only)
-subsamp = 1 # subsampling of color snapshots
-sizemax_arrows = 1.d0 # maximum size of arrows on vector plots in cm
+vecttype = 1 # display 1=displ 2=veloc 3=accel
+cutvect = 1.d0 # amplitude min en % pour vector plots
+meshvect = .false. # display mesh on vector plots or not
+modelvect = .true. # display velocity model on vector plots
+boundvect = .true. # display boundary conditions on plots
+interpol = .true. # interpolation of the display or not
+pointsdisp = 1 # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp = 2 # subsampling of color snapshots
+sizemax_arrows = 0.4d0 # maximum size of arrows on vector plots in cm
gnuplot = .false. # generate a GNUPLOT file for the grid
outputgrid = .false. # save the grid in a text file or not
# velocity and density model (nx,nz)
-nbmodels = 2 # nb de modeles differents (0,rho,vp,vs,0,0)
-1 0 2700.d0 3000.d0 0.d0 0 0
-2 0 1800.d0 2000.d0 0.d0 0 0
-nbzone = 2 # nb of zones and model number for each
-1 300 1 150 1
-1 300 151 300 2
+nbmodels = 1 # nb de modeles differents (0,rho,vp,vs,0,0)
+1 0 2700.d0 3000.d0 1732.051d0 0 0
+#2 0 1800.d0 2000.d0 1100.845d0 0 0
+#3 0 2200.d0 2500.d0 1443.375d0 0 0
+nbzone = 1 # nb of zones and model number for each
+1 250 1 150 1
+#1 80 41 60 3
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_CercleMod50_50_10Hz_Paul
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_CercleMod50_50_10Hz_Paul 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_CercleMod50_50_10Hz_Paul 2007-12-07 23:48:49 UTC (rev 8471)
@@ -34,7 +34,7 @@
xs = 2000. # source location x in meters
zs = 2850. # source location z in meters
source_type = 1 # force = 1 or explosion = 2
-time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
+time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
angleforce = 0. # angle of the source (for a force only)
Mxx = 1. # Mxx component (for a moment tensor source only)
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_acoustic
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_acoustic 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_acoustic 2007-12-07 23:48:49 UTC (rev 8471)
@@ -34,7 +34,7 @@
xs = 2000. # source location x in meters
zs = 1200. # source location z in meters
source_type = 1 # force = 1 or explosion = 2
-time_function_type = 2 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
+time_function_type = 2 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
angleforce = 0. # angle of the source (for a force only)
Mxx = 1. # Mxx component (for a moment tensor source only)
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_elastic
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_elastic 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_elastic 2007-12-07 23:48:49 UTC (rev 8471)
@@ -34,7 +34,7 @@
xs = 2000. # source location x in meters
zs = 1200. # source location z in meters
source_type = 1 # force = 1 or explosion = 2
-time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
+time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
angleforce = 0. # angle of the source (for a force only)
Mxx = 1. # Mxx component (for a moment tensor source only)
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Thomas_Lacq
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Thomas_Lacq 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Thomas_Lacq 2007-12-07 23:48:49 UTC (rev 8471)
@@ -34,7 +34,7 @@
xs = 12500. # source location x in meters
zs = 10000. # source location z in meters
source_type = 1 # force = 1 or explosion = 2
-time_function_type = 3 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
+time_function_type = 5 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
angleforce = 0. # angle of the source (for a force only)
Mxx = 1. # Mxx component (for a moment tensor source only)
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA 2007-12-07 23:48:49 UTC (rev 8471)
@@ -39,7 +39,7 @@
xs = 2000. # source location x in meters
zs = 1500. # source location z in meters
source_type = 1 # force = 1 or explosion = 2
-time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
+time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
angleforce = 0. # angle of the source (for a force only)
Mxx = 1. # Mxx component (for a moment tensor source only)
Modified: seismo/2D/SPECFEM2D/trunk/DATA/STATIONS
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/STATIONS 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/DATA/STATIONS 2007-12-07 23:48:49 UTC (rev 8471)
@@ -1,102 +1,5 @@
- 101
-S001 AA 500.0000000 2750.0000000 0.0 0.0
-S002 AA 530.0000000 2750.0000000 0.0 0.0
-S003 AA 560.0000000 2750.0000000 0.0 0.0
-S004 AA 590.0000000 2750.0000000 0.0 0.0
-S005 AA 620.0000000 2750.0000000 0.0 0.0
-S006 AA 650.0000000 2750.0000000 0.0 0.0
-S007 AA 680.0000000 2750.0000000 0.0 0.0
-S008 AA 710.0000000 2750.0000000 0.0 0.0
-S009 AA 740.0000000 2750.0000000 0.0 0.0
-S010 AA 770.0000000 2750.0000000 0.0 0.0
-S011 AA 800.0000000 2750.0000000 0.0 0.0
-S012 AA 830.0000000 2750.0000000 0.0 0.0
-S013 AA 860.0000000 2750.0000000 0.0 0.0
-S014 AA 890.0000000 2750.0000000 0.0 0.0
-S015 AA 920.0000000 2750.0000000 0.0 0.0
-S016 AA 950.0000000 2750.0000000 0.0 0.0
-S017 AA 980.0000000 2750.0000000 0.0 0.0
-S018 AA 1010.0000000 2750.0000000 0.0 0.0
-S019 AA 1040.0000000 2750.0000000 0.0 0.0
-S020 AA 1070.0000000 2750.0000000 0.0 0.0
-S021 AA 1100.0000000 2750.0000000 0.0 0.0
-S022 AA 1130.0000000 2750.0000000 0.0 0.0
-S023 AA 1160.0000000 2750.0000000 0.0 0.0
-S024 AA 1190.0000000 2750.0000000 0.0 0.0
-S025 AA 1220.0000000 2750.0000000 0.0 0.0
-S026 AA 1250.0000000 2750.0000000 0.0 0.0
-S027 AA 1280.0000000 2750.0000000 0.0 0.0
-S028 AA 1310.0000000 2750.0000000 0.0 0.0
-S029 AA 1340.0000000 2750.0000000 0.0 0.0
-S030 AA 1370.0000000 2750.0000000 0.0 0.0
-S031 AA 1400.0000000 2750.0000000 0.0 0.0
-S032 AA 1430.0000000 2750.0000000 0.0 0.0
-S033 AA 1460.0000000 2750.0000000 0.0 0.0
-S034 AA 1490.0000000 2750.0000000 0.0 0.0
-S035 AA 1520.0000000 2750.0000000 0.0 0.0
-S036 AA 1550.0000000 2750.0000000 0.0 0.0
-S037 AA 1580.0000000 2750.0000000 0.0 0.0
-S038 AA 1610.0000000 2750.0000000 0.0 0.0
-S039 AA 1640.0000000 2750.0000000 0.0 0.0
-S040 AA 1670.0000000 2750.0000000 0.0 0.0
-S041 AA 1700.0000000 2750.0000000 0.0 0.0
-S042 AA 1730.0000000 2750.0000000 0.0 0.0
-S043 AA 1760.0000000 2750.0000000 0.0 0.0
-S044 AA 1790.0000000 2750.0000000 0.0 0.0
-S045 AA 1820.0000000 2750.0000000 0.0 0.0
-S046 AA 1850.0000000 2750.0000000 0.0 0.0
-S047 AA 1880.0000000 2750.0000000 0.0 0.0
-S048 AA 1910.0000000 2750.0000000 0.0 0.0
-S049 AA 1940.0000000 2750.0000000 0.0 0.0
-S050 AA 1970.0000000 2750.0000000 0.0 0.0
-S051 AA 2000.0000000 2750.0000000 0.0 0.0
-S052 AA 2030.0000000 2750.0000000 0.0 0.0
-S053 AA 2060.0000000 2750.0000000 0.0 0.0
-S054 AA 2090.0000000 2750.0000000 0.0 0.0
-S055 AA 2120.0000000 2750.0000000 0.0 0.0
-S056 AA 2150.0000000 2750.0000000 0.0 0.0
-S057 AA 2180.0000000 2750.0000000 0.0 0.0
-S058 AA 2210.0000000 2750.0000000 0.0 0.0
-S059 AA 2240.0000000 2750.0000000 0.0 0.0
-S060 AA 2270.0000000 2750.0000000 0.0 0.0
-S061 AA 2300.0000000 2750.0000000 0.0 0.0
-S062 AA 2330.0000000 2750.0000000 0.0 0.0
-S063 AA 2360.0000000 2750.0000000 0.0 0.0
-S064 AA 2390.0000000 2750.0000000 0.0 0.0
-S065 AA 2420.0000000 2750.0000000 0.0 0.0
-S066 AA 2450.0000000 2750.0000000 0.0 0.0
-S067 AA 2480.0000000 2750.0000000 0.0 0.0
-S068 AA 2510.0000000 2750.0000000 0.0 0.0
-S069 AA 2540.0000000 2750.0000000 0.0 0.0
-S070 AA 2570.0000000 2750.0000000 0.0 0.0
-S071 AA 2600.0000000 2750.0000000 0.0 0.0
-S072 AA 2630.0000000 2750.0000000 0.0 0.0
-S073 AA 2660.0000000 2750.0000000 0.0 0.0
-S074 AA 2690.0000000 2750.0000000 0.0 0.0
-S075 AA 2720.0000000 2750.0000000 0.0 0.0
-S076 AA 2750.0000000 2750.0000000 0.0 0.0
-S077 AA 2780.0000000 2750.0000000 0.0 0.0
-S078 AA 2810.0000000 2750.0000000 0.0 0.0
-S079 AA 2840.0000000 2750.0000000 0.0 0.0
-S080 AA 2870.0000000 2750.0000000 0.0 0.0
-S081 AA 2900.0000000 2750.0000000 0.0 0.0
-S082 AA 2930.0000000 2750.0000000 0.0 0.0
-S083 AA 2960.0000000 2750.0000000 0.0 0.0
-S084 AA 2990.0000000 2750.0000000 0.0 0.0
-S085 AA 3020.0000000 2750.0000000 0.0 0.0
-S086 AA 3050.0000000 2750.0000000 0.0 0.0
-S087 AA 3080.0000000 2750.0000000 0.0 0.0
-S088 AA 3110.0000000 2750.0000000 0.0 0.0
-S089 AA 3140.0000000 2750.0000000 0.0 0.0
-S090 AA 3170.0000000 2750.0000000 0.0 0.0
-S091 AA 3200.0000000 2750.0000000 0.0 0.0
-S092 AA 3230.0000000 2750.0000000 0.0 0.0
-S093 AA 3260.0000000 2750.0000000 0.0 0.0
-S094 AA 3290.0000000 2750.0000000 0.0 0.0
-S095 AA 3320.0000000 2750.0000000 0.0 0.0
-S096 AA 3350.0000000 2750.0000000 0.0 0.0
-S097 AA 3380.0000000 2750.0000000 0.0 0.0
-S098 AA 3410.0000000 2750.0000000 0.0 0.0
-S099 AA 3440.0000000 2750.0000000 0.0 0.0
-S100 AA 3470.0000000 2750.0000000 0.0 0.0
-S101 AA 3500.0000000 2750.0000000 0.0 0.0
+ 4
+S001 AA 5000.0000000 15000.0000000 0.0 0.0
+S002 AA 10000.0000000 15000.0000000 0.0 0.0
+S003 AA 15000.0000000 15000.0000000 0.0 0.0
+S004 AA 20000.0000000 15000.0000000 0.0 0.0
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-07 23:48:49 UTC (rev 8471)
@@ -39,7 +39,7 @@
$O/plotpost.o $O/locate_receivers.o $O/locate_source_force.o $O/compute_gradient_attenuation.o\
$O/specfem2D.o $O/write_seismograms.o $O/createnum_fast.o $O/createnum_slow.o\
$O/define_shape_functions.o $O/cree_image_PNM.o $O/compute_gradient_fluid.o\
- $O/recompute_jacobian.o $O/compute_arrays_source.o $O/locate_source_moment_tensor.o
+ $O/recompute_jacobian.o $O/compute_arrays_source.o $O/locate_source_moment_tensor.o $O/numerical_recipes.o
default: meshfem2D specfem2D convolve_source_timefunction
@@ -137,6 +137,9 @@
$O/cree_image_PNM.o: cree_image_PNM.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/cree_image_PNM.o cree_image_PNM.f90
+$O/numerical_recipes.o: numerical_recipes.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/numerical_recipes.o numerical_recipes.f90
+
$O/write_seismograms.o: write_seismograms.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/write_seismograms.o write_seismograms.f90
Modified: seismo/2D/SPECFEM2D/trunk/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/constants.h 2007-12-07 23:48:49 UTC (rev 8471)
@@ -55,6 +55,9 @@
! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 4
+! error function source decay rate for Heaviside
+ double precision, parameter :: SOURCE_DECAY_RATE = 2.628d0
+
! display non lineaire pour rehausser les faibles amplitudes sur les images PNM
double precision, parameter :: POWER_DISPLAY_PNM = 0.30d0
Modified: seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90 2007-12-07 23:48:49 UTC (rev 8471)
@@ -21,9 +21,6 @@
include "constants.h"
-! source decay rate
- double precision, parameter :: SOURCE_DECAY_RATE = 2.628d0
-
integer i,j,N_j
integer number_remove
double precision, dimension(:), allocatable :: time,sem,sem_fil
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2007-12-07 23:48:49 UTC (rev 8471)
@@ -201,10 +201,15 @@
call read_value_double_precision(IIN,IGNORE_JUNK,factor)
! if Dirac source time function, use a very thin Gaussian instead
- if(time_function_type == 4) f0 = 1.d0 / (5.d0 * dt)
+! if Heaviside source time function, use a very thin error function instead
+ if(time_function_type == 4 .or. time_function_type == 5) f0 = 1.d0 / (10.d0 * dt)
-! time delay of the source in seconds, use a 20 % security margin
- t0 = 1.20d0 / f0
+! time delay of the source in seconds, use a 20 % security margin (use 2 / f0 if error function)
+ if(time_function_type == 5) then
+ t0 = 2.0d0 / f0
+ else
+ t0 = 1.20d0 / f0
+ endif
print *
print *,'Source:'
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2005-02-23 16:40:10 UTC (rev 8470)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2007-12-07 23:48:49 UTC (rev 8471)
@@ -61,7 +61,7 @@
character(len=80) datlin
integer source_type,time_function_type
- double precision x_source,z_source,xi_source,gamma_source,Mxx,Mzz,Mxz,f0,t0,factor,angleforce
+ double precision x_source,z_source,xi_source,gamma_source,Mxx,Mzz,Mxz,f0,t0,factor,angleforce,hdur,hdur_gauss
double precision, dimension(NDIM,NGLLX,NGLLZ) :: sourcearray
double precision, dimension(:,:), allocatable :: coorg
@@ -903,11 +903,18 @@
else if(time_function_type == 3 .or. time_function_type == 4) then
source_time_function(it) = factor * exp(-a*(time-t0)**2)
+! Heaviside source time function (we use a very thin error function instead)
+ else if(time_function_type == 5) then
+ hdur = 1.d0 / f0
+ hdur_gauss = hdur * 5.d0 / 3.d0
+ source_time_function(it) = factor * 0.5d0*(1.0d0+erf(SOURCE_DECAY_RATE*(time-t0)/hdur_gauss))
+
else
stop 'unknown source time function'
endif
- write(55,*) sngl(time-t0),sngl(source_time_function(it))
+! output absolute time in third column, in case user wants to check it as well
+ write(55,*) sngl(time),sngl(source_time_function(it)),sngl(time-t0)
enddo
More information about the cig-commits
mailing list