[cig-commits] r15864 - in seismo/3D/ADJOINT_TOMO/measure_adj: . PLOTS USER_MANUAL scripts_tomo
carltape at geodynamics.org
carltape at geodynamics.org
Wed Oct 21 16:40:14 PDT 2009
Author: carltape
Date: 2009-10-21 16:40:13 -0700 (Wed, 21 Oct 2009)
New Revision: 15864
Added:
seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/
seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/figures/
seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.pdf
seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.tex
Modified:
seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/9818433_T006_T030_MPM_CI_m16_iker07_win_adj.pdf
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/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:
Added a parameter to make it optional to compute adjoint sources. Started a User Manual.
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR 2009-10-21 22:34:28 UTC (rev 15863)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR 2009-10-21 23:40:13 UTC (rev 15864)
@@ -1,14 +1,15 @@
-0.585 0.0110 18200 # tstart, DT, npts: time vector for simulations
- 7 # iker (0-8) -- 0 no adj src; 1-8 see manual;
+ 7 # imeas (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
.true. # DISPLAY_DETAILS
.true. # OUTPUT_MEASUREMENT_FILES
+ .true. # COMPUTE_ADJOINT_SOURCE
-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 # 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
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/9818433_T006_T030_MPM_CI_m16_iker07_win_adj.pdf
===================================================================
(Binary files differ)
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-21 22:34:28 UTC (rev 15863)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2009-10-21 23:40:13 UTC (rev 15864)
@@ -12,7 +12,7 @@
# (Search for "USER".)
#
# EXAMPLE:
-# 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 -n MPM/CI/BH -b 0 -l -10/200 -k 7/1 -a STATIONS_ADJOINT -d DATA -s SYN -c RECON -w MEASUREMENT.WINDOWS -i m16 -j 6/30
#
#-----------------------------------
@@ -28,7 +28,7 @@
-m -- CMTSOLUTION file to plot event focal mechanism (beach ball) on the map
-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)
+ -k -- type of measurement (1-8) / plot adjoint source(0 or 1)
-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
@@ -46,8 +46,7 @@
if (!getopts('k:m:l:n:b:ra:d:s:c:w:i:j:')) {die('Check input arguments\n');}
print "STATIONS file is $opt_a\n";
if ($opt_m) {@cmtvals = &get_cmt($opt_m);}
-#if ($opt_k) {$iker = $opt_k;} # does not work for iker = 0
-$iker = $opt_k; print "\n $iker $opt_k";
+if ($opt_k) {($imeas,$iplot_adj) = split(/\//,$opt_k);}
if ($opt_l) {($tmin,$tmax) = split(/\//,$opt_l);}
if ($opt_w) {$Winfile=$opt_w;}
if ($opt_n) {($sta,$net,$chan)=split(/\//,$opt_n);}
@@ -60,17 +59,17 @@
$smodel = $opt_i;
# suffix for adjoint sources
-$klab = sprintf('iker%2.2i',$iker);
-$klab2 = sprintf('iker%2.2i',$iker-2);
+$klab = sprintf('iker%2.2i',$imeas);
+$klab2 = sprintf('iker%2.2i',$imeas-2);
@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];
+$ktitle = $ktitles[$imeas-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");}
+#if($imeas == 0) {$klab = "ik"; $ktitle = "waveform";}
+#elsif($imeas == 1) {$klab = "mtm"; $ktitle = "multitaper TT";}
+#elsif($imeas == 2) {$klab = "cc"; $ktitle = "cross-correlation TT";}
+#elsif($imeas == 3) {$klab = "bdcc"; $ktitle = "banana-doughnut CC-TT";}
+#else {die("ERROR: imeas ($imeas) must be 0,1, 2, or 3\n");}
$sTmin = sprintf("T%3.3i",$Tmin);
$sTmax = sprintf("T%3.3i",$Tmax);
@@ -90,20 +89,19 @@
# with iboth = 1, you need to have made both the CC and MT adjoint sources
if($iboth==1) {
- if( ($iker != 7) && ($iker != 8) ) {die("\nFor plotting CC and MT adjoint sources, use iker = 7 or 8.\n")}
+ if( ($imeas != 7) && ($imeas != 8) ) {die("\nFor plotting CC and MT adjoint sources, use imeas = 7 or 8.\n")}
$ktitle = "TT -- MT = black ; CC = red";
}
# plotting preferences
-$iplot_adj = 1; # plot adjoint sources
+#$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;}
+if ($imeas < 5) {$iplot_recon = 0; $iplot_CClabel = 0;}
# adjust plot size based on whether you have adjoint sources
if($iplot_adj==1) {$swid = 3.5; $fsizetitle = 10;}
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-21 22:34:28 UTC (rev 15863)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl 2009-10-21 23:40:13 UTC (rev 15864)
@@ -5,7 +5,7 @@
# plot_all_win_adj.pl
#
# EXAMPLE:
-# 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
+# plot_win_adj_all.pl -l -10/200 -m ../CMTSOLUTION_9818433 -n BH -b 0 -k 7/1 -a STATIONS_ADJOINT -d DATA -s SYN -c RECON -w MEASUREMENT.WINDOWS -i m16 -j 6/30
#
#----------------------------------------------------
@@ -15,7 +15,7 @@
sub Usage{
print STDERR <<END;
- 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
+ plot_all_win_adj.pl -m CMTFILE -a STATION_FILE -n chan -b iboth -l tmin/tmax -k imeas/iadj -d data_dir -s syn_dir -c recon_dir -w winfile -i smodel -j Tmin/Tmax
END
exit(1);
@@ -25,7 +25,7 @@
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_k) {$imeas = $opt_k;} # does not work for imeas = 1
if($opt_l) {$tcuts=$opt_l;}
if($opt_a) {$station_list=$opt_a;}
if($opt_n) {$chan=$opt_n;}
Added: seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.pdf
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.pdf
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.tex
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.tex (rev 0)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.tex 2009-10-21 23:40:13 UTC (rev 15864)
@@ -0,0 +1,213 @@
+\documentclass[11pt,titlepage,fleqn]{article}
+
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{latexsym}
+\usepackage[round]{natbib}
+\usepackage{xspace}
+%\usepackage{epsfig}
+\usepackage{graphicx}
+\usepackage{bm}
+
+%\usepackage{fancybox} # do not use when showing Table of Contents
+%\usepackage{fancyhdr}
+%\pagestyle{fancy}
+
+%-----------------------------------------------------------
+% SPACING COMMANDS (Latex Companion, p. 52)
+%-----------------------------------------------------------
+
+\usepackage{setspace}
+
+\renewcommand{\baselinestretch}{1.2}
+
+\textwidth 460pt
+\textheight 690pt
+\oddsidemargin 0pt
+\evensidemargin 0pt
+
+% see Latex Companion, p. 85
+\voffset -50pt
+\topmargin 0pt
+\headsep 20pt
+\headheight 15pt
+\headheight 0pt
+\footskip 30pt
+\hoffset 0pt
+
+\include{NEWCOMMANDS}
+
+\graphicspath{
+ {./figures/}
+}
+
+%-----------------------------------------------------------
+\begin{document}
+%-----------------------------------------------------------
+
+\begin{center}
+{\bf \Large MEAS--ADJ: A software package for making measurements and computing adjoint sources for seisic tomography}
+\end{center}
+
+\vspace{0.4cm}
+\begin{spacing}{1.0}
+\noindent Primary Developers: Qinya Liu, Carl Tape, Ying Zhou
+
+\noindent Author of Manual: Carl Tape
+
+\noindent \today
+
+\noindent
+\begin{verbatim}
+This file : /seismo/3D/measure_adj/USER_MANUAL/measure_adj_manual.tex
+\end{verbatim}
+
+\tableofcontents
+
+\end{spacing}
+
+%-------------------------------------------------------------------
+
+\pagebreak
+\section{Introduction}
+
+This is only a template manual at this point --- more to come!
+
+The open-source software package MEAS--ADJ is designed for making measurements for seismic tomography. Several misfit functions are available, such as a simple waveform difference, a simple cross-correlation traveltime difference, or a frequency dependent multitaper traveltime difference. For each misfit function, the code also provides adjoint sources, which are needed for tomographic inversions using adjoint methods \cite[\eg][]{Tromp2005,Tape2007}. This code was used to make measurements and adjoint sources for the (large-scale) tomographic inversion of \citet{Tape2009} for the southern California crust.
+
+MEAS--ADJ is designed to complement FLEXWIN \citep{Maggi2009}, which is designed to pick time windows for measurement. In particular, FLEXWIN has the option of providing a \verb+MEASUREMENT.WINDOWS+ file that can be directly read into MEAS--ADJ. However, FLEXWIN is not necessary to run MEAS--ADJ.
+
+Please email Carl Tape (\verb+carltape at fas.harvard.edu+) with suggestions or corrections.
+
+%-------------------------------------------------------------------
+
+\section{Getting started}
+
+%-------------
+
+\subsection{System requirements}
+
+In order to install and run, MEAS--ADJ requires:
+%
+\begin{itemize}
+\item UNIX operating system (Linux, Solaris, MacOS \ldots)
+\item GNU make
+\item a fortran compiler (gfortran, ifort, etc...)
+\item other packages : SAC (Seismic Analysis Code, available from IRIS); GMT (Generic Mapping Tools) for the plotting scripts
+\end{itemize}
+
+MEAS--ADJ requires the following libraries external to the package in order to compile and run: {\tt libsacio.a} and {\tt libsac.a}. Both libraries are distributed by IRIS as part of the SAC package (version 101.2 and above). The IRIS download site (as of 21-Oct-2009) is here: \verb+http://www.iris.edu/software/sac/sac.request.htm+. (To check your version, type sac.)
+
+%-------------
+
+\section{Obtaining the code}
+
+%The code is available as a gzipped tarball from CIG (Computational Infrastructure for Geodynamics, \url{http://www.geodynamics.org}). The tarball is unpacked by typing {\tt tar xvzf flexwin.tgz}.
+
+(Use SVN for now.)
+
+The package contains the MEAS--ADJ code and documentation, as well as a set of test data, examples of user files for different scenarios, and a set of utility scripts that may be useful for running MEAS--ADJ on large datasets.
+
+%-------------
+
+\section{Compilation}
+
+If your compiler of choice is gfortran, then you should be able to use the \verb+Makefile+ with only minor modifications (notably you may need to change the search path for the {\tt libsacio.a} library). If you prefer another compiler, you should modify the OPT and FC lines in the makefiles accordingly. We tested the code using gfortran version 4.3.3. (To check your version, type \verb+gfortran --version+.)
+
+Steps to compile the MEAS--ADJ package:
+%
+\begin{enumerate}
+\item \verb+make clean+
+\item \verb+make+
+\end{enumerate}
+
+You should end up with the \verb+measure_adj+ executable.
+
+%-------------
+
+\subsection{Running the test case}
+
+\begin{verbatim}
+mt_measure_adj
+csh -f run_mt_measure_adj.csh
+\end{verbatim}
+%
+If all goes well, you should reproduce the figure
+%
+\begin{verbatim}
+PLOTS/9818433_T006_T030_MPM_CI_m16_mtm_win_adj.pdf
+\end{verbatim}
+
+%-------------------------------------------------------------------
+
+\section{Running MEAS--ADJ}
+
+\subsection{Input files}
+
+\subsection{Output files}
+
+\subsection{Scripts}
+
+
+%-------------------------------------------------------------------
+
+\section{Measurement options}
+
+\begin{enumerate}
+\item No adjoint source\verb+iker = 0+
+\end{enumerate}
+
+%-------------------------------------------------------------------
+
+\section{Miscellaneous}
+
+\subsection{Bug reports and suggestions for improvements}
+
+To report bugs or suggest improvements to the code, please send an email to the CIG Computational Seismology Mailing List (cig-seismo at geodynamics.org) or Carl Tape (carltape at fas.harvard.edu), and/or use our online bug tracking system Roundup (www.geodynamics.org/roundup).
+
+%-------------
+
+\subsection{Notes and Acknowledgments}
+
+The main developers of the MEAS-ADJ source code are Qinya Liu, Carl Tape, and Ying Zhou.
+The multitaper measurement capability originated from codes used by Ying Zhou \citep[][]{YZhou2004,YZhou2005}.
+The multitaper adjoint sources were implemented by Carl Tape \citep[][Appendix~C]{Tape2009phd}.
+The organizational structure of the package was made by Qinya Liu, with adaptations by Carl Tape.
+
+The following individuals have also contributed to the development of the source code: Vala Hjorleifsdottir.
+The following individuals contributed to this manual: Carl Tape.
+
+The MEAS-ADJ code makes use of filtering and enveloping algorithms that are part of SAC (Seismic Analysis Code, Lawerence Livermore National Laboratory) provided for free to IRIS members. We thank Brian Savage for adding interfaces to these algorithms in recent SAC distributions.
+
+We acknowledge support by the National Science Foundation under grant EAR-0711177.
+
+%-------------
+
+\subsection{License}
+
+(See FLEXWIN manual for possible options.)
+
+%-------------------------------------------------------------------
+%\pagebreak
+%\small
+\begin{spacing}{1}
+\addcontentsline{toc}{section}{References}
+\bibliographystyle{agu08}
+\bibliography{preamble,REFERENCES,refs_carl}
+\end{spacing}
+%\normalsize
+%-------------------------------------------------------------------
+
+%\begin{figure}
+%\hspace{-1cm}
+%\includegraphics[width=18cm]{Hollywood_CMT.eps}
+%\caption[]
+%{{
+%CMT example figure for the {\tt SPECFEM3D$\_$BASIN} manual.
+%\label{fig:CMT_basin}
+%}}
+%\end{figure}
+
+%-------------------------------------------------------------------
+\end{document}
+%-------------------------------------------------------------------
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-10-21 22:34:28 UTC (rev 15863)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-10-21 23:40:13 UTC (rev 15864)
@@ -23,9 +23,6 @@
integer :: yr,jda,ho,mi
double precision :: sec,dist,az,baz,slat,slon
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, 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
@@ -110,18 +107,18 @@
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
+ write(adj_file_prefix,'(a,i2.2)') trim(file_prefix2)//'.iker', imeas0
+!!$ if (imeas == 0) then
!!$ adj_file_prefix = trim(file_prefix2) // '.ik'
-!!$ else if (iker == 1) then
+!!$ else if (imeas == 1) then
!!$ adj_file_prefix = trim(file_prefix2) // '.mtm'
-!!$ else if (iker == 2) then
+!!$ else if (imeas == 2) then
!!$ adj_file_prefix = trim(file_prefix2) // '.cc'
-!!$ else if (iker == 3) then
+!!$ else if (imeas == 3) then
!!$ adj_file_prefix = trim(file_prefix2) // '.bdcc'
!!$ else
-!!$ print *, 'iker = ', iker
-!!$ print *, 'Error: iker must be 0, 1, 2, or 3'
+!!$ print *, 'imeas = ', imeas
+!!$ print *, 'Error: imeas must be 0, 1, 2, or 3'
!!$ stop
!!$ endif
@@ -195,7 +192,7 @@
if(.not. use_trace) then
!stop 'Check why this MT measurement was rejected'
print *, 'Reverting from multitaper measurement to cross-correlation measurement'
- iker = iker0 - 2
+ imeas = imeas0 - 2
is_mtm = 3
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,&
@@ -208,7 +205,7 @@
endif
! check that the CC measurements are within the specified input range
- if (iker >= 5) call cc_measure_select(tshift,dlnA,cc_max)
+ if (imeas >= 5) call cc_measure_select(tshift,dlnA,cc_max)
! write frequency limits to file
if (OUTPUT_MEASUREMENT_FILES) then
@@ -235,16 +232,16 @@
! 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_syn,j,iker,&
+ file_prefix0,sta,net,chan_syn,j,imeas,&
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
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
+ if (mod(imeas,2)==1) then
+ adj_syn_all(:) = adj_syn_all(:) + tr_adj_src(:) ! imeas = 1,3,5,7
else
- adj_syn_all(:) = adj_syn_all(:) + am_adj_src(:) ! iker = 2,4,6,8
+ adj_syn_all(:) = adj_syn_all(:) + am_adj_src(:) ! imeas = 2,4,6,8
endif
endif
@@ -255,7 +252,7 @@
! CHT: (re-)set to multitaper parameters, if originally specified
if (is_mtm0 == 1) then
- iker = iker0
+ imeas = imeas0
is_mtm = is_mtm0
endif
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-10-21 22:34:28 UTC (rev 15863)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-10-21 23:40:13 UTC (rev 15864)
@@ -174,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,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc)
+ call write_average_meas(filename,imeas,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc)
!========================================
@@ -357,7 +357,7 @@
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, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
+ !call write_average_meas(filename, imeas, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
endif
enddo ! ictapers
@@ -403,7 +403,7 @@
sigma_dt = sigma_dt_cc ; sigma_dlnA = sigma_dlnA_cc
! write average multitaper measurement to file
- call write_average_meas(file_prefix, iker, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
+ call write_average_meas(file_prefix, imeas, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
!-------------------------------------------------------------------------------
! multitaper error estimation
@@ -579,7 +579,7 @@
! adjoint sources by assimulate the measurements passed from mt_measure()
!
! Input:
- ! iker -- adjoint source type: 0 for waveform, 1 for multitaper, 2 for cc, 3 for cc banana-doughnut
+ ! imeas -- adjoint source type: 0 for waveform, 1 for multitaper, 2 for cc, 3 for cc banana-doughnut
! istart -- starting index of the windowed portion of original trace, used to generate adjoint
! source that correspond to the original synthetics
! dat_dtw(:), syn_dtw(:), nlen, dt -- windowed data and synthetics
@@ -624,20 +624,20 @@
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 (imeas == 0 .and. (present(am_adj_src) .or. present(am_chi))) stop &
+ ! 'am_adj_src and am_chi are not needed for imeas = 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
+ if(imeas==1 .or. imeas==2) then
print *, 'computing waveform adjoint source'
- elseif(iker==3 .or. iker==4) then
+ elseif(imeas==3 .or. imeas==4) then
print *, 'computing banana-doughtnut adjoint source'
- elseif(iker==5 .or. iker==6) then
+ elseif(imeas==5 .or. imeas==6) then
print *, 'computing cross-correlation adjoint source'
- elseif(iker==7 .or. iker==8) then
+ elseif(imeas==7 .or. imeas==8) then
print *, 'computing multitaper adjoint source'
endif
@@ -716,7 +716,7 @@
! ----------------------------------
! CROSS CORRELATION ADJOINT SOURCES
! ----------------------------------
- if( (iker >= 3).and.(iker <= 6) ) then
+ if( (imeas >= 3).and.(imeas <= 6) ) then
! compute synthetic velocity
do i = 2, nlen-1
@@ -820,7 +820,7 @@
tr_adj_src = 0.
am_adj_src = 0.
- ! CHT: Compute the integrated waveform difference squared
+ ! integrated waveform difference squared
waveform_temp1 = 0. ; waveform_temp2 = 0. ; waveform_temp3 = 0.
do i = 1,nlen
waveform_temp1 = waveform_temp1 + ( dat_dtw(i) * time_window(i) )**2
@@ -832,34 +832,39 @@
waveform_s2 = waveform_temp2
waveform_chi = waveform_temp3
- ! compute traveltime and amplitude adjoint sources
- do i = 1,nlen
+ if (COMPUTE_ADJOINT_SOURCE) then
- i1 = istart + i
+ ! compute traveltime and amplitude adjoint sources
+ do i = 1,nlen
- ! 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
+ i1 = istart + 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)
+ ! waveform
+ if(imeas==1 .or. imeas==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 this by waveform_d2
- ! 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(imeas==3 .or. imeas==4) 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)
+ ! cross-correlation
+ elseif(imeas==5 .or. imeas==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)
- endif
- enddo
+ ! multitaper
+ elseif(imeas==7 .or. imeas==8) then
+ tr_adj_src(i1) = fp(i) * time_window(i)
+ am_adj_src(i1) = fq(i) * time_window(i)
+ endif
+ enddo
+
+ endif
+
! -------------------------------------
! COMPUTE MISFIT FUNCTION VALUE
! -------------------------------------
@@ -902,15 +907,15 @@
!!$ enddo
!!$ close(88)
- if(iker <= 2) then ! waveform
+ if(imeas <= 2) then ! waveform
tr_chi = 0.5 * waveform_chi
am_chi = 0.5 * waveform_chi
- elseif( (iker >= 3).and.(iker <= 6) ) then ! cross_correlation
+ elseif( (imeas >= 3).and.(imeas <= 6) ) then ! cross_correlation
tr_chi = window_chi(3)
am_chi = window_chi(4)
- elseif( (iker==7).or.(iker==8) ) then ! multitaper
+ elseif( (imeas==7).or.(imeas==8) ) then ! multitaper
tr_chi = window_chi(1)
am_chi = window_chi(2)
@@ -1292,33 +1297,30 @@
enddo
cc = cc/norm
- !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) <= TSHIFT_MAX) ) then
- cc_max = cc
- ishift = i
+ if (cc > cc_max) then
+ ! CHT: do not allow time shifts larger than the specified input range
+ ! This is an important criterion, since it may pick TSHIFT_MIN or TSHIFT_MAX
+ ! if cc_max within the interval occurs on the boundary.
+ if( (i*dt >= TSHIFT_MIN).and.(i*dt <= TSHIFT_MAX) ) then
+ cc_max = cc
+ ishift = i
+ endif
endif
enddo
tshift = ishift*dt
- ! amplitude MEASUREMENT (definition of Dahlen and Baig (2002), Eq. 3,17,18 : dlnA = Aobs/Asyn - 1)
- ! note that the records are UN-SHIFTED
- !dlnA = sqrt( (dt * sum( dat_dtw0(i1:i2) * dat_dtw0(i1:i2) )) &
- ! / (dt * sum( syn_dtw(i1:i2) * syn_dtw(i1:i2) )) ) - 1.
-
- ! CHT revision 15-Oct-2008
- ! The previously used expression for dlnA is a first-order perturbation of ln(A1/A2) = (A1-A2)/A2
- ! The new expression is better suited to getting values between -1 and 1,
- ! with dlnA = 0 indicating perfect fit, as before.
+ ! The previously used expression for dlnA of Dahlen and Baig (2002),
+ ! is a first-order perturbation of ln(A1/A2) = (A1-A2)/A2 .
+ ! The new expression is better suited to getting Gaussian-distributed
+ ! values between -1 and 1, with dlnA = 0 indicating perfect fit, as before.
dlnA = 0.5 * log( sum(data(i1:i2) * data(i1:i2)) / sum(syn(i1:i2) * syn(i1:i2)) )
end subroutine compute_cc
! ---------------------------------------------------------------------------
- subroutine compute_average_error(data_dtw, syn_dtw_cc, syn_dtw_cc_dt, &
- nlen, dt, sigma_dt, sigma_dlnA)
+ subroutine compute_average_error(data_dtw,syn_dtw_cc,syn_dtw_cc_dt,nlen,dt,sigma_dt,sigma_dlnA)
! CHT: Estimate the uncertainty in the CC measurement
! based on the integrated waveform difference between the data
@@ -1379,20 +1381,20 @@
! ---------------------------------------------------------------------------
- subroutine write_average_meas(filename, iker, dtau_meas, dlnA_meas, dtau_sigma, dlnA_sigma)
+ subroutine write_average_meas(filename, imeas, dtau_meas, dlnA_meas, dtau_sigma, dlnA_sigma)
character(len=*), intent(in) :: filename
double precision, intent(in) :: dtau_meas, dlnA_meas, dtau_sigma, dlnA_sigma
- integer, intent(in) :: iker
+ integer, intent(in) :: imeas
character(len=40) :: stlab, suffix
- if ( iker == 7 .or. iker == 8 ) then
+ if ( imeas == 7 .or. imeas == 8 ) then
stlab = 'Multitaper' ; suffix = 'average'
else
stlab = 'Cross-correlation' ; suffix = 'cc'
endif
- if ( iker .ge. 3 ) then
+ if ( imeas .ge. 3 ) then
if (DISPLAY_DETAILS) then
print *, trim(stlab)//' average measurements:'
print *, ' traveltime :', sngl(dtau_meas), ' +/- ', sngl(dtau_sigma)
@@ -1736,12 +1738,13 @@
open(10,file='MEASUREMENT.PAR',status='old',iostat=ios)
read(10,*) tt,dtt,nn
- read(10,*) iker0
+ read(10,*) imeas0
read(10,*) chan
read(10,*) TLONG, TSHORT
read(10,*) RUN_BANDPASS
read(10,*) DISPLAY_DETAILS
read(10,*) OUTPUT_MEASUREMENT_FILES
+ read(10,*) COMPUTE_ADJOINT_SOURCE
read(10,*) TSHIFT_MIN, TSHIFT_MAX
read(10,*) DLNA_MIN, DLNA_MAX
read(10,*) CC_MIN
@@ -1756,13 +1759,13 @@
read(10,*) NCYCLE_IN_WINDOW
close(10)
- iker = iker0
+ imeas = imeas0
! check the read-in values
print *, 'INPUTS FROM MEASUREMENT.PAR :'
print *, ' tt, dtt, nn : ',tt,dtt,nn
print *, ' chan :',chan
- print *, ' iker : ',iker
+ print *, ' imeas : ',imeas
print *, ' TLONG, TSHORT : ',TLONG, TSHORT
fstart0 = 1./TLONG ; fend0 = 1./TSHORT
print *, ' fstart, fend :', fstart0, fend0
@@ -1795,28 +1798,23 @@
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
+ if ( (itaper==1) .and. ((imeas.ge.3).and.(imeas.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
+ if ( (imeas==1).or.(imeas==2) ) then
is_mtm0 = 0
- elseif ( (iker.ge.3) .and. (iker.le.6) ) then
+ elseif ( (imeas.ge.3).and.(imeas.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
+ elseif ( (imeas==7).or.(imeas==8) ) then
+ is_mtm0 = 1 ! multitaper required for MT adjoint source
else
- print *, 'Error: iker must by 0-8'
+ print *, 'Error: imeas must by 1-8'
stop
endif
is_mtm = is_mtm0
-
- print *, ' COMPUTE_ADJOINT_SOURCE :',COMPUTE_ADJOINT_SOURCE
print *, ' is_mtm :',is_mtm
end subroutine read_par_file
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_variables.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_variables.f90 2009-10-21 22:34:28 UTC (rev 15863)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_variables.f90 2009-10-21 23:40:13 UTC (rev 15864)
@@ -21,7 +21,7 @@
double precision :: DT_SIGMA_MIN, DLNA_SIGMA_MIN
integer :: ntaper, ipwr_t, ipwr_w, ERROR_TYPE
- integer :: iker0, iker, itaper, is_mtm0, is_mtm
+ integer :: imeas0, imeas, itaper, is_mtm0, is_mtm
logical :: DISPLAY_DETAILS,OUTPUT_MEASUREMENT_FILES,RUN_BANDPASS,COMPUTE_ADJOINT_SOURCE
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-21 22:34:28 UTC (rev 15863)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl 2009-10-21 23:40:13 UTC (rev 15864)
@@ -12,12 +12,12 @@
# We do not need both for the adjoint tomography inversion.
#
# EXAMPLE:
-# 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
+# run_mt_cc_plot.pl m16 -10/200 0 -0.585/0.011/18200 7 BH 6/30 0/1/1/1 -5.0/5.0/-1.5/1.5/0.7 1/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 xxx\n")}
-($smodel,$lcut,$itest,$tvec,$iker,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
+($smodel,$lcut,$itest,$tvec,$imeas,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
#($ibp,$Ts,$lcut,$tvec,$itest,$iparbools,$par1,$par2,$par3) = @ARGV;
$pwd = $ENV{PWD};
@@ -41,7 +41,7 @@
`rm PLOTS/*ps PLOTS/*jpg PLOTS/*pdf`;
# first run in cross-correlation TT mode
-$iker = 5;
+$imeas = 5;
$itaper = 3;
# replace itaper in par3 with this value
@@ -72,8 +72,8 @@
# 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";
-print CSH "write_par_file.pl $tvec $iker $chan $Ts $iparbools $par1 $par2 $par3\n";
+#print CSH "write_par_file.pl OUTPUT_FILES $ibp $Ts $tvec $itaper $imeas $wtr/$npi $iparbools $par1 $par2 $par3\n";
+print CSH "write_par_file.pl $tvec $imeas $chan $Ts $iparbools $par1 $par2 $par3\n";
#---------------------------------------------
@@ -90,7 +90,7 @@
#////////////////////////////////////////////////////////
# second run in multitaper TT mode
-$iker = 7;
+$imeas = 7;
$itaper = 1;
# replace itaper in par3 with this value
@@ -104,7 +104,7 @@
print CSH "cp window_chi PLOTS\n";
# create MEASUREMENT.PAR file
-print CSH "write_par_file.pl $tvec $iker $chan $Ts $iparbools $par1 $par2 $par3\n";
+print CSH "write_par_file.pl $tvec $imeas $chan $Ts $iparbools $par1 $par2 $par3\n";
# run the measurement code and save output file
# (re-compilation is not necessary)
@@ -120,7 +120,7 @@
print CSH "\\mv STATIONS_ADJOINT PLOTS\n";
print CSH "cd PLOTS\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 "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -n $chan -b 1 -k $imeas -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-21 22:34:28 UTC (rev 15863)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl 2009-10-21 23:40:13 UTC (rev 15864)
@@ -21,12 +21,12 @@
# ... see write_par_file.pl
#
# 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 m16 -10/200 0 1 0 -0.585/0.011/18200 7 BH 6/30 0/1/1/1 -5.0/5.0/-1.5/1.5/0.7 1/1.0/0.5 1/0.02/2.5/2.0/2.5/3.5/1.5
#
#==========================================================
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;
+($smodel,$lcut,$itest,$iplot,$iboth,$tvec,$imeas,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
$pwd = $ENV{PWD};
@@ -85,7 +85,6 @@
if($CC_MIN < 0) {$CC_MIN = 0}
$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";
@@ -95,8 +94,8 @@
# create MEASUREMENT.PAR file
#$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 $tvec $iker $chan $Ts $iparbools $par1 $par2 $par3\n";
+#print CSH "write_par_file.pl OUTPUT_FILES $ibp $Ts $tvec $itaper $imeas $wtr/$npi $iparbools $par1 $par2 $par3\n";
+print CSH "write_par_file.pl $tvec $imeas $chan $Ts $iparbools $par1 $par2 $par3\n";
#---------------------------------------------
@@ -133,7 +132,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 -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";
+ print CSH "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -n $chan -b $iboth -k $imeas -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-21 22:34:28 UTC (rev 15863)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2009-10-21 23:40:13 UTC (rev 15864)
@@ -12,10 +12,10 @@
# NOTE: The command run_mt_measure_adj.pl below contains additional input parameters.
#
# 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
+# run_tomo.pl m16 0.011 1.0 0 BH 6/30 0/1/1/1 -5.0/5.0/-1.5/1.5/0.7 1/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/1 -5.0/5.0/-1.5/1.5/0.7 1/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/1 -3.0/3.0/-1.5/1.5/0.7 1/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/1 -2.0/2.0/-1.5/1.5/0.7 1/1.0/0.5 1/0.02/2.5/0.25/0.5/3.5/1.5
#
# socal m00 - m04, 0.2/0.2
# socal m05 - m08, 0.4/0.4
@@ -27,7 +27,7 @@
#==========================================================
if (@ARGV < 10) {die("Usage: run_tomo.pl xxx\n")}
-($smodel,$dt,$lcut_frac,$iker,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
+($smodel,$dt,$lcut_frac,$imeas,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
#($smodel,$ibp,$Ts,$dt,$iparbools,$par1,$par2,$par3,$lcut_frac) = @ARGV;
$pwd = $ENV{PWD};
@@ -158,9 +158,9 @@
$iplot = 1;
$iboth = 0;
- print CSH "run_mt_measure_adj.pl $smodel $lcut $itest $iplot $iboth $tvec $iker $chan $Ts $iparbools $par1 $par2 $par3\n";
+ print CSH "run_mt_measure_adj.pl $smodel $lcut $itest $iplot $iboth $tvec $imeas $chan $Ts $iparbools $par1 $par2 $par3\n";
- #print CSH "run_mt_cc_plot.pl $smodel $lcut $itest $tvec $iker $chan $Ts $iparbools $par1 $par2 $par3\n";
+ #print CSH "run_mt_cc_plot.pl $smodel $lcut $itest $tvec $imeas $chan $Ts $iparbools $par1 $par2 $par3\n";
#------------------------------------
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/write_par_file.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/write_par_file.pl 2009-10-21 22:34:28 UTC (rev 15863)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/write_par_file.pl 2009-10-21 23:40:13 UTC (rev 15864)
@@ -11,27 +11,27 @@
#
# INPUT:
# tstart/dt/ntime
-# iker
+# imeas
# channel (BH or LH)
# TSHORT/TLONG
-# RUN_BANDPASS/DISPLAY_DETAILS/OUTPUT_MEASUREMENT_FILES
+# RUN_BANDPASS/DISPLAY_DETAILS/OUTPUT_MEASUREMENT_FILES/COMPUTE_ADJOINT_SOURCE
# 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 -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
+# write_par_file.pl -0.585/0.011/18200 1 BH 6/30 0/1/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 < 8) {die("Usage: write_par_file.pl xxx\n");}
-($tvec,$iker,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
+($tvec,$imeas,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
# extract variables
($tstart,$dt,$ntime) = split("/",$tvec);
($TSHORT,$TLONG) = split("/",$Ts);
-($ibool1,$ibool2,$ibool3) = split("/",$iparbools);
+($ibool1,$ibool2,$ibool3,$ibool4) = 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."}
@@ -42,26 +42,27 @@
# comments for the PAR_FILE
#$line00 = "output directory"; # cannot have a string here
$line01 = "tstart, DT, npts: time vector for simulations";
-$line02 = "iker (0-8) -- 0 no adj src; 1-8 see manual;";
+$line02 = "imeas (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 = "COMPUTE_ADJOINT_SOURCE";
-$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";
+$line09 = "TSHIFT_MIN; TSHIFT_MAX";
+$line10 = "DLNA_MIN; DLNA_MAX";
+$line11 = "CC_MIN";
+$line12 = "ERROR_TYPE -- 0 none; 1 CC, MT-CC; 2 MT-jack-knife";
+$line13 = "DT_SIGMA_MIN";
+$line14 = "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";
+$line15 = "itaper -- taper type: 1 multi-taper; 2 cosine; 3 boxcar";
+$line16 = "wtr, npi (ntaper = 2*npi)";
+$line17 = "DT_FAC";
+$line18 = "ERR_FAC";
+$line19 = "DT_MAX_SCALE";
+$line20 = "NCYCLE_IN_WINDOW";
#=============================================
$ofile = "./MEASUREMENT.PAR";
@@ -71,31 +72,32 @@
open(OUT,">$ofile");
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("%23i # %s\n",$imeas,$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("%23s # %s\n",$COMPUTE_ADJOINT_SOURCE,$line08);
-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("%12.4f%11.4f # %s\n",$TSHIFT_MIN,$TSHIFT_MAX,$line09);
+print OUT sprintf("%12.4f%11.4f # %s\n",$DLNA_MIN,$DLNA_MAX,$line10);
+print OUT sprintf("%23.3f # %s\n",$CC_MIN,$line11);
+print OUT sprintf("%23i # %s\n",$ERROR_TYPE,$line12);
+print OUT sprintf("%23.3f # %s\n",$DT_SIGMA_MIN,$line13);
+print OUT sprintf("%23.3f # %s\n",$DLNA_SIGMA_MIN,$line14);
-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("%23i # %s\n",$itaper,$line15);
+print OUT sprintf("%17.3f%6.2f # %s\n",$wtr,$npi,$line16);
+print OUT sprintf("%23.3f # %s\n",$DT_FAC,$line17);
+print OUT sprintf("%23.3f # %s\n",$ERR_FAC,$line18);
+print OUT sprintf("%23.3f # %s\n",$DT_MAX_SCALE,$line19);
+print OUT sprintf("%23.3f # %s\n",$NCYCLE_IN_WINDOW,$line20);
#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("%23i # %s\n",$imeas,$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);
More information about the CIG-COMMITS
mailing list