[cig-commits] r15852 - in seismo/3D/ADJOINT_TOMO/measure_adj: . DATA PLOTS SYN scripts_tomo
carltape at geodynamics.org
carltape at geodynamics.org
Tue Oct 20 14:54:54 PDT 2009
Author: carltape
Date: 2009-10-20 14:54:52 -0700 (Tue, 20 Oct 2009)
New Revision: 15852
Added:
seismo/3D/ADJOINT_TOMO/measure_adj/DATA/9818433.CI.MPM.BHR.sac.d.T006_T030
seismo/3D/ADJOINT_TOMO/measure_adj/DATA/9818433.CI.MPM.BHT.sac.d.T006_T030
seismo/3D/ADJOINT_TOMO/measure_adj/DATA/9818433.CI.MPM.BHZ.sac.d.T006_T030
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/9818433_T006_T030_MPM_CI_m16_ikerXX_win_adj.pdf
seismo/3D/ADJOINT_TOMO/measure_adj/SYN/MPM.CI.BHR.semd.sac.m16.T006_T030
seismo/3D/ADJOINT_TOMO/measure_adj/SYN/MPM.CI.BHT.semd.sac.m16.T006_T030
seismo/3D/ADJOINT_TOMO/measure_adj/SYN/MPM.CI.BHZ.semd.sac.m16.T006_T030
Removed:
seismo/3D/ADJOINT_TOMO/measure_adj/DATA_TEST/
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/9818433_T006_T030_MPM_CI_m16_mtm_win_adj.pdf
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/EXAMPLE/
seismo/3D/ADJOINT_TOMO/measure_adj/SYN_TEST/
Modified:
seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR
seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.WINDOWS
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/make_pdf_by_event.pl
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl
seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
seismo/3D/ADJOINT_TOMO/measure_adj/mt_variables.f90
seismo/3D/ADJOINT_TOMO/measure_adj/prepare_adj_src.pl
seismo/3D/ADJOINT_TOMO/measure_adj/run_mt_measure_adj.csh
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/prepare_mt_measure_adj.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl
seismo/3D/ADJOINT_TOMO/measure_adj/write_par_file.pl
Log:
Major organizational changes, including modifications to the MEASUREMENT.PAR format. The program allows for iker = 0-8, where 0 returns no adjoint source and 1-8 returns an adjoint source: 1=data only, 2=waveform difference, 3=CC-TT banana-doughnut, 4=dlnA banana-doughnut, 5=CC-TT, 6=dlnA, 7=multitaper-TT, 8=multitaper amplitude. A PDF illustrating all nine examples is included.
Copied: seismo/3D/ADJOINT_TOMO/measure_adj/DATA/9818433.CI.MPM.BHR.sac.d.T006_T030 (from rev 15839, seismo/3D/ADJOINT_TOMO/measure_adj/DATA_TEST/9818433.CI.MPM.BHR.sac.d.T006_T030)
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/measure_adj/DATA/9818433.CI.MPM.BHR.sac.d.T006_T030
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/measure_adj/DATA/9818433.CI.MPM.BHT.sac.d.T006_T030 (from rev 15839, seismo/3D/ADJOINT_TOMO/measure_adj/DATA_TEST/9818433.CI.MPM.BHT.sac.d.T006_T030)
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/measure_adj/DATA/9818433.CI.MPM.BHT.sac.d.T006_T030
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/measure_adj/DATA/9818433.CI.MPM.BHZ.sac.d.T006_T030 (from rev 15839, seismo/3D/ADJOINT_TOMO/measure_adj/DATA_TEST/9818433.CI.MPM.BHZ.sac.d.T006_T030)
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/measure_adj/DATA/9818433.CI.MPM.BHZ.sac.d.T006_T030
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Name: svn:mergeinfo
+
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR 2009-10-20 21:54:52 UTC (rev 15852)
@@ -1,17 +1,19 @@
-OUTPUT_FILES
- 1 # is_mtm -- taper type: 0 none; 1 multi-taper; 2 cosine; 3 boxcar
- 0.020 2.50 # wtr, npi (ntaper = 2*npi)
- 1 # iker -- 0 waveform; 1 multi-taper; 2 cross-correlation;
+ -0.585 0.0110 18200 # tstart, DT, npts: time vector for simulations
+ 7 # iker (0-8) -- 0 no adj src; 1-8 see manual;
+ BH # channel: BH or LH
+ 30.000 6.000 # TLONG and TSHORT: band-pass periods for records
.false. # RUN_BANDPASS: use band-pass on records
- 30.000 6.000 # TLONG and TSHORT: band-pass periods for records
- -0.585 0.0110 18200 # tstart, DT, npts: time vector for simulations
.true. # DISPLAY_DETAILS
.true. # OUTPUT_MEASUREMENT_FILES
- 1 # INCLUDE_ERROR -- 0 none; 1 CC, MT; 2 CC, MT-jack-knife
+ -4.5000 4.5000 # TSHIFT_MIN; TSHIFT_MAX
+ -1.5000 1.5000 # DLNA_MIN; DLNA_MAX
+ 0.690 # CC_MIN
+ 0 # ERROR_TYPE -- 0 none; 1 CC, MT-CC; 2 MT-jack-knife
+ 1.000 # DT_SIGMA_MIN
+ 0.500 # DLNA_SIGMA_MIN
+ 1 # itaper -- taper type: 1 multi-taper; 2 cosine; 3 boxcar
+ 0.020 2.50 # wtr, npi (ntaper = 2*npi)
2.000 # DT_FAC
2.500 # ERR_FAC
3.500 # DT_MAX_SCALE
1.500 # NCYCLE_IN_WINDOW
- 0.6900 0.6900 # BEFORE_QUALITY; AFTER_QUALITY
- 4.5000 1.0000 # BEFORE_TSHIFT; AFTER_TSHIFT
- 1.0000 0.5000 # DT_SIGMA_MIN; DLNA_SIGMA_MIN
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.WINDOWS
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.WINDOWS 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.WINDOWS 2009-10-20 21:54:52 UTC (rev 15852)
@@ -1,15 +1,15 @@
3
-DATA_TEST/9818433.CI.MPM.BHR.sac.d.T006_T030
-SYN_TEST/MPM.CI.BHR.semd.sac.m16.T006_T030
+DATA/9818433.CI.MPM.BHR.sac.d.T006_T030
+SYN/MPM.CI.BHR.semd.sac.m16.T006_T030
3
11.2770 58.4770
58.4770 81.3770
81.3770 101.6770
-DATA_TEST/9818433.CI.MPM.BHT.sac.d.T006_T030
-SYN_TEST/MPM.CI.BHT.semd.sac.m16.T006_T030
+DATA/9818433.CI.MPM.BHT.sac.d.T006_T030
+SYN/MPM.CI.BHT.semd.sac.m16.T006_T030
1
49.3270 95.3270
-DATA_TEST/9818433.CI.MPM.BHZ.sac.d.T006_T030
-SYN_TEST/MPM.CI.BHZ.semd.sac.m16.T006_T030
+DATA/9818433.CI.MPM.BHZ.sac.d.T006_T030
+SYN/MPM.CI.BHZ.semd.sac.m16.T006_T030
1
60.9770 99.3770
Added: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/9818433_T006_T030_MPM_CI_m16_ikerXX_win_adj.pdf
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/9818433_T006_T030_MPM_CI_m16_ikerXX_win_adj.pdf
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Deleted: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/9818433_T006_T030_MPM_CI_m16_mtm_win_adj.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/make_pdf_by_event.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/make_pdf_by_event.pl 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/make_pdf_by_event.pl 2009-10-20 21:54:52 UTC (rev 15852)
@@ -1,10 +1,9 @@
#!/usr/bin/perl -w
#==========================================================
-#
# make_pdf_by_event.pl
# Carl Tape
-# 01-Oct-2008
+# 20-Oct-2009
#
# This sorts the output PDF files in PLOTS into order by
# 1 distance
@@ -31,15 +30,15 @@
($stanet,$stalon,$stlat,undef,undef,undef,undef,undef) = split(" ",$stalines[$ik-1]);
($station,$network) = split("\\.",$stanet);
- print "$ik out of $nrec : $station $network\n";
+ #print "$ik out of $nrec : $station $network\n";
@files = glob("*${station}_${network}*pdf");
$numf = @files;
if ($numf == 0) {
- print "--> no pdf file exists\n";
+ #print "$ik out of $nrec : $station $network --> no pdf file exists\n";
} elsif ($numf == 1) {
- print "--> one pdf file exists\n";
+ print "$ik out of $nrec : $station $network --> one pdf file exists\n";
$pdffile = $files[0]; chomp($pdffile);
$pdcat[$k] = $pdffile; $k = $k+1;
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2009-10-20 21:54:52 UTC (rev 15852)
@@ -2,7 +2,7 @@
#-----------------------------------
# Min Chen, 09-July-2007
-# Carl Tape, 02-Aug-2008
+# Carl Tape, 20-Oct-2009
# plot_win_adj.pl
#
# This script plots a single figure for a single source-receiver pair
@@ -12,16 +12,8 @@
# (Search for "USER".)
#
# EXAMPLE:
-# plot_win_adj.pl -m CMTSOLUTION_9818433 -k 2 -n HEC/CI -a STATIONS_ADJOINT -d DATA -s SYN -r RECON -w MEASUREMENT.WINDOWS -i m00
-# plot_win_adj.pl -m CMTSOLUTION_9818433 -k 2 -l 0/200 -n HEC/CI -a STATIONS_ADJOINT -d DATA -s SYN -r RECON -w MEASUREMENT.WINDOWS -i m00
-# plot_win_adj.pl -m CMTSOLUTION_9818433 -k 2 -r -l 0/200 -n HEC/CI -a STATIONS_ADJOINT -d DATA -s SYN -r RECON -w MEASUREMENT.WINDOWS -i m00
+# plot_win_adj.pl -m ../CMTSOLUTION_9818433 -n MPM/CI/BH -b 0 -l -10/200 -k 7 -a STATIONS_ADJOINT -d DATA -s SYN -c RECON -w MEASUREMENT.WINDOWS -i m16 -j 6/30
#
-# plot_win_adj.pl -m ../CMTSOLUTION_9818433 -k 1 -l 0/150 -n CLC/CI -a STATIONS_ADJOINT -d DATA_TEST -s SYN -r RECON_TEST -w MEASUREMENT.WINDOWS -i m00
-# plot_win_adj.pl -m ../CMTSOLUTION_14186612 -k 2 -r -l 0/200 -n SCI2/CI -a STATIONS_ADJOINT -d DATA_TEST -s SYN -r RECON_TEST -w MEASUREMENT.WINDOWS -i m00
-#
-# From the "EXAMPLE" directory:
-# plot_win_adj.pl -m CMTSOLUTION_9818433 -k 2 -l 0/200 -n CLC/CI -a STATIONS_ADJOINT -d DATA -s SYN -w MEASUREMENT.WINDOWS -i m00
-#
#-----------------------------------
use Getopt::Std;
@@ -37,7 +29,7 @@
-l -- Start/End of trace to be cut
-b -- plot MT and CC adjoint sources (=1)
-k -- type of adjoint source (0 = waveform; 1 = multi-taper; 2 = cross-correlation)
- -n -- station/network
+ -n -- station/network/chan
-r -- plot 3-compoment seismograms and adjoint sources in relatively scaled amplitude
-a -- station file to plot all the stations on the map
-d -- data directory
@@ -58,7 +50,7 @@
$iker = $opt_k; print "\n $iker $opt_k";
if ($opt_l) {($tmin,$tmax) = split(/\//,$opt_l);}
if ($opt_w) {$Winfile=$opt_w;}
-if ($opt_n) {($sta,$net)=split(/\//,$opt_n);}
+if ($opt_n) {($sta,$net,$chan)=split(/\//,$opt_n);}
if(!$opt_b) {$iboth = 0} else {$iboth = $opt_b}
if (!$opt_d) {$opt_d = "DATA"}
if (!$opt_s) {$opt_s = "SYN"}
@@ -68,12 +60,18 @@
$smodel = $opt_i;
# suffix for adjoint sources
-if($iker == 0) {$klab = "ik"; $ktitle = "waveform";}
-elsif($iker == 1) {$klab = "mtm"; $ktitle = "multitaper TT";}
-elsif($iker == 2) {$klab = "cc"; $ktitle = "cross-correlation TT";}
-elsif($iker == 3) {$klab = "bdcc"; $ktitle = "banana-doughnut CC-TT";}
-else {die("ERROR: iker ($iker) must be 0,1, 2, or 3\n");}
+$klab = sprintf('iker%2.2i',$iker);
+$klab2 = sprintf('iker%2.2i',$iker-2);
+ at ktitles = ("waveform, -d","waveform, s-d","banana-doughnut CC-TT","banana-doughnut dlnA",
+ "cross-correlation TT","amplitude dlnA","multitaper TT","multitaper dlnA");
+$ktitle = $ktitles[$iker-1];
+#if($iker == 0) {$klab = "ik"; $ktitle = "waveform";}
+#elsif($iker == 1) {$klab = "mtm"; $ktitle = "multitaper TT";}
+#elsif($iker == 2) {$klab = "cc"; $ktitle = "cross-correlation TT";}
+#elsif($iker == 3) {$klab = "bdcc"; $ktitle = "banana-doughnut CC-TT";}
+#else {die("ERROR: iker ($iker) must be 0,1, 2, or 3\n");}
+
$sTmin = sprintf("T%3.3i",$Tmin);
$sTmax = sprintf("T%3.3i",$Tmax);
$Ttag = "${sTmin}_${sTmax}";
@@ -92,16 +90,22 @@
# with iboth = 1, you need to have made both the CC and MT adjoint sources
if($iboth==1) {
- if($iker != 1) {die("\nFor plotting CC and MT adjoint sources, use iker = 1.\n")}
+ if( ($iker != 7) && ($iker != 8) ) {die("\nFor plotting CC and MT adjoint sources, use iker = 7 or 8.\n")}
$ktitle = "TT -- MT = black ; CC = red";
}
-# width of seismograms depends on whether you plot adjoint sources
+# plotting preferences
$iplot_adj = 1; # plot adjoint sources
$iplot_win = 1; # plot windows
$iplot_CClabel = 1; # plot CC measurement labels for each window
$iplot_recon = 1; # plot reconstructed synthetics
if($iplot_CClabel==1) {$iplot_win = 1;}
+
+# no reconstructed records or labels for waveform difference misfit
+if ($iker < 5) {$iplot_recon = 0; $iplot_CClabel = 0;}
+if ($iker == 0) {$iplot_adj = 0;}
+
+# adjust plot size based on whether you have adjoint sources
if($iplot_adj==1) {$swid = 3.5; $fsizetitle = 10;}
else {$swid = 7.4; $fsizetitle = 12;};
@@ -192,12 +196,12 @@
if(-f $Zdat) {$Zflag = 1} else {$Zflag = 0}
#print "\n Tflag : $Tflag, Rflag : $Rflag, Zflag : $Zflag \n";
-# SYNTHETICS -- component label always has the form BH_,
+# SYNTHETICS -- component label always has the form BH_ or LH_,
# even if the data file is something else (HH_, LH_).
-# We assume that ALL componenets exist.
-$Tlab = "$sta.$net.BHT";
-$Rlab = "$sta.$net.BHR";
-$Zlab = "$sta.$net.BHZ";
+# We assume that ALL components exist.
+$Tlab = "$sta.$net.${chan}T";
+$Rlab = "$sta.$net.${chan}R";
+$Zlab = "$sta.$net.${chan}Z";
$Tsyn = `ls ${opt_s}/${Tlab}*`; chomp($Tsyn);
$Rsyn = `ls ${opt_s}/${Rlab}*`; chomp($Rsyn);
$Zsyn = `ls ${opt_s}/${Zlab}*`; chomp($Zsyn);
@@ -565,10 +569,10 @@
#$red_pen = "-W1p,250/0/0";
#$black_pen = "-W1p,0/0/0";
- # adjoint sources (always have the BH_ channel, even if the data do not)
- $Zadj="ADJOINT_SOURCES/$sta.$net.BHZ.${klab}.adj";
- $Radj="ADJOINT_SOURCES/$sta.$net.BHR.${klab}.adj";
- $Tadj="ADJOINT_SOURCES/$sta.$net.BHT.${klab}.adj";
+ # adjoint sources (always have the BH_ or LH_ channel, even if the data do not)
+ $Zadj="ADJOINT_SOURCES/${sta}.${net}.${chan}Z.${klab}.adj";
+ $Radj="ADJOINT_SOURCES/${sta}.${net}.${chan}R.${klab}.adj";
+ $Tadj="ADJOINT_SOURCES/${sta}.${net}.${chan}T.${klab}.adj";
#$maxZ=1.0; $maxR=1.0; $maxT=1.0;
@@ -617,9 +621,10 @@
# adjoint sources for cross-correlation plot in red
if ($iboth==1) {
- $Zadj2="ADJOINT_SOURCES/$sta.$net.BHZ.cc.adj";
- $Radj2="ADJOINT_SOURCES/$sta.$net.BHR.cc.adj";
- $Tadj2="ADJOINT_SOURCES/$sta.$net.BHT.cc.adj";
+ $Zadj2="ADJOINT_SOURCES/$sta.$net.${chan}Z.${klab2}.adj";
+ $Radj2="ADJOINT_SOURCES/$sta.$net.${chan}R.${klab2}.adj";
+ $Tadj2="ADJOINT_SOURCES/$sta.$net.${chan}T.${klab2}.adj";
+
}
# VERTICAL component: adjoint source and windows
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl 2009-10-20 21:54:52 UTC (rev 15852)
@@ -5,7 +5,8 @@
# plot_all_win_adj.pl
#
# EXAMPLE:
-# plot_win_adj_all.pl -l 0/200 -m CMTSOLUTION_9818433 -b 0 -k 2 -a STATIONS_ADJOINT -d DATA -s SYN -w MEASUREMENT.WINDOWS -i m00 -j 3/30
+# plot_win_adj_all.pl -l -10/200 -m ../CMTSOLUTION_9818433 -n BH -b 0 -k 7 -a STATIONS_ADJOINT -d DATA -s SYN -c RECON -w MEASUREMENT.WINDOWS -i m16 -j 6/30
+#
#----------------------------------------------------
use Getopt::Std;
@@ -14,19 +15,20 @@
sub Usage{
print STDERR <<END;
- plot_all_win_adj.pl -m CMTFILE -a STATION_FILE -b iboth -l tmin/tmax -k iker -d data_dir -s syn_dir -c recon_dir -w winfile -i smodel -j Tmin/Tmax
+ plot_all_win_adj.pl -m CMTFILE -a STATION_FILE -n chan -b iboth -l tmin/tmax -k iker -d data_dir -s syn_dir -c recon_dir -w winfile -i smodel -j Tmin/Tmax
END
exit(1);
}
if (@ARGV == 0) { Usage(); }
-if (!getopts('m:a:b:d:k:l:s:c:w:i:j:')) {die('Check input arguments\n');}
+if (!getopts('m:a:n:b:d:k:l:s:c:w:i:j:')) {die('Check input arguments\n');}
if($opt_m) {$cmtfile=$opt_m;}
if(!$opt_b) {$iboth = 0} else {$iboth = $opt_b}
#if($opt_k) {$iker = $opt_k;} # does not work for iker = 1
if($opt_l) {$tcuts=$opt_l;}
if($opt_a) {$station_list=$opt_a;}
+if($opt_n) {$chan=$opt_n;}
if($opt_w) {$winfile=$opt_w;}
if($opt_d) {$data_dir=$opt_d;}
if($opt_s) {$syn_dir=$opt_s;}
@@ -45,8 +47,8 @@
print "$sta $net\n";
# include the -r command here if you want the traces scaled by the max Z-R-T value
- print "\n plot_win_adj.pl -m $cmtfile -n $sta/$net -b $iboth -l $tcuts -k $opt_k -a $station_list -d $data_dir -s $syn_dir -c $recon_dir -w $winfile -i $smodel -j $Ts\n";
- `plot_win_adj.pl -m $cmtfile -n $sta/$net -b $iboth -l $tcuts -k $opt_k -a $station_list -d $data_dir -s $syn_dir -c $recon_dir -w $winfile -i $smodel -j $Ts`;
+ print "\n plot_win_adj.pl -m $cmtfile -n $sta/$net/$chan -b $iboth -l $tcuts -k $opt_k -a $station_list -d $data_dir -s $syn_dir -c $recon_dir -w $winfile -i $smodel -j $Ts\n";
+ `plot_win_adj.pl -m $cmtfile -n $sta/$net/$chan -b $iboth -l $tcuts -k $opt_k -a $station_list -d $data_dir -s $syn_dir -c $recon_dir -w $winfile -i $smodel -j $Ts`;
}
Copied: seismo/3D/ADJOINT_TOMO/measure_adj/SYN/MPM.CI.BHR.semd.sac.m16.T006_T030 (from rev 15839, seismo/3D/ADJOINT_TOMO/measure_adj/SYN_TEST/MPM.CI.BHR.semd.sac.m16.T006_T030)
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/measure_adj/SYN/MPM.CI.BHR.semd.sac.m16.T006_T030
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/measure_adj/SYN/MPM.CI.BHT.semd.sac.m16.T006_T030 (from rev 15839, seismo/3D/ADJOINT_TOMO/measure_adj/SYN_TEST/MPM.CI.BHT.semd.sac.m16.T006_T030)
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/measure_adj/SYN/MPM.CI.BHT.semd.sac.m16.T006_T030
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/measure_adj/SYN/MPM.CI.BHZ.semd.sac.m16.T006_T030 (from rev 15839, seismo/3D/ADJOINT_TOMO/measure_adj/SYN_TEST/MPM.CI.BHZ.semd.sac.m16.T006_T030)
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/measure_adj/SYN/MPM.CI.BHZ.semd.sac.m16.T006_T030
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Name: svn:mergeinfo
+
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-10-20 21:54:52 UTC (rev 15852)
@@ -11,11 +11,9 @@
implicit none
- character(len=150) :: out_dir,datafile,synfile,file_prefix,file_prefix0,file_prefix2,measure_file_prefix,adj_file_prefix
-
- integer :: num_files, num_meas, i, j, k, iker, iker0, &
- ios, is_mtm, is_mtm0, npt1, npt2, npts, nn, nerr
- double precision, dimension(NDIM) :: data, syn, adj_syn_all, tr_adj_src, am_adj_src, recon_cc_all, tseis_recon_cc
+ character(len=150) :: datafile,synfile,file_prefix,file_prefix0,file_prefix2,measure_file_prefix,adj_file_prefix
+ integer :: num_files, num_meas, i, j, k, ios, npt1, npt2, npts, nn, nerr
+ double precision, dimension(NDIM) :: data, syn, adj_syn_all, tr_adj_src, am_adj_src, recon_cc_all, syn_dtw_cc
double precision :: t01, dt1, t02, dt2, t0, dt, tstart, tend, chi, tt, dtt, df
double precision, dimension(NCHI) :: window_chi
double precision :: fend0, fstart0, fend, fstart
@@ -24,11 +22,11 @@
! sac header information
integer :: yr,jda,ho,mi
double precision :: sec,dist,az,baz,slat,slon
- character(len=10) :: net,sta,chan_dat,chan,cmp
+ character(len=10) :: net,sta,chan_dat,chan,cmp,chan_syn
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 :: tshift, sigma_dt_cc, dlnA, sigma_dlnA_cc, sigma_dt, sigma_dlnA
double precision :: tr_chi, am_chi, cc_max, T_pmax_dat, T_pmax_syn
!double precision :: tshift_f1f2, cc_max_f1f2
double precision, dimension(NPT) :: dtau_w, dlnA_w, err_dt, err_dlnA, syn_dtw, data_dtw
@@ -39,81 +37,12 @@
!double precision :: trbdndw, a
!integer :: iord, passes
-
!********* PROGRAM STARTS HERE *********************
- ! input file MEASUREMENT.PAR
- ! is_mtm -- measurement type: 1 multi-taper; 2 cosine; 3 boxcar
- ! iker -- adjoint source type:
- ! 0 waveform adjoint sources
- ! 1 multi-taper adjoint sources
- ! 2 xcorr adjoint sources
- ! 3 banana doughnut adjoint sources
- ! TLONG, TSHORT -- starting and ending period of the signals
- ! tt, dtt, nn -- adj source interpolation scheme
- ! out_dir -- seismograms output directory
- ! DISPLAY_DETAILS -- controls display of detailed calculation process
- ! OUTPUT_MEASUREMENT_FILES -- controls if output dtau_w, error_dt into files, etc.
+ ! read in MEASUREMENT.PAR (see mt_sub.f90 and write_par_file.pl)
+ ! most variables are global (see mt_variables.f90)
+ call read_par_file(fstart0,fend0,tt,dtt,nn,chan)
- open(10,file='MEASUREMENT.PAR',status='old',iostat=ios)
- read(10,'(a)') out_dir
- read(10,*) is_mtm0
- read(10,*) wtr,npi
- read(10,*) iker0
- read(10,*) RUN_BANDPASS
- read(10,*) TLONG, TSHORT
- read(10,*) tt,dtt,nn
- read(10,*) DISPLAY_DETAILS
- read(10,*) OUTPUT_MEASUREMENT_FILES
- read(10,*) INCLUDE_ERROR
- read(10,*) DT_FAC
- read(10,*) ERR_FAC
- read(10,*) DT_MAX_SCALE
- read(10,*) NCYCLE_IN_WINDOW
- read(10,*) BEFORE_QUALITY, AFTER_QUALITY
- read(10,*) BEFORE_TSHIFT, AFTER_TSHIFT
- read(10,*) DT_SIGMA_MIN, DLNA_SIGMA_MIN
- close(10)
-
- out_dir = adjustl(out_dir)
- iker = iker0
- is_mtm = is_mtm0
-
- ! check the read-in values
- print *, 'INPUTS FROM MEASUREMENT.PAR :'
- print *, ' is_mtm : ',is_mtm
- print *, ' wtr, npi : ',wtr,npi
- print *, ' iker : ',iker
- print *, ' RUN_BANDPASS :',RUN_BANDPASS
- print *, ' TLONG, TSHORT : ',TLONG, TSHORT
- fstart0 = 1./TLONG ; fend0 = 1./TSHORT
- print *, ' fstart, fend :', fstart0, fend0
- print *, ' tt, dtt, nn : ',tt,dtt,nn
- print *, ' out_dir : ',trim(out_dir)
- print *, ' DISPLAY_DETAILS :',DISPLAY_DETAILS
- print *, ' OUTPUT_MEASUREMENT_FILES :',OUTPUT_MEASUREMENT_FILES
- print *, ' INCLUDE_ERROR :',INCLUDE_ERROR
- print *, ' DT_FAC :',DT_FAC
- print *, ' ERR_FAC :',ERR_FAC
- print *, ' DT_MAX_SCALE :',DT_MAX_SCALE
- print *, ' NCYCLE_IN_WINDOW :',NCYCLE_IN_WINDOW
- print *, ' BEFORE_QUALITY, AFTER_QUALITY :',BEFORE_QUALITY, AFTER_QUALITY
- print *, ' BEFORE_TSHIFT, AFTER_TSHIFT :',BEFORE_TSHIFT, AFTER_TSHIFT
- print *, ' DT_SIGMA_MIN, DLNA_SIGMA_MIN :',DT_SIGMA_MIN, DLNA_SIGMA_MIN
- !stop 'checking PAR file input'
-
- if (fstart0 >= fend0) then ; print *, 'Check input frequency range of the signal' ; stop ; endif
- if (iker == 1 .and. is_mtm /= 1) then ; print *, 'Error: make MT measurements before compute MT adjoint source' ; stop ; endif
- if (iker == 0 .and. is_mtm /= 0) then ; print *, 'Error: iker=0 and is_mtm=0 for waveform difference' ; stop ; endif
- if (iker /= 0 .and. is_mtm == 0) then ; print *, 'Error: iker=0 and is_mtm=0 for waveform difference' ; stop ; endif
- if (nn > NDIM) then ; print *, 'Error: Change interpolation nn or NDIM' ; stop ; endif
-
- ! apply filter (this should EXACTLY match the filter used in the windowing code)
- !trbdndw = 0.3
- !a = 30.
- !iord = 4
- !passes = 2
-
! input file: MEASUREMENT.WINDOWS
open(11,file='MEASUREMENT.WINDOWS',status='old',iostat=ios)
if (ios /= 0) then ; print *, 'Error opening input file: MEASUREMENT WINDOWS' ; stop ; endif
@@ -139,14 +68,6 @@
read(11,'(a)',iostat=ios) synfile
if (ios /= 0) then ; print *, 'Error reading windows file' ; stop ; endif
-!!$ ! read in data (SAC format)
-!!$ !call dread_ascfile_f(datafile,t01,dt1,npt1,data) ! ASCII format
-!!$ call dread_sacdata_f(datafile,t01,dt1,npt1,data)
-!!$
-!!$ ! read in synthetics (SAC format)
-!!$ !call dread_ascfile_f(synfile,t02,dt2,npt2,syn) ! ASCII format
-!!$ call dread_sacdata_f(synfile,t02,dt2,npt2,syn)
-
! read data and syn
call drsac1(datafile,data,npt1,t01,dt1)
call drsac1(synfile,syn,npt2,t02,dt2)
@@ -176,46 +97,34 @@
call bandpass(data,npts,dt,fstart0,fend0)
endif
-!!$ ! figure out station name, network name, component name
-!!$ call saclst_kheader_f(trim(datafile),'kstnm',sta,sta_len)
-!!$ sta(sta_len+1:) = ' '
-!!$ call saclst_kheader_f(trim(datafile),'knetwk',net,net_len)
-!!$ net(net_len+1:) = ' '
-!!$ call saclst_kheader_f(trim(datafile),'kcmpnm',chan,chan_len)
-!!$ chan(chan_len+1:) = ' '
-
! figure out station name, network name, component name, etc
call get_sacfile_header(trim(datafile),yr,jda,ho,mi,sec,net,sta,chan_dat,dist,az,baz,slat,slon)
- chan = chan_dat
- ! We assume that the synthetics always have the form BH_, but the data may not (HH_, LH_, BL_, etc).
- ! Here we adjust the channel to ALWAYS be of the form BH_ for simplicity;
- ! for flexibility, we might want to put the leading two characters in the Par_file.
- if(1==1) then
- cmp = chan_dat(3:3)
- chan = 'BH' // cmp
- !print *, chan_dat, cmp, chan
- endif
+ ! synthetics always have the form BH_ or LH_, but the data may not (HH_, LH_, BL_, etc).
+ cmp = chan_dat(3:3)
+ chan_syn = trim(chan)//trim(cmp)
- file_prefix0 = trim(sta)//'.'//trim(net)//'.'//trim(chan)
- file_prefix2 = trim(out_dir)//'/'//trim(file_prefix0)
-
- if (iker == 0) then
- adj_file_prefix = trim(file_prefix2) // '.ik'
- else if (iker == 1) then
- adj_file_prefix = trim(file_prefix2) // '.mtm'
- else if (iker == 2) then
- adj_file_prefix = trim(file_prefix2) // '.cc'
- else if (iker == 3) then
- adj_file_prefix = trim(file_prefix2) // '.bdcc'
- else
- print *, 'iker = ', iker
- print *, 'Error: iker must be 0, 1, 2, or 3'
- stop
- endif
+ file_prefix0 = trim(sta)//'.'//trim(net)//'.'//trim(chan_syn)
+ file_prefix2 = trim(OUT_DIR)//'/'//trim(file_prefix0)
print *
print *, trim(file_prefix2), ' --- '
+ ! note: MT measurement could revert to CC, but still keep the MT suffix
+ write(adj_file_prefix,'(a,i2.2)') trim(file_prefix2)//'.iker', iker0
+!!$ if (iker == 0) then
+!!$ adj_file_prefix = trim(file_prefix2) // '.ik'
+!!$ else if (iker == 1) then
+!!$ adj_file_prefix = trim(file_prefix2) // '.mtm'
+!!$ else if (iker == 2) then
+!!$ adj_file_prefix = trim(file_prefix2) // '.cc'
+!!$ else if (iker == 3) then
+!!$ adj_file_prefix = trim(file_prefix2) // '.bdcc'
+!!$ else
+!!$ print *, 'iker = ', iker
+!!$ print *, 'Error: iker must be 0, 1, 2, or 3'
+!!$ stop
+!!$ endif
+
read(11,*,iostat=ios) num_meas
if (ios /= 0) then ; print *, 'Error reading num_meas' ; stop ; endif
@@ -230,13 +139,13 @@
! write values to output file
nwin = nwin + 1 ! overall window counter
- write(12,'(a3,a8,a5,a5,3i5,2f12.3)') net,sta,chan,chan_dat,nwin,ipair,j,tstart,tend
+ write(12,'(a3,a8,a5,a5,3i5,2f12.3)') net,sta,chan_syn,chan_dat,nwin,ipair,j,tstart,tend
+ ! add taper type to file prefix
write(file_prefix,'(a,i2.2)') trim(file_prefix2)//'.', j
- ! add measurement type to file prefix
if (is_mtm == 1) then
measure_file_prefix = trim(file_prefix) // '.mtm' ! multitaper taper
- else if (is_mtm == 2) then
+ elseif (is_mtm == 2) then
measure_file_prefix = trim(file_prefix) // '.ctp' ! cosine taper
else
measure_file_prefix = trim(file_prefix) // '.btp' ! boxcar taper
@@ -264,8 +173,8 @@
! make measurements
! CHT: also compute reconstructed synthetics for CC (and MT, if specified) measurements
- call mt_measure(is_mtm,iker,datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,&
- istart,data_dtw,syn_dtw,nlen,tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tseis_recon_cc,&
+ call mt_measure(datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,&
+ istart,data_dtw,syn_dtw,nlen,tshift,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,syn_dtw_cc,&
i_pmax_dat,i_pmax_syn,i_right0,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,err_dt,err_dlnA)
i_right = i_right0
i_left = 1
@@ -274,36 +183,22 @@
T_pmax_dat = (dt*NPT) / dble(i_pmax_dat)
T_pmax_syn = (dt*NPT) / dble(i_pmax_syn)
- ! adjust measurements for adjoint source
- if (iker == 1) then ! CHT: is_mtm --> iker
+ ! adjust measurements for MT adjoint source
+ if (is_mtm == 1) then
fstart = fstart0 ; fend = fend0
- call mt_measure_select(nlen,tshift_xc,i_pmax_syn,dtau_w,dlnA_w,err_dt,err_dlnA,dt,i_left,i_right,fstart,fend,use_trace)
-
+ call mt_measure_select(nlen,tshift,i_pmax_syn,dtau_w,dlnA_w,err_dt,err_dlnA,dt,i_left,i_right,fstart,fend,use_trace)
print *, 'fstart0/fend0 :', fstart0, fend0
print *, 'fstart/fend :', fstart, fend
print *, 'Tpmax :', T_pmax_dat, T_pmax_syn
-!!$ 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, 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
! CHT: If MT measurement window is rejected by mt_measure_select, then use a CC measurement.
- ! --> May also want to consider quality checks for CC as well.
if(.not. use_trace) then
- !stop 'Ceck why this MT measurement was rejected'
+ !stop 'Check why this MT measurement was rejected'
print *, 'Reverting from multitaper measurement to cross-correlation measurement'
- iker = 2
+ iker = iker0 - 2
is_mtm = 3
- call mt_measure(is_mtm,iker,datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,istart,data_dtw,syn_dtw,nlen,&
- tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tseis_recon_cc,&
+ call mt_measure(datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,istart,data_dtw,syn_dtw,nlen,&
+ tshift,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,syn_dtw_cc,&
i_pmax_dat,i_pmax_syn,i_right,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA)
use_trace = .true.
endif
@@ -312,19 +207,8 @@
use_trace = .true.
endif
- ! CHT: If the CC timeshift is for some reason larger than the allowable max,
- ! then effectively eliminate the window by zeroing the
- ! cross-correlation traveltime and amplitude measurements.
- ! See BEFORE_TSHIFT in subroutine compute_cc in mt_sub.f90.
- ! NOTE: BEFORE_DLNA should be included in the PARFILE
- if ( abs(tshift_xc) > BEFORE_TSHIFT ) then
- if (DISPLAY_DETAILS) then
- print *, ' effectively eliminating this time window -->'
- print *, ' tshift_xc : ', tshift_xc, BEFORE_TSHIFT, abs(tshift_xc) > BEFORE_TSHIFT
- endif
- tshift_xc = 0.0
- dlnA = 0.0
- endif
+ ! check that the CC measurements are within the specified input range
+ if (iker >= 5) call cc_measure_select(tshift,dlnA,cc_max)
! write frequency limits to file
if (OUTPUT_MEASUREMENT_FILES) then
@@ -338,13 +222,9 @@
if (use_trace) then
tr_chi = 0.0 ; am_chi = 0.0 ! must be initialized
- if (iker==0) then
- call mt_adj(iker,istart,data_dtw,syn_dtw,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)
- else
- call mt_adj(iker,istart,data_dtw,syn_dtw,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,am_adj_src,am_chi)
- endif
+ call mt_adj(istart,data_dtw,syn_dtw,nlen,dt,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc,&
+ dtau_w,dlnA_w,err_dt,err_dlnA,sigma_dt,sigma_dlnA,i_left,i_right,&
+ window_chi,tr_adj_src,tr_chi,am_adj_src,am_chi)
! KEY: write misfit function values to file (two for each window)
! Here are the 20 columns of the vector window_chi
@@ -354,20 +234,27 @@
! WINDOW : 13: data power, 14: syn power, 15: (data-syn) power, 16: window duration
! FULL RECORD: 17: data power, 18: syn power, 19: (data-syn) power, 20: record duration
! Example of a reduced file: awk '{print $2,$3,$4,$5,$6,$31,$32}' window_chi > window_chi_sub
- write(13,'(a14,a8,a3,a5,i4,i4,2e14.6,20e14.6,2e14.6,2f14.6)') file_prefix0,sta,net,chan,j,iker,&
+ write(13,'(a14,a8,a3,a5,i4,i4,2e14.6,20e14.6,2e14.6,2f14.6)') &
+ file_prefix0,sta,net,chan_syn,j,iker,&
tstart,tend,window_chi(:),tr_chi,am_chi,T_pmax_dat,T_pmax_syn
print *, ' tr_chi = ', sngl(tr_chi), ' am_chi = ', sngl(am_chi)
! combine adjoint sources from different measurement windows
- adj_syn_all(:) = adj_syn_all(:) + tr_adj_src(:)
+ if (COMPUTE_ADJOINT_SOURCE) then
+ if (mod(iker,2)==1) then
+ adj_syn_all(:) = adj_syn_all(:) + tr_adj_src(:) ! iker = 1,3,5,7
+ else
+ adj_syn_all(:) = adj_syn_all(:) + am_adj_src(:) ! iker = 2,4,6,8
+ endif
+ endif
! combine CC-reconstructed records
- recon_cc_all(istart:istart+nlen-1) = recon_cc_all(istart:istart+nlen-1) + tseis_recon_cc(1:nlen)
+ recon_cc_all(istart:istart+nlen-1) = recon_cc_all(istart:istart+nlen-1) + syn_dtw_cc(1:nlen)
endif
! CHT: (re-)set to multitaper parameters, if originally specified
- if (iker0 == 1) then
+ if (is_mtm0 == 1) then
iker = iker0
is_mtm = is_mtm0
endif
@@ -377,26 +264,32 @@
!----------------------------
! write out the adjoint source
- ! filter the adjoint source (higher frequencies could enter from the tapering operations)
- ! Note: time_window in mt_adj.f90 tapers the windows.
- call bandpass(adj_syn_all,npts,dt,fstart0,fend0)
+ if (COMPUTE_ADJOINT_SOURCE) then
- ! cut and interpolate to match time-stepping for SEM
- ! NOTE: This can leave a non-zero value to start the record,
- ! which is NOT GOOD for the SEM simulation.
- call interpolate_syn(adj_syn_all,t0,dt,npts,tt,dtt,nn)
+ ! OPTIONAL: A conservative choice is to filter the adjoint source,
+ ! since higher frequencies could enter from the tapering operations.
+ ! Note: time_window in mt_adj.f90 tapers the windows.
+ !call bandpass(adj_syn_all,npts,dt,fstart0,fend0)
- ! Taper the start of the adjoint source, since the cutting of the record
- ! may have left a non-zero value to start the record, which is NOT GOOD for the SEM simulation.
- itmax = int(TSHORT/dtt)
- call taper_start(adj_syn_all,nn,itmax)
+ ! cut and interpolate to match time-stepping for SEM
+ ! NOTE: This can leave a non-zero value to start the record,
+ ! which is NOT GOOD for the SEM simulation.
+ call interpolate_syn(adj_syn_all,t0,dt,npts,tt,dtt,nn)
- ! we can have choices of outputing the adjoint source as ASCII or SAC format
- call dwascii(trim(adj_file_prefix)//'.adj',adj_syn_all,nn,tt,dtt)
- !call dwsac1(trim(adj_file_prefix)//'.adj.sac',adj_syn_all,nn,tt,dtt)
+ ! Taper the start of the adjoint source, since cutting the record
+ ! may have left a non-zero value to start the record,
+ ! which is not good for the SEM simulation.
+ itmax = int(TSHORT/dtt)
+ call taper_start(adj_syn_all,nn,itmax)
+
+ ! we can have choices of outputing the adjoint source as ASCII or SAC format
+ call dwascii(trim(adj_file_prefix)//'.adj',adj_syn_all,nn,tt,dtt)
+ !call dwsac1(trim(adj_file_prefix)//'.adj.sac',adj_syn_all,nn,tt,dtt)
!!$ call dwrite_ascfile_f(trim(adj_file_prefix)//'.adj',tt,dtt,nn,adj_syn_all)
!!$ !call dwrite_sacfile_f(trim(datafile),trim(adj_file_prefix)//'.adj',tt,nn,adj_syn_all)
+ endif
+
!----------------------------
! write out the CC-reconstructed data from synthetics
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-10-20 21:54:52 UTC (rev 15852)
@@ -32,19 +32,13 @@
! upgraded to Fortran90 by Alessia Maggi
! organized into package form by Qinya Liu
! modifications by Carl Tape and Vala Hjorleifsdottir
- !
- ! Evolution of data and syn arrays:
- ! data -- window + time taper --> dataw -- time shift --> dataw0 -- complex FFT --> wseis_dat
- ! dataw0 -- single-tapered --> dataw -- complex FFT --> wseis2
- ! syn -- window + time taper --> synw = dataw0 -- complex FFT --> wseis_syn
- ! synw0 -- single-tapered --> synw -- complex FFT --> wseis
! =================================================================================================
- subroutine mt_measure(is_mtm,iker,datafile,file_prefix,dat_dt,syn_dt,t0,dt,npts,tstart,tend, &
+ subroutine mt_measure(datafile,file_prefix,dat_dt,syn_dt,t0,dt,npts,tstart,tend, &
istart,dat_dtw,syn_dtw,nlen,tshift,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,syn_dtw_cc, &
i_pmax_dat,i_pmax_syn,i_right,trans_w,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,err_dt,err_dlnA)
- integer, intent(in) :: is_mtm,iker,npts
+ integer, intent(in) :: npts
double precision, dimension(:), intent(in) :: dat_dt,syn_dt
double precision, intent(in) :: tstart,tend,t0,dt
character(len=150), intent(in) :: file_prefix,datafile
@@ -60,7 +54,7 @@
syn_dtw_cc_dt, dat_dtw_cc, syn_dtw_h, dat_dtw_h
double precision :: sfac1,fac,omega,f0,df,df_new,dw,&
ampmax_unw,wtr_use_unw,ampmax,wtr_use,wtr_mtm,dtau_wa,dlnA_wa
- integer :: ishift,i,ictaper,j,k,fnum,i_amp_max_unw,i_amp_max,i_right_stop,idf_new,iom,iker_temp
+ integer :: ishift,i,ictaper,j,k,fnum,i_amp_max_unw,i_amp_max,i_right_stop,idf_new,iom
complex*16, dimension(NPT) :: syn_dtwo, dat_dtwo, syn_dtw_ho, dat_dtw_ho, &
top_mtm, bot_mtm, trans_mtm, wseis_recon
@@ -139,7 +133,7 @@
! cross-correlation traveltime and amplitude measurements
!------------------------------------------------------------------
- ! compute cross-correltaion time shift and also amplitude measurmement
+ ! compute cross-correlation time shift and also amplitude measurmement
! NOTE: records have already been windowed, so no information outside windows is considered
! LQY: Ying suggested to align them at relatively long periods
call compute_cc(syn_dtw, dat_dtw, nlen, dt, ishift, tshift, dlnA, cc_max)
@@ -162,7 +156,7 @@
call deconstruct_dat_cc(filename,dat_dtw,tstart,dt,nlen, &
ishift,tshift,dlnA,dat_dtw_cc)
- ! reconstruct synthetics using cross-correlation measurments
+ ! reconstruct synthetics using cross-correlation measurments (plotting purposes only)
call reconstruct_syn_cc(syn_dtw,tstart,dt,nlen,ishift,tshift,dlnA,syn_dtw_cc,syn_dtw_cc_dt)
if (OUTPUT_MEASUREMENT_FILES) then
@@ -180,7 +174,7 @@
call compute_average_error(dat_dtw,syn_dtw_cc,syn_dtw_cc_dt,nlen,dt,sigma_dt_cc,sigma_dlnA_cc)
! write cross-correlation measurement to file
- call write_average_meas(filename,IKER_CC,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc)
+ call write_average_meas(filename,iker,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc)
!========================================
@@ -296,6 +290,13 @@
elseif (is_mtm == 3) then
call boxcar(nlen, NPT, tas)
endif
+!!$ if (is_mtm == 1) then
+!!$ call staper(nlen, NPI, NTAPER, tas, NPT, ey1, ey2)
+!!$ elseif (is_mtm == 2) then
+!!$ call costaper(nlen, NPT, tas)
+!!$ elseif (is_mtm == 3) then
+!!$ call boxcar(nlen, NPT, tas)
+!!$ endif
! initialize transfer function terms
top_mtm(:) = cmplx(0.,0.)
@@ -346,7 +347,7 @@
! for cosine or boxcar tapers only -- SEE COMMENTS BELOW for the multitaper case
! NOTE 1: here we are using trans_w, not trans_mtm
- ! NOTE 2: The single-taper transfer function should give you a PERFECT FIT,
+ ! NOTE 2: The single-taper transfer function should give you a perfect fit,
! but it is not relevant from the perspective of obtaining a measurement.
if (is_mtm /= 1) then
! phase, abs(trans), travel-time and amplitude as a func of freq for single-tapered measurements
@@ -356,13 +357,13 @@
i_right,tstart,dt,nlen,syn_dtw_mt, syn_dtw_mt_dt)
!call check_recon_quality(filename,dat_dtw_cc,syn_dtw,dat_dtw,syn_dtw_mt,nlen,dt,tshift,tshift_f1f2,cc_max_f1f2,cc_max)
!call compute_average_error(dat_dtw,syn_dtw_mt,syn_dtw_mt_dt,nlen,dt,sigma_dt,sigma_dlnA)
- !call write_average_meas(filename, IKER_MT, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
+ !call write_average_meas(filename, iker, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
endif
enddo ! ictapers
! for single taper, pass back the transfer function
- if (is_mtm /= 1) return
+ if (is_mtm /= 1) return
!-------------------------------------------------------------------------------
! multitaper estimation of transfer function
@@ -595,20 +596,20 @@
! original coding by Carl Tape, finalized by Qinya Liu
! ======================================================================================================
- subroutine mt_adj(iker,istart,dat_dtw,syn_dtw,nlen,dt,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc,&
+ subroutine mt_adj(istart,dat_dtw,syn_dtw,nlen,dt,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc, &
dtau_w,dlnA_w,err_dtau,err_dlnA,sigma_dt,sigma_dlnA,i_left,i_right, &
- tr_adj_src,tr_chi,window_chi,am_adj_src,am_chi)
+ window_chi,tr_adj_src,tr_chi,am_adj_src,am_chi)
- integer, intent(in) :: iker, istart, nlen, i_left, i_right
+ integer, intent(in) :: istart, nlen, i_left, i_right
double precision, dimension(:), intent(in) :: dat_dtw, syn_dtw
double precision, intent(in) :: dt, tshift, dlnA, sigma_dt_cc, sigma_dlnA_cc, sigma_dt, sigma_dlnA
double precision, dimension(:), intent(in) :: dtau_w, dlnA_w, err_dtau, err_dlnA
- double precision, dimension(:), intent(out) :: tr_adj_src
- double precision, intent(out) :: tr_chi
+ double precision, dimension(:), intent(out) :: tr_adj_src, am_adj_src
+ double precision, intent(out) :: tr_chi, am_chi
double precision, dimension(NCHI), intent(inout) :: window_chi
- double precision, dimension(:), intent(out), optional :: am_adj_src
- double precision, intent(out), optional :: am_chi
+ !double precision, dimension(:), intent(out), optional :: am_adj_src
+ !double precision, intent(out), optional :: am_chi
double precision, dimension(NPT) :: syn_vtw, syn_vtw_h, syn_dtw_h, ey1, ey2
double precision, dimension(NPT) :: ft_bar_t, fa_bar_t, fp, fq, wp_taper, wq_taper
@@ -623,16 +624,27 @@
double precision :: waveform_chi, waveform_d2, waveform_s2, waveform_temp1, waveform_temp2, waveform_temp3
! waveform adjoint source is passed by tr_adj_src and tr_chi
- if (iker == 0 .and. (present(am_adj_src) .or. present(am_chi))) stop &
- 'am_adj_src and am_chi are not needed for iker = 0 (waveform adjoint source case)'
+ !if (iker == 0 .and. (present(am_adj_src) .or. present(am_chi))) stop &
+ ! 'am_adj_src and am_chi are not needed for iker = 0 (waveform adjoint source case)'
! check the window length
if (istart + nlen > NDIM) stop 'Check istart + nlen and NPT'
+ ! waveform
+ if(iker==1 .or. iker==2) then
+ print *, 'computing waveform adjoint source'
+ elseif(iker==3 .or. iker==4) then
+ print *, 'computing banana-doughtnut adjoint source'
+ elseif(iker==5 .or. iker==6) then
+ print *, 'computing cross-correlation adjoint source'
+ elseif(iker==7 .or. iker==8) then
+ print *, 'computing multitaper adjoint source'
+ endif
+
! ----------------------
! TAPERS
! ----------------------
- if( iker == 1 ) then
+ if( is_mtm == 1 ) then
! frequency-domain tapers
! THIS CHOICE WILL HAVE AN EFFECT ON THE ADJOINT SOURCES
ipwr_w = 10
@@ -659,15 +671,15 @@
do i = i_left, i_right ! CHT: 1 --> i_left
- if (INCLUDE_ERROR == 0) then ! no error estimate
+ if (ERROR_TYPE == 0) then ! no error estimate
wp_taper(i) = w_taper(i) / ffac
wq_taper(i) = w_taper(i) / ffac
- elseif (INCLUDE_ERROR == 1) then ! MT error estimate is assigned the CC error estimate
+ elseif (ERROR_TYPE == 1) then ! MT error estimate is assigned the CC error estimate
wp_taper(i) = w_taper(i) / ffac / (sigma_dt ** 2)
wq_taper(i) = w_taper(i) / ffac / (sigma_dlnA ** 2)
- elseif (INCLUDE_ERROR == 2) then ! MT jack-knife error estimate
+ elseif (ERROR_TYPE == 2) then ! MT jack-knife error estimate
err_t = err_dtau(i)
if (err_dtau(i) < dtau_wtr) err_t = err_t + dtau_wtr
err_A = err_dlnA(i)
@@ -683,7 +695,7 @@
!!$ enddo
!!$ close(88)
- endif ! iker == 1
+ endif ! is_mtm == 1
! post-processing time-domain taper
@@ -704,7 +716,7 @@
! ----------------------------------
! CROSS CORRELATION ADJOINT SOURCES
! ----------------------------------
- if( iker == 2 .or. iker == 3 ) then
+ if( (iker >= 3).and.(iker <= 6) ) then
! compute synthetic velocity
do i = 2, nlen-1
@@ -728,7 +740,7 @@
! MULTITAPER ADJOINT SOURCES
! -------------------------------
- if (iker == 1) then
+ if (is_mtm == 1) then
! allocate MT variables
ntaper = int(NPI * 2.0)
@@ -802,13 +814,11 @@
endif ! MT adjoint source
! -------------------------------------
- ! ASSIGN ADJOINT SOURCE and CHI
+ ! COMPUTE ADJOINT SOURCE
! -------------------------------------
- ! compute adjoint source
-
tr_adj_src = 0.
- if(iker /= 0) am_adj_src = 0.
+ am_adj_src = 0.
! CHT: Compute the integrated waveform difference squared
waveform_temp1 = 0. ; waveform_temp2 = 0. ; waveform_temp3 = 0.
@@ -818,34 +828,42 @@
waveform_temp3 = waveform_temp3 + (( dat_dtw(i) - syn_dtw(i) ) * time_window(i) )**2
enddo
! NOTE: does not include DT factor or normalization by duration of window
- waveform_d2 = 0.5 * waveform_temp1
- waveform_s2 = 0.5 * waveform_temp2
- waveform_chi = 0.5 * waveform_temp3
+ waveform_d2 = waveform_temp1
+ waveform_s2 = waveform_temp2
+ waveform_chi = waveform_temp3
! compute traveltime and amplitude adjoint sources
do i = 1,nlen
i1 = istart + i
- if(iker==0) then ! waveform
- tr_adj_src(i1) = ( syn_dtw(i) - dat_dtw(i) ) * time_window(i)
+ ! waveform
+ if(iker==1 .or. iker==2) then
+ tr_adj_src(i1) = -dat_dtw(i)/waveform_d2 * time_window(i)
+ am_adj_src(i1) = ( syn_dtw(i) - dat_dtw(i) ) * time_window(i) ! consider normalizing by waveform_d2
- elseif(iker==1) then ! multitaper
- tr_adj_src(i1) = fp(i) * time_window(i)
- am_adj_src(i1) = fq(i) * time_window(i)
+ ! banana-doughnut kernel adjoint source (no measurement)
+ elseif(iker==3 .or. iker==4) then
+ tr_adj_src(i1) = ft_bar_t(i) * time_window(i)
+ am_adj_src(i1) = fa_bar_t(i) * time_window(i)
- elseif(iker==2) then ! cross-correlation
+ ! cross-correlation
+ elseif(iker==5 .or. iker==6) then
tr_adj_src(i1) = -(tshift / sigma_dt_cc**2 ) * ft_bar_t(i) * time_window(i)
am_adj_src(i1) = -(dlnA / sigma_dlnA_cc**2 ) * fa_bar_t(i) * time_window(i)
- ! banana-doughnut kernel adjoint source (no measurement)
- elseif(iker==3) then
- tr_adj_src(i1) = ft_bar_t(i) * time_window(i)
- am_adj_src(i1) = fa_bar_t(i) * time_window(i)
+ ! multitaper
+ elseif(iker==7 .or. iker==8) then
+ tr_adj_src(i1) = fp(i) * time_window(i)
+ am_adj_src(i1) = fq(i) * time_window(i)
endif
enddo
+ ! -------------------------------------
+ ! COMPUTE MISFIT FUNCTION VALUE
+ ! -------------------------------------
+
! CHT: compute misfit function value and measurement value
! Note: The taper functions for MT may include error estimates.
! 1: multitaper, TT
@@ -855,27 +873,27 @@
!window_chi(:) = 0.
! misfit function value
- if(iker==1) window_chi(1) = 0.5 * 2.0 * df * sum( (dtau_w(1:i_right))**2 * wp_taper(1:i_right) )
- if(iker==1) window_chi(2) = 0.5 * 2.0 * df * sum( (dlnA_w(1:i_right))**2 * wq_taper(1:i_right) )
+ if(is_mtm==1) window_chi(1) = 0.5 * 2.0 * df * sum( (dtau_w(1:i_right))**2 * wp_taper(1:i_right) )
+ if(is_mtm==1) window_chi(2) = 0.5 * 2.0 * df * sum( (dlnA_w(1:i_right))**2 * wq_taper(1:i_right) )
window_chi(3) = 0.5 * (tshift/sigma_dt_cc)**2
window_chi(4) = 0.5 * (dlnA/sigma_dlnA_cc)**2
! measurement (no uncertainty estimates)
- if(iker==1) window_chi(5) = sum( dtau_w(1:i_right) * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
- if(iker==1) window_chi(6) = sum( dlnA_w(1:i_right) * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
+ if(is_mtm==1) window_chi(5) = sum( dtau_w(1:i_right) * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
+ if(is_mtm==1) window_chi(6) = sum( dlnA_w(1:i_right) * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
window_chi(7) = tshift
window_chi(8) = dlnA
! estimated mesurement uncertainties
- if(iker==1) window_chi(9) = sigma_dt
- if(iker==1) window_chi(10) = sigma_dlnA
+ if(is_mtm==1) window_chi(9) = sigma_dt
+ if(is_mtm==1) window_chi(10) = sigma_dlnA
window_chi(11) = sigma_dt_cc
window_chi(12) = sigma_dlnA_cc
! for normalization, divide by duration of window
- window_chi(13) = waveform_d2
- window_chi(14) = waveform_s2
- window_chi(15) = waveform_chi
+ window_chi(13) = 0.5 * waveform_d2
+ window_chi(14) = 0.5 * waveform_s2
+ window_chi(15) = 0.5 * waveform_chi
window_chi(16) = nlen*dt
!!$ open(88,file='testing.dat')
@@ -884,17 +902,18 @@
!!$ enddo
!!$ close(88)
- if(iker==0) then ! waveform
- tr_chi = waveform_chi
+ if(iker <= 2) then ! waveform
+ tr_chi = 0.5 * waveform_chi
+ am_chi = 0.5 * waveform_chi
- elseif(iker==1) then ! multitaper
+ elseif( (iker >= 3).and.(iker <= 6) ) then ! cross_correlation
+ tr_chi = window_chi(3)
+ am_chi = window_chi(4)
+
+ elseif( (iker==7).or.(iker==8) ) then ! multitaper
tr_chi = window_chi(1)
am_chi = window_chi(2)
- elseif(iker==2 .or. iker==3) then ! cross_correlation
- tr_chi = window_chi(3)
- am_chi = window_chi(4)
-
endif
end subroutine mt_adj
@@ -1009,6 +1028,34 @@
!-----------------------------------------------------------------------------
+ subroutine cc_measure_select(tshift,dlnA,cc_max)
+
+ ! CHT: If the CC timeshift is for some reason larger than the allowable max,
+ ! then effectively eliminate the window by zeroing the
+ ! cross-correlation traveltime and amplitude measurements.
+ ! See subroutine compute_cc in mt_sub.f90.
+
+ double precision, intent(inout) :: tshift, dlnA, cc_max
+
+ if( (cc_max < CC_MIN) .or. (tshift < TSHIFT_MIN) .or. (tshift > TSHIFT_MAX) &
+ .or. (dlnA < DLNA_MIN) .or. (dlnA > DLNA_MAX) ) then
+ ! zero the CC measurments
+ tshift = 0.0
+ dlnA = 0.0
+ if (DISPLAY_DETAILS) then
+ print *, 'Fail if ANY of these is true :'
+ print *, ' cc_max : ', cc_max, CC_MIN, cc_max < CC_MIN
+ print *, ' tshift : ', tshift, TSHIFT_MIN, tshift < TSHIFT_MIN
+ print *, ' tshift : ', tshift, TSHIFT_MAX, tshift > TSHIFT_MAX
+ print *, ' dlnA : ', dlnA, DLNA_MIN, dlnA < DLNA_MIN
+ print *, ' dlnA : ', dlnA, DLNA_MAX, dlnA > DLNA_MAX
+ endif
+ endif
+
+ end subroutine cc_measure_select
+
+ !-----------------------------------------------------------------------------
+
subroutine mt_measure_select(nlen,tshift,i_pmax_syn,dtau_w,dlnA_w,err_dt,err_dlnA, &
dt,i_left,i_right,fstart,fend,use_trace)
@@ -1247,7 +1294,7 @@
!if(cc > cc_max) then
! CHT, 07-Sept-2008: Do not allow time shifts larger than the specified input
- if(cc .gt. cc_max .and. abs(i*dt) <= BEFORE_TSHIFT ) then
+ if( (cc .gt. cc_max) .and. (abs(i*dt) <= TSHIFT_MAX) ) then
cc_max = cc
ishift = i
endif
@@ -1314,7 +1361,7 @@
endif
! make final adjustments to uncertainty estimate
- if (INCLUDE_ERROR == 0) then
+ if (ERROR_TYPE == 0) then
! set uncertainty factors to 1 if you do not want to incorporate them
! into the adjoint sources and the misfit function values
sigma_dt = 1.0
@@ -1339,14 +1386,13 @@
integer, intent(in) :: iker
character(len=40) :: stlab, suffix
- if (iker == 1) then
+ if ( iker == 7 .or. iker == 8 ) then
stlab = 'Multitaper' ; suffix = 'average'
- elseif (iker == 2) then
+ else
stlab = 'Cross-correlation' ; suffix = 'cc'
endif
- if (iker == 1 .or. iker == 2) then
-
+ if ( iker .ge. 3 ) then
if (DISPLAY_DETAILS) then
print *, trim(stlab)//' average measurements:'
print *, ' traveltime :', sngl(dtau_meas), ' +/- ', sngl(dtau_sigma)
@@ -1362,7 +1408,6 @@
write(72,*) dlnA_meas, dlnA_sigma
close(72)
endif
-
endif
end subroutine write_average_meas
@@ -1678,6 +1723,106 @@
!-------------------------------------------------------------------
+ subroutine read_par_file(fstart0,fend0,tt,dtt,nn,chan)
+
+ double precision, intent(out) :: fstart0,fend0,tt,dtt
+ integer, intent(out) :: nn
+ character(len=10), intent(out) :: chan
+ integer :: ios
+
+ ! input file MEASUREMENT.PAR -- see write_par_file.pl for details
+
+ OUT_DIR = 'OUTPUT_FILES' ! default
+
+ open(10,file='MEASUREMENT.PAR',status='old',iostat=ios)
+ read(10,*) tt,dtt,nn
+ read(10,*) iker0
+ read(10,*) chan
+ read(10,*) TLONG, TSHORT
+ read(10,*) RUN_BANDPASS
+ read(10,*) DISPLAY_DETAILS
+ read(10,*) OUTPUT_MEASUREMENT_FILES
+ read(10,*) TSHIFT_MIN, TSHIFT_MAX
+ read(10,*) DLNA_MIN, DLNA_MAX
+ read(10,*) CC_MIN
+ read(10,*) ERROR_TYPE
+ read(10,*) DT_SIGMA_MIN
+ read(10,*) DLNA_SIGMA_MIN
+ read(10,*) itaper
+ read(10,*) wtr,npi
+ read(10,*) DT_FAC
+ read(10,*) ERR_FAC
+ read(10,*) DT_MAX_SCALE
+ read(10,*) NCYCLE_IN_WINDOW
+ close(10)
+
+ iker = iker0
+
+ ! check the read-in values
+ print *, 'INPUTS FROM MEASUREMENT.PAR :'
+ print *, ' tt, dtt, nn : ',tt,dtt,nn
+ print *, ' chan :',chan
+ print *, ' iker : ',iker
+ print *, ' TLONG, TSHORT : ',TLONG, TSHORT
+ fstart0 = 1./TLONG ; fend0 = 1./TSHORT
+ print *, ' fstart, fend :', fstart0, fend0
+ print *, ' RUN_BANDPASS :',RUN_BANDPASS
+ print *, ' DISPLAY_DETAILS :',DISPLAY_DETAILS
+ print *, ' OUTPUT_MEASUREMENT_FILES :',OUTPUT_MEASUREMENT_FILES
+ print *, ' TSHIFT_MIN, TSHIFT_MAX :',TSHIFT_MIN, TSHIFT_MAX
+ print *, ' DLNA_MIN, DLNA_MAX :',DLNA_MIN, DLNA_MAX
+ print *, ' CC_MIN :'
+ print *, ' ERROR_TYPE :',ERROR_TYPE
+ print *, ' DT_SIGMA_MIN :',DT_SIGMA_MIN
+ print *, ' DLNA_SIGMA_MIN :',DLNA_SIGMA_MIN
+ print *, ' itaper : ',itaper
+ print *, ' wtr, npi : ',wtr,npi
+ print *, ' DT_FAC :',DT_FAC
+ print *, ' ERR_FAC :',ERR_FAC
+ print *, ' DT_MAX_SCALE :',DT_MAX_SCALE
+ print *, ' NCYCLE_IN_WINDOW :',NCYCLE_IN_WINDOW
+ !stop 'checking PAR file input'
+
+ ! assign additional parameters and stop for certain inconsistencies
+ if (fstart0.ge.fend0) then
+ print *, 'Check input frequency range of the signal'
+ stop
+ endif
+
+ if (nn > NDIM) then
+ print *, 'Error: Change interpolation nn or NDIM'
+ stop
+ endif
+
+ ! for CC kernels, itaper must be a single taper (2 or 3)
+ if ( (itaper==1) .and. ((iker.ge.3).and.(iker.le.6)) ) then
+ print *, 'Error: Change itaper to 2 or 3'
+ stop
+ endif
+
+ COMPUTE_ADJOINT_SOURCE = .true.
+ if (iker==0) COMPUTE_ADJOINT_SOURCE = .false.
+
+ if (iker.le.2) then
+ is_mtm0 = 0
+ elseif ( (iker.ge.3) .and. (iker.le.6) ) then
+ is_mtm0 = itaper ! 2 or 3
+ elseif (iker==7 .or. iker==8) then
+ is_mtm0 = 1 ! MT taper for MT adjoint source
+ else
+ print *, 'Error: iker must by 0-8'
+ stop
+ endif
+
+ is_mtm = is_mtm0
+
+ print *, ' COMPUTE_ADJOINT_SOURCE :',COMPUTE_ADJOINT_SOURCE
+ print *, ' is_mtm :',is_mtm
+
+ end subroutine read_par_file
+
+!-------------------------------------------------------------------
+
subroutine get_sacfile_header(data_file,yr,jda,ho,mi,sec,ntw,sta,comp,dist,az,baz,slat,slon)
character(len=*),intent(in) :: data_file
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_variables.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_variables.f90 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_variables.f90 2009-10-20 21:54:52 UTC (rev 15852)
@@ -12,11 +12,17 @@
!
! See write_par_file.pl and mt_measure_adj.f90
+ character(len=150) :: OUT_DIR
+
double precision :: TLONG, TSHORT
double precision :: WTR, NPI, DT_FAC, ERR_FAC, DT_MAX_SCALE, NCYCLE_IN_WINDOW
- double precision :: BEFORE_QUALITY, AFTER_QUALITY, BEFORE_TSHIFT, AFTER_TSHIFT
+ !double precision :: BEFORE_QUALITY, AFTER_QUALITY, BEFORE_TSHIFT, AFTER_TSHIFT
+ double precision :: TSHIFT_MIN, TSHIFT_MAX, DLNA_MIN, DLNA_MAX, CC_MIN
double precision :: DT_SIGMA_MIN, DLNA_SIGMA_MIN
- integer :: ntaper, ipwr_t, ipwr_w, INCLUDE_ERROR
- logical :: DISPLAY_DETAILS, OUTPUT_MEASUREMENT_FILES, RUN_BANDPASS
+ integer :: ntaper, ipwr_t, ipwr_w, ERROR_TYPE
+ integer :: iker0, iker, itaper, is_mtm0, is_mtm
+
+ logical :: DISPLAY_DETAILS,OUTPUT_MEASUREMENT_FILES,RUN_BANDPASS,COMPUTE_ADJOINT_SOURCE
+
end module mt_variables
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/prepare_adj_src.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/prepare_adj_src.pl 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/prepare_adj_src.pl 2009-10-20 21:54:52 UTC (rev 15852)
@@ -10,7 +10,7 @@
# a particular component (say, Z), but IS made on another component (say, R or T).
#
# EXAMPLE:
-# prepare_adj_src.pl -m PLOTS/CMTSOLUTION_9818433 -s STATIONS_CI -o ADJOINT_SOURCES OUTPUT_FILES/*adj
+# prepare_adj_src.pl -m PLOTS/CMTSOLUTION_9818433 -z BH -s STATIONS_CI -o ADJOINT_SOURCES OUTPUT_FILES/*adj
#
#--------------------------------------------------------------------
@@ -21,11 +21,12 @@
use DELAZ5;
#Reading input arguments:
- at ARGV>0 || die("prepare_adj_src.pl -m CMT -s STATION -o OUTDIR all_adj_files\n");
+ at ARGV>0 || die("prepare_adj_src.pl -m CMT -z BH -s STATION -o OUTDIR all_adj_files\n");
-if (!getopts('m:s:o:')) {die("Check arguments\n");}
+if (!getopts('m:z:s:o:')) {die("Check arguments\n");}
if ($opt_m) {$cmt = $opt_m;} else {$cmt = "CMTSOLUTION";}
+if ($opt_z) {$chan = $opt_z;} else {$chan = "BH";}
if ($opt_s) {$stafile = $opt_s;} else {$stafile = "STATIONS";}
if ($opt_o) {$outdir = $opt_o;} else {$outdir = "ADJ_OUT";}
@@ -48,28 +49,32 @@
#die("testing");
-# NOTE: Based on the output from mt_measure_adj.f90, the channels are
-# all of the form BH_, no matter what the channel is for the data records.
+$cmpZ = "${chan}Z";
+$cmpR = "${chan}R";
+$cmpT = "${chan}T";
+$cmpE = "${chan}E";
+$cmpN = "${chan}N";
+# NOTE: channel is either BH or LH (from synthetics)
$stafile_temp = "station.tmp";
system("rm -f ${stafile_temp}");
$nstation = 0;
print "Copying stations adjoint source file: $outdir ...\n";
foreach $stanet (keys(%adj)) {
($sta,$net) = split(/\./,$stanet);
- if (defined $adj{$stanet}{BHT} or defined $adj{$stanet}{BHR} or defined $adj{$stanet}{BHZ}) {
- if (defined $adj{$stanet}{BHT} ) {
- $tcomp = $adj{$stanet}{BHT}; $rcomp = $tcomp; $zcomp = $tcomp;
- $ncomp = $rcomp; $ecomp = $tcomp; $zcomp=~s/BHT/BHZ/;
- $rcomp =~s/BHT/BHR/; $ncomp=~s/BHT/BHN/; $ecomp=~s/BHT/BHE/;}
- elsif (defined $adj{$stanet}{BHR}) {
- $rcomp = $adj{$stanet}{BHR}; $tcomp = $rcomp; $zcomp = $rcomp;
- $ncomp = $rcomp; $ecomp = $rcomp; $zcomp=~s/BHR/BHZ/;
- $tcomp =~s/BHR/BHT/; $ncomp=~s/BHR/BHN/; $ecomp=~s/BHR/BHE/;}
+ if (defined $adj{$stanet}{$cmpT} or defined $adj{$stanet}{$cmpR} or defined $adj{$stanet}{$cmpZ}) {
+ if (defined $adj{$stanet}{$cmpT} ) {
+ $tcomp = $adj{$stanet}{$cmpT}; $rcomp = $tcomp; $zcomp = $tcomp;
+ $ncomp = $rcomp; $ecomp = $tcomp; $zcomp=~s/$cmpT/$cmpZ/;
+ $rcomp =~s/$cmpT/$cmpR/; $ncomp=~s/$cmpT/$cmpN/; $ecomp=~s/$cmpT/$cmpE/;}
+ elsif (defined $adj{$stanet}{$cmpR}) {
+ $rcomp = $adj{$stanet}{$cmpR}; $tcomp = $rcomp; $zcomp = $rcomp;
+ $ncomp = $rcomp; $ecomp = $rcomp; $zcomp=~s/$cmpR/$cmpZ/;
+ $tcomp =~s/$cmpR/$cmpT/; $ncomp=~s/$cmpR/$cmpN/; $ecomp=~s/$cmpR/$cmpE/;}
else {
- $zcomp = $adj{$stanet}{BHZ}; $tcomp = $zcomp; $rcomp = $zcomp;
- $ncomp = $rcomp; $ecomp = $rcomp; $rcomp=~s/BHZ/BHR/;
- $tcomp =~s/BHZ/BHT/; $ncomp=~s/BHZ/BHN/; $ecomp=~s/BHZ/BHE/;}
+ $zcomp = $adj{$stanet}{$cmpZ}; $tcomp = $zcomp; $rcomp = $zcomp;
+ $ncomp = $rcomp; $ecomp = $rcomp; $rcomp=~s/$cmpZ/$cmpR/;
+ $tcomp =~s/$cmpZ/$cmpT/; $ncomp=~s/$cmpZ/$cmpN/; $ecomp=~s/$cmpZ/$cmpE/;}
# get station latitude and longitude from stations file
# NOTE: There could be stations in the list with the same name
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/run_mt_measure_adj.csh
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/run_mt_measure_adj.csh 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/run_mt_measure_adj.csh 2009-10-20 21:54:52 UTC (rev 15852)
@@ -1,10 +1,7 @@
\rm -rf PLOTS/ADJOINT_SOURCES
-mkdir PLOTS/ADJOINT_SOURCES
+ mkdir PLOTS/ADJOINT_SOURCES
cp OUTPUT_FILES/*adj PLOTS/ADJOINT_SOURCES
cp OUTPUT_FILES/*recon.cc* PLOTS/RECON
-cp DATA_TEST/* PLOTS/DATA_TEST/
-cp SYN_TEST/* PLOTS/SYN_TEST/
-cp MEASUREMENT.WINDOWS PLOTS/
cp window_chi PLOTS
\rm -rf ADJOINT_SOURCES
mkdir ADJOINT_SOURCES
@@ -12,5 +9,5 @@
cp STATIONS_ADJOINT ADJOINT_SOURCES
\mv STATIONS_ADJOINT PLOTS
cd PLOTS
-plot_win_adj_all.pl -l -10/200 -m ../CMTSOLUTION_9818433 -b 0 -k 1 -a STATIONS_ADJOINT -d DATA_TEST -s SYN_TEST -c RECON -w MEASUREMENT.WINDOWS -i m16 -j 6/30
-cd -
\ No newline at end of file
+plot_win_adj_all.pl -l -10/200 -m ../CMTSOLUTION_9818433 -n BH -b 0 -k 7 -a STATIONS_ADJOINT -d DATA -s SYN -c RECON -w MEASUREMENT.WINDOWS -i m16 -j 6/30
+cd /data1/cig/seismo/3D/ADJOINT_TOMO/measure_adj_work
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 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/prepare_mt_measure_adj.pl 2009-10-20 21:54:52 UTC (rev 15852)
@@ -4,7 +4,7 @@
#
# prepare_mt_measure_adj.pl
# Carl Tape
-# 01-Oct-2008
+# 19-Oct-2009
#
# This script prepares the data and synthetics for running in the measurement code.
# At the moment, it is somewhat tailored to the SoCal dataset.
@@ -16,20 +16,11 @@
# Tmin/Tmax bandpass periods for filtering data and synthetics
#
# EXAMPLE:
-# prepare_mt_measure_adj.pl m03 0 1 0 6/40 9818433
+# prepare_mt_measure_adj.pl m16 0 1 0 6/30 9818433
#
-# prepare_mt_measure_adj.pl m00 1 1 1 6/40 9818433 (MEASURE_TEST -- window dir)
-# prepare_mt_measure_adj.pl m00 0 1 1 6/40 9818433 (MEASURE -- window dir)
-# prepare_mt_measure_adj.pl m00 0 1 0 6/40 9818433 (MEASURE)
-#
-# prepare_mt_measure_adj.pl m00 0 1 1 6/40 14095540 (MEASURE -- window dir)
-#
-# prepare_mt_measure_adj.pl m00 0 1 0 2/40 9818433
-# prepare_mt_measure_adj.pl m00 1 1 1 2/40 9818433
-#
#==========================================================
-if (@ARGV < 6) {die("Usage: prepare_mt_measure_adj.pl itest iplot iwindow Tmin/Tmax eid\n");}
+if (@ARGV < 6) {die("Usage: prepare_mt_measure_adj.pl smodel itest iplot iwindow Tmin/Tmax eid\n");}
($smodel,$itest,$iplot,$iwindow,$Ts,$eid) = @ARGV;
# iwindow: controls where the data, synthetics, and windows files are
@@ -50,6 +41,14 @@
if($iwindow==0) {die("Exit: must have iwindow = 1 for TEST option.");}
}
+#============================================
+# USER INPUT
+
+# CMTSOLUTION file
+$dir_cmt = "/home/carltape/results/SOURCES/socal_16/v16_files";
+$cmtfile = "${dir_cmt}/CMTSOLUTION_${eid}";
+if (not -f $cmtfile) {die("check if cmtfile $cmtfile exist or not\n")}
+
# directories for data, synthetics, and window files : MUST BE CHANGED FOR EACH USER
if ($iwindow==1) {
# windowing code run directory
@@ -75,6 +74,8 @@
#$dir_win_out = "/net/sierra/raid1/carltape/results/WINDOWS/model_m0/";
}
+#============================================
+
if (not -e ${dir_win_run_syn}) {die("check if ${dir_win_run_syn} exist or not\n");}
if (not -e ${dir_win_run_data}) {die("check if ${dir_win_run_data} exist or not\n");}
if (not -e ${win_in}) {die("check if ${win_in} exist or not\n");}
@@ -97,7 +98,7 @@
# copy ALL data and synthetics into measurement directory (BHZ,BHR,BHT)
# NOTE: copy the whole directory, because the list might be too long for using cp
-# It might be simpler to LINK the data and syn directories to the local run directory.
+# Another option is to LINK the data and syn directories to the local run directory.
print CSH "echo copying files into the measurement directory...\n";
print CSH "cp -r ${dir_win_run_syn} ${dir_meas_syn}\n sleep 2s\n";
print CSH "cp -r ${dir_win_run_data} ${dir_meas_data}\n sleep 2s\n";
@@ -105,9 +106,6 @@
#print CSH "cp ${dir_win_run_data}/* ${dir_meas_data}\n";
# CMTSOLUTION file
-$dir_cmt = "/home/carltape/results/SOURCES/socal_16/v16_files";
-$cmtfile = "${dir_cmt}/CMTSOLUTION_${eid}";
-if (not -f $cmtfile) {die("check if cmtfile $cmtfile exist or not\n")}
print CSH "\\rm CMTSOLUTION*\n";
print CSH "echo $cmtfile\n";
print CSH "cp $cmtfile .\n";
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl 2009-10-20 21:54:52 UTC (rev 15852)
@@ -4,7 +4,7 @@
#
# run_mt_cc_plot.pl
# Carl Tape
-# 27-July-2009
+# 20-Oct-2009
#
# This runs mt_measure_adj.f90 first in cross-correlation mode,
# and then in multitaper mode. The purpose is to have
@@ -12,12 +12,13 @@
# We do not need both for the adjoint tomography inversion.
#
# EXAMPLE:
-# run_mt_cc_plot.pl m16 0 6/30 -10/200 -0.585/0.011/18200 0 1/1/1 2.0/2.5/3.5/1.5 0.7/0.7/5.0/1.0 1.0/0.5 (DATA/SYN)
+# run_mt_cc_plot.pl m16 -10/200 0 -0.585/0.011/18200 7 BH 6/30 0/1/1 -5.0/5.0/-1.5/1.5/0.7 0/1.0/0.5 1/0.02/2.5/2.0/2.5/3.5/1.5
#
#==========================================================
-if (@ARGV < 10) {die("Usage: run_mt_cc_plot.pl smodel Tmin/Tmax tstart/dt/npt Ldur itest iparbools par1 par2 par3\n")}
-($smodel,$ibp,$Ts,$lcut,$tvec,$itest,$iparbools,$par1,$par2,$par3) = @ARGV;
+if (@ARGV < 10) {die("Usage: run_mt_cc_plot.pl xxx\n")}
+($smodel,$lcut,$itest,$tvec,$iker,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
+#($ibp,$Ts,$lcut,$tvec,$itest,$iparbools,$par1,$par2,$par3) = @ARGV;
$pwd = $ENV{PWD};
@@ -39,10 +40,14 @@
# delete all figures in PLOTS
`rm PLOTS/*ps PLOTS/*jpg PLOTS/*pdf`;
-# first run in cross-correlation mode
-$iker = 2;
+# first run in cross-correlation TT mode
+$iker = 5;
$itaper = 3;
+# replace itaper in par3 with this value
+(undef,$v2,$v3,$v4,$v5,$v6,$v7) = split("/",$par3);
+$par3 = "$itaper/$v2/$v3/$v4/$v5/$v6/$v7";
+
$plot_adj_dir = "PLOTS/ADJOINT_SOURCES";
#$cmtfile = "CMTSOLUTION_9818433";
@@ -66,8 +71,9 @@
#print CSH "\\rm -rf OUTPUT_FILES ; mkdir OUTPUT_FILES ; make clean\n make\n";
# create MEASUREMENT.PAR file
-$wtr = 0.02; $npi = 2.5; # "default" values
-print CSH "write_par_file.pl OUTPUT_FILES $ibp $Ts $tvec $itaper $iker $wtr/$npi $iparbools $par1 $par2 $par3\n";
+#$wtr = 0.02; $npi = 2.5; # "default" values
+#print CSH "write_par_file.pl OUTPUT_FILES $ibp $Ts $tvec $itaper $iker $wtr/$npi $iparbools $par1 $par2 $par3\n";
+print CSH "write_par_file.pl $tvec $iker $chan $Ts $iparbools $par1 $par2 $par3\n";
#---------------------------------------------
@@ -83,10 +89,14 @@
#////////////////////////////////////////////////////////
-# second run in multitaper mode
-$iker = 1;
+# second run in multitaper TT mode
+$iker = 7;
$itaper = 1;
+# replace itaper in par3 with this value
+(undef,$v2,$v3,$v4,$v5,$v6,$v7) = split("/",$par3);
+$par3 = "$itaper/$v2/$v3/$v4/$v5/$v6/$v7";
+
# empty the output files directory
print CSH "\\rm -rf OUTPUT_FILES ; mkdir OUTPUT_FILES\n";
@@ -94,7 +104,7 @@
print CSH "cp window_chi PLOTS\n";
# create MEASUREMENT.PAR file
-print CSH "write_par_file.pl OUTPUT_FILES $ibp $Ts $tvec $itaper $iker $wtr/$npi $iparbools $par1 $par2 $par3\n";
+print CSH "write_par_file.pl $tvec $iker $chan $Ts $iparbools $par1 $par2 $par3\n";
# run the measurement code and save output file
# (re-compilation is not necessary)
@@ -105,12 +115,12 @@
# prepare_adj_src.pl dumps the ZEN adjoint sources into $adj_dir
$adj_dir = "ADJOINT_SOURCES";
print CSH "\\rm -rf ${adj_dir}\n mkdir ${adj_dir}\n";
-print CSH "prepare_adj_src.pl -m $cmtfile -s PLOTS/$stafile1 -o ${adj_dir} OUTPUT_FILES/*adj\n";
+print CSH "prepare_adj_src.pl -m $cmtfile -z $chan -s PLOTS/$stafile1 -o ${adj_dir} OUTPUT_FILES/*adj\n";
print CSH "cp STATIONS_ADJOINT ${adj_dir}\n";
print CSH "\\mv STATIONS_ADJOINT PLOTS\n";
print CSH "cd PLOTS\n";
-print CSH "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -b 1 -k 1 -a $stafile2 -d $dir_data -s $dir_syn -c $dir_recon -w MEASUREMENT.WINDOWS -i $smodel -j $Ts\n";
+print CSH "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -n $chan -b 1 -k $iker -a $stafile2 -d $dir_data -s $dir_syn -c $dir_recon -w MEASUREMENT.WINDOWS -i $smodel -j $Ts\n";
print CSH "\\rm test*.pdf\n";
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl 2009-10-20 21:54:52 UTC (rev 15852)
@@ -4,43 +4,29 @@
#
# run_mt_measure_adj.pl
# Carl Tape
-# 11-Oct-2007
+# 19-Oct-2009
#
-# This script runs after prepare_mt_measure.pl.
+# This script runs following prepare_mt_measure.pl.
# It executes mt_measure_adj and then ouputs a complete set of
# un-rotated adjoint sources for SPECFEM3D as well as STATIONS_ADJOINT.
# It also provides the option of making a set of plots by
# calling plot_win_adj_all.pl.
#
# INPUT:
-# iker type of measurement (and adjoint source)
-# Tmin/Tmax bandpass periods for filtering data and synthetics
-# tstart/dt/npt time vector -- must match SEM time vector
-# -9.00/0.011/18200
-# -0.0825/0.011/18200
+# smodel index label for model
+# lcut cut times for plotting records
# itest DATA_TEST/SYN_TEST or DATA/SYN
# iplot make plots of data, syn, windows, and adjoint sources
# iboth plot both MT and CC adjoint sources (=1) or not (=0)
+# ... see write_par_file.pl
#
-# EXAMPLE (T = [6s, 40s]):
-# run_mt_measure_adj.pl m00 0 0 0 6/40 -10/180 -0.6/0.011/18200 1 1 1 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.1 (DATA_TEST/SYN_TEST)
-# run_mt_measure_adj.pl m00 2 3 0 6/40 -10/180 -0.6/0.011/18200 1 1 1 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.1 (DATA_TEST/SYN_TEST)
-# run_mt_measure_adj.pl m00 1 1 0 6/40 -10/180 -0.6/0.011/18200 1 1 1 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.1 (DATA_TEST/SYN_TEST)
+# EXAMPLE (T = [6s, 30s]):
+# run_mt_measure_adj.pl m16 -10/200 0 1 0 -0.585/0.011/18200 7 BH 6/30 0/1/1 -5.0/5.0/-1.5/1.5/0.7 0/1.0/0.5 1/0.02/2.5/2.0/2.5/3.5/1.5
#
-# run_mt_measure_adj.pl m00 2 3 0 6/40 -10/180 -0.6/0.011/18200 0 1 0 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.1 (DATA/SYN, TT-xcorr, errors)
-# run_mt_measure_adj.pl m00 2 3 0 6/40 -10/180 -0.6/0.011/18200 0 1 0 1/1/0 2.0/2.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.1 (DATA/SYN, TT-xcorr, no errors)
-#
-# No errors included:
-# run_mt_measure_adj.pl m00 2 3 0 6/40 -10/180 -0.6/0.011/18200 1 1 1 1/1/0 2.0/2.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.1 (DATA_TEST/SYN_TEST)
-#
-# EXAMPLE (T = [2s, 40s]):
-# run_mt_measure_adj.pl m00 2 3 0 2/40 -10/120 -0.6/0.011/18200 1 1 0 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.1 (DATA_TEST/SYN_TEST)
-# run_mt_measure_adj.pl m00 2 3 0 2/40 -10/120 -0.6/0.011/18200 0 1 0 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.1 (DATA/SYN)
-#
#==========================================================
-if (@ARGV < 14) {die("Usage: run_mt_measure_adj.pl smodel iker itaper ibp Tmin/Tmax tstart/dt/npt itest iplot iboth iparbools par1 par2 par3\n")}
-($smodel,$iker,$itaper,$ibp,$Ts,$lcut,$tvec,$itest,$iplot,$iboth,$iparbools,$par1,$par2,$par3) = @ARGV;
+if (@ARGV < 13) {die("Usage: run_mt_measure_adj.pl xxx\n")}
+($smodel,$lcut,$itest,$iplot,$iboth,$tvec,$iker,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
$pwd = $ENV{PWD};
@@ -77,8 +63,8 @@
#print CSH "\\rm -rf OUTPUT_FILES ; mkdir OUTPUT_FILES ; make clean\n make\n";
#----------------------------
-# choose ipar2 based on the values used in the windowing code PAR_FILE (Feb 2009)
-# This is a more sensible option, since you guarantee consistency.
+# choose ipar1 based on the values used in the windowing code PAR_FILE (Feb 2009)
+# This is a more sensible option, since you guarantee consistency with FLEXWIN.
if(1==1) {
# Allow for some "loosening" of the parameters, since the exact nature
# of the measurement in the measurement code may exclude some input windows.
@@ -90,25 +76,27 @@
(undef,undef,$dlna_base) = split(" ",`grep DLNA_BASE PAR_FILE`);
(undef,undef,$dlna_ref) = split(" ",`grep DLNA_REFERENCE PAR_FILE`);
(undef,undef,$cc_base) = split(" ",`grep CC_BASE PAR_FILE`);
- print "Values from PAR_FILE:\n";
- print "-- $tshift_base -- $tshift_ref -- $dlna_base -- $dlna_ref -- $cc_base --\n";
- $BEFORE_QUALITY = $cc_base - $cc_fac;
- $AFTER_QUALITY = $cc_base - $cc_fac;
- $BEFORE_TSHIFT = abs($tshift_ref) + $tshift_base + $dt_fac;
- $BEFORE_DLNA = abs($dlna_ref) + $dlna_base;
+
+ $TSHIFT_MIN = $tshift_ref - $tshift_base - $dt_fac;
+ $TSHIFT_MAX = $tshift_ref + $tshift_base + $dt_fac;
+ $DLNA_MIN = $dlna_ref - $dlna_base;
+ $DLNA_MAX = $dlna_ref + $dlna_base;
+ $CC_MIN = $cc_base - $cc_fac;
+ if($CC_MIN < 0) {$CC_MIN = 0}
- $par2_old = $par2;
- ($BEFORE_QUALITY_OLD,$AFTER_QUALITY_OLD,$BEFORE_TSHIFT_OLD,$AFTER_TSHIFT_OLD) = split("/",$par2);
- $AFTER_TSHIFT = $AFTER_TSHIFT_OLD;
- $par2 = "${BEFORE_QUALITY}/${AFTER_QUALITY}/${BEFORE_TSHIFT}/${AFTER_TSHIFT}";
- print "Updating par2:\n old -- $par2_old --\n new -- $par2 --\n";
+ $par1_old = $par1;
+ ($TSHIFT_MIN_OLD,$TSHIFT_MAX_OLD,$DLNA_MIN_OLD,$DLNA_MAX_OLD,$CC_MIN_OLD) = split("/",$par1);
+ $par1 = "${TSHIFT_MIN}/${TSHIFT_MAX}/${DLNA_MIN}/${DLNA_MAX}/${CC_MIN}";
+ print "Updating values from PAR_FILE:\n";
+ print " old -- $par1_old --\n new -- $par1 --\n";
}
#----------------------------
# create MEASUREMENT.PAR file
-$wtr = 0.02; $npi = 2.5; # "default" values
+#$wtr = 0.02; $npi = 2.5; # "default" values
#$itaper = 3;
-print CSH "write_par_file.pl OUTPUT_FILES $ibp $Ts $tvec $itaper $iker $wtr/$npi $iparbools $par1 $par2 $par3\n";
+#print CSH "write_par_file.pl OUTPUT_FILES $ibp $Ts $tvec $itaper $iker $wtr/$npi $iparbools $par1 $par2 $par3\n";
+print CSH "write_par_file.pl $tvec $iker $chan $Ts $iparbools $par1 $par2 $par3\n";
#---------------------------------------------
@@ -145,7 +133,7 @@
#print CSH "\\rm PLOTS/*pdf PLOTS/*jpg PLOTS/*ps\n";
print CSH "cd PLOTS\n";
- print CSH "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -b $iboth -k $iker -a $stafile2 -d $dir_data -s $dir_syn -c $dir_recon -w MEASUREMENT.WINDOWS -i $smodel -j $Ts\n";
+ print CSH "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -n $chan -b $iboth -k $iker -a $stafile2 -d $dir_data -s $dir_syn -c $dir_recon -w MEASUREMENT.WINDOWS -i $smodel -j $Ts\n";
# # make a single PDF file
# $isort = 3;
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2009-10-20 21:54:52 UTC (rev 15852)
@@ -4,27 +4,19 @@
#
# run_tomo.pl
# Carl Tape
-# 27-March-2009
+# 20-Oct-2009
#
# This is the master program for running the measurement code on all the
# events for which the windowing code has been run.
#
# NOTE: The command run_mt_measure_adj.pl below contains additional input parameters.
#
-# EXAMPLE (from working measure_adj directory) :
-# run_tomo.pl m16 0 6/30 0.011 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/5.0/1.0 1.0/0.5 1.0
-# run_tomo.pl m16 0 3/30 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/3.0/1.0 1.0/0.5 1.0
-# run_tomo.pl m16 0 2/30 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/2.0/1.0 1.0/0.5 1.0
+# EXAMPLES (19-Oct-2009):
+# run_tomo.pl m16 0.011 1.0 0 BH 6/30 0/1/1 -5.0/5.0/-1.5/1.5/0.7 0/1.0/0.5 3/0.02/2.5/2.0/2.5/3.5/1.5
+# run_tomo.pl m16 0.011 1.0 7 BH 6/30 0/1/1 -5.0/5.0/-1.5/1.5/0.7 0/1.0/0.5 1/0.02/2.5/2.0/2.5/3.5/1.5
+# run_tomo.pl m16 0.011 1.0 7 BH 3/30 0/1/1 -3.0/3.0/-1.5/1.5/0.7 0/1.0/0.5 1/0.02/2.5/0.25/0.5/3.5/1.5
+# run_tomo.pl m16 0.011 1.0 7 BH 2/30 0/1/1 -2.0/2.0/-1.5/1.5/0.7 0/1.0/0.5 1/0.02/2.5/0.25/0.5/3.5/1.5
#
-# socal m00 - m04
-# 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
-#
-# Multitaper (m00, m01, m02, m03, m15, m16):
-# run_tomo.pl m16 0 6/30 0.011 1/1/1 2.0/2.5/3.5/1.5 0.7/0.7/5.0/1.0 1.0/0.5 1.0 # 0.8
-# run_tomo.pl m16 0 3/30 0.011 1/1/1 0.25/0.5/3.5/1.5 0.7/0.7/3.0/1.0 1.0/0.5 1.0 # 0.6
-# run_tomo.pl m16 0 2/30 0.011 1/1/1 0.25/0.5/3.5/1.5 0.7/0.7/2.0/1.0 1.0/0.5 1.0 # 0.4
-#
# socal m00 - m04, 0.2/0.2
# socal m05 - m08, 0.4/0.4
# socal m09, MAX: 4.0 pm 0.6 (2-30s), 4.0 pm 0.8 (3-30s), 5.0 pm 0.8 (6-30s)
@@ -34,8 +26,9 @@
#
#==========================================================
-if (@ARGV < 8) {die("Usage: run_tomo.pl smodel ibandpass Tmin/Tmax dt iparbools par1 par2 par3 lcut_frac\n")}
-($smodel,$ibp,$Ts,$dt,$iparbools,$par1,$par2,$par3,$lcut_frac) = @ARGV;
+if (@ARGV < 10) {die("Usage: run_tomo.pl xxx\n")}
+($smodel,$dt,$lcut_frac,$iker,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
+#($smodel,$ibp,$Ts,$dt,$iparbools,$par1,$par2,$par3,$lcut_frac) = @ARGV;
$pwd = $ENV{PWD};
@@ -64,7 +57,7 @@
if (not -f ${file_eid}) {die("check if ${file_eid} exist or not\n")}
open(IN,"${file_eid}"); @elines = <IN>; close(IN);
$nevent = @elines;
-print "\n $nevent events in the list \n";
+print "\n $nevent events in the list\n";
# full stations file (copy into PLOTS directory)
$file_stations = "/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem";
@@ -84,8 +77,8 @@
#die("testing");
$imin = 1; $imax = $nevent;
-#$imin = 100; $imax = $nevent;
-#$imin = 183; $imax = $imin; # 183 Yorba Linda
+#$imin = 1; $imax = 265;
+#$imin = 261; $imax = $imin; # 183 Yorba Linda
# EXAMPLE lines from EIDs_simulation_duration
# 2 9967025 200.000 18200 0.29 4.084
@@ -160,23 +153,15 @@
print CSH "prepare_mt_measure_adj.pl $smodel 0 1 0 $Ts $eid\n";
#------------------------------------
- # EXAMPLES for socal dataset
- # MT and CC, with both adjoint source plotted
- #print CSH "run_mt_cc_plot.pl $smodel $ibp $Ts $lcut $tvec 0 $iparbools $par1 $par2 $par3\n";
+ $itest = 0;
+ $iplot = 1;
+ $iboth = 0;
- # CC only
- #print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+ print CSH "run_mt_measure_adj.pl $smodel $lcut $itest $iplot $iboth $tvec $iker $chan $Ts $iparbools $par1 $par2 $par3\n";
- # multitaper only
- print CSH "run_mt_measure_adj.pl $smodel 1 1 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+ #print CSH "run_mt_cc_plot.pl $smodel $lcut $itest $tvec $iker $chan $Ts $iparbools $par1 $par2 $par3\n";
- # coverage kernels for m16 (CC banana-doughnut adjoint sources)
- #print CSH "run_mt_measure_adj.pl $smodel 3 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
-
- # waveform difference
- #print CSH "run_mt_measure_adj.pl $smodel 0 0 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
-
#------------------------------------
# copy run files into RUN directory
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/write_par_file.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/write_par_file.pl 2009-10-20 17:17:34 UTC (rev 15851)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/write_par_file.pl 2009-10-20 21:54:52 UTC (rev 15852)
@@ -4,66 +4,65 @@
#
# write_par_file.pl
# Carl Tape
-# 24-July-2009
+# 19-Oct-2009
#
# This script writes a parameters file to be used in the measurement code.
# See the reading-in of this file in mt_measure_adj.f90.
#
# INPUT:
-# odir
-# RUN_BANDPASS
-# TSHORT/TLONG
# tstart/dt/ntime
-# itaper
# iker
-# wtr/npi
-# DISPLAY_DETAILS/OUTPUT_MEASUREMENT_FILES/INCLUDE_ERROR
-# DT_FAC/ERR_FAC/DT_MAX_SCALE/NCYCLE_IN_WINDOW
-# BEFORE_QUALITY/AFTER_QUALITY/BEFORE_TSHIFT/AFTER_TSHIFT
-# DT_SIGMA_MIN/DLNA_SIGMA_MIN
+# channel (BH or LH)
+# TSHORT/TLONG
+# RUN_BANDPASS/DISPLAY_DETAILS/OUTPUT_MEASUREMENT_FILES
+# par1: TSHIFT_MIN/TSHIFT_MAX/DLNA_MIN/DLNA_MAX/CC_MIN
+# par2: ERROR_TYPE/DT_SIGMA_MIN/DLNA_SIGMA_MIN
+# par3: itaper/wtr/npi/DT_FAC/ERR_FAC/DT_MAX_SCALE/NCYCLE_IN_WINDOW
#
# EXAMPLE:
-# write_par_file.pl OUTPUT_FILES 0 6/30 -0.585/0.011/18200 3 2 0.02/2.5 1/1/1 2.0/2.5/3.5/3.0 0.69/0.69/4.5/1.0 1.0/0.5
+# write_par_file.pl -0.585/0.011/18200 1 BH 6/30 0/1/1 -4.5/4.5/-1.5/1.5/0.69 1/1.0/0.5 1/0.02/2.5/2.0/2.5/3.5/1.5
#
#
#==========================================================
-if (@ARGV < 11) {die("Usage: write_par_file.pl odir TSHORT/TLONG tstart/dt/ntime itaper iker ntaper/wtr/npi iparbools par1 par2 par3\n");}
-($odir,$ibp,$Ts,$tvec,$itaper,$iker,$mtvals,$iparbools,$par1,$par2,$par3) = @ARGV;
+if (@ARGV < 8) {die("Usage: write_par_file.pl xxx\n");}
+($tvec,$iker,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
# extract variables
+($tstart,$dt,$ntime) = split("/",$tvec);
($TSHORT,$TLONG) = split("/",$Ts);
-($tstart,$dt,$ntime) = split("/",$tvec);
-($wtr,$npi) = split("/",$mtvals);
-($ibool1,$ibool2,$INCLUDE_ERROR) = split("/",$iparbools);
-($DT_FAC,$ERR_FAC,$DT_MAX_SCALE,$NCYCLE_IN_WINDOW) = split("/",$par1);
-($BEFORE_QUALITY,$AFTER_QUALITY,$BEFORE_TSHIFT,$AFTER_TSHIFT) = split("/",$par2);
-($DT_SIGMA_MIN,$DLNA_SIGMA_MIN) = split("/",$par3);
+($ibool1,$ibool2,$ibool3) = split("/",$iparbools);
+if($ibool1==1) {$RUN_BANDPASS = ".true."} else {$RUN_BANDPASS = ".false."}
+if($ibool2==1) {$DISPLAY_DETAILS = ".true."} else {$DISPLAY_DETAILS = ".false."}
+if($ibool3==1) {$OUTPUT_MEASUREMENT_FILES = ".true."} else {$OUTPUT_MEASUREMENT_FILES = ".false."}
+($TSHIFT_MIN,$TSHIFT_MAX,$DLNA_MIN,$DLNA_MAX,$CC_MIN) = split("/",$par1);
+($ERROR_TYPE,$DT_SIGMA_MIN,$DLNA_SIGMA_MIN) = split("/",$par2);
+($itaper,$wtr,$npi,$DT_FAC,$ERR_FAC,$DT_MAX_SCALE,$NCYCLE_IN_WINDOW) = split("/",$par3);
-if($ibp==1) {$RUN_BANDPASS = ".true."} else {$RUN_BANDPASS = ".false."}
-if($ibool1==1) {$DISPLAY_DETAILS = ".true."} else {$DISPLAY_DETAILS = ".false."}
-if($ibool2==1) {$OUTPUT_MEASUREMENT_FILES = ".true."} else {$OUTPUT_MEASUREMENT_FILES = ".false."}
-#if($ibool3==1) {$INCLUDE_ERROR = ".true."} else {$INCLUDE_ERROR = ".false."}
-
# comments for the PAR_FILE
#$line00 = "output directory"; # cannot have a string here
-$line01 = "is_mtm -- taper type: 0 none; 1 multi-taper; 2 cosine; 3 boxcar";
-$line02 = "wtr, npi (ntaper = 2*npi)";
-$line03 = "iker -- 0 waveform; 1 multi-taper; 2 cross-correlation;";
-$line04 = "RUN_BANDPASS: use band-pass on records";
-$line05 = "TLONG and TSHORT: band-pass periods for records";
-$line06 = "tstart, DT, npts: time vector for simulations";
-$line07 = "DISPLAY_DETAILS";
-$line08 = "OUTPUT_MEASUREMENT_FILES";
-$line09 = "INCLUDE_ERROR -- 0 none; 1 CC, MT; 2 CC, MT-jack-knife";
-$line10 = "DT_FAC";
-$line11 = "ERR_FAC";
-$line12 = "DT_MAX_SCALE";
-$line13 = "NCYCLE_IN_WINDOW";
-$line14 = "BEFORE_QUALITY; AFTER_QUALITY";
-$line15 = "BEFORE_TSHIFT; AFTER_TSHIFT";
-$line16 = "DT_SIGMA_MIN; DLNA_SIGMA_MIN";
+$line01 = "tstart, DT, npts: time vector for simulations";
+$line02 = "iker (0-8) -- 0 no adj src; 1-8 see manual;";
+$line03 = "channel: BH or LH";
+$line04 = "TLONG and TSHORT: band-pass periods for records";
+$line05 = "RUN_BANDPASS: use band-pass on records";
+$line06 = "DISPLAY_DETAILS";
+$line07 = "OUTPUT_MEASUREMENT_FILES";
+$line08 = "TSHIFT_MIN; TSHIFT_MAX";
+$line09 = "DLNA_MIN; DLNA_MAX";
+$line10 = "CC_MIN";
+$line11 = "ERROR_TYPE -- 0 none; 1 CC, MT-CC; 2 MT-jack-knife";
+$line12 = "DT_SIGMA_MIN";
+$line13 = "DLNA_SIGMA_MIN";
+
+$line14 = "itaper -- taper type: 1 multi-taper; 2 cosine; 3 boxcar";
+$line15 = "wtr, npi (ntaper = 2*npi)";
+$line16 = "DT_FAC";
+$line17 = "ERR_FAC";
+$line18 = "DT_MAX_SCALE";
+$line19 = "NCYCLE_IN_WINDOW";
+
#=============================================
$ofile = "./MEASUREMENT.PAR";
`\\rm $ofile`; # remove file if it exists
@@ -71,23 +70,45 @@
print "\nWriting to $ofile ...\n";
open(OUT,">$ofile");
-print OUT sprintf("%s\n",$odir);
-print OUT sprintf("%23i # %s\n",$itaper,$line01);
-print OUT sprintf("%17.3f%6.2f # %s\n",$wtr,$npi,$line02);
-print OUT sprintf("%23i # %s\n",$iker,$line03);
-print OUT sprintf("%23s # %s\n",$RUN_BANDPASS,$line04);
-print OUT sprintf("%12.3f%11.3f # %s\n",$TLONG,$TSHORT,$line05);
-print OUT sprintf("%8.3f%7.4f%8i # %s\n",$tstart,$dt,$ntime,$line06);
-print OUT sprintf("%23s # %s\n",$DISPLAY_DETAILS,$line07);
-print OUT sprintf("%23s # %s\n",$OUTPUT_MEASUREMENT_FILES,$line08);
-print OUT sprintf("%23i # %s\n",$INCLUDE_ERROR,$line09);
-print OUT sprintf("%23.3f # %s\n",$DT_FAC,$line10);
-print OUT sprintf("%23.3f # %s\n",$ERR_FAC,$line11);
-print OUT sprintf("%23.3f # %s\n",$DT_MAX_SCALE,$line12);
-print OUT sprintf("%23.3f # %s\n",$NCYCLE_IN_WINDOW,$line13);
-print OUT sprintf("%12.4f%11.4f # %s\n",$BEFORE_QUALITY,$AFTER_QUALITY,$line14);
-print OUT sprintf("%12.4f%11.4f # %s\n",$BEFORE_TSHIFT,$AFTER_TSHIFT,$line15);
-print OUT sprintf("%12.4f%11.4f # %s\n",$DT_SIGMA_MIN,$DLNA_SIGMA_MIN,$line16);
+print OUT sprintf("%8.3f%7.4f%8i # %s\n",$tstart,$dt,$ntime,$line01);
+print OUT sprintf("%23i # %s\n",$iker,$line02);
+print OUT sprintf("%23s # %s\n",$chan,$line03);
+print OUT sprintf("%12.3f%11.3f # %s\n",$TLONG,$TSHORT,$line04);
+print OUT sprintf("%23s # %s\n",$RUN_BANDPASS,$line05);
+print OUT sprintf("%23s # %s\n",$DISPLAY_DETAILS,$line06);
+print OUT sprintf("%23s # %s\n",$OUTPUT_MEASUREMENT_FILES,$line07);
+
+print OUT sprintf("%12.4f%11.4f # %s\n",$TSHIFT_MIN,$TSHIFT_MAX,$line08);
+print OUT sprintf("%12.4f%11.4f # %s\n",$DLNA_MIN,$DLNA_MAX,$line09);
+print OUT sprintf("%23.3f # %s\n",$CC_MIN,$line10);
+print OUT sprintf("%23i # %s\n",$ERROR_TYPE,$line11);
+print OUT sprintf("%23.3f # %s\n",$DT_SIGMA_MIN,$line12);
+print OUT sprintf("%23.3f # %s\n",$DLNA_SIGMA_MIN,$line13);
+
+print OUT sprintf("%23i # %s\n",$itaper,$line14);
+print OUT sprintf("%17.3f%6.2f # %s\n",$wtr,$npi,$line15);
+print OUT sprintf("%23.3f # %s\n",$DT_FAC,$line16);
+print OUT sprintf("%23.3f # %s\n",$ERR_FAC,$line17);
+print OUT sprintf("%23.3f # %s\n",$DT_MAX_SCALE,$line18);
+print OUT sprintf("%23.3f # %s\n",$NCYCLE_IN_WINDOW,$line19);
+
+#print OUT sprintf("%s\n",$odir);
+#print OUT sprintf("%23i # %s\n",$itaper,$line01);
+#print OUT sprintf("%17.3f%6.2f # %s\n",$wtr,$npi,$line02);
+#print OUT sprintf("%23i # %s\n",$iker,$line03);
+#print OUT sprintf("%23s # %s\n",$RUN_BANDPASS,$line04);
+#print OUT sprintf("%12.3f%11.3f # %s\n",$TLONG,$TSHORT,$line05);
+#print OUT sprintf("%8.3f%7.4f%8i # %s\n",$tstart,$dt,$ntime,$line06);
+#print OUT sprintf("%23s # %s\n",$DISPLAY_DETAILS,$line07);
+#print OUT sprintf("%23s # %s\n",$OUTPUT_MEASUREMENT_FILES,$line08);
+#print OUT sprintf("%23i # %s\n",$INCLUDE_ERROR,$line09);
+#print OUT sprintf("%23.3f # %s\n",$DT_FAC,$line10);
+#print OUT sprintf("%23.3f # %s\n",$ERR_FAC,$line11);
+#print OUT sprintf("%23.3f # %s\n",$DT_MAX_SCALE,$line12);
+#print OUT sprintf("%23.3f # %s\n",$NCYCLE_IN_WINDOW,$line13);
+#print OUT sprintf("%12.4f%11.4f # %s\n",$BEFORE_QUALITY,$AFTER_QUALITY,$line14);
+#print OUT sprintf("%12.4f%11.4f # %s\n",$BEFORE_TSHIFT,$AFTER_TSHIFT,$line15);
+#print OUT sprintf("%12.4f%11.4f # %s\n",$DT_SIGMA_MIN,$DLNA_SIGMA_MIN,$line16);
close(OUT);
#=================================================================
More information about the CIG-COMMITS
mailing list