[cig-commits] [commit] devel, master: adds noise example to EXAMPLES/ directory (4634584)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Jun 18 15:19:13 PDT 2014
Repository : https://github.com/geodynamics/specfem2d
On branches: devel,master
Link : https://github.com/geodynamics/specfem2d/compare/fc67e6fd7ad890705b2b72b4b3c509accb22249e...e9ca46c40131588d89d7b0883250bc6584ce6b4c
>---------------------------------------------------------------
commit 46345846ebec960cafda71caf4669d50401354bc
Author: Ryan Modrak <rmodrak at princeton.edu>
Date: Fri Jul 29 19:34:21 2011 +0000
adds noise example to EXAMPLES/ directory
>---------------------------------------------------------------
46345846ebec960cafda71caf4669d50401354bc
.../Par_file_noise_1 | 41 +-
.../Par_file_noise_2 | 43 +--
.../Par_file_noise_3 | 43 +--
.../SOURCE_noise | 16 +-
noise_uniform/STATIONS_target_noise | 3 +
noise_uniform/adj_cc.f90 | 415 +++++++++++++++++++++
noise_uniform/process.sh | 67 ++++
.../uniform.dat | 11 +-
8 files changed, 555 insertions(+), 84 deletions(-)
diff --git a/M2_UPPA/Par_file_M2_UPPA b/noise_uniform/Par_file_noise_1
similarity index 83%
copy from M2_UPPA/Par_file_M2_UPPA
copy to noise_uniform/Par_file_noise_1
index d76a229..266fa90 100644
--- a/M2_UPPA/Par_file_M2_UPPA
+++ b/noise_uniform/Par_file_noise_1
@@ -1,8 +1,9 @@
# title of job
-title = Test for M2 UPPA
+title = Noise_2D
# forward or adjoint simulation
SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
+NOISE_TOMOGRAPHY = 1 # 0 = earthquake simulation, 1/2/3 = noise simulation
SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
# parameters concerning partitioning
@@ -18,10 +19,10 @@ TURN_ATTENUATION_ON = .false. # turn attenuation on or off fo
TURN_VISCATTENUATION_ON = .false. # turn viscous attenuation on or off
Q0 = 1 # quality factor for viscous attenuation
freq0 = 10 # frequency for viscous attenuation
-p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
+p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 1600 # total number of time steps
+nt = 6000 # total number of time steps
deltat = 1.d-3 # duration of a time step
# source parameters
@@ -34,18 +35,18 @@ f0_attenuation = 5.196152422706633 # (Hz) relevant only if
# receiver line parameters for seismograms
seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
-generate_STATIONS = .true. # creates a STATION file in ./DATA
+generate_STATIONS = .false. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
# first receiver line (repeat these 6 lines and adjust nreceiverlines accordingly)
-nrec = 11 # number of receivers
-xdeb = 300. # first receiver x in meters
-zdeb = 2200. # first receiver z in meters
-xfin = 3700. # last receiver x in meters (ignored if onlyone receiver)
-zfin = 2200. # last receiver z in meters (ignored if onlyone receiver)
-enreg_surf_same_vertical = .true. # receivers inside the medium or at the surface
+nrec = 1 # number of receivers
+xdeb = 3000. # first receiver x in meters
+zdeb = 1500. # first receiver z in meters
+xfin = 3000. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 1500. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
# display parameters
NTSTEP_BETWEEN_OUTPUT_INFO = 100 # display frequency in time steps
@@ -63,20 +64,17 @@ sizemax_arrows = 1.d0 # maximum size of arrows on vec
gnuplot = .false. # generate a GNUPLOT file for the grid
output_grid = .false. # save the grid in a text file or not
output_energy = .false. # compute and output acoustic and elastic energy (slows down the code significantly)
-output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
+output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
# velocity and density models
-nbmodels = 4 # nb of different models
+nbmodels = 1 # nb of different models
# define models as
# I: (model_number 1 rho Vp Vs 0 0 Qp Qs 0 0 0 0 0 0) or
# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qs).
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
-1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
-2 1 2500.d0 2700.d0 0 0 0 9999 9999 0 0 0 0 0 0
-3 1 2200.d0 2500.d0 1443.375d0 0 0 9999 9999 0 0 0 0 0 0
-4 1 2200.d0 2200.d0 1343.375d0 0 0 9999 9999 0 0 0 0 0 0
+1 1 2300.d0 2450.d0 1425.d0 0 0 9999 9999 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
@@ -100,7 +98,7 @@ tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the
# PARAMETERS FOR INTERNAL MESHING
# file containing interfaces for internal mesh
-interfacesfile = ../interfaces_M2_UPPA_curved.dat
+interfacesfile = ../EXAMPLES/noise_uniform/uniform.dat
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
@@ -110,12 +108,9 @@ nx = 80 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
absorbright = .true.
-absorbtop = .false.
+absorbtop = .true.
absorbleft = .true.
# define the different regions of the model in the (nx,nz) spectral element mesh
-nbregions = 4 # nb of regions and model number for each
-1 80 1 20 1
-1 80 21 40 2
-1 80 41 60 3
-60 70 21 40 4
+nbregions = 1 # nb of regions and model number for each
+1 80 1 60 1
diff --git a/M2_UPPA/Par_file_M2_UPPA b/noise_uniform/Par_file_noise_2
similarity index 82%
copy from M2_UPPA/Par_file_M2_UPPA
copy to noise_uniform/Par_file_noise_2
index d76a229..d942f0c 100644
--- a/M2_UPPA/Par_file_M2_UPPA
+++ b/noise_uniform/Par_file_noise_2
@@ -1,9 +1,10 @@
# title of job
-title = Test for M2 UPPA
+title = Noise_2D
# forward or adjoint simulation
SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
-SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
+NOISE_TOMOGRAPHY = 2 # 0 = earthquake simulation, 1/2/3 = noise simulation
+SAVE_FORWARD = .true. # save the last frame, needed for adjoint simulation
# parameters concerning partitioning
nproc = 1 # number of processes
@@ -18,10 +19,10 @@ TURN_ATTENUATION_ON = .false. # turn attenuation on or off fo
TURN_VISCATTENUATION_ON = .false. # turn viscous attenuation on or off
Q0 = 1 # quality factor for viscous attenuation
freq0 = 10 # frequency for viscous attenuation
-p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
+p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 1600 # total number of time steps
+nt = 6000 # total number of time steps
deltat = 1.d-3 # duration of a time step
# source parameters
@@ -34,18 +35,18 @@ f0_attenuation = 5.196152422706633 # (Hz) relevant only if
# receiver line parameters for seismograms
seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
-generate_STATIONS = .true. # creates a STATION file in ./DATA
+generate_STATIONS = .false. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
# first receiver line (repeat these 6 lines and adjust nreceiverlines accordingly)
-nrec = 11 # number of receivers
-xdeb = 300. # first receiver x in meters
-zdeb = 2200. # first receiver z in meters
-xfin = 3700. # last receiver x in meters (ignored if onlyone receiver)
-zfin = 2200. # last receiver z in meters (ignored if onlyone receiver)
-enreg_surf_same_vertical = .true. # receivers inside the medium or at the surface
+nrec = 1 # number of receivers
+xdeb = 3000. # first receiver x in meters
+zdeb = 1500. # first receiver z in meters
+xfin = 3000. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 1500. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
# display parameters
NTSTEP_BETWEEN_OUTPUT_INFO = 100 # display frequency in time steps
@@ -63,20 +64,17 @@ sizemax_arrows = 1.d0 # maximum size of arrows on vec
gnuplot = .false. # generate a GNUPLOT file for the grid
output_grid = .false. # save the grid in a text file or not
output_energy = .false. # compute and output acoustic and elastic energy (slows down the code significantly)
-output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
+output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
# velocity and density models
-nbmodels = 4 # nb of different models
+nbmodels = 1 # nb of different models
# define models as
# I: (model_number 1 rho Vp Vs 0 0 Qp Qs 0 0 0 0 0 0) or
# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qs).
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
-1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
-2 1 2500.d0 2700.d0 0 0 0 9999 9999 0 0 0 0 0 0
-3 1 2200.d0 2500.d0 1443.375d0 0 0 9999 9999 0 0 0 0 0 0
-4 1 2200.d0 2200.d0 1343.375d0 0 0 9999 9999 0 0 0 0 0 0
+1 1 2300.d0 2450.d0 1425.d0 0 0 9999 9999 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
@@ -100,7 +98,7 @@ tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the
# PARAMETERS FOR INTERNAL MESHING
# file containing interfaces for internal mesh
-interfacesfile = ../interfaces_M2_UPPA_curved.dat
+interfacesfile = ../EXAMPLES/noise_uniform/uniform.dat
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
@@ -110,12 +108,9 @@ nx = 80 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
absorbright = .true.
-absorbtop = .false.
+absorbtop = .true.
absorbleft = .true.
# define the different regions of the model in the (nx,nz) spectral element mesh
-nbregions = 4 # nb of regions and model number for each
-1 80 1 20 1
-1 80 21 40 2
-1 80 41 60 3
-60 70 21 40 4
+nbregions = 1 # nb of regions and model number for each
+1 80 1 60 1
diff --git a/M2_UPPA/Par_file_M2_UPPA b/noise_uniform/Par_file_noise_3
similarity index 82%
copy from M2_UPPA/Par_file_M2_UPPA
copy to noise_uniform/Par_file_noise_3
index d76a229..263df60 100644
--- a/M2_UPPA/Par_file_M2_UPPA
+++ b/noise_uniform/Par_file_noise_3
@@ -1,8 +1,9 @@
# title of job
-title = Test for M2 UPPA
+title = Noise_2D
# forward or adjoint simulation
-SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
+SIMULATION_TYPE = 2 # 1 = forward, 2 = adjoint + kernels
+NOISE_TOMOGRAPHY = 3 # 0 = earthquake simulation, 1/2/3 = noise simulation
SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
# parameters concerning partitioning
@@ -18,10 +19,10 @@ TURN_ATTENUATION_ON = .false. # turn attenuation on or off fo
TURN_VISCATTENUATION_ON = .false. # turn viscous attenuation on or off
Q0 = 1 # quality factor for viscous attenuation
freq0 = 10 # frequency for viscous attenuation
-p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
+p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 1600 # total number of time steps
+nt = 6000 # total number of time steps
deltat = 1.d-3 # duration of a time step
# source parameters
@@ -34,18 +35,18 @@ f0_attenuation = 5.196152422706633 # (Hz) relevant only if
# receiver line parameters for seismograms
seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
-generate_STATIONS = .true. # creates a STATION file in ./DATA
+generate_STATIONS = .false. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
# first receiver line (repeat these 6 lines and adjust nreceiverlines accordingly)
-nrec = 11 # number of receivers
-xdeb = 300. # first receiver x in meters
-zdeb = 2200. # first receiver z in meters
-xfin = 3700. # last receiver x in meters (ignored if onlyone receiver)
-zfin = 2200. # last receiver z in meters (ignored if onlyone receiver)
-enreg_surf_same_vertical = .true. # receivers inside the medium or at the surface
+nrec = 1 # number of receivers
+xdeb = 3000. # first receiver x in meters
+zdeb = 1500. # first receiver z in meters
+xfin = 3000. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 1500. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
# display parameters
NTSTEP_BETWEEN_OUTPUT_INFO = 100 # display frequency in time steps
@@ -63,20 +64,17 @@ sizemax_arrows = 1.d0 # maximum size of arrows on vec
gnuplot = .false. # generate a GNUPLOT file for the grid
output_grid = .false. # save the grid in a text file or not
output_energy = .false. # compute and output acoustic and elastic energy (slows down the code significantly)
-output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
+output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
# velocity and density models
-nbmodels = 4 # nb of different models
+nbmodels = 1 # nb of different models
# define models as
# I: (model_number 1 rho Vp Vs 0 0 Qp Qs 0 0 0 0 0 0) or
# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qs).
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
-1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
-2 1 2500.d0 2700.d0 0 0 0 9999 9999 0 0 0 0 0 0
-3 1 2200.d0 2500.d0 1443.375d0 0 0 9999 9999 0 0 0 0 0 0
-4 1 2200.d0 2200.d0 1343.375d0 0 0 9999 9999 0 0 0 0 0 0
+1 1 2300.d0 2450.d0 1425.d0 0 0 9999 9999 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
@@ -100,7 +98,7 @@ tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the
# PARAMETERS FOR INTERNAL MESHING
# file containing interfaces for internal mesh
-interfacesfile = ../interfaces_M2_UPPA_curved.dat
+interfacesfile = ../EXAMPLES/noise_uniform/uniform.dat
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
@@ -110,12 +108,9 @@ nx = 80 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
absorbright = .true.
-absorbtop = .false.
+absorbtop = .true.
absorbleft = .true.
# define the different regions of the model in the (nx,nz) spectral element mesh
-nbregions = 4 # nb of regions and model number for each
-1 80 1 20 1
-1 80 21 40 2
-1 80 41 60 3
-60 70 21 40 4
+nbregions = 1 # nb of regions and model number for each
+1 80 1 60 1
diff --git a/DATA_to_sort_older_examples/SOURCE_attenuation_2D b/noise_uniform/SOURCE_noise
similarity index 60%
copy from DATA_to_sort_older_examples/SOURCE_attenuation_2D
copy to noise_uniform/SOURCE_noise
index 5e26462..e4e23e3 100644
--- a/DATA_to_sort_older_examples/SOURCE_attenuation_2D
+++ b/noise_uniform/SOURCE_noise
@@ -1,13 +1,13 @@
#source 1. 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.
source_surf = .false. # source inside the medium or at the surface
-xs = 1000. # source location x in meters
-zs = 1000. # source location z in meters
+xs = 0. # source location x in meters
+zs = 0. # source location z in meters
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
+time_function_type = 3 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+f0 = 10.0 # dominant source frequency (Hz) if not Dirac or Heaviside
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)
-Mxz = 0. # Mxz component (for a moment tensor source only)
-factor = 1.d10 # amplification factor
+Mxx = 0.d0 # Mxx component (for a moment tensor source only)
+Mzz = 0.d0 # Mzz component (for a moment tensor source only)
+Mxz = 0.d0 # Mxz component (for a moment tensor source only)
+factor = 0.d0 # amplification factor
diff --git a/noise_uniform/STATIONS_target_noise b/noise_uniform/STATIONS_target_noise
new file mode 100644
index 0000000..ea64c3f
--- /dev/null
+++ b/noise_uniform/STATIONS_target_noise
@@ -0,0 +1,3 @@
+S0001 AA 1000.0000000 1500.0000000 0.0 0.0
+S0002 AA 2000.0000000 1500.0000000 0.0 0.0
+S0003 AA 3000.0000000 1500.0000000 0.0 0.0
diff --git a/noise_uniform/adj_cc.f90 b/noise_uniform/adj_cc.f90
new file mode 100644
index 0000000..9e9e1a1
--- /dev/null
+++ b/noise_uniform/adj_cc.f90
@@ -0,0 +1,415 @@
+program adj_cc
+
+implicit none
+
+! flags
+logical, parameter :: use_filtering = .true.
+logical, parameter :: use_positive_branch = .false.
+
+! FILTERING PARAMETERS
+real freq_low,freq_high
+data freq_low / 1d-2 /
+data freq_high / 1d1 /
+
+! WINDOW PARAMETERS
+real :: w_delay, w_width, w_tukey
+data w_delay / 1.25 /
+data w_width / 1.00 /
+data w_tukey / 0.4 /
+!see explanation of
+!window parameters, below
+
+! time variables
+integer :: it, nt, nthalf
+double precision :: dt
+
+! data variables
+double precision, dimension(:), allocatable :: seismo_1, seismo_2, seismo_3, seismo_4, &
+ seismo_adj, t, w
+
+! input/ output
+character(len=64) :: file_in
+integer :: ios
+
+! miscellaneous
+double precision, parameter :: PI = 3.141592653589793
+integer :: it_off, it_wdt, it_begin, it_end, k
+integer :: ifreq, nfreq
+real :: F1,F2,D(8),G,DELT
+real :: alpha, beta
+
+
+
+
+! EXPLANATION OF WINDOW PARAMETERS
+
+!To select the desired branch of the cross-correlogram, we employ a Tukey window. A Tukey taper is just a variant of a cosine taper. We use three control parameters
+
+!W_DELAY controls the time offset of the window
+!W_WIDTH controls the width of the window (i.e., the total time range over which the window has non-zero support)
+!W_TUKEY controls the sharpness of the drop-off
+
+!In noise tomography applications, W_DELAY should be roughly equal to the surface wave travel time from the one receiver to the other.
+
+!Checks on W_WIDTH are carried out to make sure that the window makes sense and lies within a single branch of the cross-correlogram. If the the window falls outside these bounds, it will be adjusted.
+
+!W_TUKEY is a number between 0 and 1, 0 being pure boxcar and 1 being pure cosine
+
+
+
+! Get file info
+call getarg(1,file_in)
+call getlen(file_in,nt)
+call getinc(file_in,nt,dt)
+nthalf = (nt+1)/2
+
+write(*,*) ''
+write(*,*) 'This routine works only for evenly sampled cross-correlograms.'
+write(*,*) 'Reading from file: '//trim(file_in)
+write(*,'(a,i10)') ' nt: ', nt
+write(*,'(a,f10.3)') ' dt: ', dt
+
+! Allocate, initialize
+allocate(t(nt))
+allocate(w(nt))
+allocate(seismo_1(nt))
+allocate(seismo_2(nt))
+allocate(seismo_3(nt))
+allocate(seismo_4(nt))
+allocate(seismo_adj(nt))
+w(:) = 0.0d0
+seismo_1(:) = 0.0d0
+seismo_2(:) = 0.0d0
+seismo_3(:) = 0.0d0
+seismo_4(:) = 0.0d0
+seismo_adj(:) = 0.0d0
+
+
+!!!!!!!!!! READ INPUT !!!!!!!!!!!!!!!!!!!!
+open(unit=1001,file=trim(file_in),status='old',action='read')
+do it = 1, nt
+ read(1001,*) t(it), seismo_1(nt-it+1)
+end do
+close(1001)
+
+
+!!!!!!!!!! DIFFERENTIATE !!!!!!!!!!!!!!!!!!!
+seismo_1(1) = 0.0
+seismo_1(nt) = 0.0
+do it = 2, nt-1
+ seismo_2(it) = ( seismo_2(it+1) - seismo_1(it-1) ) / (2*dt)
+end do
+
+
+!!!!!!!!!! FILTER !!!!!!!!!!!!!!!!!!!!
+seismo_3 = seismo_2
+if (use_filtering) then
+! THIS SECTION CALCULATES THE FILTER AND MUST BE CALLED BEFORE
+! FILTER IS CALLED
+DELT = 1.0d3 * dt
+F1=freq_low
+F2=freq_high
+call BNDPAS(F1,F2,DELT,D,G,nt)
+! F1 = LOW FREQUENCY CUTOFF (6 DB DOWN)
+! F2 = HIGH FREQUENCY CUTOFF (6 DB DOWN)
+! DELT = SAMPLE INTERVAL IN MILLISECONDS
+! D = WILL CONTAIN 8 Z DOMAIN COEFICIENTS OF RECURSIVE FILTER
+! G = WILL CONTAIN THE GAIN OF THE FILTER,
+call FILTER(seismo_3,nt,D,G,2)
+! X = DATA VECTOR OF LENGTH N CONTAINING DATA TO BE FILTERED
+! D = FILTER COEFFICIENTS CALCULATED BY BNDPAS
+! G = FILTER GAIN
+! IG = 1 one pass
+! IG = 2 two passes
+end if
+
+
+!!!!!!!!!! WINDOW !!!!!!!!!!!!!!!!!!!!
+it_off = floor(w_delay/dt)
+it_wdt = 2*floor(w_width/(2.*dt))
+alpha = w_tukey
+
+if (use_positive_branch) then
+ write(*,*) 'Choosing positive branch'
+ it_begin = nthalf + it_off - it_wdt/2
+ it_end = nthalf + it_off + it_wdt/2
+ if (it_begin < nthalf) it_begin = nthalf
+ if (it_end > nt) it_end = nt
+else
+ write(*,*) 'Choosing negative branch'
+ it_begin = nthalf - it_off - it_wdt/2
+ it_end = nthalf - it_off + it_wdt/2
+ if (it_begin < 1) it_begin = 1
+ if (it_end > nthalf) it_end = nthalf
+endif
+
+write(*,'(a,2f10.3)') ' Time range: ', t(1), t(nt)
+write(*,'(a,2f10.3)') ' Window: ', t(it_begin), t(it_end)
+write(*,'(a,f10.3,f10.3)') ' Filtering: ', 1./freq_high, 1./freq_low
+
+!! Tukey taper
+k=0
+do it = it_begin,it_end
+ k=k+1
+ beta = real(k-1)/(it_end-it_begin)
+
+ if (beta<alpha/2.) then
+ w(it) = 0.5*(1.+cos(2.*pi/alpha*(beta-alpha/2.)))
+
+ elseif (beta>alpha/2. .and. beta<1.-alpha/2.) then
+ w(it) = 1.0
+
+ else
+ w(it) = 0.5*(1.+cos(2*pi/w_tukey*(beta-1.+alpha/2.)))
+
+ endif
+end do
+seismo_4 = w * seismo_3
+
+
+!!!!!!!!!! NORMALIZE !!!!!!!!!!!!!!!!!!!!
+seismo_adj = - seismo_4/(DOT_PRODUCT(seismo_4,seismo_4)*dt)
+
+
+!!!!!!!!!! WRITE ADJOINT SOURCE !!!!!!!!!!!!!!!!!!!!
+open(unit=1002,file=trim(file_in)//'.adj',status='unknown',iostat=ios)
+if (ios /= 0) write(*,*) 'Error opening output file.'
+
+write(*,*) ''
+write(*,*) 'Writing to file: '//trim(file_in)//'.adj'
+
+do it = 1,nt
+ write(1002,*), t(it), seismo_adj(it)
+end do
+close(1002)
+
+write(*,*) 'Finished writing to file.'
+write(*,*) ''
+
+
+end program adj_cc
+
+
+
+!=====================================================================
+subroutine getlen(filename,len)
+
+implicit none
+
+!input
+character(len=64) :: filename
+
+!output
+integer :: len
+
+!local
+integer, parameter :: IMAX = 1000000
+integer :: i,ios
+real :: dummy1, dummy2
+
+open(unit=1001,file=trim(filename),status='old',action='read')
+len=0
+do i=1,IMAX
+ read(1001,*,iostat=ios) dummy1, dummy2
+ if (ios==-1) exit
+ len=len+1
+enddo
+close(1001)
+
+end subroutine getlen
+
+
+
+!=====================================================================
+subroutine getinc(filename,len,inc)
+
+implicit none
+
+!input
+character(len=64) :: filename
+integer :: len
+
+!output
+double precision :: inc
+
+!local
+integer :: it
+double precision, dimension(len) :: t
+double precision :: sumdt
+real :: dummy
+
+open(unit=1001,file=trim(filename),status='old',action='read')
+do it=1,len
+ read(1001,*) t(it), dummy
+enddo
+close(1001)
+
+sumdt = 0.0d0
+do it=1,len-1
+ sumdt = sumdt + t(it+1) - t(it)
+enddo
+inc=sumdt/(len-1)
+
+end subroutine getinc
+
+
+!=====================================================================
+SUBROUTINE BNDPAS(F1,F2,DELT,D,G,N)
+! RECURSIVE BUTTERWORTH BAND PASS FILTER (KANASEWICH, TIME SERIES
+! ANALYSIS IN GEOPHYSICS, UNIVERSITY OF ALBERTA PRESS, 1975; SHANKS,
+! JOHN L, RECURSION FILTERS FOR DIGITAL PROCESSING, GEOPHYSICS, V32,
+! FILTER. THE FILTER WILL HAVE 8 POLES IN THE S PLANE AND IS
+! APPLIED IN FORWARD AND REVERSE DIRECTIONS SO AS TO HAVE ZERO
+! PHASE SHIFT. THE GAIN AT THE TWO FREQUENCIES SPECIFIED AS
+! CUTOFF FREQUENCIES WILL BE -6DB AND THE ROLLOFF WILL BE ABOUT
+! THE FILTER TO PREVENT ALIASING PROBLEMS.
+ COMPLEX P(4),S(8),Z1,Z2
+ real D(8),XC(3),XD(3),XE(3)
+ double precision :: X(N)
+ DATA ISW/0/,TWOPI/6.2831853/
+! THIS SECTION CALCULATES THE FILTER AND MUST BE CALLED BEFORE
+! FILTER IS CALLED
+
+! F1 = LOW FREQUENCY CUTOFF (6 DB DOWN)
+! F2 = HIGH FREQUENCY CUTOFF (6 DB DOWN)
+! DELT = SAMPLE INTERVAL IN MILLISECONDS
+! D = WILL CONTAIN 8 Z DOMAIN COEFICIENTS OF RECURSIVE FILTER
+! G = WILL CONTAIN THE GAIN OF THE FILTER,
+
+ DT=DELT/1000.0
+ TDT=2.0/DT
+ FDT=4.0/DT
+ ISW=1
+ P(1)=CMPLX(-.3826834,.9238795)
+ P(2)=CMPLX(-.3826834,-.9238795)
+ P(3)=CMPLX(-.9238795,.3826834)
+ P(4)=CMPLX(-.9238795,-.3826834)
+ W1=TWOPI*F1
+ W2=TWOPI*F2
+ W1=TDT*TAN(W1/TDT)
+ W2=TDT*TAN(W2/TDT)
+ HWID=(W2-W1)/2.0
+ WW=W1*W2
+ DO 19 I=1,4
+ Z1=P(I)*HWID
+ Z2=Z1*Z1-WW
+ Z2=CSQRT(Z2)
+ S(I)=Z1+Z2
+ 19 S(I+4)=Z1-Z2
+ G=.5/HWID
+ G=G*G
+ G=G*G
+ DO 29 I=1,7,2
+ B=-2.0*REAL(S(I))
+ Z1=S(I)*S(I+1)
+ C=REAL(Z1)
+ A=TDT+B+C/TDT
+ G=G*A
+ D(I)=(C*DT-FDT)/A
+ 29 D(I+1)=(A-2.0*B)/A
+ G=G*G
+ 5 FORMAT ('-FILTER GAIN IS ', 9E12.6)
+ RETURN
+
+ ENTRY FILTER(X,N,D,G,IG)
+
+! X = DATA VECTOR OF LENGTH N CONTAINING DATA TO BE FILTERED
+! D = FILTER COEFFICIENTS CALCULATED BY BNDPAS
+! G = FILTER GAIN
+! IG = 1 one pass
+! ig = 2 two passes
+
+ IF (ISW.EQ.1) GO TO 31
+ WRITE (6,6)
+ 6 FORMAT ('1BNDPAS MUST BE CALLED BEFORE FILTER')
+ return
+
+! APPLY FILTER IN FORWARD DIRECTION
+
+ 31 XM2=X(1)
+ XM1=X(2)
+ XM=X(3)
+ XC(1)=XM2
+ XC(2)=XM1-D(1)*XC(1)
+ XC(3)=XM-XM2-D(1)*XC(2)-D(2)*XC(1)
+ XD(1)=XC(1)
+ XD(2)=XC(2)-D(3)*XD(1)
+ XD(3)=XC(3)-XC(1)-D(3)*XD(2)-D(4)*XD(1)
+ XE(1)=XD(1)
+ XE(2)=XD(2)-D(5)*XE(1)
+ XE(3)=XD(3)-XD(1)-D(5)*XE(2)-D(6)*XE(1)
+ X(1)=XE(1)
+ X(2)=XE(2)-D(7)*X(1)
+ X(3)=XE(3)-XE(1)-D(7)*X(2)-D(8)*X(1)
+ DO 39 I=4,N
+ XM2=XM1
+ XM1=XM
+ XM=X(I)
+ K=I-((I-1)/3)*3
+ GO TO (34,35,36),K
+ 34 M=1
+ M1=3
+ M2=2
+ GO TO 37
+ 35 M=2
+ M1=1
+ M2=3
+ GO TO 37
+ 36 M=3
+ M1=2
+ M2=1
+ 37 XC(M)=XM-XM2-D(1)*XC(M1)-D(2)*XC(M2)
+ XD(M)=XC(M)-XC(M2)-D(3)*XD(M1)-D(4)*XD(M2)
+ XE(M)=XD(M)-XD(M2)-D(5)*XE(M1)-D(6)*XE(M2)
+ 39 X(I)=XE(M)-XE(M2)-D(7)*X(I-1)-D(8)*X(I-2)
+!
+!
+ if(ig.eq.1) goto 3333
+ XM2=X(N)
+ XM1=X(N-1)
+ XM=X(N-2)
+ XC(1)=XM2
+ XC(2)=XM1-D(1)*XC(1)
+ XC(3)=XM-XM2-D(1)*XC(2)-D(2)*XC(1)
+ XD(1)=XC(1)
+ XD(2)=XC(2)-D(3)*XD(1)
+ XD(3)=XC(3)-XC(1)-D(3)*XD(2)-D(4)*XD(1)
+ XE(1)=XD(1)
+ XE(2)=XD(2)-D(5)*XE(1)
+ XE(3)=XD(3)-XD(1)-D(5)*XE(2)-D(6)*XE(1)
+ X(N)=XE(1)
+ X(N-1)=XE(2)-D(7)*X(1)
+ X(N-2)=XE(3)-XE(1)-D(7)*X(2)-D(8)*X(1)
+ DO 49 I=4,N
+ XM2=XM1
+ XM1=XM
+ J=N-I+1
+ XM=X(J)
+ K=I-((I-1)/3)*3
+ GO TO (44,45,46),K
+ 44 M=1
+ M1=3
+ M2=2
+ GO TO 47
+ 45 M=2
+ M1=1
+ M2=3
+ GO TO 47
+ 46 M=3
+ M1=2
+ M2=1
+ 47 XC(M)=XM-XM2-D(1)*XC(M1)-D(2)*XC(M2)
+ XD(M)=XC(M)-XC(M2)-D(3)*XD(M1)-D(4)*XD(M2)
+ XE(M)=XD(M)-XD(M2)-D(5)*XE(M1)-D(6)*XE(M2)
+ 49 X(J)=XE(M)-XE(M2)-D(7)*X(J+1)-D(8)*X(J+2)
+ 3333 continue
+ if (ig.eq.1) then
+ gg=sqrt(g) ! if only pass once, modify gain
+ else
+ gg=g
+ endif
+ DO 59 I=1,N
+ 59 X(I)=X(I)/gg
+ RETURN
+END
+
diff --git a/noise_uniform/process.sh b/noise_uniform/process.sh
new file mode 100755
index 0000000..8726c3d
--- /dev/null
+++ b/noise_uniform/process.sh
@@ -0,0 +1,67 @@
+#!/bin/sh
+
+USAGE="USAGE: submit par_file_directory"
+if [ $# -eq 0 ]; then PAR_DIR=uniform; fi
+if [ $# -eq 1 ]; then PAR_DIR=$1; fi
+if [ $# -ne 0 ] && [ $# -ne 1 ]; then echo "$USAGE"; exit 1; fi
+
+
+PAR_DIR_FULL=$PWD/$PAR_DIR
+RUN_DIR=../..
+cd $RUN_DIR
+
+
+# prepare directories
+rm -rf SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL
+mkdir -p SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL
+
+
+# prepare files
+cp $PAR_DIR_FULL/SOURCE_noise DATA/SOURCE
+cp $PAR_DIR_FULL/STATIONS_target_noise DATA/STATIONS_target
+#cp $PAR_DIR_FULL/S_squared NOISE_TOMOGRAPHY
+echo 1 > NOISE_TOMOGRAPHY/irec_master
+
+
+#simulation 1
+cp $PAR_DIR_FULL/Par_file_noise_1 DATA/Par_file
+make; bin/xmeshfem2D; bin/xspecfem2D
+mkdir OUTPUT_ALL/step_1
+mv OUTPUT_FILES/image* OUTPUT_ALL/step_1
+mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_1
+mv DATA/Par_file OUTPUT_ALL/step_1
+
+
+#simulation 2
+cp $PAR_DIR_FULL/Par_file_noise_2 DATA/Par_file
+bin/xmeshfem2D; bin/xspecfem2D
+mkdir OUTPUT_ALL/step_2
+
+ADJ_CODE=$PAR_DIR_FULL/adj_cc.f90
+gfortran $ADJ_CODE -o adj_cc
+cp OUTPUT_FILES/*.semd SEM
+./adj_cc SEM/S0003.AA.BXY.semd
+
+cd SEM
+rename '.semd' '' *.adj
+awk '{printf(" %20.10f %20.10f\n",$1,0.)}' < S0003.AA.BXY.adj > ZEROS
+cp ZEROS S0001.AA.BXX.adj; cp ZEROS S0001.AA.BXY.adj; cp ZEROS S0001.AA.BXZ.adj
+cp ZEROS S0002.AA.BXX.adj; cp ZEROS S0002.AA.BXY.adj; cp ZEROS S0002.AA.BXZ.adj
+cp ZEROS S0003.AA.BXX.adj; cp ZEROS S0003.AA.BXZ.adj
+
+cd ..
+
+mv OUTPUT_FILES/image* OUTPUT_ALL/step_2
+mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_2
+mv DATA/Par_file OUTPUT_ALL/step_2
+
+
+#simulation 3
+cp $PAR_DIR_FULL/Par_file_noise_3 DATA/Par_file
+bin/xmeshfem2D; bin/xspecfem2D
+mkdir OUTPUT_ALL/step_3
+mv OUTPUT_FILES/image* OUTPUT_ALL/step_3
+mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_3
+mv OUTPUT_FILES/proc* OUTPUT_ALL/step_3
+mv DATA/Par_file OUTPUT_ALL/step_3
+
diff --git a/DATA_to_sort_older_examples/interfaces_attenuation_analytic.dat b/noise_uniform/uniform.dat
similarity index 81%
copy from DATA_to_sort_older_examples/interfaces_attenuation_analytic.dat
copy to noise_uniform/uniform.dat
index 45d1cbf..f444c9b 100644
--- a/DATA_to_sort_older_examples/interfaces_attenuation_analytic.dat
+++ b/noise_uniform/uniform.dat
@@ -10,16 +10,17 @@
#
2
0 0
- 5000 0
+ 4000 0
#
# interface number 2
#
2
- 0 2000
- 5000 2000
+ 0 3000
+ 4000 3000
+#
# for each layer, we give the number of spectral elements in the vertical direction
#
#
-# layer number 1 (bottom layer)
+# layer number 1
#
- 44
+ 60
More information about the CIG-COMMITS
mailing list