[cig-commits] r21081 - in seismo/3D/ADJOINT_TOMO/measure_adj: . PLOTS USER_MANUAL
liuqy at geodynamics.org
liuqy at geodynamics.org
Tue Nov 27 15:55:19 PST 2012
Author: liuqy
Date: 2012-11-27 15:55:19 -0800 (Tue, 27 Nov 2012)
New Revision: 21081
Modified:
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl
seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.pdf
seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.tex
seismo/3D/ADJOINT_TOMO/measure_adj/ma_sub.f90
seismo/3D/ADJOINT_TOMO/measure_adj/measure_adj.f90
seismo/3D/ADJOINT_TOMO/measure_adj/prepare_adj_src.pl
seismo/3D/ADJOINT_TOMO/measure_adj/rotate_adj_src.f90
seismo/3D/ADJOINT_TOMO/measure_adj/run_measure_adj.csh
Log:
STA.NT.CMP.recon.sac is now reconstructed from both CC and MTM
tweak scripts to run for different datasets.
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2012-11-27 20:16:42 UTC (rev 21080)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2012-11-27 23:55:19 UTC (rev 21081)
@@ -16,6 +16,7 @@
use Getopt::Std;
use POSIX;
+use List::Util qw[min max];
sub Usage{
print STDERR <<END;
@@ -206,7 +207,7 @@
#$Zsyn = "${opt_s}/${Zlab}.semd.sac${suffix}";
# RECONSTRUCTED SYNTHETICS
-$suffix = ".recon.cc.sac"; # USER MUST CHANGE
+$suffix = ".recon.sac"; # USER MUST CHANGE
$Trecon = "${opt_c}/${Tlab}${suffix}";
$Rrecon = "${opt_c}/${Rlab}${suffix}";
$Zrecon = "${opt_c}/${Zlab}${suffix}";
@@ -302,7 +303,7 @@
for ($i=0; $i < $nT; $i++) {
$nline=$nline+1;
$win=`sed -n ${nline}p $Winfile`; chomp($win);
- (undef,$winb,$wine)=split(/ +/,$win);
+ ($winb,$wine)=split(" ",$win);
push @Twinb, "$winb";
push @Twine, "$wine";
@@ -326,7 +327,7 @@
$win=`sed -n ${nline}p $Winfile`;
chomp($win);
# print "$win\n";
- (undef,$winb,$wine)=split(/ +/,$win);
+ ($winb,$wine)=split(" ",$win);
# print "$winb\n";
# print "$wine\n";
push @Rwinb, "$winb";
@@ -341,20 +342,20 @@
# windows and measurements: vertical component
$nzline=`grep -n \"$Zsyn\" $Winfile | awk -F: '{ print \$1 }'`;
chomp($nzline);
-#print "$nzline\n";
if ($nzline) {
$nline=$nzline;
$nline=$nline+1;
$nZ=`sed -n ${nline}p $Winfile`;
@mlines = `grep $Zlab ${meas_file}`; # measurements
+
for ($i=0; $i < $nZ; $i++) {
$nline=$nline+1;
$win=`sed -n ${nline}p $Winfile`;
chomp($win);
- # print "$win\n";
- (undef,$winb,$wine)=split(/ +/,$win);
- # print "$winb\n";
- # print "$wine\n";
+ # print "$win\n";
+ ($winb,$wine)=split(" ",$win);
+ # print "$winb\n";
+ # print "$wine\n";
push @Zwinb, "$winb";
push @Zwine, "$wine";
@@ -392,6 +393,10 @@
$maxD = max($maxT,max($maxR,$maxZ));
$maxS = max($maxT2,max($maxR2,$maxZ2));
$max = max($maxD,$maxS);
+if ($maxTDS == 0) {$maxTDS=$maxD;}
+if ($maxRDS == 0) {$maxRDS=$maxD;}
+if ($maxZDS == 0) {$maxZDS=$maxD;}
+
#print "\n-- $Tsyn -- $Tdat\n -- $Rsyn -- $Rdat\n -- $Zsyn -- $Zdat\n";
#print "\n Overall max values: $maxD (data), $maxS (syn) $max (overall)\n";
@@ -417,9 +422,18 @@
$sizeZ = $size*$maxZ/$maxZDS;
$sizeZ2 = $size*$maxZ2/$maxZDS;
}
+$maxsize=max($sizeT,$sizeR,$sizeZ);
+$maxsize2=max($sizeT2,$sizeR2,$sizeZ2);
+if ($sizeT == 0) {$sizeT = $maxsize;}
+if ($sizeR == 0) {$sizeR = $maxsize;}
+if ($sizeZ == 0) {$sizeZ = $maxsize;}
+if ($sizeT2 == 0) {$sizeT2 = $maxsize2;}
+if ($sizeR2 == 0) {$sizeR2 = $maxsize2;}
+if ($sizeZ2 == 0) {$sizeZ2 = $maxsize2;}
+
# limits for reconstructed records
-if ($ntline) {
+if ($Tflag && $ntline) {
(undef,$minT3,$maxT3) = split(" ",`$saclst depmin depmax f $Trecon`);
$maxT3 = max(abs($minT3),abs($maxT3));
if ($opt_r) {
@@ -428,7 +442,7 @@
$sizeT3 = $size*$maxT3/$maxTDS;
}
}
-if ($nrline) {
+if ($Rflag && $nrline) {
(undef,$minR3,$maxR3) = split(" ",`$saclst depmin depmax f $Rrecon`);
$maxR3 = max(abs($minR3),abs($maxR3));
if ($opt_r) {
@@ -437,7 +451,7 @@
$sizeR3 = $size*$maxR3/$maxRDS;
}
}
-if ($nzline) {
+if ($Zflag && $nzline) {
(undef,$minZ3,$maxZ3) = split(" ",`$saclst depmin depmax f $Zrecon`);
$maxZ3 = max(abs($minZ3),abs($maxZ3));
if ($opt_r) {
@@ -456,6 +470,7 @@
# TRANSVERSE component: data, synthetics, and windows
#print "pssac2 $Tsyn -X-2 -Y4.5 $proj $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O $tick \n";
+if (-f $Tsyn){
`pssac2 $Tsyn -X-2 -Y4.5 $proj $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O $tick >> $ps_file`; # synthetics
if ($ntline) {
if ($iplot_win==1) {
@@ -482,15 +497,17 @@
`echo \"$xtext $ytextpos3 7 0 1 $just $stlabT_dA[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
}
}
-}
+}}
if ($Tflag==1) {`pssac2 $Tdat $proj $bounds -Ent-3 -M${sizeT} $black_pen -N -K -O >> $ps_file`;} # data
-`pssac2 $Tsyn $proj $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O >> $ps_file`; # synthetics
-if ($ntline && ($iplot_recon==1)) {`pssac2 $Trecon $proj $bounds -Ent-3 -M${sizeT3} $recon_pen -N -K -O >> $ps_file`;} # reconstructed
+if (-f $Tsyn){
+`pssac2 $Tsyn $proj $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O >> $ps_file`;} # synthetics
+if ($Tflag && $ntline && ($iplot_recon==1)) {`pssac2 $Trecon $proj $bounds -Ent-3 -M${sizeT3} $recon_pen -N -K -O >> $ps_file`;} # reconstructed
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC T\" | pstext $proj $bounds -N -O -K >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strT\" | pstext $proj $bounds -N -O -K >> $ps_file`;
# RADIAL component: data, synthetics, and windows
-`pssac2 $Rsyn -Y$dY $proj $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O $tick2 >> $ps_file`;
+if (-f $Rsyn) {
+`pssac2 $Rsyn -Y$dY $proj $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O $tick2 >> $ps_file`;}
if ($nrline) {
if ($iplot_win==1) {
for ($i=0; $i<$nR; $i++) {
@@ -514,12 +531,14 @@
}
}
if ($Rflag==1) {`pssac2 $Rdat $proj $bounds -Ent-3 -M${sizeR} $black_pen -N -K -O >> $ps_file`;}
-`pssac2 $Rsyn $proj $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O >> $ps_file`;
-if ($nrline && ($iplot_recon==1)) {`pssac2 $Rrecon $proj $bounds -Ent-3 -M${sizeR3} $recon_pen -N -K -O >> $ps_file`;}
+if (-f $Rsyn){
+`pssac2 $Rsyn $proj $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O >> $ps_file`;}
+if ($Rflag && $nrline && ($iplot_recon==1)) {`pssac2 $Rrecon $proj $bounds -Ent-3 -M${sizeR3} $recon_pen -N -K -O >> $ps_file`;}
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC R\" | pstext $proj $bounds -N -O -K >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strR\" | pstext $proj $bounds -N -O -K >> $ps_file`;
# VERTICAL component: data, synthetics, and windows
+if (-f $Zsyn) {
`pssac2 $Zsyn -Y$dY $proj $bounds -Ent-3 -M${sizeZ2} $red_pen -N -K -O $tick2 >> $ps_file`;
if ($nzline) {
if ($iplot_win==1) {
@@ -527,6 +546,7 @@
open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
print GMT "$Zwinb[$i] 0\n $Zwine[$i] 0\n $Zwine[$i] 2\n $Zwinb[$i] 2\n";
close(GMT);
+ #print "*** $Zwinb[$i] $Zwine[$i]\n";
}
}
if ($iplot_CClabel==1) {
@@ -542,10 +562,11 @@
`echo \"$xtext $ytextpos3 7 0 1 $just $stlabZ_dA[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
}
}
-}
+}}
if ($Zflag==1) {`pssac2 $Zdat $proj $bounds -Ent-3 -M${sizeZ} $black_pen -N -K -O $tick2 >> $ps_file`;}
-`pssac2 $Zsyn $proj $bounds -Ent-3 -M${sizeZ2} $red_pen -N -O -K >> $ps_file`;
-if ($nzline && ($iplot_recon==1)) {`pssac2 $Zrecon $proj $bounds -Ent-3 -M${sizeZ3} $recon_pen -N -K -O >> $ps_file`;}
+if (-f $Zsyn){
+`pssac2 $Zsyn $proj $bounds -Ent-3 -M${sizeZ2} $red_pen -N -O -K >> $ps_file`;}
+if ($Zflag && $nzline && ($iplot_recon==1)) {`pssac2 $Zrecon $proj $bounds -Ent-3 -M${sizeZ3} $recon_pen -N -K -O >> $ps_file`;}
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC Z\" | pstext $proj $bounds -N -K -O >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strZ\" | pstext $proj $bounds -N -O -K >> $ps_file`;
@@ -718,9 +739,9 @@
# subroutine name: max
# Input: number1, number2
# returns greater of 2 numbers
-sub max {
- if ($_[0]<$_[1]) {return $_[1]} else {return $_[0]};
-}
+#sub max {
+# if ($_[0]<$_[1]) {return $_[1]} else {return $_[0]};
+#}
sub get_cmt {
my ($cmt_file)=@_;
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.tex
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.tex 2012-11-27 20:16:42 UTC (rev 21080)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.tex 2012-11-27 23:55:19 UTC (rev 21081)
@@ -199,7 +199,7 @@
\begin{verbatim}
-0.585 0.0110 18200 # tstart, DT, npts: time vector for simulations
7 # imeas (1-8; see manual)
- BH # channel of synthetics: BH or LH
+ BH # channel of synthetics: BH or LH
30.000 6.000 # TLONG and TSHORT: band-pass periods for records
.false. # RUN_BANDPASS: use band-pass on records
.true. # DISPLAY_DETAILS
@@ -398,7 +398,7 @@
prefix.recon_syn_dt.sac # reconstructed syn with only dtau(om) applied
prefix.err_ph/abs/dt/dlnA[_full] # err_phi/abs/dt/dlnA_mtm(idf_new:idf_new:i_right)
prefix.freq_limits # fstart/fend (input),df, fstart/fend(adjusted for given window)
-prefix0.recon.cc.sac # recon. syn for the entire data trace
+prefix0.recon.sac # recon. syn after applying CC or MTM for the entire data trace
\end{verbatim}
IF \verb+DISPLAY_MORE_DETAILS = .true.+ in the parameter file, following files will be written for every window pair
@@ -427,6 +427,7 @@
prepare_adj_src.pl -m CMTSOLUTION_9818433 -s PLOTS/STATIONS_TOMO -o ADJOINT_SOURCES \
OUTPUT_FILES/*adj
\end{verbatim}
+It loops over all the stations with adjoint sources present in \verb+OUTPUT_FILES+, and write an associate \verb+STATION_ADJOINT+ to local dir. It also rotates adjoint sources to E-N-Z, and deposits them into \verb+ADJOINT_SOURCES+ dir.
The key is that it will create an all-zeros record if no measurement is made on a particular component (say, Z), but IS made on another component (say, R or T). It uses the following fortran
program:
\begin{verbatim}
@@ -463,7 +464,11 @@
\end{verbatim}
gives plots of detailed information on data, syn, reconstructed syn and corresponding adjoint sources, as well as
source/station distributions on a map (Figure \ref{fig:iker07} for MTM, \verb+imeas=7+).
+Some tweaking is probably necessary before successful application to a different dataset, including: setting
+evid in CMT file (correspond to data name evid.net.sta), setting corresponding kstnam,knetwk,kcmpnm for both
+data and synthetics, setting o header properly, and all synthetics need to be present (data/recon-syn not needed).
+
\item plot the statistics of window files (check file structure before use)
\begin{verbatim}
plot_win_stats_all.pl Tmin/Tmax model iplot [plot histogram or not]
@@ -604,6 +609,7 @@
\en
For instance, weights can be given to windows in \verb+P_SV/SH/Rayleigh/Love+ - \verb+Z/R/T+, 6 different categories of measurements separately. They are only used to change the traveltime adjoint sources (and corresponding $\phi$ values) at this
point.
+\red{How about weighting related to the noise-level of the data trace itself?}
\subsection{Time-domain taper}
Time domain tapers are first applied to windowed data and synthetics (no matter what other options are)
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/ma_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/ma_sub.f90 2012-11-27 20:16:42 UTC (rev 21080)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/ma_sub.f90 2012-11-27 23:55:19 UTC (rev 21081)
@@ -11,7 +11,8 @@
subroutine mt_measure(datafile,filename,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)
+ i_pmax_dat,i_pmax_syn,i_right,trans_w,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,syn_dtw_mt, &
+ err_dt,err_dlnA)
! ======================================================================
! subroutine mt_measure(): making measurements on windowed data and synthetics
@@ -54,13 +55,13 @@
double precision, intent(out) :: tshift,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,sigma_dt,sigma_dlnA
integer, intent(out) :: istart,nlen,i_pmax_dat,i_pmax_syn,i_right
- double precision, dimension(:), intent(out) :: syn_dtw,dat_dtw,syn_dtw_cc,dtau_w,dlnA_w
+ double precision, dimension(:), intent(out) :: syn_dtw,dat_dtw,syn_dtw_cc,syn_dtw_mt,dtau_w,dlnA_w
complex*16, dimension(:), intent(out) :: trans_w
double precision, dimension(:), intent(out), optional :: err_dt,err_dlnA
! local variables
- double precision, dimension(NPT) :: syn_dtw_mt, syn_dtw_mt_dt, &
+ double precision, dimension(NPT) :: syn_dtw_mt_dt, &
syn_dtw_cc_dt, dat_dtw_cc, syn_dtw_h, dat_dtw_h
double precision :: sfac1,fac,f0,df,df_new,dw, &
ampmax_unw,wtr_use_unw,ampmax,wtr_use,wtr_mtm,dtau_wa,dlnA_wa !omega
@@ -663,7 +664,7 @@
! note: 2 is needed for the integration from -inf to inf
df = 1. /(NPT*dt)
ffac = 2.0 * df * sum(w_taper(i_left:i_right) ) ! CHT: 1 --> i_left
- if (DISPLAY_DETAILS) print *, ' f-dom taper normalization factor, ffac = ', ffac
+ if (DISPLAY_DETAILS) print *, ' f-dom taper normalization factor, ffac = ', sngl(ffac)
! wp_taper and wq_taper are modified frequency-domain tapers
! Notice the option to include the frequency-dependent error.
@@ -1958,7 +1959,7 @@
! for data_file trace.
- ! integer header variables
+ ! integer header variables (these values are not really used!)
call getnhv('nzyear',yr,nerr)
call getnhv('nzjday',jda,nerr)
call getnhv('nzhour',ho,nerr)
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/measure_adj.f90 2012-11-27 20:16:42 UTC (rev 21080)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/measure_adj.f90 2012-11-27 23:55:19 UTC (rev 21081)
@@ -32,7 +32,7 @@
character(len=150) :: datafile,synfile,file_prefix,file_prefix0,file_prefix2,measure_file_prefix,adj_file_prefix
integer :: num_meas, j, ios, npt1, npt2, npts, nn
- double precision, dimension(NDIM) :: data, syn, adj_syn_all, tr_adj_src, am_adj_src, recon_cc_all, syn_dtw_cc
+ double precision, dimension(NDIM) :: data, syn, adj_syn_all, tr_adj_src, am_adj_src, recon_cc_all, syn_dtw_cc, syn_dtw_mt
double precision :: t01, dt1, t02, dt2, t0, dt, tstart, tend, tt, dtt, df
double precision, dimension(NCHI) :: window_chi
double precision :: fend0, fstart0, fend, fstart
@@ -199,7 +199,7 @@
! also compute reconstructed synthetics for CC (and MT, if specified) measurements
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_pmax_dat,i_pmax_syn,i_right0,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,syn_dtw_mt,err_dt,err_dlnA)
i_right = i_right0
i_left = 1 ! LQY: is it feasible that i_left is not 1? mt_adj() inherently assumes it.
@@ -226,7 +226,7 @@
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)
+ i_pmax_dat,i_pmax_syn,i_right,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,syn_dtw_mt)
else
print *, ' using this MTM. '
endif
@@ -330,7 +330,11 @@
endif
! combine CC-reconstructed records
- recon_cc_all(istart:istart+nlen-1) = recon_cc_all(istart:istart+nlen-1) + syn_dtw_cc(1:nlen)
+ if (imeas >= 7) then
+ recon_cc_all(istart:istart+nlen-1) = recon_cc_all(istart:istart+nlen-1) + syn_dtw_mt(1:nlen)
+ else
+ recon_cc_all(istart:istart+nlen-1) = recon_cc_all(istart:istart+nlen-1) + syn_dtw_cc(1:nlen)
+ endif
endif ! COMPUTE_ADJOINT_SOURCE
@@ -349,7 +353,7 @@
! write out the CC-reconstructed data from synthetics
if (OUTPUT_MEASUREMENT_FILES) &
- call dwsac1(trim(file_prefix2)//'.recon.cc.sac',recon_cc_all,npts,t0,dt)
+ call dwsac1(trim(file_prefix2)//'.recon.sac',recon_cc_all,npts,t0,dt)
! OPTIONAL: A conservative choice is to filter the adjoint source,
! since higher frequencies could enter from the tapering operations.
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/prepare_adj_src.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/prepare_adj_src.pl 2012-11-27 20:16:42 UTC (rev 21080)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/prepare_adj_src.pl 2012-11-27 23:55:19 UTC (rev 21081)
@@ -24,7 +24,7 @@
#Reading input arguments:
@ARGV>0 || die("prepare_adj_src.pl -m CMT -z BH -s STATION -o OUTDIR -i all_adj_RTZ_files\n");
-if (!getopts('m:z:s:o:i:')) {die("Check arguments\n");}
+if (!getopts('m:z:s:o:i')) {die("Check arguments\n");}
if ($opt_m) {$cmt = $opt_m;} else {$cmt = "CMTSOLUTION";}
if ($opt_z) {$chan = $opt_z;} else {$chan = "BH";}
@@ -38,6 +38,7 @@
($elat,$elon) = get_cmt_location($cmt);
# fill up the hash table based on adjoint sources
+print "Processing adjoint source files ...\n";
foreach $file (@ARGV) {
if (not -f $file) {die("Check file $file\n");}
($basefile) = basename($file);
@@ -47,7 +48,7 @@
$adj{$stanet}{$cmp} = $file;
if (defined $adj{$stanet}{all}) {$adj{$stanet}{all} ++ ;}
else {$adj{$stanet}{all} = 1;}
- print "$basefile -- $stanet\n";
+ print " $basefile -- $stanet $cmp\n";
}
#die("testing");
@@ -60,11 +61,12 @@
# NOTE: channel is either BH or LH (from synthetics)
$stafile_temp = "station.tmp";
-system("rm -f ${stafile_temp}");
+if (-f ${stafile_temp}) {system("rm -f ${stafile_temp}");}
$nstation = 0;
-print "Copying stations adjoint source file: $outdir ...\n";
+print "Copying station adjoint source files to $outdir ...\n";
foreach $stanet (keys(%adj)) {
- ($sta,$net) = split(/\./,$stanet);
+ ($sta,$net) = split(/\./,$stanet);
+
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;
@@ -89,7 +91,7 @@
#print "\n -- $sta1,$net1,$slat1,$slon1 -- \n";
if("$sta1.$net1" eq $stanet) {
$slon = $slon1; $slat = $slat1;
- print "*** $stanet ***\n";
+ print " *** $stanet ***\n";
$text = sprintf("%8s %4s %.6f %.6f %.6f %.6f",$sta1,$net1,$slat1,$slon1,$selev,$sburial);
`echo $text >> ${stafile_temp}`;
$nstation ++ ;
@@ -105,7 +107,7 @@
# get back-azimuth and apply rotation using rotate_adj_src.f90,
# which will WRITE OUT the adjoint sources on the east and north components
(undef,undef,undef,$baz) = delaz5($slat,$slon,$elat,$elon,0); # station - event loc in radians
- # print "***$tcomp***\nbin/rotate_adj_src $baz $zcomp $tcomp $rcomp $ecomp $ncomp\n\n";
+ #print "***$tcomp***\nbin/rotate_adj_src $baz $zcomp $tcomp $rcomp $ecomp $ncomp\n\n";
system("rotate_adj_src $baz $zcomp $tcomp $rcomp $ecomp $ncomp");
if ($? != 0) {
die("Error: rotate_adj_src $baz $tcomp $rcomp $ecomp $ncomp\n");
@@ -113,7 +115,7 @@
system("cp -f $ecomp $ncomp $zcomp $outdir");
}
- }
+ } else {die("check hash table for adj{$stanet}{$cmpT,$cmpR,$cmpZ}\n");}
}
# write out STATIONS_ADJOINT file
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/rotate_adj_src.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/rotate_adj_src.f90 2012-11-27 20:16:42 UTC (rev 21080)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/rotate_adj_src.f90 2012-11-27 23:55:19 UTC (rev 21081)
@@ -37,8 +37,8 @@
! at least one file (T,R,Z) should be present
if (.not. t_exist .and. .not. r_exist) then
if (.not. z_exist) stop 'At least one file should exist: zfile, tfile, rfile'
- ! no need to read Z comp adjoint source (keep the same)
- ! call drascii(zfile,tdata,npts,t0,dt)
+ ! need to read Z comp adjoint source for [to,dt,npts]
+ call drascii(zfile,tdata,npts,t0,dt)
endif
! read in T file
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/run_measure_adj.csh
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/run_measure_adj.csh 2012-11-27 20:16:42 UTC (rev 21080)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/run_measure_adj.csh 2012-11-27 23:55:19 UTC (rev 21081)
@@ -11,7 +11,7 @@
cp window_chi PLOTS
\rm -rf ADJOINT_SOURCES
mkdir ADJOINT_SOURCES
-prepare_adj_src.pl -m CMTSOLUTION_9818433 -s PLOTS/STATIONS_TOMO -o ADJOINT_SOURCES OUTPUT_FILES/*adj
+prepare_adj_src.pl -m CMTSOLUTION_9818433 -s PLOTS/STATIONS_TOMO -o ADJOINT_SOURCES -i OUTPUT_FILES/*adj
cp STATIONS_ADJOINT ADJOINT_SOURCES
\mv STATIONS_ADJOINT PLOTS
cd PLOTS
More information about the CIG-COMMITS
mailing list