[cig-commits] r12507 - in seismo/3D/ADJOINT_TOMO/measure_adj: . scripts_tomo
carltape at geodynamics.org
carltape at geodynamics.org
Fri Aug 1 13:12:01 PDT 2008
Author: carltape
Date: 2008-08-01 13:12:01 -0700 (Fri, 01 Aug 2008)
New Revision: 12507
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:
Update in main code to effectively reject windows within which the synthetics are much larger in amplitude than the data. The cutoff of dlnA < -0.80 was emperically based on the socal records, but will probably be about right for any records. Eventually, this should probably be included as an input parameter in the MEASUREMENT.PAR file. For now, it saves the user from manually rejecting the window, since it zeros out the measurement (and, thus, adjoint source).
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2008-08-01 20:04:13 UTC (rev 12506)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2008-08-01 20:12:01 UTC (rev 12507)
@@ -24,6 +24,8 @@
double precision :: sec,dist,az,baz,slat,slon
character(len=10) :: net,sta,chan_dat,chan,cmp
+ double precision :: BEFORE_DLNA ! temporary: this should be an input parameter
+
double precision :: tshift_xc, sigma_dt_cc, dlnA, sigma_dlnA_cc, sigma_dt, sigma_dlnA
double precision :: tr_chi, am_chi, tshift_f1f2, cc_max_f1f2, cc_max
double precision, dimension(NPT) :: dtau_w, dlnA_w, err_dt, err_dlnA, dzr0, dzr20
@@ -289,13 +291,16 @@
! 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
+ ! NOTE: BEFORE_DLNA should be included in the PARFILE
+ BEFORE_DLNA = -0.80
+ if ( abs(tshift_xc) > BEFORE_TSHIFT .or. dlnA < BEFORE_DLNA ) then
if (DISPLAY_DETAILS) then
print *, ' effectively eliminating this time window -->'
print *, ' tshift_xc : ', tshift_xc, BEFORE_TSHIFT, abs(tshift_xc) > BEFORE_TSHIFT
+ print *, ' dlnA : ', dlnA, BEFORE_DLNA, dlnA < BEFORE_DLNA
endif
+ tshift_xc = 0.0
+ dlnA = 0.0
endif
! write frequency limits to file
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-08-01 20:04:13 UTC (rev 12506)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/combine_adj_src_all.pl 2008-08-01 20:12:01 UTC (rev 12507)
@@ -1,14 +1,14 @@
#!/usr/bin/perl -w
#-----------------------------------
-# Carl Tape, 15-June-2008
+# Carl Tape, 25-July-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 cc cc m05
+# combine_adj_src_all.pl 2 6 cc cc m06
#
#-----------------------------------
@@ -23,7 +23,7 @@
# list of event IDs
#$file_eids = "/net/sierra/raid1/carltape/results/SOURCES/eid_run"; # favorites
-$file_eids = "/net/sierra/raid1/carltape/results/SOURCES/socal_7/SOCAL_FINAL_CMT_v7_eid";
+$file_eids = "/net/sierra/raid1/carltape/results/SOURCES/socal_9/EIDs_only_eid";
if (not -f $file_eids) {die("\n check if $file_eids exists\n")}
open(IN,$file_eids); @eids = <IN>; $nevent = @eids;
@@ -31,6 +31,14 @@
$dir_run = "/net/sierra/raid1/carltape/socal/socal_3D/RUNS";
if (not -e $dir_run) {die("check if $dir_run exist or not\n")}
+if (0==1) {
+ for ($i = 1; $i <= $nevent; $i = $i+1) {
+ $eid = $eids[$i-1]; chomp($eid);
+ print "$i -- $eid\n";
+ }
+ die("testing");
+}
+
#----------------------------------------------------
#$nevent = 5;
@@ -44,7 +52,8 @@
$k = 0;
$imin = 1; $imax = $nevent; # default
#$imin = 1; $imax = 100;
-#$imin = 189; $imax = $imin; # 133 (9818433)
+#$imin = 156; $imax = $imin;
+
for ($i = $imin; $i <= $imax; $i = $i+1) {
$eid = $eids[$i-1]; chomp($eid);
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-08-01 20:04:13 UTC (rev 12506)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/prepare_mt_measure_adj.pl 2008-08-01 20:12:01 UTC (rev 12507)
@@ -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_7/CMT_files_pre_inverted";
+$dir_cmt = "/net/sierra/raid1/carltape/results/SOURCES/socal_9/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-08-01 20:04:13 UTC (rev 12506)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2008-08-01 20:12:01 UTC (rev 12507)
@@ -10,13 +10,13 @@
# events for which the windowing code has been run.
#
# EXAMPLE (from working measure_adj directory) :
-# 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 m06 -->
+# run_tomo.pl m07 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.4/0.4 1.0
+# run_tomo.pl m07 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
+# run_tomo.pl m04 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 m04 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
#
#==========================================================
@@ -26,7 +26,7 @@
($Tmin,$Tmax) = split("/",$Ts);
$sTmin = sprintf("T%3.3i",$Tmin);
-# for plotting purposes only
+# cut times for plotting purposes only
($lcut_min,$lcut_max) = split("/",$lcut);
$dir0 = "/net/sierra/raid1/carltape/socal/socal_3D/RUNS";
@@ -35,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_7/EIDs_simulation_duration";
+$file_eid = "/net/sierra/raid1/carltape/results/SOURCES/socal_9/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;
@@ -54,7 +54,7 @@
$imin = 1; $imax = $nevent;
#$imin = 1; $imax = 160;
-#$imin = 168; $imax = $imin;
+$imin = 204; $imax = $imin;
# EXAMPLE lines from EIDs_simulation_duration
# 5 10222753 300.0 27300 0.26 332.0 4.010
@@ -106,7 +106,7 @@
print "--> DO NOT RUN: no window file exists ($window_file)\n";
} else {
- print "--> lets make adjoint sources\n";
+ print "--> run the measurement code and 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
More information about the cig-commits
mailing list