[cig-commits] r12405 - in seismo/3D/ADJOINT_TOMO/flexwin: . scripts scripts/prepare_scripts/socal user_files/socal_3D
carltape at geodynamics.org
carltape at geodynamics.org
Sun Jul 13 16:15:26 PDT 2008
Author: carltape
Date: 2008-07-13 16:15:26 -0700 (Sun, 13 Jul 2008)
New Revision: 12405
Added:
seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T002_m05
seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T006_m05
seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m05.f90
Modified:
seismo/3D/ADJOINT_TOMO/flexwin/Makefile.in
seismo/3D/ADJOINT_TOMO/flexwin/configure
seismo/3D/ADJOINT_TOMO/flexwin/configure.ac
seismo/3D/ADJOINT_TOMO/flexwin/make_intel
seismo/3D/ADJOINT_TOMO/flexwin/make_intel_caltech
seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows.pl
seismo/3D/ADJOINT_TOMO/flexwin/scripts/prepare_scripts/socal/process_data_and_syn.pl
Log:
Changed the names in the make files and configure files from test2 to flexwin. The name of the main program had already been changed from test2.f90 to flexwin.f90. Also added parameter files for m05 socal.
Modified: seismo/3D/ADJOINT_TOMO/flexwin/Makefile.in
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/Makefile.in 2008-07-09 22:36:02 UTC (rev 12404)
+++ seismo/3D/ADJOINT_TOMO/flexwin/Makefile.in 2008-07-13 23:15:26 UTC (rev 12405)
@@ -5,10 +5,10 @@
LDFLAGS = @LDFLAGS@
LIBS = @LIBS@
-PROGS = test2
+PROGS = flexwin
-test2_LIBS = -lSacLib -lsac -ltau
-test2_OBJS = \
+flexwin_LIBS = -lSacLib -lsac -ltau
+flexwin_OBJS = \
distaz.o \
io_subs.o \
maxima.o \
@@ -38,8 +38,8 @@
$(FC) $(FCFLAGS) -c -o $@ $<
-test2: test2.f90 $(test2_OBJS)
- $(FC) $(FCFLAGS) test2.f90 -o test2 $(test2_OBJS) $(LDFLAGS) $(test2_LIBS) $(LIBS)
+flexwin: flexwin.f90 $(flexwin_OBJS)
+ $(FC) $(FCFLAGS) flexwin.f90 -o flexwin $(flexwin_OBJS) $(LDFLAGS) $(flexwin_LIBS) $(LIBS)
clean:
rm -f $(PROGS) *.o *.mod
Modified: seismo/3D/ADJOINT_TOMO/flexwin/configure
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/configure 2008-07-09 22:36:02 UTC (rev 12404)
+++ seismo/3D/ADJOINT_TOMO/flexwin/configure 2008-07-13 23:15:26 UTC (rev 12405)
@@ -578,7 +578,7 @@
PACKAGE_STRING='flexwin 1.0'
PACKAGE_BUGREPORT='cig-seismo at geodynamics.org'
-ac_unique_file="test2.f90"
+ac_unique_file="flexwin.f90"
ac_subst_vars='SHELL
PATH_SEPARATOR
PACKAGE_NAME
Modified: seismo/3D/ADJOINT_TOMO/flexwin/configure.ac
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/configure.ac 2008-07-09 22:36:02 UTC (rev 12404)
+++ seismo/3D/ADJOINT_TOMO/flexwin/configure.ac 2008-07-13 23:15:26 UTC (rev 12405)
@@ -11,7 +11,7 @@
AC_PREREQ(2.61)
AC_INIT([flexwin], [1.0], [cig-seismo at geodynamics.org])
-AC_CONFIG_SRCDIR([test2.f90])
+AC_CONFIG_SRCDIR([flexwin.f90])
# Checks for programs.
Modified: seismo/3D/ADJOINT_TOMO/flexwin/make_intel
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/make_intel 2008-07-09 22:36:02 UTC (rev 12404)
+++ seismo/3D/ADJOINT_TOMO/flexwin/make_intel 2008-07-13 23:15:26 UTC (rev 12405)
@@ -13,7 +13,7 @@
SHARE=${A_CODE}/share
LIBS= -lSacLib -lsac -lutil -ltau -lnumrec -lm
-PROGS= test2
+PROGS= flexwin
default:${PROGS}
LIBSRC2=user_parameters.f90 seismo_subs.f90 user_functions.f90 maxima.f90 travel_times.f90 measurement_module.f90 xcorr-measure.f90 measure_windows_xcorr.f90 io_subs.f90 select_windows_stalta2.f90
@@ -36,8 +36,8 @@
$(FC) ${FFLAGS} -c -o $@ $< -I${INCDIR}
-test2:test2.f90 ${LIB2}
- $(FC) ${FFLAGS} test2.f90 -o test2 ${LIB2} -L${LIBDIR} -L${OPTLIBDIR} ${LIBS}
+flexwin:flexwin.f90 ${LIB2}
+ $(FC) ${FFLAGS} flexwin.f90 -o flexwin ${LIB2} -L${LIBDIR} -L${OPTLIBDIR} ${LIBS}
if [ ! -e iasp91.hed ] ; then ( ln -s ${SHARE}/iasp91.hed . ; ln -s ${SHARE}/iasp91.tbl . ) ; fi
clean:
Modified: seismo/3D/ADJOINT_TOMO/flexwin/make_intel_caltech
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/make_intel_caltech 2008-07-09 22:36:02 UTC (rev 12404)
+++ seismo/3D/ADJOINT_TOMO/flexwin/make_intel_caltech 2008-07-13 23:15:26 UTC (rev 12405)
@@ -15,7 +15,7 @@
SHARE=${A_CODE}/share
LIBS= -lSacLib -lsac -lutil -ltau -lnumrec -lm
-PROGS= test2
+PROGS= flexwin
default:${PROGS}
LIBSRC2=user_parameters.f90 seismo_subs.f90 user_functions.f90 select_windows_stalta2.f90 maxima.f90 travel_times.f90 measurement_module.f90 xcorr-measure.f90 measure_windows_xcorr.f90 io_subs.f90
@@ -38,8 +38,8 @@
$(FC) ${FFLAGS} -c -o $@ $< -I${INCDIR}
-test2:test2.f90 ${LIB2}
- $(FC) ${FFLAGS} test2.f90 -o test2 ${LIB2} -L${LIBDIR} -L${OPTLIBDIR} ${LIBS}
+flexwin:flexwin.f90 ${LIB2}
+ $(FC) ${FFLAGS} flexwin.f90 -o flexwin ${LIB2} -L${LIBDIR} -L${OPTLIBDIR} ${LIBS}
if [ ! -e iasp91.hed ] ; then ( ln -s ${SHARE}/iasp91.hed . ; ln -s ${SHARE}/iasp91.tbl . ) ; fi
clean:
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows.pl 2008-07-09 22:36:02 UTC (rev 12404)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows.pl 2008-07-13 23:15:26 UTC (rev 12405)
@@ -38,8 +38,8 @@
# prepare_meas_all.pl
#
# EXAMPLES:
-# pick_all_windows.pl m04 0 6/40 1/203 1/1/1/0 1 0 # both, T = 6s
-# pick_all_windows.pl m04 0 2/40 1/203 1/1/1/1 1 0 # both, T = 2s
+# pick_all_windows.pl m05 0 6/40 1/203 1/1/1/0 1 0 # both, T = 6s
+# pick_all_windows.pl m05 0 2/40 1/203 1/1/1/1 1 0 # both, T = 2s
#
# pick_all_windows.pl m00 0 6/40 179/179 1/1/0/0 1 0 # plots only
# pick_all_windows.pl m00 0 6/40 179/179 0/0/1/0 1 0 # measurement code only
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 2008-07-09 22:36:02 UTC (rev 12404)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/prepare_scripts/socal/process_data_and_syn.pl 2008-07-13 23:15:26 UTC (rev 12405)
@@ -3,7 +3,7 @@
#==========================================================
#
# Carl Tape
-# 03-July-2007
+# 12-July-2007
# process_data_and_syn.pl
#
# This script processes data and 3D-Socal-SEM synthetics for Southern California.
@@ -19,30 +19,34 @@
# dat_ext extension for processed data
# Trange band-pass filter
#
-# EXAMPLES:
-# ~/UTILS/process_data_and_syn.pl 0 1 d d 20 6/40 PROCESSED # synthetics
-# ~/UTILS/process_data_and_syn.pl 1 0 d d 20 6/40 PROCESSED # data
-# ~/UTILS/process_data_and_syn.pl 1 1 d d 20 6/40 PROCESSED # both
+# ORDER OF OPERATIONS:
+# ~/UTILS/process_data_and_syn.pl 1 1 0 d d 20 6/40 PROCESSED # data, initial processing
+# ~/UTILS/process_data_and_syn.pl 1 0 1 d d 20 6/40 PROCESSED # syn, initial processing
+# ~/UTILS/process_data_and_syn.pl 2 1 1 d d 20 6/40 PROCESSED # both, create cut files
+# ~/UTILS/process_data_and_syn.pl 3 1 0 d d 20 6/40 PROCESSED # data, execute cut file
+# ~/UTILS/process_data_and_syn.pl 3 0 1 d d 20 6/40 PROCESSED # syn, execute cut file
#
-# ~/UTILS/process_data_and_syn.pl 0 1 d d 20 0.2/100 PROCESSED # synthetics
-# ~/UTILS/process_data_and_syn.pl 1 0 d d 20 0.2/100 PROCESSED # data
-# ~/UTILS/process_data_and_syn.pl 1 1 d d 20 0.2/100 PROCESSED # both
+# ~/UTILS/process_data_and_syn.pl 4 1 0 d d 20 6/40 PROCESSED # data, bandpass T=6-40
+# ~/UTILS/process_data_and_syn.pl 4 0 1 d d 20 6/40 PROCESSED # syn, bandpass T=6-40
#
+# ~/UTILS/process_data_and_syn.pl 4 1 0 d d 20 2/40 PROCESSED # data, bandpass T=2-40
+# ~/UTILS/process_data_and_syn.pl 4 0 1 d d 20 2/40 PROCESSED # syn, bandpass T=2-40
+#
#==========================================================
-if (@ARGV < 7) {die("Usage: process_data_and_syn.pl idata isyn syn_ext dat_ext sps Tmin/Tmax \n")}
-($idata,$isyn,$syn_ext,$dat_ext,$sps,$Trange,$pdir) = @ARGV;
+if (@ARGV < 8) {die("Usage: process_data_and_syn.pl idata isyn syn_ext dat_ext sps Tmin/Tmax pdir\n")}
+($iprocess,$idata,$isyn,$syn_ext,$dat_ext,$sps,$Trange,$pdir) = @ARGV;
-$smodel = "m04"; # KEY: model iteration index
+$smodel = "m05"; # KEY: model iteration index
$iexecute = 0;
$ilist = 1;
$iprocess0 = 0;
# THESE MUST BE DONE IN SEQUENTIAL ORDER
-$iprocess1 = 0; # only done ONCE, no matter how many band-pass filters you use
-$iprocess2 = 1; # cut records and pad zeros (and bandpass, if iprocess3=1)
-$iprocess3 = 1; # bandpass filter (iprocess2=1 also)
+#$iprocess1 = 0; # only done ONCE, no matter how many band-pass filters you use
+#$iprocess2 = 1; # cut records and pad zeros (and bandpass, if iprocess3=1)
+#$iprocess3 = 0; # bandpass filter (iprocess2=1 also)
#$Lrange = "-l -40/220";
#$Lrange = " ";
@@ -58,8 +62,8 @@
# directories
#$CMT_list = "/net/sierra/raid1/carltape/socal/socal_3D/SYN/model_${smodel}";
-$dir_source = "/net/sierra/raid1/carltape/results/SOURCES/socal_6";
-$dirCMT = "${dir_source}/v6_files";
+$dir_source = "/net/sierra/raid1/carltape/results/SOURCES/socal_7";
+$dirCMT = "${dir_source}/v7_files";
$CMT_list = "${dir_source}/EIDs_only_eid";
$dirdat0 = "/net/sierra/raid1/carltape/socal/socal_3D/DATA/FINAL";
$dirsyn0 = "/net/sierra/raid1/carltape/socal/socal_3D/SYN/model_${smodel}";
@@ -114,14 +118,14 @@
print "\nWriting to $cshfile ...\n";
open(CSH,">$cshfile");
-if($iprocess0==1) {
+if($iprocess==0) {
$sfile = "nprocess_syn";
open(SYN,">$sfile");
}
$imin = 1; $imax = $ncmt; # default
#$imin = 1; $imax = 96;
-#$imin = 174; $imax = $imin;
+#$imin = 178; $imax = $imin;
#----------------------------------------------------------------------
@@ -149,9 +153,11 @@
$dirdat_pro_2 = "${dirdat_pro_1}/$pdirbpass";
$dirsyn_pro_2 = "${dirsyn_pro_1}/$pdirbpass";
- # cut times file
- $cutfile = "$dirdat/${eid}_cut";
- $cutfile_done = "${cutfile}_done";
+ # cut times files
+ $cutfile_dat = "$dirdat/${eid}_dat_cut";
+ $cutfile_dat_done = "${cutfile_dat}_done";
+ $cutfile_syn = "$dirdat/${eid}_syn_cut"; # note: data directory
+ $cutfile_syn_done = "$dirsyn/${eid}_syn_cut_done"; # note: syn directory
# optional -- delete pre-processed directories
#print CSH "rm -rf $dirsyn/PRO*\n";
@@ -160,7 +166,7 @@
#----------------------------------------------------------------------
# PROCESSING PART 0: check the number of processed synthetic files for each event
- if ($iprocess0 == 1) {
+ if ($iprocess == 0) {
if(-e ${dirsyn_pro_1}) {
($nfile,undef,undef) = split(" ",`ls -1 ${dirsyn_pro_1}/* | wc`);
print SYN "$eid $nfile\n";
@@ -170,7 +176,7 @@
#----------------------------------------------------------------------
# PROCESSING PART 1: assign SAC headers, interpolate, and pick P and S arrivals (based on a 1D socal model)
- if ($iprocess1 == 1) {
+ if ($iprocess == 1) {
# synthetics -- this will convolve with the source half-duration (prior to interpolating)
if ($isyn == 1) {
@@ -204,409 +210,304 @@
}
} # idata
- # # process synthetics
- # if ($isyn == 1) {
- # $dirsyn = "${dirsyn0}/${eid}";
- # #if (not -e $dirsyn) {die("check if $dirsyn exist or not\n")}
- # if (-e $dirsyn) {
- # if (not -e "$dirsyn/$pdir") {
- # print CSH "cd $dirsyn\n";
- # print CSH "\\rm -rf $pdir\n";
- # #print CSH "~/UTILS/interp_records.pl -S -t 0.011 *.semd.sac \n";
- # #print CSH "process_trinet_syn_new.pl -S -m $cmtfile -a $stafile -p -t $Trange -s 20 -d $pdir -x ${syn_ext} *.semd.sac \n";
- # print CSH "process_trinet_syn_new.pl -S -m $cmtfile -h -a $stafile -p -d $pdir -x ${syn_ext} *.semd.sac \n";
- # print CSH "cd $pdir\n";
- # print CSH "rotate.pl *E.semd.sac.${syn_ext} \n";
+ } # iprocess = 1
- # #print CSH "process_trinet_syn_new.pl -m $cmtfile -a $stafile *.semd \n";
- # #print CSH "\\rm *sac.d *BHR* *BHT* \n";
- # #print CSH "process_trinet_syn_new.pl -S -s 20 *.sac \n";
- # }
- # } else {
- # print "\n $dirsyn does not exist";
- # }
- # } # isyn
-
- # # process data and make a list of files for which we are missing the PZ
- # # At the end, from CHECK_DIR : cat *pz* > no_PoleZero_files
- # if ($idata == 1) {
- # $dirdat = "${dirdat0}/${eid}";
- # if (-e $dirdat) {
- # if (not -e "$dirdat/$pdir") {
-
- # $ofile1 = "${eid}_ofile";
- # $ofile2 = "${eid}_no_pz_files";
- # $odir = "${dir_dat}/CHECK_DIR";
- # print CSH "mkdir -p $odir\n";
- # print CSH "cd $dirdat\n";
- # print CSH "\\rm -rf $pdir; mkdir $pdir\n";
-
- # # note -I to convert to displacement but NOT subsequently bandpass
- # #print CSH "process_cal_data.pl -m $cmtfile -i none -p -t $Trange -s 20 -d $pdir -x ${dat_ext} *.sac > $ofile1 \n";
- # print CSH "process_cal_data.pl -m $cmtfile -I none -p -t $Trange -d $pdir -x ${dat_ext} *.sac > $ofile1 \n";
-
- # print CSH "grep skip $ofile1 | awk '{print \$7}' > $ofile2 \n";
- # print CSH "\\cp $ofile1 $ofile2 $odir\n";
- # print CSH "cd $pdir\n";
- # print CSH "rotate.pl *E.sac.${dat_ext} \n";
- # }
- # } else {
- # print "\n $dirdat does not exist";
- # }
- # } # idata
-
- }
-
#----------------------------------------------------------------------
- # PROCESSING PART 2: cutting records, padding zeros, and bandpass
+ # PROCESSING PART 2: getting the cut times for the records
- if ($iprocess2 == 1) {
+ if ($iprocess == 2) {
- # BOTH the (initially) processed synthetics and data directories must exist
+ # BOTH the INITIALLY processed synthetics and data directories must exist,
+ # even if you only want to process synthetics.
if ( (not -e ${dirdat_pro_1}) || (not -e ${dirsyn_pro_1}) ) {
print "--> dirdat ${dirdat_pro_1} and dirsyn ${dirsyn_pro_1} do not both exist\n";
- } else {
+ } elsif ( ((-f $cutfile_syn) || (-f $cutfile_syn_done)) || ((-f $cutfile_dat) || (-f $cutfile_dat_done)) ) {
- if ( (-e ${dirdat_pro_2} && $idata==1) || (-e ${dirsyn_pro_2} && $isyn==1) ) {
- if ($idata == 1) {
- print "${dirdat_pro_2} already exists\n";
- }
- if ($isyn == 1) {
- print "${dirsyn_pro_2} already exists\n";
- }
+ if (-f $cutfile_syn) {
+ print "cutfile_syn ${cutfile_syn} already exists\n";
+ }
+ if (-f $cutfile_syn_done) {
+ print "cutfile_syn_done ${cutfile_syn_done} already exists\n";
+ }
+ if (-f $cutfile_dat) {
+ print "cutfile_dat ${cutfile_dat} already exists\n";
+ }
+ if (-f $cutfile_dat_done) {
+ print "cutfile_dat_done ${cutfile_dat_done} already exists\n";
+ }
+ print "--> you are probably ready for cutting or bandpassing...\n";
- } else {
+ } else {
- if (-f $cutfile_done) {
- print "cutfile ${cutfile_done} already exists -- now bandpass\n";
+ print "\nWriting to cutfiles ...\n";
+ open(CUTDAT,">${cutfile_dat}");
+ open(CUTSYN,">${cutfile_syn}");
- if ($iprocess3==0) {
- print "set iprocess3=1 if you want to bandpass\n";
+ # grab all DATA files
+ @files = glob("${dirdat_pro_1}/*");
+ $nfile = @files;
+ print "\n $nfile data files to line up with synthetics\n";
- } else {
- print "iprocess3=1 -- now bandpass\n";
+ foreach $datfile (@files) {
+ # read the sac headers -- network, station, component
+ (undef,$net,$sta,$chan) = split(" ",`saclst knetwk kstnm kcmpnm f $datfile`);
+ $comp = `echo $chan | awk '{print substr(\$1,3,1)}'`;
+ chomp($comp);
- if ($isyn == 1) {
- if (-e ${dirsyn_pro_1}) {
- if (not -e ${dirsyn_pro_2}) {
- print CSH "cd ${dirsyn_pro_1}\n";
- #print CSH "\\rm -rf $pdirbpass\n";
- print CSH "process_trinet_syn_new.pl -S -t $Trange -d $pdirbpass -x $Ttag *.${syn_suffix} \n";
- print CSH "cd $pdirbpass\n";
- print CSH "rotate.pl *E.${syn_suffix}.${Ttag}\n";
+ # synthetics are always BH_ component
+ $synfile = "${dirsyn_pro_1}/${sta}.${net}.BH${comp}.${syn_suffix}";
- } else {
- print "$pdirbpass already exists\n";
- }
- } else {
- print "${dirsyn_pro_1} does not exist\n";
- }
- } # isyn
+ # if the synthetic file exists, then go on
+ if (-f $synfile) {
+ print "$datfile $synfile\n"; # only list files that are pairs
- # data
- if ($idata == 1) {
- if (-e ${dirdat_pro_1}) {
- if (not -e ${dirdat_pro_2}) {
+ # get info on data and synthetics
+ (undef,$bd,$ed,$deltad,$nptd) = split(" ",`saclst b e delta npts f $datfile`);
+ (undef,$bs,$es,$deltas,$npts) = split(" ",`saclst b e delta npts f $synfile`);
+ $tlend = $ed - $bd;
+ $tlens = $es - $bs;
+
+ # dt should be the same for both records ALREADY
+ if (log($deltad/$deltas) > 0.01) {
+ print "$datfile $synfile\n";
+ print "DT values are not close enough: $deltad, $deltas\n";
+ die("fix the DT values\n");
+ } else {
+ $dt = $deltad;
+ }
- $ofile1 = "${eid}_ofile";
- $ofile2 = "${eid}_no_pz_files";
- $odir = "${dirdat0}/CHECK_DIR";
- print CSH "mkdir -p $odir\n";
+ # determine the cut for the records
+ # b = earliest start time
+ # e = earliest end time, multiplied by some factor
+ if ($bd < $bs) {
+ $b0 = $bd;
+ } else {
+ $b0 = $bs;
+ }
+ if ($b0 < $bmin) {
+ $b0 = $bmin;
+ }
+ if ($ed < $es) {
+ $e0 = $ed;
+ } else {
+ $e0 = $es;
+ }
- print CSH "cd ${dirdat_pro_1}\n";
- #print CSH "\\rm -rf $pdirbpass\n";
- print CSH "process_cal_data.pl -i none -t $Trange -d $pdirbpass -x $Ttag *.${dat_suffix} > $ofile1\n";
- print CSH "grep skip $ofile1 | awk '{print \$7}' > $ofile2\n";
- print CSH "\\cp $ofile1 $ofile2 $odir\n";
+ $b = $b0;
+ $tlen0 = $e0 - $b;
+ $tlen = $tfac * $tlen0; # extend record length (if desired)
+ $e = $b0 + $tlen;
+ $npt = int( ($e-$b)/$dt );
- print CSH "cd $pdirbpass\n";
- print CSH "rotate.pl *E.${dat_suffix}.${Ttag}\n";
+ #print CUT "$datfile $synfile $b $e $npt $dt\n";
+ print CUTDAT "$datfile $b $e $npt $dt\n";
+ print CUTSYN "$synfile $b $e $npt $dt\n";
- } else {
- print "$pdirbpass already exists\n";
- }
- } else {
- print "${dirdat_pro_1} does not exist\n";
- }
- } # idata
+ if (0==1) {
+ print "\n Data : $bd $ed $deltad $nptd -- $tlend";
+ print "\n Syn : $bs $es $deltas $npts -- $tlens";
+ print "\n b0, e0, tlen0 : $b0, $e0, $tlen0 ";
+ print "\n b : $bd, $bs, $b ";
+ print "\n e : $ed, $es, $e ";
+ print "\n npt : $nptd, $npts, $npt ";
+ print "\n $tlen = $tfac * ($e0 - $b)";
+ print "\n $e = $b0 + $tlen \n";
+ }
+ } # if synfile exists
+ } # for all data files
+ print "\n Done making cutfile $cutfile_dat\n";
+ print "\n Done making cutfile $cutfile_syn\n";
+ close(CUTDAT);
+ close(CUTSYN);
+
+ } # dirdat_pro_1 and dirsyn_pro_1 exist
+ } # iprocess=2
- } # iprocess3==1
+ #----------------------------------------------------------------------
+ # PROCESSING PART 3: cutting records and padding zeros
- } elsif (-f $cutfile) {
- print "cutfile $cutfile exists -- now read it in sac and execute\n";
+ if ($iprocess == 3) {
- # read cut file
- open(IN,"$cutfile"); @lines = <IN>; close(IN); $nlines = @lines;
+ if ($idata == 1) {
+ if (-e ${dirdat_pro_2}) {
+ print "--> dirdat ${dirdat_pro_2} already exists\n";
- $sacfile = "sac.mac";
- `echo echo on > $sacfile`;
- `echo readerr badfile fatal >> $sacfile`;
+ } else {
+ if (-f $cutfile_dat_done) {
+ print "cutfile_dat_done ${cutfile_dat_done} already exists\n";
- for ($j = 1; $j <= $nlines; $j++) {
-
- $line = $lines[$j-1]; chomp($line);
- ($datfile,$synfile,$b,$e,$npt,$dt) = split(" ",$line);
- print "$j out of $nlines -- $datfile\n";
- #print "-- $datfile -- $synfile -- $b -- $e -- $npt -- $dt -- \n";
-
- # cut records and fill zeros
- `echo r $datfile $synfile >> $sacfile`;
- `echo cuterr fillz >> $sacfile`;
- `echo "cut $b n $npt" >> $sacfile`;
- `echo r $datfile $synfile >> $sacfile`; # cut both
- `echo cut off >> $sacfile`;
- `echo w over >> $sacfile`;
-
- }
- `echo quit >> $sacfile`;
-
- # KEY: execute SAC command
- `sac $sacfile`;
- `sleep 5s`;
- `rm $sacfile`;
- print "\n Done with pre-processing, part 2\n";
- `mv $cutfile ${cutfile_done}`;
-
} else {
- print "\nWriting to $cutfile ...\n";
- open(CUT,">$cutfile");
+ if (not -f $cutfile_dat) {
+ print "cutfile_dat $cutfile_dat does not exist -- try iprocess = 2\n";
- # grab all DATA files
- @files = glob("${dirdat_pro_1}/*");
- $nfile = @files;
- print "\n $nfile data files to line up with synthetics\n";
+ } else {
- foreach $datfile (@files) {
- # read the sac headers -- network, station, component
- (undef,$net,$sta,$chan) = split(" ",`saclst knetwk kstnm kcmpnm f $datfile`);
- $comp = `echo $chan | awk '{print substr(\$1,3,1)}'`;
- chomp($comp);
+ # read cut file
+ open(IN,"${cutfile_dat}"); @lines = <IN>; close(IN); $nlines = @lines;
- # synthetics are always BH_ component
- $synfile = "${dirsyn_pro_1}/${sta}.${net}.BH${comp}.${syn_suffix}";
+ $sacfile = "sacdat.mac";
+ `echo echo on > $sacfile`;
+ `echo readerr badfile fatal >> $sacfile`;
- # if the synthetic file exists, then go on
- if (-f $synfile) {
- print "$datfile $synfile\n"; # only list files that are pairs
+ for ($j = 1; $j <= $nlines; $j++) {
- # get info on data and synthetics
- (undef,$bd,$ed,$deltad,$nptd) = split(" ",`saclst b e delta npts f $datfile`);
- (undef,$bs,$es,$deltas,$npts) = split(" ",`saclst b e delta npts f $synfile`);
- $tlend = $ed - $bd;
- $tlens = $es - $bs;
-
- # dt should be the same for both records ALREADY
- if (log($deltad/$deltas) > 0.01) {
- print "$datfile $synfile\n";
- print "DT values are not close enough: $deltad, $deltas\n";
- die("fix the DT values\n");
- } else {
- $dt = $deltad;
- }
+ $line = $lines[$j-1]; chomp($line);
+ ($datfile,$b,$e,$npt,$dt) = split(" ",$line);
+ print "$j out of $nlines -- $datfile\n";
+ #print "-- $datfile -- $synfile -- $b -- $e -- $npt -- $dt -- \n";
- # determine the cut for the records
- # b = earliest start time
- # e = earliest end time, multiplied by some factor
- if ($bd < $bs) {
- $b0 = $bd;
- } else {
- $b0 = $bs;
- }
- if ($b0 < $bmin) {
- $b0 = $bmin;
- }
- if ($ed < $es) {
- $e0 = $ed;
- } else {
- $e0 = $es;
- }
+ # cut records and fill zeros
+ `echo r $datfile >> $sacfile`;
+ `echo cuterr fillz >> $sacfile`;
+ `echo "cut $b n $npt" >> $sacfile`;
+ `echo r $datfile >> $sacfile`;
+ `echo cut off >> $sacfile`;
+ `echo w over >> $sacfile`;
+ }
+ `echo quit >> $sacfile`;
- $b = $b0;
- $tlen0 = $e0 - $b;
- $tlen = $tfac * $tlen0; # extend record length (if desired)
- $e = $b0 + $tlen;
- $npt = int( ($e-$b)/$dt );
+ # KEY: execute SAC command
+ `sac $sacfile`;
+ `sleep 5s`;
+ `rm $sacfile`;
+ print "\n Done cutting pre-processed data files\n";
- print CUT "$datfile $synfile $b $e $npt $dt\n";
+ # rename cut file in data directory
+ `mv ${cutfile_dat} ${cutfile_dat_done}`;
- if (0==1) {
- print "\n Data : $bd $ed $deltad $nptd -- $tlend";
- print "\n Syn : $bs $es $deltas $npts -- $tlens";
- print "\n b0, e0, tlen0 : $b0, $e0, $tlen0 ";
- print "\n b : $bd, $bs, $b ";
- print "\n e : $ed, $es, $e ";
- print "\n npt : $nptd, $npts, $npt ";
- print "\n $tlen = $tfac * ($e0 - $b)";
- print "\n $e = $b0 + $tlen \n";
- }
- } # if synfile exists
- } # for all data files
- print "\n Done making cutfile $cutfile\n";
- close(CUT);
+ } # cutfile exist
+ } # cutfile_done exist
+ } # bandpass dir exist
+ } # idata
- }
- } # if data (or syn) have already been bandpass filtered
- } # dirdat_pro_1 and dirsyn_pro_1 exist
- } # iprocess2=1
+ #------------------
+ if ($isyn == 1) {
+ if (-e ${dirsyn_pro_2}) {
+ print "--> dirsyn ${dirsyn_pro_2} already exists\n";
-# if ($iprocess2 == 1) {
-# # grab all DATA files
-# @files = glob("${dirdat_pro_1}/*");
-# $nfile = @files;
-# print "\n $nfile data files to line up with synthetics\n";
+ } else {
+ if (-f $cutfile_syn_done) {
+ print "cutfile_syn_done ${cutfile_syn_done} already exists\n";
-# #------------------------
+ } else {
-# $sacfile = "sac.mac";
-# `echo echo on > $sacfile`;
-# `echo readerr badfile fatal >> $sacfile`;
+ if (not -f $cutfile_syn) {
+ print "cutfile_syn $cutfile_syn does not exist -- try iprocess = 2\n";
-# foreach $file (@files) {
-# # read the sac headers -- network, station, component
-# (undef,$net,$sta,$chan) = split(" ",`saclst knetwk kstnm kcmpnm f $file`);
-# $comp = `echo $chan | awk '{print substr(\$1,3,1)}'`;
-# chomp($comp);
-# #print "\n $net $sta $chan $comp\n";
+ } else {
-# # synthetics are always BH_ component
-# $synt = "${dirsyn_pro_1}/${sta}.${net}.BH${comp}.${syn_suffix}";
-# #print "\n $file $synt ";
+ # read cut file
+ open(IN,"${cutfile_syn}"); @lines = <IN>; close(IN); $nlines = @lines;
-# # if the synthetic file exists, then go on
-# if (-f $synt) {
-# print "\n $file $synt"; # only list files that are pairs
+ $sacfile = "sacsyn.mac";
+ `echo echo on > $sacfile`;
+ `echo readerr badfile fatal >> $sacfile`;
-# # get info on data and synthetics
-# `echo r $file >> $sacfile`;
-# `echo r $synt >> $sacfile`;
-# (undef,$bd,$ed,$deltad,$nptd) = split(" ",`saclst b e delta npts f $file`);
-# (undef,$bs,$es,$deltas,$npts) = split(" ",`saclst b e delta npts f $synt`);
-# $tlend = $ed - $bd;
-# $tlens = $es - $bs;
-
-# # dt should be the same for both records ALREADY
-# if (log($deltad/$deltas) > 0.01) {
-# print "$file $synt\n";
-# print "DT values are not close enough: $deltad, $deltas\n";
-# die("fix the DT values\n");
-# } else {
-# $dt = $deltad;
-# }
+ for ($j = 1; $j <= $nlines; $j++) {
-# # determine the cut for the records
-# # b = earliest start time
-# # e = earliest end time, multiplied by some factor
-# if ($bd < $bs) {
-# $b0 = $bd;
-# } else {
-# $b0 = $bs;
-# }
-# if ($b0 < $bmin) {
-# $b0 = $bmin;
-# }
+ $line = $lines[$j-1]; chomp($line);
+ ($synfile,$b,$e,$npt,$dt) = split(" ",$line);
+ print "$j out of $nlines -- $synfile\n";
+ #print "-- $datfile -- $synfile -- $b -- $e -- $npt -- $dt -- \n";
-# if ($ed < $es) {
-# $e0 = $ed;
-# } else {
-# $e0 = $es;
-# }
+ # cut records and fill zeros
+ `echo r $synfile >> $sacfile`;
+ `echo cuterr fillz >> $sacfile`;
+ `echo "cut $b n $npt" >> $sacfile`;
+ `echo r $synfile >> $sacfile`;
+ `echo cut off >> $sacfile`;
+ `echo w over >> $sacfile`;
+ }
+ `echo quit >> $sacfile`;
-# $b = $b0;
-# $tlen0 = $e0 - $b;
-# $tlen = $tfac * $tlen0; # extend record length (if desired)
-# $e = $b0 + $tlen;
+ # KEY: execute SAC command
+ `sac $sacfile`;
+ `sleep 5s`;
+ `rm $sacfile`;
+ print "\n Done cutting pre-processed syn files\n";
-# $npt = int( ($e-$b)/$dt );
+ # copy syn cut file into syn directory
+ `cp ${cutfile_syn} ${cutfile_syn_done}`;
-# if (0==1) {
-# print "\n Data : $bd $ed $deltad $nptd -- $tlend";
-# print "\n Syn : $bs $es $deltas $npts -- $tlens";
-# print "\n b0, e0, tlen0 : $b0, $e0, $tlen0 ";
-# print "\n b : $bd, $bs, $b ";
-# print "\n e : $ed, $es, $e ";
-# print "\n npt : $nptd, $npts, $npt ";
-# print "\n $tlen = $tfac * ($e0 - $b)";
-# print "\n $e = $b0 + $tlen \n";
-# }
+ } # cutfile exist
+ } # cutfile_done exist
+ } # bandpass dir exist
+ } # isyn
+
+ } # iprocess = 3
-# # cut records and fill zeros
-# `echo r $file $synt >> $sacfile`;
-# `echo cuterr fillz >> $sacfile`;
-# `echo "cut $b n $npt" >> $sacfile`;
-# `echo r $file $synt >> $sacfile`; # cut both
-# `echo cut off >> $sacfile`;
-# `echo w over >> $sacfile`;
-# }
-# }
-# `echo quit >> $sacfile`;
-# `sac $sacfile`; # KEY: EXECUTE the SAC commands
-# `sleep 5s`;
-# `rm $sacfile`;
+ #----------------------------------------------------------------------
+ # PROCESSING PART 4: bandpass
-# print "\n Done with pre-processing, part 2\n";
+ if ($iprocess == 4) {
-# }
+ if ($isyn == 1) {
+ if (-e ${dirsyn_pro_2}) {
+ print "--> dirsyn ${dirsyn_pro_2} already exists\n";
-# #----------------------------------------------------------------------
-# # PROCESSING PART 3: band-pass filter and rotate
+ } else {
+ if (not -f ${cutfile_syn_done}) {
+ print "cutfile_syn_done ${cutfile_syn_done} does not exist\n";
-# if ($iprocess3 == 1) {
+ } else {
+ if (not -e ${dirsyn_pro_1}) {
+ print "${dirsyn_pro_1} does not exist\n";
-# if ($isyn == 1) {
-# if (-e ${dirsyn_pro_1}) {
-# if (not -e ${dirsyn_pro_2}) {
-# print CSH "cd ${dirsyn_pro_1}\n";
-# #print CSH "\\rm -rf $pdirbpass\n";
-# print CSH "process_trinet_syn_new.pl -S -t $Trange -d $pdirbpass -x $Ttag *.${syn_suffix} \n";
-# print CSH "cd $pdirbpass\n";
-# print CSH "rotate.pl *E.${syn_suffix}.${Ttag}\n";
+ } else {
+ print CSH "cd ${dirsyn_pro_1}\n";
+ #print CSH "\\rm -rf $pdirbpass\n";
+ print CSH "process_trinet_syn_new.pl -S -t $Trange -d $pdirbpass -x $Ttag *.${syn_suffix} \n";
+ print CSH "cd $pdirbpass\n";
+ print CSH "rotate.pl *E.${syn_suffix}.${Ttag}\n";
+ }
+ }
+ }
+ } # isyn
-# } else {
-# print "$pdirbpass already exists\n";
-# }
-# } else {
-# print "${dirsyn_pro_1} does not exist\n";
-# }
-# } # isyn
+ #-----------
+ if ($idata == 1) {
+ if (-e ${dirdat_pro_2}) {
+ print "--> dirdat ${dirdat_pro_2} already exists\n";
-# # data
-# if ($idata == 1) {
-# if (-e ${dirdat_pro_1}) {
-# if (not -e ${dirdat_pro_2}) {
+ } else {
+ if (not -f ${cutfile_dat_done}) {
+ print "cutfile_dat_done ${cutfile_dat_done} does not exist\n";
-# $ofile1 = "${eid}_ofile";
-# $ofile2 = "${eid}_no_pz_files";
-# $odir = "${dirdat0}/CHECK_DIR";
-# print CSH "mkdir -p $odir\n";
+ } else {
+ if (not -e ${dirdat_pro_1}) {
+ print "${dirdat_pro_1} does not exist\n";
-# print CSH "cd ${dirdat_pro_1}\n";
-# #print CSH "\\rm -rf $pdirbpass\n";
-# print CSH "process_cal_data.pl -i none -t $Trange -d $pdirbpass -x $Ttag *.${dat_suffix} > $ofile1\n";
-# print CSH "grep skip $ofile1 | awk '{print \$7}' > $ofile2\n";
-# print CSH "\\cp $ofile1 $ofile2 $odir\n";
+ } else {
+ $ofile1 = "${eid}_ofile";
+ $ofile2 = "${eid}_no_pz_files";
+ $odir = "${dirdat0}/CHECK_DIR";
+ print CSH "mkdir -p $odir\n";
-# print CSH "cd $pdirbpass\n";
-# print CSH "rotate.pl *E.${dat_suffix}.${Ttag}\n";
+ print CSH "cd ${dirdat_pro_1}\n";
+ #print CSH "\\rm -rf $pdirbpass\n";
+ print CSH "process_cal_data.pl -i none -t $Trange -d $pdirbpass -x $Ttag *.${dat_suffix} > $ofile1\n";
+ print CSH "grep skip $ofile1 | awk '{print \$7}' > $ofile2\n";
+ print CSH "\\cp $ofile1 $ofile2 $odir\n";
-# } else {
-# print "$pdirbpass already exists\n";
-# }
-# } else {
-# print "${dirdat_pro_1} does not exist\n";
-# }
-# } # idata
-# } # iprocess3=1
+ print CSH "cd $pdirbpass\n";
+ print CSH "rotate.pl *E.${dat_suffix}.${Ttag}\n";
+ }
+ }
+ }
+ } # idata
-}
+ } # iprocess==4
-if($iprocess0==1) {close(SYN);}
+} # END OF LOOP OVER EVENTS
+if($iprocess==0) {close(SYN);}
+
#======================
close(CSH);
print "closing $cshfile\n";
Added: seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T002_m05
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T002_m05 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T002_m05 2008-07-13 23:15:26 UTC (rev 12405)
@@ -0,0 +1,69 @@
+# -------------------------------------------------------------
+#
+# 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 = 2.00
+WIN_MAX_PERIOD = 40.00
+
+# -------------------------------------------------------------
+# E(t) water level
+STALTA_BASE = 0.07
+
+# -------------------------------------------------------------
+# TSHIFT
+TSHIFT_BASE = 3.0
+
+# -------------------------------------------------------------
+# limit on CC for window acceptance
+CC_BASE = 0.85
+
+# -------------------------------------------------------------
+# limit on dlnA (dA/A) for window acceptance
+DLNA_BASE = 1.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 = 2.5
+SNR_MAX_BASE = 3.5
+
+# -------------------------------------------------------------
+# limit on signal to noise ratio in a particular window.
+WINDOW_AMP_BASE = 4.0
+
+# -------------------------------------------------------------
+# Fine tuning constants
+C_0 (internal minima) = 1.0
+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 = 1.0
+WEIGHT_AVERAGE_CC = 0.0
+WEIGHT_N_WINDOWS = 5.0
Added: seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T006_m05
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T006_m05 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T006_m05 2008-07-13 23:15:26 UTC (rev 12405)
@@ -0,0 +1,69 @@
+# -------------------------------------------------------------
+#
+# 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 = .false.
+
+# -------------------------------------------------------------
+# period min/max for filtering
+RUN_BANDPASS = .false.
+WIN_MIN_PERIOD = 6.00
+WIN_MAX_PERIOD = 40.00
+
+# -------------------------------------------------------------
+# E(t) water level (0.23)
+STALTA_BASE = 0.22
+
+# -------------------------------------------------------------
+# TSHIFT
+TSHIFT_BASE = 4.0
+
+# -------------------------------------------------------------
+# limit on CC for window acceptance
+CC_BASE = 0.74
+
+# -------------------------------------------------------------
+# limit on dlnA (dA/A) for window acceptance
+DLNA_BASE = 1.5
+
+# -------------------------------------------------------------
+# 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 = 3.0
+SNR_MAX_BASE = 3.5
+
+# -------------------------------------------------------------
+# limit on signal to noise ratio in a particular window.
+WINDOW_AMP_BASE = 2.5
+
+# -------------------------------------------------------------
+# Fine tuning constants
+C_0 (internal minima) = 0.7
+C_1 (small windows) = 2.0
+C_2 (prominence) = 0.0
+C_3a (separation height) = 3.0
+C_3b (separation time) = 2.0
+C_4a (curtail on left) = 2.5
+C_4b (curtail on right) = 12.0
+
+WEIGHT_SPACE_COVERAGE = 1.0
+WEIGHT_AVERAGE_CC = 0.0
+WEIGHT_N_WINDOWS = 0.0
Added: seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m05.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m05.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m05.f90 2008-07-13 23:15:26 UTC (rev 12405)
@@ -0,0 +1,135 @@
+! -------------------------------------------------------------
+! 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
+
+ ! for qinya's scsn picking
+ double precision :: Pnl_start, S_end, Sw_start, Sw_end
+
+!===========================
+
+! -----------------------------------------------------------------
+! 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_AMP_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.
+ if (BODY_WAVE_ONLY) then
+ !Pnl_start = P_pick - 5.0
+ !S_end = S_pick + 5.0
+ Pnl_start = P_pick - 2.5*WIN_MIN_PERIOD
+ S_end = S_pick + 2.5*WIN_MIN_PERIOD
+ Sw_start = -15.0 + dist_km/3.5
+ Sw_end = 35.0 + dist_km/3.1
+
+ else
+ Pnl_start = -5.0 + dist_km/7.8
+ Sw_start = -15.0 + dist_km/3.5
+ Sw_end = 35.0 + dist_km/3.1
+ S_end = Sw_start
+ endif
+
+ ! variables for signal to noise ratio criteria.
+ signal_end = Sw_end
+ noise_end = Pnl_start
+ if(DEBUG) then
+ write(*,*) 'DEBUG : P_pick = ', sngl(P_pick)
+ write(*,*) 'DEBUG : signal_end = ', sngl(sigmal_end)
+ write(*,*) 'DEBUG : noise_end = ', sngl(noise_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.Pnl_start) then
+ STALTA_W_LEVEL(i) = 10.*STALTA_BASE
+ endif
+
+ ! raises STA/LTA water level after surface wave arrives
+ if (BODY_WAVE_ONLY) then
+ if(time.gt.S_end) then
+ STALTA_W_LEVEL(i) = 10.*STALTA_BASE
+ endif
+
+ else
+!!$ ! set time- and distance-specific Tshift and DlnA to mimic Qinya's criteria
+!!$ ! (see Liu et al., 2004, p. 1755; note comment above)
+!!$ if(time.ge.Pnl_start .and. time.lt.Sw_start) then
+!!$ !DLNA_LIMIT(i) = 1.5 ! ratio is 2.5, and dlna is ratio-1
+!!$ TSHIFT_LIMIT(i) = 3.0 + dist_km/80.0
+!!$ endif
+!!$ if(time.ge.Sw_start .and. time.le.Sw_end) then
+!!$ !DLNA_LIMIT(i) = 1.5 ! ratio is 2.5, and dlna is ratio-1
+!!$ TSHIFT_LIMIT(i) = 3.0 + dist_km/50.0
+!!$ endif
+
+ ! double the STA/LTA water level after the surface waves
+ if(time.gt.Sw_end) then
+ STALTA_W_LEVEL(i) = 2.0*STALTA_BASE
+ endif
+
+ endif
+
+ enddo
+
+! 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_AMP_BASE
+! endif
+! enddo
+!
+! End of user-dependent portion
+! -----------------------------------------------------------------
+
+end subroutine set_up_criteria_arrays
+! -------------------------------------------------------------
More information about the cig-commits
mailing list