[cig-commits] r12397 - in seismo/3D/ADJOINT_TOMO/flexwin: scripts scripts/prepare_scripts/socal user_files/socal_3D
carltape at geodynamics.org
carltape at geodynamics.org
Tue Jul 8 07:28:38 PDT 2008
Author: carltape
Date: 2008-07-08 07:28:37 -0700 (Tue, 08 Jul 2008)
New Revision: 12397
Added:
seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T002_m04
seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T006_m04
seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m04.f90
Modified:
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/prepare_scripts/socal/process_data_and_syn.pl
Log:
Updated run scripts to reflect renaming of test2.f90 to flexwin.f90. Added new user files for socal data.
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows.pl 2008-07-08 12:27:37 UTC (rev 12396)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows.pl 2008-07-08 14:28:37 UTC (rev 12397)
@@ -38,8 +38,8 @@
# prepare_meas_all.pl
#
# EXAMPLES:
-# pick_all_windows.pl m03 0 6/40 1/203 1/1/1/0 1 0 # both, T = 6s
-# pick_all_windows.pl m03 0 2/40 1/203 1/1/1/1 1 0 # both, T = 2s
+# 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 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
@@ -94,9 +94,9 @@
# specify various directories (MUST BE MODIFIED FOR EACH USER)
$dir0 = "/net/denali/scratch1/carltape/svn/cig/seismo/3D";
-$dir_win_code = "$dir0/ADJOINT_TOMO_WORK/flexwin";
+$dir_win_code = "$dir0/ADJOINT_TOMO/flexwin_work";
$dir_win_run = "$dir0/flexwin_run";
-#$dir_win_code = "$dir0/ADJOINT_TOMO_WORK/flexwin_copy2";
+#$dir_win_code = "$dir0/ADJOINT_TOMO/flexwin_work_copy2";
#$dir_win_run = "$dir0/flexwin_run_copy2";
$dir_scripts = "${dir_win_code}/scripts";
@@ -191,7 +191,7 @@
`cp ${userfun_file_in} ${userfun_file}`;
# name of executable windowing code (in $dir_win_code)
-$win_execute = "test2";
+$win_execute = "flexwin";
# make command for windowing code
$make = "make -f make_intel_caltech";
Modified: seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows_local.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows_local.pl 2008-07-08 12:27:37 UTC (rev 12396)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/pick_all_windows_local.pl 2008-07-08 14:28:37 UTC (rev 12397)
@@ -24,7 +24,7 @@
# prepare_meas_all.pl
#
# EXAMPLES:
-# /net/denali/scratch1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO_WORK/flexwin/scripts/pick_all_windows_local.pl
+# /net/denali/scratch1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO/flexwin_work/scripts/pick_all_windows_local.pl
#
# pick_all_windows_local.pl 2/40 1/1/1/1 input MEASURE
# pick_all_windows_local.pl 6/40 1/1/1/0 input MEASURE
@@ -39,7 +39,7 @@
($Ts,$ibools,$input,$dir) = @ARGV;
# specify source code directory (MUST BE MODIFIED FOR EACH USER)
-$dir_win_code = "/net/denali/scratch1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO_WORK/flexwin";
+$dir_win_code = "/net/denali/scratch1/carltape/svn/cig/seismo/3D/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")}
@@ -68,7 +68,7 @@
$par_file = "${dir_win_code}/PAR_FILE";
# name of executable windowing code (in $dir_win_code)
-$win_execute = "test2";
+$win_execute = "flexwin";
# make command for windowing code
$make = "make -f make_intel_caltech";
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-08 12:27:37 UTC (rev 12396)
+++ seismo/3D/ADJOINT_TOMO/flexwin/scripts/prepare_scripts/socal/process_data_and_syn.pl 2008-07-08 14:28:37 UTC (rev 12397)
@@ -3,7 +3,7 @@
#==========================================================
#
# Carl Tape
-# 19-June-2007
+# 03-July-2007
# process_data_and_syn.pl
#
# This script processes data and 3D-Socal-SEM synthetics for Southern California.
@@ -33,16 +33,17 @@
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;
-$smodel = "m03"; # KEY: model iteration index
+$smodel = "m04"; # KEY: model iteration index
$iexecute = 0;
$ilist = 1;
-# THESE MUST BE DONE IN SEQUENTIAL ORDER
$iprocess0 = 0;
-$iprocess1 = 1; # only done ONCE, no matter how many band-pass filters you use
-$iprocess2 = 0; # only done ONCE, no matter how many band-pass filters you use
-$iprocess3 = 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)
+
#$Lrange = "-l -40/220";
#$Lrange = " ";
@@ -101,7 +102,15 @@
#die("testing");
# write the C-shell script to file
-$cshfile = "process_data_and_syn.csh";
+if($idata==1 && $isyn==1) {
+ $cshfile = "process_data_and_syn.csh";
+} elsif($idata==1 && $isyn==0) {
+ $cshfile = "process_data.csh";
+} elsif($idata==0 && $isyn==1) {
+ $cshfile = "process_syn.csh";
+} else {
+ die("check idata and isyn\n");
+}
print "\nWriting to $cshfile ...\n";
open(CSH,">$cshfile");
@@ -111,8 +120,8 @@
}
$imin = 1; $imax = $ncmt; # default
-#$imin = 1; $imax = $imin;
-#$imin = 1; $imax = 160;
+#$imin = 1; $imax = 96;
+#$imin = 174; $imax = $imin;
#----------------------------------------------------------------------
@@ -132,34 +141,47 @@
print "$ievent, $imax, Event $eid\n";
print CSH "echo $ievent, $imax, Event $eid\n";
+ # data and synthetics directories
+ $dirsyn = "${dirsyn0}/${eid}";
+ $dirdat = "${dirdat0}/${eid}";
+ $dirdat_pro_1 = "${dirdat0}/${eid}/$pdir";
+ $dirsyn_pro_1 = "${dirsyn0}/${eid}/$pdir";
+ $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";
+
+ # optional -- delete pre-processed directories
+ #print CSH "rm -rf $dirsyn/PRO*\n";
+ #print CSH "rm -rf $dirdat/PRO*\n";
+
#----------------------------------------------------------------------
# PROCESSING PART 0: check the number of processed synthetic files for each event
if ($iprocess0 == 1) {
- $dirsyn = "${dirsyn0}/${eid}";
- $dirsyn_pro = "$dirsyn/$pdir";
- if(-e $dirsyn_pro) {
- ($nfile,undef,undef) = split(" ",`ls -1 $dirsyn/$pdir/* | wc`);
+ if(-e ${dirsyn_pro_1}) {
+ ($nfile,undef,undef) = split(" ",`ls -1 ${dirsyn_pro_1}/* | wc`);
print SYN "$eid $nfile\n";
}
}
#----------------------------------------------------------------------
- # PROCESSING PART 1: assigning SAC headers and interpolating
+ # PROCESSING PART 1: assign SAC headers, interpolate, and pick P and S arrivals (based on a 1D socal model)
if ($iprocess1 == 1) {
- # synthetics -- NOTE: convolve with half-duration (prior to interpolating)
+ # synthetics -- this will convolve with the source half-duration (prior to interpolating)
if ($isyn == 1) {
- $dirsyn = "${dirsyn0}/${eid}";
if (-e $dirsyn) {
- if (not -e "$dirsyn/$pdir") {
+ if (not -e ${dirsyn_pro_1}) {
print CSH "cd $dirsyn\n";
#print CSH "mv $pdir ${pdir}_OLD\n";
#print CSH "\\rm -rf $pdir\n";
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 "$pdir already exists\n";
+ print "dir ${dirsyn_pro_1} already exists\n";
}
} else {
print "$dirsyn does not exist\n";
@@ -168,15 +190,14 @@
# data
if ($idata == 1) {
- $dirdat = "${dirdat0}/${eid}";
if (-e $dirdat) {
- if (not -e "$dirdat/$pdir") {
+ if (not -e ${dirdat_pro_1}) {
print CSH "cd $dirdat\n";
#print CSH "mv $pdir ${pdir}_OLD\n";
#print CSH "\\rm -rf $pdir; mkdir $pdir\n";
print CSH "process_cal_data.pl -m $cmtfile -p -s $sps -d $pdir -x ${dat_ext} *.${dat_suffix0}\n";
} else {
- print "$pdir already exists\n";
+ print "dir ${dirdat_pro_1} already exists\n";
}
} else {
print "$dirdat does not exist\n";
@@ -237,73 +258,149 @@
}
#----------------------------------------------------------------------
- # PROCESSING PART 2: cutting records and padding zeros
+ # PROCESSING PART 2: cutting records, padding zeros, and bandpass
- if ($iprocess2 == 1 || $iprocess3 == 1) {
+ if ($iprocess2 == 1) {
- # both the processed synthetics and processed data directories must exist
- # base directories containing pre-processed data and synthetics
- $dirdat = "${dirdat0}/${eid}/$pdir";
- $dirsyn = "${dirsyn0}/${eid}/$pdir";
- $dirdat_pro = "$dirdat/$pdirbpass";
- $dirsyn_pro = "$dirsyn/$pdirbpass";
+ # BOTH the (initially) processed synthetics and data directories must exist
- if ( (not -e $dirdat) || (not -e $dirsyn) ) {
- print "--> dirdat $dirdat and dirsyn $dirsyn do not both exist\n";
+ 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 {
- if ( (-e $dirdat_pro && $idata==1) || (-e $dirsyn_pro && $isyn==1) ) {
- if ($idata == 1) {
- print "$dirdat_pro already exists\n";
+
+ 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 already exists\n";
+ print "${dirsyn_pro_2} already exists\n";
}
} else {
- if ($iprocess2 == 1) {
-
- # base directories containing pre-processed data and synthetics
- $dirdat = "${dirdat0}/${eid}/$pdir";
- $dirsyn = "${dirsyn0}/${eid}/$pdir";
+ if (-f $cutfile_done) {
+ print "cutfile ${cutfile_done} already exists -- now bandpass\n";
+ if ($iprocess3==0) {
+ print "set iprocess3=1 if you want to bandpass\n";
+
+ } else {
+ print "iprocess3=1 -- now bandpass\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 "$pdirbpass already exists\n";
+ }
+ } else {
+ print "${dirsyn_pro_1} does not exist\n";
+ }
+ } # isyn
+
+ # data
+ if ($idata == 1) {
+ if (-e ${dirdat_pro_1}) {
+ if (not -e ${dirdat_pro_2}) {
+
+ $ofile1 = "${eid}_ofile";
+ $ofile2 = "${eid}_no_pz_files";
+ $odir = "${dirdat0}/CHECK_DIR";
+ print CSH "mkdir -p $odir\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";
+
+ print CSH "cd $pdirbpass\n";
+ print CSH "rotate.pl *E.${dat_suffix}.${Ttag}\n";
+
+ } else {
+ print "$pdirbpass already exists\n";
+ }
+ } else {
+ print "${dirdat_pro_1} does not exist\n";
+ }
+ } # idata
+
+ } # iprocess3==1
+
+ } elsif (-f $cutfile) {
+ print "cutfile $cutfile exists -- now read it in sac and execute\n";
+
+ # read cut file
+ open(IN,"$cutfile"); @lines = <IN>; close(IN); $nlines = @lines;
+
+ $sacfile = "sac.mac";
+ `echo echo on > $sacfile`;
+ `echo readerr badfile fatal >> $sacfile`;
+
+ 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");
+
# grab all DATA files
- @files = glob("${dirdat}/*");
+ @files = glob("${dirdat_pro_1}/*");
$nfile = @files;
print "\n $nfile data files to line up with synthetics\n";
- #------------------------
-
- `echo echo on > sac.mac`;
- `echo readerr badfile fatal >> sac.mac`;
-
- foreach $file (@files) {
- # read the sac headers
- (undef,$net,$sta,$chan) = split(" ",`saclst knetwk kstnm kcmpnm f $file`);
+ 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);
- #print "\n $net $sta $chan $comp\n";
# synthetics are always BH_ component
- $synt = "${dirsyn}/${sta}.${net}.BH${comp}.${syn_suffix}";
- #print "\n $file $synt ";
+ $synfile = "${dirsyn_pro_1}/${sta}.${net}.BH${comp}.${syn_suffix}";
# if the synthetic file exists, then go on
- if (-f $synt) {
- print "\n $file $synt"; # only list files that are pairs
+ if (-f $synfile) {
+ print "$datfile $synfile\n"; # only list files that are pairs
# get info on data and synthetics
- `echo r $file >> sac.mac`;
- `echo r $synt >> sac.mac`;
- (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`);
+ (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 "$file $synt\n";
+ print "$datfile $synfile\n";
print "DT values are not close enough: $deltad, $deltas\n";
die("fix the DT values\n");
} else {
@@ -321,7 +418,6 @@
if ($b0 < $bmin) {
$b0 = $bmin;
}
-
if ($ed < $es) {
$e0 = $ed;
} else {
@@ -330,11 +426,12 @@
$b = $b0;
$tlen0 = $e0 - $b;
- $tlen = $tfac * $tlen0; # extend record length
+ $tlen = $tfac * $tlen0; # extend record length (if desired)
$e = $b0 + $tlen;
-
$npt = int( ($e-$b)/$dt );
+ print CUT "$datfile $synfile $b $e $npt $dt\n";
+
if (0==1) {
print "\n Data : $bd $ed $deltad $nptd -- $tlend";
print "\n Syn : $bs $es $deltas $npts -- $tlens";
@@ -345,87 +442,174 @@
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);
- # cut records and fill zeros
- `echo r $file $synt >> sac.mac`;
- `echo cuterr fillz >> sac.mac`;
- `echo "cut $b n $npt" >> sac.mac`;
- `echo r $file $synt >> sac.mac`; # cut both
- `echo cut off >> sac.mac`;
- `echo w over >> sac.mac`;
- }
+ }
+ } # if data (or syn) have already been bandpass filtered
+ } # dirdat_pro_1 and dirsyn_pro_1 exist
+ } # iprocess2=1
- }
- `echo quit >> sac.mac`;
- `sac sac.mac`; # KEY: EXECUTE the SAC commands
- `rm sac.mac`;
- print "\n Done with pre-processing, part 2\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";
- }
+# #------------------------
- #----------------------------------------------------------------------
- # PROCESSING PART 3: band-pass filter and rotate
+# $sacfile = "sac.mac";
+# `echo echo on > $sacfile`;
+# `echo readerr badfile fatal >> $sacfile`;
- if ($iprocess3 == 1) {
+# 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";
- if ($isyn == 1) {
- $dirsyn = "${dirsyn0}/${eid}/$pdir";
- if (-e $dirsyn) {
- if (not -e "$dirsyn/$pdirbpass") {
- print CSH "cd $dirsyn\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
+# $synt = "${dirsyn_pro_1}/${sta}.${net}.BH${comp}.${syn_suffix}";
+# #print "\n $file $synt ";
- } else {
- print "$pdirbpass already exists\n";
- }
- } else {
- print "$dirsyn does not exist\n";
- }
- } # isyn
+# # if the synthetic file exists, then go on
+# if (-f $synt) {
+# print "\n $file $synt"; # only list files that are pairs
+# # 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;
+# }
- # data
- if ($idata == 1) {
- $dirdat = "${dirdat0}/${eid}/$pdir";
- if (-e $dirdat) {
- if (not -e "$dirdat/$pdirbpass") {
+# # 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;
+# }
- $ofile1 = "${eid}_ofile";
- $ofile2 = "${eid}_no_pz_files";
- $odir = "${dirdat0}/CHECK_DIR";
- print CSH "mkdir -p $odir\n";
+# if ($ed < $es) {
+# $e0 = $ed;
+# } else {
+# $e0 = $es;
+# }
- print CSH "cd $dirdat\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;
- print CSH "cd $pdirbpass\n";
- print CSH "rotate.pl *E.${dat_suffix}.${Ttag}\n";
+# $npt = int( ($e-$b)/$dt );
- } else {
- print "$pdirbpass already exists\n";
- }
- } else {
- print "$dirdat does not exist\n";
- }
- } # idata
- } # iprocess3=1
- } # if data (or syn) have already been bandpass filtered
- } # dirdat and dirsyn exist
- } # iprocess2=1 or iprocess3=1
+# 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";
+# }
+# # 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`;
+
+# print "\n Done with pre-processing, part 2\n";
+
+# }
+
+# #----------------------------------------------------------------------
+# # PROCESSING PART 3: band-pass filter and rotate
+
+# if ($iprocess3 == 1) {
+
+# 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 "$pdirbpass already exists\n";
+# }
+# } else {
+# print "${dirsyn_pro_1} does not exist\n";
+# }
+# } # isyn
+
+
+# # data
+# if ($idata == 1) {
+# if (-e ${dirdat_pro_1}) {
+# if (not -e ${dirdat_pro_2}) {
+
+# $ofile1 = "${eid}_ofile";
+# $ofile2 = "${eid}_no_pz_files";
+# $odir = "${dirdat0}/CHECK_DIR";
+# print CSH "mkdir -p $odir\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";
+
+# print CSH "cd $pdirbpass\n";
+# print CSH "rotate.pl *E.${dat_suffix}.${Ttag}\n";
+
+# } else {
+# print "$pdirbpass already exists\n";
+# }
+# } else {
+# print "${dirdat_pro_1} does not exist\n";
+# }
+# } # idata
+# } # iprocess3=1
+
}
if($iprocess0==1) {close(SYN);}
#======================
close(CSH);
+print "closing $cshfile\n";
if($iexecute==1) {system("csh -f $cshfile");}
print "\n ";
Added: seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T002_m04
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T002_m04 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T002_m04 2008-07-08 14:28:37 UTC (rev 12397)
@@ -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_m04
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T006_m04 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/PAR_FILE_T006_m04 2008-07-08 14:28:37 UTC (rev 12397)
@@ -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_m04.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m04.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/flexwin/user_files/socal_3D/user_functions_m04.f90 2008-07-08 14:28:37 UTC (rev 12397)
@@ -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