[cig-commits] r20007 - seismo/1D/SPECFEM1D/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Apr 30 09:04:52 PDT 2012
Author: dkomati1
Date: 2012-04-30 09:04:52 -0700 (Mon, 30 Apr 2012)
New Revision: 20007
Modified:
seismo/1D/SPECFEM1D/trunk/Makefile
seismo/1D/SPECFEM1D/trunk/constants.h
seismo/1D/SPECFEM1D/trunk/diffusion.f90
seismo/1D/SPECFEM1D/trunk/plot_all_snapshots.gnu
seismo/1D/SPECFEM1D/trunk/source_time_function.f90
seismo/1D/SPECFEM1D/trunk/wave.f90
Log:
added a point source, as suggested by Martin van Driel
Modified: seismo/1D/SPECFEM1D/trunk/Makefile
===================================================================
--- seismo/1D/SPECFEM1D/trunk/Makefile 2012-04-30 12:31:24 UTC (rev 20006)
+++ seismo/1D/SPECFEM1D/trunk/Makefile 2012-04-30 16:04:52 UTC (rev 20007)
@@ -32,11 +32,11 @@
# Intel
#F90 = ifort
-#FLAGS =-O3 -xP -vec-report0 -e03 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe3 -ftz
+#FLAGS = -O1 -vec-report0 -std95 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check all -align sequence -assume byterecl -ftz -traceback -ftrapuv
# GNU gfortran
F90 = gfortran
-FLAGS = -std=gnu -fimplicit-none -frange-check -O2 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow
+FLAGS = -std=gnu -fimplicit-none -frange-check -O2 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -fbounds-check
wave: constants.h \
$O/wave.o \
Modified: seismo/1D/SPECFEM1D/trunk/constants.h
===================================================================
--- seismo/1D/SPECFEM1D/trunk/constants.h 2012-04-30 12:31:24 UTC (rev 20006)
+++ seismo/1D/SPECFEM1D/trunk/constants.h 2012-04-30 16:04:52 UTC (rev 20007)
@@ -25,14 +25,29 @@
!=====================================================================
! number of spectral elements
- integer, parameter :: NSPEC = 11
+ integer, parameter :: NSPEC = 250
! number of GLL points (polynomial degree plus one)
- integer, parameter :: NGLL = 7
+ integer, parameter :: NGLL = 5
+! Courant–Friedrichs–Lewy (CFL) stability value
+! see e.g. http://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition
+ double precision, parameter :: courant_CFL = 0.20d0
+
+! number of time steps to compute
+ integer, parameter :: NSTEP = 20000
+
+! model parameters (SI)
+ double precision, parameter :: LENGTH = 3.0d+03 ! m
+ double precision, parameter :: DENSITY = 2.5d+03 ! kg/m^3
+ double precision, parameter :: RIGIDITY = 3.0d+10 ! Pa
+
! number of global points
integer, parameter :: NGLOB = (NGLL-1)*NSPEC + 1
! for the Gauss-Lobatto-Legendre points and weights
double precision, parameter :: GAUSSALPHA = 0.d0,GAUSSBETA = 0.d0
+! pi
+ double precision, parameter :: PI = 3.141592653589793d0
+
Modified: seismo/1D/SPECFEM1D/trunk/diffusion.f90
===================================================================
--- seismo/1D/SPECFEM1D/trunk/diffusion.f90 2012-04-30 12:31:24 UTC (rev 20006)
+++ seismo/1D/SPECFEM1D/trunk/diffusion.f90 2012-04-30 16:04:52 UTC (rev 20007)
@@ -51,18 +51,10 @@
include "constants.h"
-! number of timesteps
- integer, parameter :: NSTEP = 30000
-
-! time step in seconds
- double precision, parameter :: DT = 100000000. ! s
-
! fixed boundary conditions
logical, parameter :: FIXED_BC = .true.
! model parameters (SI)
- double precision, parameter :: LENGTH = 3.0d+03 ! m
- double precision, parameter :: DENSITY = 2.5d+03 ! kg/m^3
double precision, parameter :: THERMALCONDUCTIVITY = 10.0d-01 ! cal/m/s/K
double precision, parameter :: HEATCAPACITY = 0.3d+03 ! cal/kg/K
@@ -103,7 +95,7 @@
! time marching
double precision deltat,deltatover2,deltatsqover2
- double precision dh,diffusivity,time_step
+ double precision dh,diffusivity
! end fluxes
double precision flux_1,flux_NGLOB
@@ -168,8 +160,9 @@
! estimate the time step
dh = LENGTH/dble(NGLOB-1)
diffusivity = THERMALCONDUCTIVITY/(HEATCAPACITY*DENSITY)
- time_step = 0.5*dh*dh/diffusivity
-! print *,'time step estimate: ',time_step,' seconds'
+!!!!!! deltat = 0.5*dh*dh/diffusivity
+ deltat = 800000.
+ print *,'time step estimate: ',deltat,' seconds'
if(FIXED_BC) then
! set up the temperatures at the ends
@@ -182,7 +175,6 @@
endif
! time marching parameters
- deltat = DT
deltatover2 = deltat/2.
deltatsqover2 = deltat*deltat/2.
@@ -259,7 +251,8 @@
temperature(:) = temperature(:) + deltatover2*grad_temperature(:)
! write out snapshots
- if(mod(it-1,1000) == 0) then
+ if(mod(it,100) == 0) then
+ print *,'time step ',it,' out of ',NSTEP
write(moviefile,"('snapshot',i5.5)") it
open(unit=10,file=moviefile,status='unknown')
do iglob = 1,NGLOB
Modified: seismo/1D/SPECFEM1D/trunk/plot_all_snapshots.gnu
===================================================================
--- seismo/1D/SPECFEM1D/trunk/plot_all_snapshots.gnu 2012-04-30 12:31:24 UTC (rev 20006)
+++ seismo/1D/SPECFEM1D/trunk/plot_all_snapshots.gnu 2012-04-30 16:04:52 UTC (rev 20007)
@@ -1,37 +1,404 @@
-# Gnuplot script to plot all the snapshots
-
set term x11
-
+#set term wxt
set xrange [0:3000]
-set yrange [-1.5:+1.5]
-
-plot "snapshot00001" w l lc 1
-pause -1 "hit any key..."
-
-plot "snapshot01001" w l lc 1
-pause -1 "hit any key..."
-
-plot "snapshot02001" w l lc 1
-pause -1 "hit any key..."
-
-plot "snapshot03001" w l lc 1
-pause -1 "hit any key..."
-
-plot "snapshot04001" w l lc 1
-pause -1 "hit any key..."
-
-plot "snapshot05001" w l lc 1
-pause -1 "hit any key..."
-
-plot "snapshot06001" w l lc 1
-pause -1 "hit any key..."
-
-plot "snapshot07001" w l lc 1
-pause -1 "hit any key..."
-
-plot "snapshot08001" w l lc 1
-pause -1 "hit any key..."
-
-plot "snapshot09001" w l lc 1
-pause -1 "hit any key..."
-
+#set yrange [-1.5:+1.5]
+plot "snapshot00100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot00200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot00300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot00400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot00500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot00600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot00700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot00800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot00900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot01000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot01100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot01200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot01300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot01400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot01500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot01600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot01700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot01800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot01900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot02000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot02100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot02200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot02300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot02400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot02500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot02600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot02700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot02800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot02900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot03000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot03100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot03200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot03300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot03400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot03500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot03600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot03700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot03800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot03900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot04000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot04100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot04200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot04300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot04400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot04500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot04600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot04700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot04800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot04900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot05000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot05100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot05200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot05300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot05400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot05500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot05600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot05700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot05800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot05900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot06000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot06100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot06200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot06300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot06400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot06500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot06600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot06700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot06800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot06900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot07000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot07100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot07200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot07300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot07400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot07500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot07600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot07700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot07800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot07900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot08000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot08100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot08200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot08300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot08400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot08500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot08600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot08700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot08800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot08900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot09000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot09100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot09200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot09300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot09400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot09500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot09600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot09700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot09800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot09900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot10000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot10100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot10200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot10300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot10400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot10500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot10600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot10700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot10800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot10900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot11000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot11100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot11200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot11300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot11400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot11500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot11600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot11700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot11800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot11900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot12000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot12100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot12200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot12300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot12400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot12500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot12600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot12700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot12800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot12900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot13000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot13100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot13200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot13300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot13400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot13500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot13600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot13700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot13800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot13900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot14000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot14100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot14200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot14300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot14400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot14500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot14600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot14700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot14800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot14900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot15000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot15100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot15200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot15300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot15400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot15500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot15600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot15700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot15800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot15900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot16000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot16100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot16200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot16300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot16400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot16500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot16600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot16700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot16800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot16900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot17000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot17100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot17200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot17300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot17400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot17500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot17600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot17700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot17800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot17900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot18000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot18100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot18200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot18300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot18400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot18500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot18600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot18700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot18800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot18900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot19000" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot19100" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot19200" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot19300" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot19400" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot19500" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot19600" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot19700" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot19800" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot19900" w l lc 1
+pause -1 'hit any key...'
+plot "snapshot20000" w l lc 1
+pause -1 'hit any key...'
Modified: seismo/1D/SPECFEM1D/trunk/source_time_function.f90
===================================================================
--- seismo/1D/SPECFEM1D/trunk/source_time_function.f90 2012-04-30 12:31:24 UTC (rev 20006)
+++ seismo/1D/SPECFEM1D/trunk/source_time_function.f90 2012-04-30 16:04:52 UTC (rev 20007)
@@ -39,10 +39,10 @@
alpha = decay_rate/hdur
! Gaussian
- source_time_function = alpha*exp(-alpha*alpha*t*t)/sqrt(PI)
+! source_time_function = alpha*exp(-alpha*alpha*t*t)/sqrt(PI)
! Ricker wavelet
-! source_time_function = -2.*(alpha**3)*t*exp(-alpha*alpha*t*t)/sqrt(PI)
+ source_time_function = -2.*(alpha**3)*t*exp(-alpha*alpha*t*t)/sqrt(PI)
end function source_time_function
Modified: seismo/1D/SPECFEM1D/trunk/wave.f90
===================================================================
--- seismo/1D/SPECFEM1D/trunk/wave.f90 2012-04-30 12:31:24 UTC (rev 20006)
+++ seismo/1D/SPECFEM1D/trunk/wave.f90 2012-04-30 16:04:52 UTC (rev 20007)
@@ -51,23 +51,12 @@
include "constants.h"
-! number of timesteps
- integer, parameter :: NSTEP = 10000
+! use a plane wave source or a point source (a force)
+ logical, parameter :: USE_PLANE_WAVE_SOURCE = .false.
-! time step in seconds
- double precision, parameter :: DT = 0.0002 ! s
-
! fixed boundary conditions
logical, parameter :: FIXED_BC = .true.
-! model parameters (SI)
- double precision, parameter :: LENGTH = 3.0d+03 ! m
- double precision, parameter :: DENSITY = 2.5d+03 ! kg/m^3
- double precision, parameter :: RIGIDITY = 3.0d+10 ! Pa
-
-! pi
- double precision, parameter :: PI = 3.141592653589793d0
-
integer ispec,i,j,iglob,it
! Gauss-Lobatto-Legendre points of integration
@@ -104,12 +93,12 @@
integer, dimension(NGLL,NSPEC) :: ibool
! time marching
- double precision dh,v,courant_CFL,time_step
+ double precision dh,v
double precision deltat,deltatover2,deltatsqover2
! source
- integer ispec_source,i_source
- double precision hdur,source_amp
+ integer ispec_source,i_source,iglob_source
+ double precision hdur,source_amp,stf
double precision, external :: source_time_function
! receiver
@@ -169,35 +158,35 @@
enddo
enddo
-! estimate the time step
+! estimate the time step in seconds
dh = LENGTH/dble(NGLOB-1)
v = dsqrt(RIGIDITY/DENSITY)
- courant_CFL = 0.2
- time_step = courant_CFL*dh/v
-! print *,'time step estimate: ',time_step,' seconds'
+ deltat = courant_CFL*dh/v
+ print *,'time step estimate: ',deltat,' seconds'
! set the source
- ispec_source = (NSPEC+1)/2
- i_source = (NGLL+1)/2
- hdur = 50.*DT
- source_amp = 0.001
+ ispec_source = NSPEC/2
+ i_source = 2
+ hdur = 500.*deltat
+ source_amp = 10000000.d0
! set the receiver
ireceiver = 2*NGLL-2
! time marching parameters
- deltat = DT
deltatover2 = deltat/2.
deltatsqover2 = deltat*deltat/2.
-! initialize
+! initialize the fields to zero
displ(:) = 0.
veloc(:) = 0.
accel(:) = 0.
- do iglob = 1,NGLOB
- displ(iglob) = dsin(PI*x(iglob)/LENGTH)
- enddo
+ if(USE_PLANE_WAVE_SOURCE) then
+ do iglob = 1,NGLOB
+ displ(iglob) = dsin(PI*x(iglob)/LENGTH)
+ enddo
+ endif
do it = 1,NSTEP
@@ -241,9 +230,11 @@
enddo ! end loop over all spectral elements
! add source at global level
-! iglob_source = ibool(ispec_source,i_source)
-! stf = source_time_function(dble(it-1)*DT-hdur,hdur)
-! accel(iglob_source) = accel(iglob_source) + stf*source_amp
+ if(.not. USE_PLANE_WAVE_SOURCE) then
+ iglob_source = ibool(i_source,ispec_source)
+ stf = source_time_function(dble(it-1)*deltat-hdur,hdur)
+ accel(iglob_source) = accel(iglob_source) + stf*source_amp
+ endif
! fixed boundary conditions at global level
if(FIXED_BC) then
@@ -258,7 +249,8 @@
veloc(:) = veloc(:) + deltatover2*accel(:)
! write out snapshots
- if(mod(it-1,1000) == 0) then
+ if(mod(it,100) == 0) then
+ print *,'time step ',it,' out of ',NSTEP
write(moviefile,"('snapshot',i5.5)") it
open(unit=10,file=moviefile,status='unknown')
do iglob = 1,NGLOB
@@ -273,7 +265,7 @@
open(unit=12,file='seismogram',status='unknown')
do it = 1,NSTEP
- write(12,*) sngl(dble(it-1)*DT-hdur),sngl(seismogram(it))
+ write(12,*) sngl(dble(it-1)*deltat-hdur),sngl(seismogram(it))
enddo
close(12)
More information about the CIG-COMMITS
mailing list