[cig-commits] r15876 - in seismo/3D/ADJOINT_TOMO/flexwin: . USER_MANUAL scripts scripts/prepare_scripts/socal user_files/socal_3D
carltape at geodynamics.org
carltape at geodynamics.org
Sun Oct 25 16:06:35 PDT 2009
Author: carltape
Date: 2009-10-25 16:06:34 -0700 (Sun, 25 Oct 2009)
New Revision: 15876
Added:
seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T003_T010_m16
seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m16_man.f90
Modified:
seismo/3D/ADJOINT_TOMO/flexwin/USER_MANUAL/AM-allcitations.bib
seismo/3D/ADJOINT_TOMO/flexwin/USER_MANUAL/flexwin_manual.pdf
seismo/3D/ADJOINT_TOMO/flexwin/USER_MANUAL/manual_tuning.tex
seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows.pl
seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows_local.pl
seismo/3D/ADJOINT_TOMO/flexwin/scripts/plot_windows_all.pl
seismo/3D/ADJOINT_TOMO/flexwin/scripts/prepare_scripts/socal/process_data_and_syn.pl
seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90
seismo/3D/ADJOINT_TOMO/flexwin/select_windows_stalta2.f90
seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m00.f90
seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m16.f90
seismo/3D/ADJOINT_TOMO/flexwin/user_functions.f90
Log:
Added two variables required for user_functions.f90: noise_start and signal_start, to complement noise_end and signal_end. This provides flexibility in how to specify the noise time interval and the signal time interval. This required modifying the subroutine check_data_quality in select_windows_stalta2.f90; the manual was also modified to include the change. Minor updates to socal run files.
Modified: seismo/3D/ADJOINT_TOMO/flexwin/USER_MANUAL/AM-allcitations.bib
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/USER_MANUAL/AM-allcitations.bib 2009-10-25 20:44:01 UTC (rev 15875)
+++ seismo/3D/ADJOINT_TOMO/flexwin/USER_MANUAL/AM-allcitations.bib 2009-10-25 23:06:34 UTC (rev 15876)
@@ -4185,10 +4185,10 @@
@article{MaggiEtal2009,
Author = {Maggi, A. and Tape, C. and Chen, M. and Chao, D. and Tromp, J.},
Journal = gji,
+ Pages = {257--281},
Title = {An automated time-window selection algorithm for seismic tomography},
Year = 2009,
- note = {(in press)}
- }
+ Volume = 178}
@article{MJJ99,
Author = {Mahatsente, R. and Jentzsch, G. and Jahr, T.},
Modified: seismo/3D/ADJOINT_TOMO/flexwin/USER_MANUAL/flexwin_manual.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/ADJOINT_TOMO/flexwin/USER_MANUAL/manual_tuning.tex
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/USER_MANUAL/manual_tuning.tex 2009-10-25 20:44:01 UTC (rev 15875)
+++ seismo/3D/ADJOINT_TOMO/flexwin/USER_MANUAL/manual_tuning.tex 2009-10-25 23:06:34 UTC (rev 15876)
@@ -48,7 +48,7 @@
\item[{\tt TSHIFT\_BASE}]Corresponds to $\Delta\tau_0$ in Table~\ref{tb:params}, and is the maximum cross-correlation lag (in seconds) between synthetic and observed seismogram for a window to be acceptable.
\item[{\tt DLNA\_BASE}]Corresponds to $\Delta\ln{A}_0$ in Table~\ref{tb:params}, and is the maximum amplitude ratio ($\Delta\ln{A}$ or $\Delta A/A$) between synthetic and observed seismogram for a window to be acceptable.
\item[{\tt TSHIFT\_REFERENCE}]Corresponds to $\Delta\tau_{\rm ref}$ in Table~\ref{tb:params}, and allows for a systematic traveltime bias in the synthetics.
-\item[{\tt TSHIFT\_REFERENCE}]Corresponds to $\Delta\ln{A}_{\rm ref}$ in Table~\ref{tb:params}, and allows for a systematic amplitude bias in the synthetics.
+\item[{\tt DLNA\_REFERENCE}]Corresponds to $\Delta\ln{A}_{\rm ref}$ in Table~\ref{tb:params}, and allows for a systematic amplitude bias in the synthetics.
\item[{\tt C\_0}]Corresponds to $C_0$ in Table~\ref{tb:params}, and is expressed as a multiple of $w_E$. No window may contain a local minimum in its STA:LTA waveform that falls below the local value of $C_0 w_E$. See Figure~\ref{fg:win_composite}b.
\item[{\tt C\_1}]Corresponds to $C_1$ in Table~\ref{tb:params}, and is expressed as a multiple of $T_0$. No window may be shorter than $C_1 T_0$.
\item[{\tt C\_2}]Corresponds to $C_2$ in Table~\ref{tb:params}, and is expressed as a multiple of $w_E$. A window whose seed maximum on the STA:LTA waveform rises less than $C_2 w_E$ above either of its adjacent minima is rejected. See Figure~\ref{fg:win_composite}c.
@@ -352,7 +352,24 @@
%\Delta\tau_0(t) & = \Delta\tau_0.
\end{align}
+%-----------------------
+\subsection{Time intervals for signal-to-noise calculations}
+
+If the parameter {\tt DATA$\_$QUALITY = .true.} in the {\tt PAR$\_$FILE}, then the {\tt user\_functions.f90} file requires the specification of four times that define the ``noise'' and ``signal'' time intervals. For example, the default {\tt user\_functions.f90} file contains these lines:
+%
+\begin{verbatim}
+ ! these values will be used for signal2noise calculations
+ if (DATA_QUALITY) then
+ noise_start=b
+ noise_end=ph_times(1)-WIN_MAX_PERIOD
+ signal_start=noise_end
+ signal_end=b+(npts-1)*dt
+ endif
+\end{verbatim}
+%
+The use of all four variables allows provides the flexibility of choosing these time intervals. (For example, the start of the noise does not need to be the beginning of the seismogram. Or the start of the signal does not have to coincide with the end of the noise.)
+
%-----------------------
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows.pl 2009-10-25 20:44:01 UTC (rev 15875)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows.pl 2009-10-25 23:06:34 UTC (rev 15876)
@@ -45,6 +45,8 @@
# pick_all_windows.pl m00 0 6/30 179/179 1/1/0/0 1 0 # make plots only, T = 6-30s
# pick_all_windows.pl m00 0 6/30 179/179 0/0/1/0 1 0 # make WINDOWS file, T = 6-30s
#
+# pick_all_windows.pl m16 0 3/10 198/198 1/1/1/0 1 0 # socal 14236768
+#
#==========================================================
if (@ARGV < 7) {die("Usage: pick_all_windows.pl model Tmin/Tmax imin/imax idebug/iplot/imeas/ibody idataset iexecute\n");}
@@ -73,18 +75,15 @@
if($ibp==1) {$bbp = ".true."} else {$bbp = ".false."}
# directories for data and synthetics
+# 1=socal; 2=Japan
if($idataset == 1) {
-
- $dir_data = "/home/carltape/SOCAL_ADJOINT/DATA/FINAL"; # Socal data (Carl)
- $dir_syn = "/home/carltape/SOCAL_ADJOINT/SYN/model_${smodel}"; # Socal syn
- #$dir_syn = "/home/carltape/SOCAL_ADJOINT/SYN/model_pre_${smodel}"; # Socal syn
+ $dir_data = "/home/carltape/SOCAL_ADJOINT/DATA/FINAL";
+ $dir_syn = "/home/carltape/SOCAL_ADJOINT/SYN/model_${smodel}";
$sdataset = "socal";
-
} elsif ($idataset == 2) {
- $dir_data = "/net/sierra/raid1/mchen/DATA/TEST"; # Japan data (Min)
- $dir_syn = "/net/sierra/raid1/mchen/SEM/TEST"; # Japan syn (Min)
+ $dir_data = "/home/mchen/DATA/TEST";
+ $dir_syn = "/home/mchen/SEM/TEST";
$sdataset = "japan";
-
} else {
die("\n idataset must be 1 (socal) or 2 (japan)\n");
}
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows_local.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows_local.pl 2009-10-25 20:44:01 UTC (rev 15875)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows_local.pl 2009-10-25 23:06:34 UTC (rev 15876)
@@ -4,7 +4,7 @@
#
# pick_all_windows_local.pl
# Carl Tape
-# 02-Oct-2008
+# 20-Ocit-2009
#
# This is the "local version" of pick_all_windows_local.pl
# This is run from the directory containing MEASURE, DATA, SYN,
@@ -24,7 +24,7 @@
# prepare_meas_all.pl
#
# EXAMPLES:
-# /net/denali/raid1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO/flexwin_work/scripts/pick_all_windows_local.pl
+# /data1/cig/seismo/3D/ADJOINT_TOMO/flexwin_work/scripts/pick_all_windows_local.pl
#
# pick_all_windows_local.pl 2/30 1/1/1/1 input MEASURE
# pick_all_windows_local.pl 6/30 1/1/1/0 input MEASURE
@@ -34,20 +34,27 @@
# pick_all_windows_local.pl 6/30 1/1/0/0 input_test MEASURE_TEST # plots only
# pick_all_windows_local.pl 2/30 1/1/0/1 input_test MEASURE_TEST # plots only
#
+# pick_all_windows_local.pl 3/10 1/1/1/0 input_test MEASURE_TEST # plots AND window file
+#
#==========================================================
if (@ARGV < 4) {die("Usage: pick_all_windows_local.pl Tmin/Tmax idebug/iplot/imeas/ibody input_file MEASURE_dir\n");}
($Ts,$ibools,$input,$dir) = @ARGV;
-# specify source code directory (MUST BE MODIFIED FOR EACH USER)
-$dir_win_code = "/net/denali/raid1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO/flexwin_work";
+$pwd = $ENV{PWD};
+#-------------------------------------
+# USER INPUT
+
+# source code directory
+$dir_win_code = "/home/carltape/ADJOINT_TOMO/flexwin_work";
+
+#-------------------------------------
+
$dir_scripts = "${dir_win_code}/scripts";
if (not -e ${dir_win_code}) {die("check if ${dir_win_code} exist or not\n")}
if (not -e ${dir_scripts}) {die("check if ${dir_scripts} exist or not\n")}
-#-------------------------------------
-
# split input entries
($Tmin,$Tmax) = split("/",$Ts);
($imin,$imax) = split("/",$inds);
@@ -72,8 +79,7 @@
$win_execute = "flexwin";
# make command for windowing code
-$make = "make -f make_gfortran"; # change as of 9-30-08
-#$make = "make -f make_intel_caltech";
+$make = "make -f make_gfortran";
# location of plot_windows_all.pl
$plot_windows_perl = "${dir_scripts}/plot_windows_all.pl";
@@ -118,7 +124,7 @@
# remove contents in target directory
# NOTE: This copies in the PAR_FILE from the main directory -- it is not the local version!
-print CSH "cd -\n";
+print CSH "cd $pwd\n";
print CSH "\\cp ${par_file} .\n";
print CSH "\\rm -rf ${dir}\n mkdir ${dir}\n";
@@ -136,6 +142,7 @@
print CSH "sort -g -k 2 $reclist0 > ${reclist_dist}\n";
# generate composite PDF file using plot_windows_all.pl
+ # make sure Carl's settings are on in plot_seismos_gmt.sh
$pdffile = "$dir/${eout}_all.pdf";
print CSH "${plot_windows_perl} ${dir_win_code} ${dir} ${reclist_dist} $pdffile\n";
}
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/plot_windows_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/plot_windows_all.pl 2009-10-25 20:44:01 UTC (rev 15875)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/plot_windows_all.pl 2009-10-25 23:06:34 UTC (rev 15876)
@@ -68,6 +68,7 @@
$obsfile = "${dir_win_run_meas}/${filename}.obs";
# generate output PDF files -- make sure that plot_seismos_gmt.sh is working
+ # NOTE: turn Carl's settings on in plot_seismos_gmt.sh
if (-f $obsfile) {
print "$obsfile\n";
$ftag = "${dir_win_run_meas}/${filename}";
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/prepare_scripts/socal/process_data_and_syn.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/prepare_scripts/socal/process_data_and_syn.pl 2009-10-25 20:44:01 UTC (rev 15875)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/prepare_scripts/socal/process_data_and_syn.pl 2009-10-25 23:06:34 UTC (rev 15876)
@@ -3,7 +3,7 @@
#==========================================================
#
# Carl Tape
-# 27-March-2009
+# 19-Oct-2009
# process_data_and_syn.pl
#
# This script processes data and 3D-Socal-SEM synthetics for Southern California.
@@ -25,20 +25,20 @@
# Trange band-pass filter
#
# ORDER OF OPERATIONS:
-# ~/UTILS/process_data_and_syn.pl 1 m16 1 0 d 6/30 # data, initial pre-processing
-# ~/UTILS/process_data_and_syn.pl 1 m16 0 1 d 6/30 # syn, initial pre-processing
-# ~/UTILS/process_data_and_syn.pl 2 m16 1 1 d 6/30 # both, create cut files
-# ~/UTILS/process_data_and_syn.pl 3 m16 1 0 d 6/30 # data, execute cut file
-# ~/UTILS/process_data_and_syn.pl 3 m16 0 1 d 6/30 # syn, execute cut file
+# ~/UTILS/process_data_and_syn.pl 1 m16 1 0 d 6/30 # data, initial pre-processing
+# ~/UTILS/process_data_and_syn.pl 1 m16 0 1 d 6/30 # syn, initial pre-processing (REPEAT)
+# ~/UTILS/process_data_and_syn.pl 2 m16 1 1 d 6/30 # both, create cut files
+# ~/UTILS/process_data_and_syn.pl 3 m16 1 0 d 6/30 # data, execute cut file
+# ~/UTILS/process_data_and_syn.pl 3 m16 0 1 d 6/30 # syn, execute cut file
#
-# ~/UTILS/process_data_and_syn.pl 4 m16 1 0 d 6/30 # data, bandpass T=6-30
-# ~/UTILS/process_data_and_syn.pl 4 m16 0 1 d 6/30 # syn, bandpass T=6-30
+# ~/UTILS/process_data_and_syn.pl 4 m16 1 0 d 6/30 # data, bandpass T=6-30
+# ~/UTILS/process_data_and_syn.pl 4 m16 0 1 d 6/30 # syn, bandpass T=6-30
#
-# ~/UTILS/process_data_and_syn.pl 4 m16 1 0 d 3/30 # data, bandpass T=3-30
-# ~/UTILS/process_data_and_syn.pl 4 m16 0 1 d 3/30 # syn, bandpass T=3-30
+# ~/UTILS/process_data_and_syn.pl 4 m16 1 0 d 3/30 # data, bandpass T=3-30
+# ~/UTILS/process_data_and_syn.pl 4 m16 0 1 d 3/30 # syn, bandpass T=3-30
#
-# ~/UTILS/process_data_and_syn.pl 4 m16 1 0 d 2/30 # data, bandpass T=2-30
-# ~/UTILS/process_data_and_syn.pl 4 m16 0 1 d 2/30 # syn, bandpass T=2-30
+# ~/UTILS/process_data_and_syn.pl 4 m16 1 0 d 2/30 # data, bandpass T=2-30
+# ~/UTILS/process_data_and_syn.pl 4 m16 0 1 d 2/30 # syn, bandpass T=2-30
#
# FILE-NAMING CONVENTIONS (example for event 9818433):
# DATA -- main data directory
@@ -73,6 +73,7 @@
$sps = 20; # samples per second for interpolation/sub-sampling
$tfac = 1.0; # factor to extend lengths of records (should be > 1.0)
$bmin = -40; # minimum allowable time before origin time for cut records
+$chdur = 0; # convolve synthetics with Gaussian half duration (default = 1)
$pdir = "PROCESSED"; # tag for processed directories
$syn_ext = $smodel; # KEY: include model index for comparison among synthetics from different models
@@ -85,21 +86,21 @@
$ilist = 1; # read EIDs from a list (=1) or grab all CMTSOLUTION files (=0)
# EID list
-$CMT_list = "/net/sierra/raid1/carltape/socal/socal_3D/SYN/model_${smodel}";
-$dir_source = "/net/sierra/raid1/carltape/results/SOURCES/socal_16";
+$CMT_list = "/home/carltape/SOCAL_ADJOINT/SYN/model_${smodel}";
+$dir_source = "/home/carltape/results/SOURCES/socal_16";
$dirCMT = "${dir_source}/v16_files";
$CMT_list = "${dir_source}/EIDs_only_loc";
-#$CMT_list = "/net/sierra/raid1/carltape/results/EID_LISTS/syn_run_${smodel}";
-#$CMT_list = "/net/sierra/raid1/carltape/results/EID_LISTS/syn_run_m12";
+#$CMT_list = "/home/carltape/results/EID_LISTS/syn_run_${smodel}";
+#$CMT_list = "/home/carltape/results/EID_LISTS/syn_run_m12";
# directories
-$dirdat0 = "/net/sierra/raid1/carltape/socal/socal_3D/DATA/FINAL";
-$dirsyn0 = "/net/sierra/raid1/carltape/socal/socal_3D/SYN/model_${smodel}";
-#$dirsyn0 = "/net/sierra/raid1/carltape/socal/socal_3D/SYN/model_pre_${smodel}";
+$dirdat0 = "/home/carltape/SOCAL_ADJOINT/DATA/FINAL";
+$dirsyn0 = "/home/carltape/SOCAL_ADJOINT/SYN/model_${smodel}";
+#$dirsyn0 = "/home/carltape/SOCAL_ADJOINT/SYN/model_pre_${smodel}";
# STATIONS file
-#$stafile = "/net/denali/home1/carltape/gmt/stations/seismic/Matlab_output/STATIONS";
-$stafile = "/net/denali/home1/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem";
-#$stafile = "/net/denali/home1/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_OUTER_specfem";
+#$stafile = "/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS";
+$stafile = "/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem";
+#$stafile = "/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_OUTER_specfem";
#-------------------------------------
@@ -158,7 +159,7 @@
$imin = 1; $imax = $ncmt; # default
#$imin = 1; $imax = 10;
-#$imin = 228; $imax = $imin;
+$imin = 76; $imax = $imin;
#----------------------------------------------------------------------
@@ -236,7 +237,11 @@
print CSH "process_trinet_syn_new.pl -m $cmtfile -a $stafile *semd\n";
print CSH "sleep 5s\n";
} else {
- print CSH "process_trinet_syn_new.pl -S -m $cmtfile -h -a $stafile -s $sps -p -d $pdir -x ${syn_ext} *.${syn_suffix0}\n";
+ if ($chdur == 1) {
+ print CSH "process_trinet_syn_new.pl -S -m $cmtfile -h -a $stafile -s $sps -p -d $pdir -x ${syn_ext} *.${syn_suffix0}\n";
+ } else {
+ print CSH "process_trinet_syn_new.pl -S -m $cmtfile -a $stafile -s $sps -p -d $pdir -x ${syn_ext} *.${syn_suffix0}\n";
+ }
}
} else {
print "dir ${dirsyn_pro_1} already exists\n";
Modified: seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90 2009-10-25 20:44:01 UTC (rev 15875)
+++ seismo/3D/ADJOINT_TOMO/flexwin/seismo_subs.f90 2009-10-25 23:06:34 UTC (rev 15876)
@@ -111,9 +111,12 @@
character*8, dimension(MAX_PHASES) :: ph_names
double precision, dimension(MAX_PHASES) :: ph_times
- ! end times for signal and noise
- double precision :: signal_end, noise_end
+ ! start/end times for signal and noise
+ double precision :: noise_start, noise_end, signal_start, signal_end
+ ! indices of start/end times for signal and noise
+ integer :: in1, in2, is1, is2
+
end module seismo_variables
!----------------------------------------------------------------------
Modified: seismo/3D/ADJOINT_TOMO/flexwin/select_windows_stalta2.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/select_windows_stalta2.f90 2009-10-25 20:44:01 UTC (rev 15875)
+++ seismo/3D/ADJOINT_TOMO/flexwin/select_windows_stalta2.f90 2009-10-25 23:06:34 UTC (rev 15876)
@@ -15,7 +15,7 @@
double precision :: w_level_local
logical :: data_ok
- double precision :: time_obs_noise, time_obs_signal, noise_int, signal_int, snr_int, signal_amp, noise_amp, snr_amp
+ double precision :: noise_int, signal_int, snr_int, signal_amp, noise_amp, snr_amp
! initialise global arrays
num_win = 0
@@ -37,7 +37,7 @@
if(DATA_QUALITY) then
if(DEBUG) print *, "Checking data quality"
! returns data_ok = .true. if two signal-to-noise criteria are met
- call check_data_quality(data_ok,signal_int,noise_int,snr_int,time_obs_signal,time_obs_noise,signal_amp,noise_amp,snr_amp)
+ call check_data_quality(data_ok,signal_int,noise_int,snr_int,signal_amp,noise_amp,snr_amp)
! if data quality is not ok, then return without selecting windows
if (.not. data_ok) then
num_win = 0
@@ -45,105 +45,125 @@
endif
endif
- ! find list of maxima and minima, ignoring water level
- w_level_local = 0.0
- if(DEBUG) print *, "Finding max min "
- call find_maxima(STA_LTA,npts,maxima_lp,nmax,w_level_local)
- call find_minima(STA_LTA,npts,minima_lp,nmin,w_level_local)
+!!$ if (1==1) then
+!!$ ! keep the manually selected window that has met the signal-to-noise criteria;
+!!$ ! this option is used for time-reversing coda surface waves
+!!$ nwin = 1
+!!$ iL(1) = is1
+!!$ iM(1) = int((is1+is2)/2)
+!!$ iR(1) = is2
+!!$
+!!$ ! reject window if it does not satisfy F2, tshift and dlnA criteria;
+!!$ ! here we are only interested in using the dlnA criterion
+!!$ call reject_on_fit_criteria(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
+!!$ if (nwin.eq.0) then
+!!$ write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
+!!$ num_win = 0
+!!$ return
+!!$ endif
+!!$ if(DEBUG) call display_windows_full(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
+!!$
+!!$ else
- ! find all possible windows within seismogram, with central maxima above the water level
- if(DEBUG) print *, "Making all windows "
- call setup_M_L_R(nmax,maxima_lp,nmin,minima_lp,nwin,iM,iL,iR)
- if (nwin.eq.0) then
- write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
- num_win = 0
- return
- endif
+ ! find list of maxima and minima, ignoring water level
+ w_level_local = 0.0
+ if(DEBUG) print *, "Finding max min "
+ call find_maxima(STA_LTA,npts,maxima_lp,nmax,w_level_local)
+ call find_minima(STA_LTA,npts,minima_lp,nmin,w_level_local)
- ! remove windows with internal minima below the water level
- if(DEBUG) print *, "Rejecting on water level "
- call reject_on_water_level(nwin,iM,iL,iR,nmin,minima_lp,C_0)
- if (nwin.eq.0) then
- write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
- num_win = 0
- return
- endif
- ! sort windows
- call sort_on_start_time(nwin,iM,iL,iR)
- !if(DEBUG) call display_windows(nwin,iM,iL,iR)
+ ! find all possible windows within seismogram, with central maxima above the water level
+ if(DEBUG) print *, "Making all windows "
+ call setup_M_L_R(nmax,maxima_lp,nmin,minima_lp,nwin,iM,iL,iR)
+ if (nwin.eq.0) then
+ write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
+ num_win = 0
+ return
+ endif
- ! remove small windows
- if(DEBUG) print *, "Rejecting on size "
- call reject_on_window_width(nwin,iM,iL,iR,C_1)
- if (nwin.eq.0) then
- write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
- num_win = 0
- return
- endif
- !if(DEBUG) call display_windows(nwin,iM,iL,iR)
+ ! remove windows with internal minima below the water level
+ if(DEBUG) print *, "Rejecting on water level "
+ call reject_on_water_level(nwin,iM,iL,iR,nmin,minima_lp,C_0)
+ if (nwin.eq.0) then
+ write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
+ num_win = 0
+ return
+ endif
+ ! sort windows
+ call sort_on_start_time(nwin,iM,iL,iR)
+ !if(DEBUG) call display_windows(nwin,iM,iL,iR)
- ! reject on prominence of central maximum
- if(DEBUG) print *, "Rejecting on prominence "
- call reject_on_prominence(nwin,iM,iL,iR,nmin,minima_lp,C_2)
- if (nwin.eq.0) then
- write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
- num_win = 0
- return
- endif
- !if(DEBUG) call display_windows(nwin,iM,iL,iR)
+ ! remove small windows
+ if(DEBUG) print *, "Rejecting on size "
+ call reject_on_window_width(nwin,iM,iL,iR,C_1)
+ if (nwin.eq.0) then
+ write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
+ num_win = 0
+ return
+ endif
+ !if(DEBUG) call display_windows(nwin,iM,iL,iR)
- ! reject windows containing more than one distinct phase
- call reject_on_phase_separation(nwin,iM,iL,iR,nmin,minima_lp,nmax,maxima_lp,C_3a, C_3b)
- if (nwin.eq.0) then
- write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
- num_win = 0
- return
- endif
- !if(DEBUG) call display_windows(nwin,iM,iL,iR)
+ ! reject on prominence of central maximum
+ if(DEBUG) print *, "Rejecting on prominence "
+ call reject_on_prominence(nwin,iM,iL,iR,nmin,minima_lp,C_2)
+ if (nwin.eq.0) then
+ write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
+ num_win = 0
+ return
+ endif
+ !if(DEBUG) call display_windows(nwin,iM,iL,iR)
- ! curtail window ends by time decay
- call curtail_window_length(nwin,iL,iR,nmax,maxima_lp,C_4a,C_4b)
- !if(DEBUG) call display_windows(nwin,iM,iL,iR)
+ ! reject windows containing more than one distinct phase
+ call reject_on_phase_separation(nwin,iM,iL,iR,nmin,minima_lp,nmax,maxima_lp,C_3a, C_3b)
+ if (nwin.eq.0) then
+ write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
+ num_win = 0
+ return
+ endif
+ !if(DEBUG) call display_windows(nwin,iM,iL,iR)
- ! REPEAT: remove small windows, since curtailing my have shrunk them
- ! ALESSIA: Perhaps this only needs to be done once (here)?
- if(DEBUG) print *, "Rejecting on size (REPEAT)"
- call reject_on_window_width(nwin,iM,iL,iR,C_1)
- if (nwin.eq.0) then
- write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
- num_win = 0
- return
- endif
- if(DEBUG) call display_windows(nwin,iM,iL,iR)
+ ! curtail window ends by time decay
+ call curtail_window_length(nwin,iL,iR,nmax,maxima_lp,C_4a,C_4b)
+ !if(DEBUG) call display_windows(nwin,iM,iL,iR)
- ! check window quality
- if(DATA_QUALITY) call check_window_s2n(nwin,iM,iL,iR)
- if (nwin.eq.0) then
- write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
- num_win = 0
- return
- endif
- if(DEBUG) call display_windows(nwin,iM,iL,iR)
+ ! REPEAT: remove small windows, since curtailing may have shrunk them
+ ! NOTE: perhaps this only needs to be done once (here)
+ if(DEBUG) print *, "Rejecting on size (REPEAT)"
+ call reject_on_window_width(nwin,iM,iL,iR,C_1)
+ if (nwin.eq.0) then
+ write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
+ num_win = 0
+ return
+ endif
+ if(DEBUG) call display_windows(nwin,iM,iL,iR)
- ! reject windows that do not satisfy F2, tshift and dlnA criteria
- call reject_on_fit_criteria(nwin,iM,iL,iR,CC_local,Tshift_local, dlnA_local)
- if (nwin.eq.0) then
- write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
- num_win = 0
- return
- endif
- if(DEBUG) call display_windows_full(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
+ ! check window quality
+ if(DATA_QUALITY) call check_window_s2n(nwin,iM,iL,iR)
+ if (nwin.eq.0) then
+ write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
+ num_win = 0
+ return
+ endif
+ if(DEBUG) call display_windows(nwin,iM,iL,iR)
- ! now that all criteria are satisfied, reject any duplicate windows
- call reject_on_duplicate(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
- if(DEBUG) call display_windows_full(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
+ ! reject windows that do not satisfy F2, tshift and dlnA criteria
+ call reject_on_fit_criteria(nwin,iM,iL,iR,CC_local,Tshift_local, dlnA_local)
+ if (nwin.eq.0) then
+ write(*,*) 'NO WINDOWS SELECTED FOR THIS TRACE'
+ num_win = 0
+ return
+ endif
+ if(DEBUG) call display_windows_full(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
- ! resolve the overlaps
- call resolve_overlaps(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
- if(DEBUG) call display_windows_full(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
+ ! now that all criteria are satisfied, reject any duplicate windows
+ call reject_on_duplicate(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
+ if(DEBUG) call display_windows_full(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
- ! for each window
+ ! resolve the overlaps
+ call resolve_overlaps(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
+ if(DEBUG) call display_windows_full(nwin,iM,iL,iR,CC_local,Tshift_local,dlnA_local)
+!!$ endif
+
! set global number of windows
num_win = nwin
do k = 1,num_win
@@ -209,7 +229,9 @@
pow_ratio = signal_pow / noise_pow
if(DEBUG) then
- write (*,*) 'DEBUG : iwin, amp_ratio : ', iwin, amp_ratio
+ write (*,*) 'DEBUG : amp_ratio,signal_amp,noise_amp : ',amp_ratio,signal_amp,noise_amp
+ write (*,*) 'DEBUG : pow_ratio,signal_pow,noise_pow : ',pow_ratio,signal_pow,noise_pow
+ write (*,*) 'DEBUG : iwin, amp_ratio : ', iwin, amp_ratio, S2N_LIMIT(iM(iwin))
endif
if (amp_ratio < S2N_LIMIT(iM(iwin))) then
@@ -236,7 +258,7 @@
end subroutine check_window_s2n
- subroutine check_data_quality(data_ok,signal_int,noise_int,snr_int,time_obs_signal,time_obs_noise,signal_amp,noise_amp,snr_amp)
+ subroutine check_data_quality(data_ok,signal_int,noise_int,snr_int,signal_amp,noise_amp,snr_amp)
use seismo_variables
! data_ok is a logical value to indicate whether or not
@@ -245,9 +267,10 @@
! and that particular station will be rejected.
logical, intent(out) :: data_ok
- double precision :: time_obs_noise, time_obs_signal, noise_int, signal_int, snr_int, signal_amp, noise_amp, snr_amp
+ double precision, intent(out) :: signal_int,noise_int,snr_int,signal_amp,noise_amp,snr_amp
- integer :: i, j
+ double precision :: time_obs_noise, time_obs_signal
+ integer :: i, j, k
! initialize values
data_ok = .true.
@@ -258,39 +281,43 @@
noise_amp = 0.0
snr_amp = 0.0
- if(DEBUG) write(*,*) 'DEBUG : noise_end = ', sngl(noise_end)
- if(DEBUG) write(*,*) 'DEBUG : signal_end = ', sngl(signal_end)
+ if(DEBUG) then
+ write(*,*) 'DEBUG : noise_start = ', sngl(noise_start)
+ write(*,*) 'DEBUG : noise_end = ', sngl(noise_end)
+ write(*,*) 'DEBUG : signal_start = ', sngl(signal_start)
+ write(*,*) 'DEBUG : signal_end = ', sngl(signal_end)
+ endif
- ! Integrating the signal and noise. We will take the RMS of the signal
- ! and noise, and take the ratio. Noise is integrated from start time to
- ! time noise_end, and signal is integrated from start time
- ! to time signal_end. We also reject based on the amplitude
- ! ratio of the largest data point before noise_end and the
- ! largest data point between noise_end and signal_end.
-
- ! These two loops prepare noise, noise_amp, signal, signal_amp.
- ! noise_end and signal_end are defined in user_functions.f90
- ! (set_up_criteria_arrays).
-
- ! noise loop
- do i = 1, NDIM
+ ! compute noise
+ ! see user_functions.f90 (default: noise_start = b)
+ do m = 1, npts
+ time_obs_noise = b+(m-1)*dt
+ if (time_obs_noise > noise_start) exit
+ enddo
+ in1 = m-1 ! index corresponding to noise_start
+ do i = 1, npts
time_obs_noise = b+(i-1)*dt
- if (time_obs_noise > noise_end) exit
+ if (time_obs_noise > noise_end ) exit
enddo
- ! i is the index corresponding to noise_end
- noise_int = sum( (obs_lp(1:i))**2 )/(i-1)
- noise_amp = maxval( abs(obs_lp(1:i)) )
+ in2 = i-1 ! index corresponding to noise_end
+ noise_int = sum( (obs_lp(in1:in2))**2 )/(in2-in1)
+ noise_amp = maxval( abs(obs_lp(in1:in2)) )
- ! signal loop
- do j = 1, NDIM
+ ! compute signal
+ ! see user_functions.f90 (default: signal_start = noise_end)
+ do j = 1, npts
time_obs_signal = b+(j-1)*dt
+ if (time_obs_signal > signal_start) exit
+ enddo
+ is1 = j-1 ! index corresponding to signal_start
+ do k = 1, npts
+ time_obs_signal = b+(k-1)*dt
if (time_obs_signal > signal_end) exit
enddo
- ! j is the signal corresponding to signal_end
+ is2 = k-1 ! index corresponding to signal_end
+ signal_int = sum( (obs_lp(is1:is2))**2 )/(is2-is1)
+ signal_amp = maxval( abs(obs_lp(is1:is2)) )
- signal_int = sum( (obs_lp(i:j))**2 )/(j-i)
- signal_amp = maxval( abs(obs_lp(i:j)) )
-
! Calculate signal to noise ratio and amplitude ratio.
snr_int = signal_int / noise_int
snr_amp = signal_amp / noise_amp
@@ -306,15 +333,16 @@
endif
if (DEBUG) then
+ print *, 'DEBUG : in1,in2,is1,is2,npts :', in1,in2,is1,is2,npts
print *, 'DEBUG : data_ok :', data_ok
- print *, 'DEBUG : signal_int : ', signal_int
- print *, 'DEBUG : noise_int :', noise_int
- print *, 'DEBUG : snr_int :', snr_int
- print *, 'DEBUG : signal_amp : ', signal_amp
- print *, 'DEBUG : noise_amp :', noise_amp
- print *, 'DEBUG : snr_amp :', snr_amp
- print *, 'DEBUG : SNR_INTEGRATE_BASE :', SNR_INTEGRATE_BASE
- print *, 'DEBUG : SNR_MAX_BASE :', SNR_MAX_BASE
+ print *, 'DEBUG : signal_int : ', sngl(signal_int)
+ print *, 'DEBUG : noise_int :', sngl(noise_int)
+ print *, 'DEBUG : snr_int :', sngl(snr_int)
+ print *, 'DEBUG : signal_amp : ', sngl(signal_amp)
+ print *, 'DEBUG : noise_amp :', sngl(noise_amp)
+ print *, 'DEBUG : snr_amp :', sngl(snr_amp)
+ print *, 'DEBUG : SNR_INTEGRATE_BASE :', sngl(SNR_INTEGRATE_BASE)
+ print *, 'DEBUG : SNR_MAX_BASE :', sngl(SNR_MAX_BASE)
endif
end subroutine check_data_quality
@@ -801,13 +829,11 @@
integer, dimension(*), intent(inout) :: iM, iL, iR
double precision, dimension(*), intent(out) :: CC_local, Tshift_local, dlnA_local
- integer :: iwin
- integer :: nwin_new
+ integer :: iwin, nwin_new
integer, dimension(NWINDOWS) :: iM_new, iL_new, iR_new
double precision :: CC_temp, Tshift_temp, dlnA_temp
double precision :: tshift_min, tshift_max, dlnA_min, dlnA_max
logical :: accept
-
nwin_new = 0
! for each proto-window, check that all included maxima are part of the
Added: seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T003_T010_m16
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T003_T010_m16 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T003_T010_m16 2009-10-25 23:06:34 UTC (rev 15876)
@@ -0,0 +1,71 @@
+# -------------------------------------------------------------
+#
+# This is the parameter file for FLEXWIN. It is based on the
+# same syntax as the Par_file for SPECFEM. Variable names are
+# put first, values are placed after the 34th column.
+#
+# Comment lines and blank lines are significant. If you
+# change the layout of this file or add/remove parameters
+# you must also modify the user_variables module and the
+# read_parameter_file subroutine at the start of seismo_subs.f90.
+#
+# -------------------------------------------------------------
+
+# -------------------------------------------------------------
+# boolean parameters
+DEBUG = .true.
+MAKE_SEISMO_PLOTS = .true.
+MAKE_WINDOW_FILES = .true.
+BODY_WAVE_ONLY = .true.
+
+# -------------------------------------------------------------
+# period min/max for filtering
+RUN_BANDPASS = .false.
+WIN_MIN_PERIOD = 3.00
+WIN_MAX_PERIOD = 30.00
+
+# -------------------------------------------------------------
+# E(t) water level
+STALTA_BASE = 0.11
+
+# -------------------------------------------------------------
+# maximum allowable time shift from reference TSHIFT
+TSHIFT_BASE = 100.0
+TSHIFT_REFERENCE = 0.0
+
+# -------------------------------------------------------------
+# maximum allowable amplitude measurement relative to reference DLNA
+DLNA_BASE = 2.25
+DLNA_REFERENCE = 2.00
+
+# -------------------------------------------------------------
+# limit on CC for window acceptance
+CC_BASE = 0.0
+
+# -------------------------------------------------------------
+# boolean switch for check_data_quality
+DATA_QUALITY = .true.
+
+# if DATA_QUALITY = .true. and if two different measurements of
+# signal-to-noise ratios exceeds these two base levels,
+# then the data time series (and syn) is kept
+SNR_INTEGRATE_BASE = 8.0
+SNR_MAX_BASE = 4.0
+
+# -------------------------------------------------------------
+# limit on signal to noise ratio in a particular window.
+WINDOW_SNR_BASE = 4.0
+
+# -------------------------------------------------------------
+# Fine tuning constants
+C_0 (internal minima) = 1.3
+C_1 (small windows) = 4.0
+C_2 (prominence) = 0.0
+C_3a (separation height) = 4.0
+C_3b (separation time) = 2.5
+C_4a (curtail on left) = 2.0
+C_4b (curtail on right) = 6.0
+
+WEIGHT_SPACE_COVERAGE = 0.70
+WEIGHT_AVERAGE_CC = 0.25
+WEIGHT_N_WINDOWS = 0.05
Modified: seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m00.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m00.f90 2009-10-25 20:44:01 UTC (rev 15875)
+++ seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m00.f90 2009-10-25 23:06:34 UTC (rev 15876)
@@ -74,8 +74,11 @@
endif
! variables for signal to noise ratio criteria.
- signal_end = Sw_end
- noise_end = Pnl_start
+ noise_start = b
+ noise_end = Pnl_start
+ signal_start = noise_end
+ signal_end = Sw_end
+
if(DEBUG) then
if (BODY_WAVE_ONLY) then
write(*,*) 'DEBUG : P_pick = ', P_pick
Modified: seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m16.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m16.f90 2009-10-25 20:44:01 UTC (rev 15875)
+++ seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m16.f90 2009-10-25 23:06:34 UTC (rev 15876)
@@ -74,8 +74,11 @@
endif
! variables for signal to noise ratio criteria.
- signal_end = Sw_end
- noise_end = Pnl_start
+ noise_start = b
+ noise_end = Pnl_start
+ signal_start = noise_end
+ signal_end = Sw_end
+
if(DEBUG) then
if (BODY_WAVE_ONLY) then
write(*,*) 'DEBUG : P_pick = ', P_pick
Added: seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m16_man.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m16_man.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m16_man.f90 2009-10-25 23:06:34 UTC (rev 15876)
@@ -0,0 +1,159 @@
+! -------------------------------------------------------------
+! edit here to change T0 and T1 on some condition
+! Note, this function is called AFTER the seismogram has been
+! read but before it is filtered.
+! -------------------------------------------------------------
+
+subroutine modify_T0_T1_on_condition
+ use seismo_variables
+
+ ! do nothing
+
+ ! adjust fstart and fend accordingly
+ !FSTART=1./WIN_MAX_PERIOD
+ !FEND=1./WIN_MIN_PERIOD
+
+end subroutine modify_T0_T1_on_condition
+
+! -------------------------------------------------------------
+! Edit here to change the time dependent properties of the selection criteria
+! Note, this function is called AFTER the seismogram has been read.
+! -------------------------------------------------------------
+subroutine set_up_criteria_arrays
+ use seismo_variables
+
+ integer :: i
+ double precision :: time
+
+ double precision :: Pnl_start, S_end, Sw_start, Sw_end
+ double precision :: Nlam, dtresh, vref
+
+!===========================
+
+! -----------------------------------------------------------------
+! This is the basic version of the subroutine - no variation with time
+! -----------------------------------------------------------------
+ do i = 1, npts
+ time = b+(i-1)*dt
+ DLNA_LIMIT(i) = DLNA_BASE
+ CC_LIMIT(i) = CC_BASE
+ TSHIFT_LIMIT(i) = TSHIFT_BASE ! WIN_MIN_PERIOD/2.0
+ STALTA_W_LEVEL(i) = STALTA_BASE
+ S2N_LIMIT(i) = WINDOW_S2N_BASE
+ enddo
+
+!!$ if (.not. BODY_WAVE_ONLY) then
+!!$ Pnl_start = -5.0 + dist_km/7.8
+!!$ Sw_start = -15.0 + dist_km/3.5
+!!$ Sw_end = 35.0 + dist_km/3.1
+!!$ else
+!!$ Pnl_start = P_pick - 5.0
+!!$ S_end = S_pick + 5.0
+!!$ Sw_start = -15.0 + dist_km/3.5
+!!$ Sw_end = 35.0 + dist_km/3.1
+!!$ endif
+
+ ! regional (Qinya's formulation):
+ ! -------------------------------------------------------------
+ ! see Liu et al. (2004), p. 1755, but note that the PARENTHESES
+ ! that are listed in the publication should not be there
+ ! THESE ARE PROBABLY NOT ACCURATE ENOUGH FOR LONGER PATHS.
+
+ !Sw_start = -15.0 + dist_km/3.5
+ !Sw_end = 35.0 + dist_km/3.1
+ Sw_end = 35.0 + dist_km/3.1
+
+ if (BODY_WAVE_ONLY) then
+ Pnl_start = P_pick - 2.5*WIN_MIN_PERIOD
+ else
+ !Pnl_start = -5.0 + dist_km/7.8
+ Pnl_start = -10.0 + dist_km/7.8
+ endif
+
+ ! variables for signal to noise ratio criteria.
+ !noise_start = b
+ !noise_end = Pnl_start
+ !signal_start = noise_end
+ !signal_end = Sw_end
+ noise_start = b
+ noise_end = Pnl_start
+ signal_start = Sw_end + 50.0
+ signal_end = signal_start + 1000.0 ! (end of record)
+
+ if(DEBUG) then
+ if (BODY_WAVE_ONLY) then
+ write(*,*) 'DEBUG : P_pick = ', P_pick
+ write(*,*) 'DEBUG : S_pick = ', S_pick
+ endif
+ !write(*,*) 'DEBUG : noise_start = ', sngl(noise_start)
+ !write(*,*) 'DEBUG : noise_end = ', sngl(noise_end)
+ !write(*,*) 'DEBUG : signal_start = ', sngl(signal_start)
+ !write(*,*) 'DEBUG : signal_end = ', sngl(signal_end)
+ endif
+
+ ! --------------------------------
+ ! modulate criteria in time
+ do i = 1, npts
+ time = b+(i-1)*dt ! time
+
+ ! raises STA/LTA water level before P wave arrival.
+ if(time.lt. signal_start ) then
+ STALTA_W_LEVEL(i) = 10.*STALTA_BASE
+ endif
+
+ ! raises STA/LTA water level before P wave arrival.
+ if(time.gt. signal_end ) then
+ STALTA_W_LEVEL(i) = 10.*STALTA_BASE
+ endif
+
+ enddo
+
+ ! --------------------------------
+ ! if the distance to the station is less than N wavelengths, then reject records
+ ! by reasing the entire water level
+
+ Nlam = 0.0 ! number of wavelengths
+ vref = 2.0 ! reference velocity, km/s
+ dtresh = Nlam*WIN_MIN_PERIOD*vref
+ if (dist_km .le. dtresh ) then
+ if(DEBUG) then
+ write(*,*) 'REJECT by raising water level: station is too close for this period range'
+ write(*,*) 'dist_km, dtresh = Nlam*WIN_MIN_PERIOD, Nlam, WIN_MIN_PERIOD :'
+ write(*,'(4f12.4)') dist_km, dtresh, Nlam, WIN_MIN_PERIOD
+ endif
+ do i = 1,npts
+ STALTA_W_LEVEL(i) = 10.*STALTA_BASE
+ enddo
+ endif
+
+! The following is for check_window quality_s2n
+
+! -----------------------------------------------------------------
+! Start of user-dependent portion
+
+! This is where you modulate the time dependence of the selection
+! criteria. You have access to the following parameters from the
+! seismogram itself:
+!
+! dt, b, kstnm, knetwk, kcmpnm
+! evla, evlo, stla, stlo, evdp, azimuth, backazimuth, dist_deg, dist_km
+! num_phases, ph_names, ph_times
+!
+! Example of modulation:
+!-----------------------
+! To increase s2n limit after arrival of R1 try
+!
+! R_vel=3.2
+! R_time=dist_km/R_vel
+! do i = 1, npts
+! time=b+(i-1)*dt
+! if (time.gt.R_time) then
+! S2N_LIMIT(i)=2*WINDOW_S2N_BASE
+! endif
+! enddo
+!
+! End of user-dependent portion
+! -----------------------------------------------------------------
+
+end subroutine set_up_criteria_arrays
+! -------------------------------------------------------------
Modified: seismo/3D/ADJOINT_TOMO/flexwin/user_functions.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/user_functions.f90 2009-10-25 20:44:01 UTC (rev 15875)
+++ seismo/3D/ADJOINT_TOMO/flexwin/user_functions.f90 2009-10-25 23:06:34 UTC (rev 15876)
@@ -47,7 +47,9 @@
! these values will be used for signal2noise calculations
! if DATA_QUALITY=.true.
if (DATA_QUALITY) then
+ noise_start=b
noise_end=ph_times(1)-WIN_MAX_PERIOD
+ signal_start=noise_end
signal_end=b+(npts-1)*dt
endif
More information about the CIG-COMMITS
mailing list