[cig-commits] r12413 - in seismo/3D/ADJOINT_TOMO/measure_adj: . scripts_tomo
carltape at geodynamics.org
carltape at geodynamics.org
Mon Jul 14 11:40:09 PDT 2008
Author: carltape
Date: 2008-07-14 11:40:09 -0700 (Mon, 14 Jul 2008)
New Revision: 12413
Modified:
seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/combine_adj_src_all.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/prepare_mt_measure_adj.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl
Log:
Updated run scripts. Updated mt_measure_adj.f90 to effectively ignore input windows whose cross-correlation traveltime measurement exceeds (in absolute value) the BEFORE_TSHIFT parameter in MEASUREMENT.PAR. Before, this parameter was only used to reject a multi-taper measurement possibility.
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2008-07-14 18:10:40 UTC (rev 12412)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2008-07-14 18:40:09 UTC (rev 12413)
@@ -257,14 +257,15 @@
print *, 'fstart0/fend0 :', fstart0, fend0
print *, ' fstart/fend :', fstart, fend
- if (use_trace .and. (cc_max < BEFORE_QUALITY .or. cc_max_f1f2 < AFTER_QUALITY .or. tshift_xc > BEFORE_TSHIFT .or. tshift_f1f2 > AFTER_TSHIFT)) then
+ ! CHT: tshift_xc --> abs(tshift_xc) and tshift_f1f2 --> abs(tshift_f1f2)
+ if (use_trace .and. (cc_max < BEFORE_QUALITY .or. cc_max_f1f2 < AFTER_QUALITY .or. abs(tshift_xc) > BEFORE_TSHIFT .or. abs(tshift_f1f2) > AFTER_TSHIFT)) then
use_trace = .false.
if (DISPLAY_DETAILS) then
print *, 'Fail if ANY of these is true :'
print *, ' cc_max : ', cc_max, BEFORE_QUALITY, cc_max < BEFORE_QUALITY
print *, ' cc_max_f1f2 : ', cc_max_f1f2, AFTER_QUALITY, cc_max_f1f2 < AFTER_QUALITY
- print *, ' tshift_xc : ', tshift_xc, BEFORE_TSHIFT, tshift_xc > BEFORE_TSHIFT
- print *, ' tshift_f1f2 : ', tshift_f1f2, AFTER_TSHIFT, tshift_f1f2 > AFTER_TSHIFT
+ print *, ' tshift_xc : ', tshift_xc, BEFORE_TSHIFT, abs(tshift_xc) > BEFORE_TSHIFT
+ print *, ' tshift_f1f2 : ', tshift_f1f2, AFTER_TSHIFT, abs(tshift_f1f2) > AFTER_TSHIFT
endif
endif
print *, ' Status: ', trim(measure_file_prefix), use_trace
@@ -286,6 +287,17 @@
use_trace = .true.
endif
+ ! CHT: if the CC traveltime shift is greater than the specified value, then we effectively eliminate the window
+ ! In principle, the windowing code should take care of this threshold.
+ if ( abs(tshift_xc) > BEFORE_TSHIFT ) then
+ tshift_xc = 0.0
+ dlnA = 0.0
+ if (DISPLAY_DETAILS) then
+ print *, ' effectively eliminating this time window -->'
+ print *, ' tshift_xc : ', tshift_xc, BEFORE_TSHIFT, abs(tshift_xc) > BEFORE_TSHIFT
+ endif
+ endif
+
! write frequency limits to file
if (OUTPUT_MEASUREMENT_FILES) then
df = 1./(dt*NPT)
@@ -297,7 +309,7 @@
! compute adjoint sources and misfit function values and also the CC-reconstructed records
if (use_trace) then
- tr_chi = 0. ; am_chi = 0. ! must be initialized
+ tr_chi = 0.0 ; am_chi = 0.0 ! must be initialized
if (iker==0) then
call mt_adj(iker,istart,dzr20,dzr0,nlen,dt,tshift_xc,dlnA,sigma_dt_cc,sigma_dlnA_cc,&
dtau_w,dlnA_w,err_dt,err_dlnA,sigma_dt,sigma_dlnA,i_left,i_right,tr_adj_src,tr_chi,window_chi)
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/combine_adj_src_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/combine_adj_src_all.pl 2008-07-14 18:10:40 UTC (rev 12412)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/combine_adj_src_all.pl 2008-07-14 18:40:09 UTC (rev 12413)
@@ -1,29 +1,29 @@
#!/usr/bin/perl -w
#-----------------------------------
-# Carl Tape, 15-Nov-2007
+# Carl Tape, 15-June-2008
# combine_adj_src_all.pl
#
# This script calls combine_adj_src.pl to plot make composite adjoint sources.
# for a set of events.
#
# EXAMPLE:
-# combine_adj_src_all.pl 2/6 m03
+# combine_adj_src_all.pl 2 6 cc cc m05
#
#-----------------------------------
-if (@ARGV < 2) {die("Usage: combine_adj_src_all.pl Tmin/Tmax model iplot\n")}
-($Ts,$smodel) = @ARGV;
+if (@ARGV < 5) {die("Usage: combine_adj_src_all.pl Tmin Tmax tag1 tag2 model iplot\n")}
+($Tmin1,$Tmin2,$tag1,$tag2,$smodel) = @ARGV;
$pwd = $ENV{PWD};
-($Tmin1,$Tmin2) = split("/",$Ts);
+#($Tmin1,$Tmin2) = split("/",$Ts);
$sTmin1 = sprintf("T%3.3i",$Tmin1);
$sTmin2 = sprintf("T%3.3i",$Tmin2);
# list of event IDs
#$file_eids = "/net/sierra/raid1/carltape/results/SOURCES/eid_run"; # favorites
-$file_eids = "/net/sierra/raid1/carltape/results/SOURCES/socal_6/SOCAL_FINAL_CMT_v6_eid";
+$file_eids = "/net/sierra/raid1/carltape/results/SOURCES/socal_7/SOCAL_FINAL_CMT_v7_eid";
if (not -f $file_eids) {die("\n check if $file_eids exists\n")}
open(IN,$file_eids); @eids = <IN>; $nevent = @eids;
@@ -43,7 +43,7 @@
# loop over all events
$k = 0;
$imin = 1; $imax = $nevent; # default
-$imin = 1; $imax = 100;
+#$imin = 1; $imax = 100;
#$imin = 189; $imax = $imin; # 133 (9818433)
for ($i = $imin; $i <= $imax; $i = $i+1) {
@@ -83,7 +83,7 @@
if ( -f $sta_file1 && -f $sta_file2) {
# Case 1: both STATIONS_ADJOINT files exist: execute combine_adj_src.pl
print CSH "echo running combine_adj_src.pl...\n";
- print CSH "combine_adj_src.pl $dir_adj1 $dir_adj2 $odir cc mtm\n";
+ print CSH "combine_adj_src.pl $dir_adj1 $dir_adj2 $odir $tag1 $tag2\n";
} else {
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/prepare_mt_measure_adj.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/prepare_mt_measure_adj.pl 2008-07-14 18:10:40 UTC (rev 12412)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/prepare_mt_measure_adj.pl 2008-07-14 18:40:09 UTC (rev 12413)
@@ -99,7 +99,7 @@
# CMTSOLUTION file
#$dir_cmt = "/net/sierra/raid1/carltape/results/SOURCES/socal_6/CMT_files_post_inverted";
-$dir_cmt = "/net/sierra/raid1/carltape/results/SOURCES/socal_6/CMT_files_post_inverted_mod";
+$dir_cmt = "/net/sierra/raid1/carltape/results/SOURCES/socal_7/CMT_files_pre_inverted";
$cmtfile = "${dir_cmt}/CMTSOLUTION_${eid}";
if (not -f $cmtfile) {die("check if cmtfile $cmtfile exist or not\n")}
print CSH "\\rm CMTSOLUTION*\n";
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2008-07-14 18:10:40 UTC (rev 12412)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2008-07-14 18:40:09 UTC (rev 12413)
@@ -4,14 +4,19 @@
#
# run_tomo.pl
# Carl Tape
-# 12-May-2008
+# 12-June-2008
#
# This is the master program for running the measurement code on all the
# events for which the windowing code has been run.
#
# EXAMPLE (from working measure_adj directory) :
-# run_tomo.pl m03 0 6/40 0.011 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.2 0.8
-# run_tomo.pl m03 0 2/40 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.2 0.5
+# socal m05 -->
+# run_tomo.pl m05 0 6/40 0.011 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/5.0/1.0 0.2/0.2 0.8
+# run_tomo.pl m05 0 2/40 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/4.0/1.0 0.2/0.2 0.5
+#
+# socal m00 - m04
+# run_tomo.pl m05 0 6/40 0.011 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.2 0.8
+# run_tomo.pl m05 0 2/40 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.2 0.5
#
#==========================================================
@@ -30,7 +35,7 @@
if (not -e $dir_all) {die("check if $dir_all exist or not\n")}
# file containing event IDs, half-durations for source, and number of time-steps for simulations
-$file_eid = "/net/sierra/raid1/carltape/results/SOURCES/socal_6/EIDs_simulation_duration";
+$file_eid = "/net/sierra/raid1/carltape/results/SOURCES/socal_7/EIDs_simulation_duration";
if (not -f ${file_eid}) {die("check if ${file_eid} exist or not\n")}
open(IN,"${file_eid}"); @elines = <IN>; close(IN);
$nevent = @elines;
@@ -48,8 +53,8 @@
#die("testing");
$imin = 1; $imax = $nevent;
-#$imin = 1; $imax = 20;
-#$imin = 114; $imax = $imin;
+#$imin = 1; $imax = 160;
+#$imin = 168; $imax = $imin;
# EXAMPLE lines from EIDs_simulation_duration
# 5 10222753 300.0 27300 0.26 332.0 4.010
@@ -104,14 +109,17 @@
print "--> lets make adjoint sources\n";
print CSH "prepare_mt_measure_adj.pl $smodel 0 1 0 $Ts $eid\n";
- # KEY: for T >=6, do both surface wave and body wave measurements
- if ($Tmin <= 2 || $eid == 14095540) {
- # special case for a particular event -- do cross-correlation only
- print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
- } else {
- print CSH "run_mt_cc_plot.pl $smodel $ibp $Ts $lcut $tvec 0 $iparbools $par1 $par2 $par3\n";
- }
+# # KEY: for T >=6, do both surface wave and body wave measurements
+# if ($Tmin <= 2 || $eid == 14095540) {
+# # special case for a particular event -- do cross-correlation only
+# print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+# } else {
+# print CSH "run_mt_cc_plot.pl $smodel $ibp $Ts $lcut $tvec 0 $iparbools $par1 $par2 $par3\n";
+# }
+ # SOCAL: as of model m04, we are no longer making MT measurements
+ print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+
$infile1 = "window_chi"; $outfile1 = "$dir_meas/${ftag}_${infile1}";
$infile2 = "run_file"; $outfile2 = "$dir_meas/${ftag}_${infile2}";
#$tfile = `ls PLOTS/*mt_cc*`; $infile3 = `basename $tfile`;
More information about the cig-commits
mailing list