[cig-commits] r21075 - in seismo/3D/ADJOINT_TOMO/measure_adj: . PLOTS USER_MANUAL
liuqy at geodynamics.org
liuqy at geodynamics.org
Tue Nov 27 09:05:45 PST 2012
Author: liuqy
Date: 2012-11-27 09:05:45 -0800 (Tue, 27 Nov 2012)
New Revision: 21075
Added:
seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure.bib
seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/run.latex
Modified:
seismo/3D/ADJOINT_TOMO/measure_adj/Makefile_ifort_caltech
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/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_constants.f90
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/write_par_file.pl
Log:
Mostly cosmetic changes made to main programs and scripts.
Update user manual to include detailed description of the misfit/adjoint
source for various imeas.
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/Makefile_ifort_caltech
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/Makefile_ifort_caltech 2012-11-27 02:07:13 UTC (rev 21074)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/Makefile_ifort_caltech 2012-11-27 17:05:45 UTC (rev 21075)
@@ -2,15 +2,17 @@
F90_FLAGS = -O2 -132
-LIB = -L/opt/seismo-util/lib -lDRWFiles -lf90recipes -lDSacio -lDSacLib -lSacTools -lm
+SACLIBDIR = ${SACHOME}/lib
+LIB = -lsacio -lsac
-MOD = ma_constants ma_variables ma_sub2 ma_sub
+# NOTE: order matters if modules depend on other modules
+MOD = ma_constants ma_variables ascii_rw ma_sub2 ma_sub
SRC_DIR = .
MOD_DIR = mod
OBJ_DIR = obj
BIN_DIR = .
-MAIN = measure_adj
+MAIN = measure_adj rotate_adj_src
MOD_FLAG = module
MOD_OBJ = $(patsubst %,$(OBJ_DIR)/%.o,$(MOD))
@@ -20,7 +22,7 @@
all : measure_adj rotate_adj_src
$(MAIN) : % : $(SRC_DIR)/%.f90 $(F90_OBJ) $(MOD_OBJ)
- $(F90) -o $(BIN_DIR)/$* $(F90_FLAGS) $(SRC_DIR)/$*.f90 -$(MOD_FLAG) $(MOD_DIR) $(OBJ) $(LIB)
+ $(F90) -o $(BIN_DIR)/$* $(F90_FLAGS) $(SRC_DIR)/$*.f90 -$(MOD_FLAG) $(MOD_DIR) $(OBJ) -L${SACLIBDIR} $(LIB)
$(F90_OBJ): $(OBJ_DIR)/%.o : $(SRC_DIR)/%.f90
$(F90) -o $@ $(F90_FLAGS) -c $(SRC_DIR)/$*.f90
@@ -29,7 +31,7 @@
$(F90) -o $@ $(F90_FLAGS) -c $(SRC_DIR)/$*.f90 -$(MOD_FLAG) $(MOD_DIR)
rotate_adj_src: rotate_adj_src.f90
- $(F90) -o rotate_adj_src rotate_adj_src.f90 -L/opt/seismo-util/lib -lDRWFiles
+ $(F90) -o rotate_adj_src rotate_adj_src.f90 -$(MOD_FLAG) $(MOD_DIR) $(OBJ) -L${SACLIBDIR} $(LIB)
.PHONY : clean
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 02:07:13 UTC (rev 21074)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2012-11-27 17:05:45 UTC (rev 21075)
@@ -1,4 +1,4 @@
-#!/usr/bin/perl -w
+#!/usr/bin/perl
#
#-----------------------------------
# plot_win_adj.pl
@@ -28,14 +28,14 @@
-b -- plot MT and CC adjoint sources (=1)
-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
+ -r -- plot 3-component seismograms and adjoint sources in relatively scaled amplitude
-a -- station file to plot all the stations on the map
-d -- data directory
-s -- synthetics directory
-c -- reconstructed directory
- -w -- the file with the windows
- -i -- label for the model iteration
- -j -- Tmin/Tmax
+ -w -- MEASUREMENT.WINDOWS file
+ -i -- label for the model iteration (for title)
+ -j -- Tmin/Tmax (for title)
END
exit(1);
}
@@ -136,7 +136,7 @@
$jpg_file = "$name.jpg";
#$saclst="/opt/seismo-util/bin/saclst";
-$saclst="/usr/local/sac/bin/saclst";
+$saclst="saclst";
#------------
$xmin = -122; $xmax = -114;
@@ -155,14 +155,14 @@
# plot southern California faults
if($ifault==1) {
- $dir0 = "/home/carltape/gmt";
- $fault_file = "$dir0/faults/jennings.xy";
- $kcf_file = "$dir0/faults/kcf.xy";
- $breck_file = "$dir0/faults/breck.xy";
+ $dir0 = "/opt/seismo/datalib/SC/";
+ $fault_file = "$dir0/jennings.xy";
+# $kcf_file = "$dir0/faults/kcf.xy";
+# $breck_file = "$dir0/faults/breck.xy";
$fault_infoK = "-m -W0.5p,0/0/0";
`psxy $fault_file $bounds $proj $fault_infoK -O -K >> $ps_file`;
- `psxy $kcf_file $bounds $proj $fault_infoK -O -K >> $ps_file`;
- `psxy $breck_file $bounds $proj $fault_infoK -O -K >> $ps_file`;
+# `psxy $kcf_file $bounds $proj $fault_infoK -O -K >> $ps_file`;
+# `psxy $breck_file $bounds $proj $fault_infoK -O -K >> $ps_file`;
}
# plot the stations
@@ -282,9 +282,9 @@
$tick2 = "-Ba50f10/a1f1S";
$pen = "0.7p"; # line thickness for seismograms
-$red_pen = "-W${pen},250/0/0";
-$black_pen = "-W${pen},0/0/0";
-$recon_pen = "-W${pen},0/255/255";
+$red_pen = "-W${pen},250/0/0"; # syn
+$black_pen = "-W${pen},0/0/0"; # data
+$recon_pen = "-W${pen},0/255/255"; # recon_syn
my(@Twinb); my(@Twine);
my(@Rwinb); my(@Rwine);
@@ -455,6 +455,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";
`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) {
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 2012-11-27 02:07:13 UTC (rev 21074)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl 2012-11-27 17:05:45 UTC (rev 21075)
@@ -16,7 +16,7 @@
sub Usage{
print STDERR <<END;
- 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
+ plot_win_adj_all.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);
Added: seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure.bib
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure.bib (rev 0)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure.bib 2012-11-27 17:05:45 UTC (rev 21075)
@@ -0,0 +1,123 @@
+Automatically generated by Mendeley 1.6
+Any changes to this file will be lost if it is regenerated by Mendeley.
+
+ at article{YZhou2005,
+author = {Zhou, Y. and Dahlen, F. A. and Nolet, Guust and Laske, Gabi},
+doi = {10.1111/j.1365-246X.2005.02780.x},
+file = {:data2/data3/Research\_Papers/Zhou05.pdf:pdf;:data2/data3/Research\_Papers/Ying06.pdf:pdf},
+journal = {Geophys. J. Int.},
+keywords = {1 i n t,Frechet derivatives,Surface-wave inversions,echet derivatives,fr,global seismology,global seismology of the upper mantle,has been pro-,of the upper mantle,ro d u c,sensitivity,sensitivity kernels,surface waves,t i o n,the large-scale seismic structure,tomography,upper-mantle tomogrpahy},
+mendeley-tags = {Surface-wave inversions,sensitivity kernels,upper-mantle tomogrpahy},
+number = {3},
+pages = {1087--1111},
+publisher = {John Wiley $\backslash$\& Sons},
+title = {{Finite-frequency effects in global surface-wave tomography}},
+url = {http://scholar.google.com/scholar?hl=en\&btnG=Search\&q=intitle:Finite-frequency+effects+in+global+surface-wave+tomography\#0 http://www3.interscience.wiley.com/journal/118661525/abstract},
+volume = {163},
+year = {2005}
+}
+ at phdthesis{Tape2009phd,
+author = {Tape, C. H.},
+booktitle = {Area},
+file = {:data2/data3/Research\_Papers/CarlTapePhDThesis.pdf:pdf},
+keywords = {adjoint methods,adjoint tomography,southern California Tomogrpahy},
+mendeley-tags = {adjoint methods,adjoint tomography,southern California Tomogrpahy},
+school = {Caltech},
+title = {{Seismic Tomography of Southern California Using Adjoint Methods}},
+url = {http://thesis.library.caltech.edu/1656/},
+volume = {2009},
+year = {2009}
+}
+ at article{Tape2007,
+author = {Tape, C. H. and Liu, Q. and Tromp, J.},
+file = {:data2/data3/Research\_Papers/Tape07.pdf:pdf},
+journal = {Geophys. J. Int.},
+keywords = {adjoint methods,finite-frequency tomography},
+mendeley-tags = {adjoint methods,finite-frequency tomography},
+pages = {1105--29},
+title = {{Finite-frequency tomography using adjoint methods - Methodology and examples using membrane surface waves}},
+volume = {168},
+year = {2007}
+}
+ at article{Tromp2005,
+author = {Tromp, J. and Tape, C. H. and Liu, Q.},
+file = {:data2/data3/Research\_Papers/Tromp05.pdf:pdf},
+journal = {Geophys. J. Int.},
+keywords = {adjoint methods,banana-doughnut kernels,seismic tomography,time-reversal},
+mendeley-tags = {adjoint methods,banana-doughnut kernels,seismic tomography,time-reversal},
+pages = {195--216},
+publisher = {[Oxford]: Published for the Royal Astronomical Society, the Deutsche Geophysikalische Gesellschaft, and the European Geophysical Society by Blackwell Scientific Publications, c1989-},
+title = {{Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels}},
+url = {http://scholar.google.com/scholar?hl=en\&btnG=Search\&q=intitle:Seismic+tomography,+adjoint+methods,+time+reversal+and+banana-doughnut+kernels\#0},
+volume = {160},
+year = {2005}
+}
+ at article{Tape2010,
+author = {Tape, C. H. and Liu, Q. and Maggi, A. and Tromp, J.},
+doi = {10.1111/j.1365-246X.2009.04429.x},
+file = {:data2/data3/Research\_Papers/TapeLiu10.pdf:pdf},
+issn = {0956540X},
+journal = {Geophys. J. Int.},
+keywords = {adjoint tomography,between simulated,body waves,computational seismology,crustal structure,inverse theory,or,oscillations,seismic,seismic tomography uses measurements,southern California tomography,surface waves and free oscillations,tomography},
+mendeley-tags = {adjoint tomography,southern California tomography},
+pages = {433--462},
+title = {{Seismic tomography of the southern California crust based on spectral-element and adjoint methods}},
+url = {http://blackwell-synergy.com/doi/abs/10.1111/j.1365-246X.2009.04429.x},
+volume = {180},
+year = {2010}
+}
+ at article{Maggi2009,
+author = {Maggi, A. and Tape, C. H. and Chen, M. and Chao, D. and Tromp, J.},
+doi = {10.1111/j.1365-246X.2009.04099.x},
+file = {:data2/data3/Research\_Papers/MaggiTape09.pdf:pdf},
+issn = {0956540X},
+journal = {Geophys. J. Int.},
+keywords = {3d,a 3-d reference,body waves,of,opened up the possibility,oscillations,seismic tomography,seismic tomography based upon,surface waves and free,that is,the numeric kernels have,time series analysis,tomography},
+month = jul,
+pages = {257--281},
+title = {{An automated time-window selection algorithm for seismic tomography}},
+url = {http://blackwell-synergy.com/doi/abs/10.1111/j.1365-246X.2009.04099.x},
+volume = {178},
+year = {2009}
+}
+ at article{Tape2009,
+abstract = {Using an inversion strategy based on adjoint methods, we developed a three-dimensional seismological model of the southern California crust. The resulting model involved 16 tomographic iterations, which required 6800 wavefield simulations and a total of 0.8 million central processing unit hours. The new crustal model reveals strong heterogeneity, including local changes of +/-30\% with respect to the initial three-dimensional model provided by the Southern California Earthquake Center. The model illuminates shallow features such as sedimentary basins and compositional contrasts across faults. It also reveals crustal features at depth that aid in the tectonic reconstruction of southern California, such as subduction-captured oceanic crustal fragments. The new model enables more realistic and accurate assessments of seismic hazard.},
+author = {Tape, C. H. and Liu, Q. and Maggi, Alessia and Tromp, J.},
+doi = {10.1126/science.1175298},
+file = {:data2/data3/Research\_Papers/TapeLiuMaggiTromp09.pdf:pdf;:data2/data3/Research\_Papers/TapeLiuMaggiTromp09.pdf:pdf;:data2/data3/Research\_Papers/Tape\_SOM.pdf:pdf},
+issn = {1095-9203},
+journal = {Science},
+keywords = {adjoint tomography,southern california velocity models},
+mendeley-tags = {adjoint tomography,southern california velocity models},
+month = aug,
+number = {5943},
+pages = {988--992},
+pmid = {19696349},
+title = {{Adjoint tomography of the southern California crust}},
+url = {http://www.ncbi.nlm.nih.gov/pubmed/19696349},
+volume = {325},
+year = {2009}
+}
+ at article{LiuTromp2006,
+author = {Liu, Q. and Tromp, J.},
+file = {:data2/data3/Research\_Papers/LiuTr06a.pdf:pdf},
+journal = {Bull. Seism. Soc. Am.},
+number = {6},
+pages = {2283--2397},
+title = {{Finite-frequency kernels based on adjoint methods}},
+url = {http://www.bssaonline.org/cgi/content/abstract/94/1/187},
+volume = {96},
+year = {2006}
+}
+ at article{YZhou2004,
+author = {Zhou, Y. and Nolet, G. and Dahlen, F. A.},
+file = {:data2/data3/Research\_Papers/Zhou\_2004.pdf:pdf},
+journal = {Geophys. J. Int.},
+keywords = {Surface-wave inversions,measurements,sensitivity kernels},
+mendeley-tags = {Surface-wave inversions,measurements,sensitivity kernels},
+pages = {142--168},
+title = {{Three-dimensional sensitivity kernels for surface wave observables}},
+url = {http://scholar.google.com/scholar?hl=en\&btnG=Search\&q=intitle:Three-dimensional+sensitivity+kernels+for+surface+wave+observables,+Geophys\#0},
+volume = {158},
+year = {2004}
+}
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 02:07:13 UTC (rev 21074)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/measure_adj_manual.tex 2012-11-27 17:05:45 UTC (rev 21075)
@@ -1,5 +1,4 @@
\documentclass[11pt,titlepage,fleqn]{article}
-
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{latexsym}
@@ -8,6 +7,7 @@
%\usepackage{epsfig}
\usepackage{graphicx}
\usepackage{bm}
+\usepackage[usenames]{color}
%\usepackage{fancybox} # do not use when showing Table of Contents
%\usepackage{fancyhdr}
@@ -35,7 +35,10 @@
\footskip 30pt
\hoffset 0pt
-\include{NEWCOMMANDS}
+\include{newcommands}
+\newcommand{\eg}{e.g.,}
+\newcommand{\refSec}[1]{section \ref{#1}}
+\newcommand{\vala}{Hj\"{o}rleifsd\'{o}ttir}
\graphicspath{
{./figures/}
@@ -46,14 +49,14 @@
%-----------------------------------------------------------
\begin{center}
-{\bf \Large MEAS--ADJ: A software package for making measurements and computing adjoint sources for seismic tomography}
+{\bf \Large \verb+measure_adj+: A software package for making measurements and computing adjoint sources for seismic 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 Author of Manual: Carl Tape, Qinya Liu
\noindent \today
@@ -75,9 +78,9 @@
{\bf NOTE: This is a preliminary user manual with many modifications 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,LiuTromp2006,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.
+The open-source software package \verb+measure_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,LiuTromp2006,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.
+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 \verb+measure_adj+. However, FLEXWIN is not necessary to run \verb+measure_adj+.
There are two additional sets of notes on multitaper measurements and adjoint sources, \verb+multitaper_notes.pdf+ and \verb+multitaper_vala.pdf+. Some sections of \verb+multitaper_notes.pdf+ were published as ``Multitaper measurements for adjoint tomography'' in \citet[][Appendix~C]{Tape2009phd}.
@@ -91,7 +94,7 @@
\subsection{System requirements}
-In order to install and run, MEAS--ADJ requires:
+In order to install and run, \verb+measure_adj+ program requires:
%
\begin{itemize}
\item UNIX operating system (Linux, Solaris, MacOS \ldots)
@@ -110,7 +113,7 @@
(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.
+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 \verb+measure_adj+ on large datasets.
%-------------
@@ -134,7 +137,7 @@
After compiling, execute the following two commands:
%
\begin{verbatim}
-mt_measure_adj
+measure_adj
csh -f run_mt_measure_adj.csh
\end{verbatim}
%
@@ -144,6 +147,17 @@
PLOTS/9818433_T006_T030_MPM_CI_m16_mtm_win_adj.pdf
\end{verbatim}
+A set of auxilliary programs/files are required for the successful run of
+post-processing and plotting scripts, including
+\begin{itemize}
+ \item Southern California faults xy files in \verb+PLOTS/plot_win_adj.pl+:
+\begin{verbatim}
+ $dir0 = "/opt/seismo/datalib/SC/";
+ $fault_file = "$dir0/jennings.xy";
+\end{verbatim}
+\item Dr. Lupei Zhu's \verb+saclst+ program
+\item Dr. Lupei Zhu's \verb+pssac2+ plotting program (\verb+pssac+ may work as well).
+\end{itemize}
%-------------------------------------------------------------------
\pagebreak
@@ -174,6 +188,8 @@
\end{verbatim}
%
The first line contains the number of pairs of records to be read in. Each pair of files is followed by the number of windows within which measurements will be made, followed by the time intervals for each window.
+Note that data and synthetics require to have exactly the same \verb+t0, dt, npts+ for the code to run through. Therefore
+no interpolation of \verb+dt+ is really necessary at any point.
%-------------
@@ -183,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: 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
@@ -192,12 +208,12 @@
-4.5000 4.5000 # TSHIFT_MIN; TSHIFT_MAX
-1.5000 1.5000 # DLNA_MIN; DLNA_MAX
0.690 # CC_MIN
- 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
+ 1 # ERROR_TYPE -- 0 none; 1 CC, MT-CC; 2 MT-jack-knife (all for MTM!)
+ 1.000 # DT_SIGMA_MIN: water level of traveltime uncertainty
+ 0.500 # DLNA_SIGMA_MIN: water level of dlnA uncertainty
+ 1 # ITAPER:1 multi-taper; 2 cosine; 3 boxcar taper type for freq-dep. meas.
0.020 2.50 # WTR, NPI (ntaper = 2*NPI)
- 2.000 # DT_FAC
+ 2.000 # DT_FAC (following 4 pars are for MTM --> CC check)
2.500 # ERR_FAC
3.500 # DT_MAX_SCALE
1.500 # NCYCLE_IN_WINDOW
@@ -224,7 +240,8 @@
\item Limits on the cross-correlation measurement: \verb+CC_MIN+. This can optionally be read in from the \verb+PAR_FILE+ used in the FLEXWIN algorithm.
-\item \verb+ERROR_TYPE+ = 0--2: Type of estimated errors associated with each measurement (and adjoint source). No error (0), cross-correlation error estimate (1), multitaper jack-knife error estimate (2). The cross-correlation error estimate is presented in \citet[][Appendix~A]{Tape2010}.
+\item \verb+ERROR_TYPE+ = 0--2: Type of estimated errors associated with \blue{multi-taper} measurement (and adjoint source). No error (0), cross-correlation error estimate (1), multitaper jack-knife error estimate (2). The cross-correlation error estimate is presented in \citet[][Appendix~A]{Tape2010}. For imeas=5,6 (cc), \verb+sigma_dt/dlnA_cc+
+is automatically used. For \verb+ERROR_TYPE = 0+, all \verb+sigma_tau+ and \verb+sigma_dlnA+ are set to be $1$ for multitaper adjoint sources.
\item \verb+DT_SIGMA_MIN+: water-level minimum cross-correlation traveltime difference uncertainty estimate.
@@ -232,7 +249,9 @@
\item Multitaper parameter, \verb+ITAPER+: Type of taper to use in constructing the transfer function between synthetics and data. Taper options are multitaper (1), cosine taper (2), and boxcar taper (3). For the single-taper options (2--3) the transfer function is not used, as the adjoint sources are constructed directly from the synthetic seismograms. For the multitaper option, the number of tapers is fixed to be twice \verb+NPI+ (see \verb+multitaper_vala.pdf+).
-\item Multitaper parameters, \verb+WTR+, \verb+NPI+ (see \verb+multitaper_vala.pdf+).
+\item Multitaper parameters, \verb+WTR+, \verb+NPI+ (see \verb+multitaper_vala.pdf+). Note \verb+WTR+ is also used to
+determine \verb+i_right+ corresponding to the maximum frequency of valid frequency-dependent measurements. It requires
+that the power at this frequency is above $10$ times the \verb+WTR*max_syn_power+.
\item Multitaper parameter, \verb+DT_FAC+ (see \refSec{sec:MTparm}).
@@ -248,19 +267,31 @@
\subsubsection{Multitaper parameters}
\label{sec:MTparm}
-The selection of the multitaper parameters \verb+DT_FAC+, \verb+ERR_FAC+, \verb+DT_MAX_SCALE+, and \verb+NCYCLE_IN_WINDOW+ are not easy to determine for each particular dataset. These parameters determine whether a multitaper measurement is reasonable enough to retain (otherwise revert to a cross-correlation measurement). The key subroutine is \verb+mt_measure_select.f90+.
+The selection of the multitaper parameters \verb+DT_FAC+, \verb+ERR_FAC+, \verb+DT_MAX_SCALE+, and \verb+NCYCLE_IN_WINDOW+ are not easy for each particular dataset. These parameters determine whether a multitaper measurement is reasonable enough to retain (otherwise revert to a cross-correlation measurement). The key subroutine is \verb+mt_measure_select.f90+.
-This information used to be specified in \verb+mt_constants.f90+, and two example sets of parameters are:
+MTMs are rejected (\verb+ user_trace = .false.+) based on:
+\begin{itemize}
+\item number of cycles in the window
+\item number of frequency points are given within \verb+[fstart,fend]+. Even when \verb+RUN_BANDPASS+ is turned off, the \verb+TLONG+ and \verb+TSHORT+ are still converted to \verb+fstart,fend+. They are best set to the values used in \verb+FLEXWIN+
+\item \verb+tshift_cc <= dt+: too small a time shift
+\item \verb+ abs(dtau_w(j)) > Tvec(j)/DT_FAC+: $\tau$ at a specific frequency is too high
+\item \verb+ err_dt(j) > Tvec(j)/ERR_FAC+: error for a specific frequency is too high
+\item \verb+ abs(dtau_w(j)) > DT_MAX_SCALE*abs(tshift))+ $\tau$ at specific frequency deviates too much for global \verb+tshift_cc+
+\end{itemize}
+Note that the cross-correlation measurements are then run through \verb+cc_measure_select()+ to determine if they are
+usable or not (\verb+tshift, dlnA = 0+).
+
+Two example sets of parameters are given
%
\begin{verbatim}
! ********* adjust following for global traces ************************
! ratio of current period with respect to dt and err_dt measurements
- double precision, parameter :: DT_FAC = 1.0
+ double precision, parameter :: DT_FAC = 1.0
double precision, parameter :: ERR_FAC = 8.0
! max time shift allowed at all freq should be DT_MAX_SCALE * Tshift_xc
double precision, parameter :: DT_MAX_SCALE = 5.0
- ! use 3 cycles for surfaces waves
- double precision, parameter :: NCYCLE_IN_WINDOW = 2.0
+ ! few cycles for surfaces waves
+ double precision, parameter :: NCYCLE_IN_WINDOW = 2.0 ! window length > 2.*50
!***********************************************************************
! ********* adjust following for socal 3-30s traces ********************
@@ -277,15 +308,59 @@
%-------------
\subsection{Output files}
-
-With \verb+OUTPUT_MEASUREMENT_FILES = .true.+ in \verb+MEASUREMENT.PAR+, you should find numerous output files in \verb+OUTPUT_FILES+.
-{\bf NOTE: MORE HERE.}
-
-There are two output files in the local run directory:
+There are two output files related to measurements in the local run directory:
%
\begin{enumerate}
+
+% unit 13
\item \verb+window_chi+, a comprehensive output file of misfit values, with one line per window.
-\item \verb+window_index+, an abbreviated output file with the indexing for each window. The columns are \verb+net, sta, chan_syn, chan_dat, nwin, ipair, j, tstart, tend+, where \verb+nwin+ is the overall window counter, \verb+ipair+ is the seismogram (pair) counter, and \verb+j+ is the local window counter. For example, for the test case run:
+\begin{verbatim}
+! KEY: write misfit function values to file (for each window)
+! Here are some columns of values in window_chi (add 8 for the actual column number)
+! 1: MT-TT chi (imeas=7), 2: MT-dlnA chi (imeas=8)
+! 3: XC-TT chi (imeas=5), 4: XC-dlnA chi (imeas=6)
+! 5: MT-TT meas (freq-average), 6: MT-dlnA meas, 7: XC-TT meas, 8: XC-dlnA meas
+! 9: MT-TT error, 10: MT-dlnA error, 11: XC-TT error, 12: XC-dlnA error
+! 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
+! WINDOW : 21: tr_chi, 22: am_chi
+! Example of a reduced file: awk '{print $2,$3,$4,$5,$6,$29}' window_chi > window_tr_chi
+
+write(13,'(a14,a8,a3,a5,i4,i4,2e14.6,20e14.6,2e14.6,2f14.6)') &
+ file_prefix0,sta,net,chan_syn,j,imeas,&
+ tstart,tend,window_chi(:),tr_chi,am_chi,T_pmax_dat,T_pmax_syn
+
+! misfit function value
+if(is_mtm==1) window_chi(1) = 0.5 * 2.0 * df * sum( (dtau_w(1:i_right))**2 * wp_taper(1:i_right) ) ! tr_chi
+if(is_mtm==1) window_chi(2) = 0.5 * 2.0 * df * sum( (dlnA_w(1:i_right))**2 * wq_taper(1:i_right) ) ! amp_chi
+window_chi(3) = 0.5 * (tshift/sigma_dt_cc)**2 ! tr_chi
+window_chi(4) = 0.5 * (dlnA/sigma_dlnA_cc)**2 ! amp_chi
+! cc/averaged-mt tshift/dlnA measurement
+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 measurement uncertainties
+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) = 0.5 * sum(dat_dtw(:)**2)
+window_chi(14) = 0.5 * sum(syn_dtw(:)**2)
+window_chi(15) = 0.5 * sum((dat_dtw-syn_dtw)**2) ! tr/amp_chi for imeas=1/2
+window_chi(16) = nlen*dt
+window_chi(17) = 0.5 * sum( data**2 ) ! power of entire trace
+window_chi(18) = 0.5 * sum( syn**2 )
+window_chi(19) = 0.5 * sum( (data-syn)**2 )
+window_chi(20) = npts*dt
+\end{verbatim}
+% unit 12
+\item \verb+window_index+, an abbreviated output file with the indexing for each window. The columns are
+\begin{verbatim}
+write(12,'(a3,a8,a5,a5,3i5,2f12.3)') net,sta,chan_syn,chan_dat,nwin,ipair,j,tstart,tend
+\end{verbatim}
+where \verb+nwin+ is the overall window counter, \verb+ipair+ is the seismogram (pair) counter, and \verb+j+ is the local window counter. For example, for the test case run:
%
\begin{verbatim}
CI MPM BHR BHR 1 1 1 11.277 58.477
@@ -297,6 +372,46 @@
%
\end{enumerate}
+With \verb+COMPUTE_ADJOINT_SOURCE = .true.+,
+\begin{enumerate}
+\item adjoint source files (e.g., \verb+MPM.CI.BHZ.iker07.adj+) will appear
+in \verb+OUTPUT_FILES/+. These are ascii files with the time colume defined by \verb+tstart,DT,npts+ in the parameter file.
+(if \verb+DO_RAY_DENSITY_SOURCE = .true.+, adjoint file names are \verb+MPM.CI.BHZ.iker7.density.adj+).
+
+\item \verb+window_chi_sum+ in local directory, the sum of misfit values (CC or MTM) of all windows, with weights taken into account if \verb+DO_WEIGHTING = .true.+. This scalar value literally corresponds to all the adjoint sources produced in \verb+OUTPUT_FILES/+. If input includes all windows for a particular event, this scalar value is the event misfit.
+More comprehensive sum of misfit can be found in the script
+\begin{verbatim}
+/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
+\end{verbatim}
+% added by LQY
+\end{enumerate}
+
+With \verb+OUTPUT_MEASUREMENT_FILES = .true.+ in \verb+MEASUREMENT.PAR+, you should find numerous output files in \verb+OUTPUT_FILES/+.
+Other files include (\verb+prefix=MPM.CI.BHZ.01.mtm+, \verb+prefix0=MPM.CI.BHZ+)
+\begin{verbatim}
+prefix.recon_syn_cc.sac # syn_dtw_cc
+prefix.recon_syn_cc_dt.sac # syn_dtw_cc_dt
+prefix.dt/dlnA_average/cc # dtau/dlnA_meas, dtaul/dlnA_sigma ('average' means 'mtm')
+prefix.ph/abs/dlnA/ph_cor/dt # transfer fun.: phi/abs/log(abs)/phi(corr)/dtau(idf_new:idf_new:i_right)
+ # tshift/dlnA from cc are added to dtau/dlnA_wt
+prefix.recon_syn.sac # recon. syn with both dtau(om) and dlnA(om) applied for all imeas >=3 cases.
+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
+\end{verbatim}
+
+IF \verb+DISPLAY_MORE_DETAILS = .true.+ in the parameter file, following files will be written for every window pair
+\begin{verbatim}
+prefix.obs (original data/syn)
+prefix.syn
+prefix.dat.sac (time windowed data/syn)
+prefix.syn.sac
+prefix.obs.power % abs(dat_dtwo(1:i_right)), i.e., dat_dtw_cc
+prefix.syn.power % abs(syn_dtwo(1:i_right)), i.e., syn_dtw
+\end{verbatim}
+
+
%-------------
\subsection{Scripts}
@@ -304,22 +419,63 @@
The three scripts in the run directory are
%
\begin{enumerate}
-\item \verb+write_par_file.pl+. Write the file \verb+MEASUREMENT.PAR+ from a set of parameters. This is called within the run scripts in \verb+scripts_tomo+.
+\item \verb+write_par_file.pl+. Write the file \verb+MEASUREMENT.PAR+ from a given set of parameters. This is called within the run scripts in \verb+scripts_tomo+.
-\item \verb+prepare_adj_src.pl+. This script reads in a set of adjoint sources made from Z-R-T records and outputs a set of adjoint sources in Z-E-N that can be used in SPECFEM3D. 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).
+\item \verb+prepare_adj_src.pl+. This script reads in a set of adjoint sources made from Z-R-T records and outputs a set of adjoint sources in Z-E-N that can be used in SPECFEM3D.
+\begin{verbatim}
+prepare_adj_src.pl -m CMT -z BH -s STATION -o OUTDIR -i [rotation] all_adj_ZRT_files
+prepare_adj_src.pl -m CMTSOLUTION_9818433 -s PLOTS/STATIONS_TOMO -o ADJOINT_SOURCES \
+ OUTPUT_FILES/*adj
+\end{verbatim}
+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}
+rotate_adj_src baz(radian!) zfile tfile rfile efile nfile
+\end{verbatim}
+The perl script also outputs \verb+STATIONS_ADJOINT+ file for all the stations with adjoint sources.
+%{\bf NOTE: We need to eliminate the dependence on the library /opt/seismo-util/lib/perl.} --- LQY: done
-{\bf NOTE: We need to eliminate the dependence on the library
-
-/opt/seismo-util/lib/perl.}
-
-\item \verb+combine_adj_src.pl+. This script combines two sets of adjoint sources (for example, two different bandpassed versions), and outputs the new adjoint sources, along with a new STATIONS file to be used.
+\item \verb+combine_adj_src.pl+. This script combines two sets of adjoint sources (for example, two different bandpassed versions), and outputs the new adjoint sources, along with a new \verb+STATIONS_ADJOINT+ file to be used.
+\begin{verbatim}
+Usage: combine_adj_src.pl DIR_1 DIR_2 DIR_NEW imeas1 imeas2 chan
+\end{verbatim}
+Adjoint sources from \verb+dir1+ and \verb+dir2+ with the same station/net/component names are added up; station files
+from \verb+dir1/STATIONS_ADJOINT+ and \verb+dir2/STATIONS_ADJOINT+ are combined into \verb+dir_new/STATIONS_ADJOINT+.
\end{enumerate}
%
-There are three other directories with scripts:
+
+Scripts in the \verb+PLOT/+ directory:
%
\begin{itemize}
\item The primary plotting script is \verb+PLOTS/plot_win_adj.pl+.
+\begin{verbatim}
+plot_win_adj_all.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
+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
+\end{verbatim}
+which uses another script
+\begin{verbatim}
+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
+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
+\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+).
+\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]
+plot_win_stats.pl Tmin/Tmax eid name meas_file sta_file eid_text cmtall_psmeca
+\end{verbatim}
+
+\item organize \verb+plot_win_adj_all.pl+ output by station or event (???)
+\begin{verbatim}
+make_pdf_by_event.pl
+make_pdf_by_station.pl
+\end{verbatim}
+
\item For large datasets, we have included an additional set of scripts in the directory \verb+scripts_tomo+. Each user will have to make some modifications to these scripts.
\item For plotting output files for individual measurements, we have included a set of scripts in \verb+scripts_meas+. These scripts have not been extensively tested.
@@ -331,19 +487,157 @@
\section{Measurement options}
\label{sec:meas}
-There are several choices of measurements, and each choice leads to a different adjoint source, as illustrated in \citet{Tromp2005}. The user must specify a value of \verb+imeas+ in \verb+MEASUREMENT.PAR+ with the following options:
+There are several choices of measurements (or definitions of misfit functions), and each choice leads to a different adjoint source, as illustrated in \citet{Tromp2005}. The user must specify a value of \verb+imeas+ in \verb+MEASUREMENT.PAR+ with the following options:
%
\begin{enumerate}
\item \verb+imeas = 1+, normalized waveform difference. Adjoint source is constructed from the {\em data only}, with the form $-d(t)/\|d(t)\|^2$.
+\eqa
+\phi &=& \frac{1}{2}\frac{\int [d(t)-s(t)]^2\,dt }{\int d^2(t)\,dt} \quad \text{at}\quad s(t)=0 \nn\\
+f^\dagger(t) &=& -\frac{d(t)}{\|d(t)\|^2} \nn
+\ena
+If \verb+NO_WAVEFORM_DIFFERENCE = .true.+, this becomes (why???)
+\eq
+f^\dagger(t) = d(t) \nn
+\en
+
\item \verb+imeas = 2+, waveform difference, $s(t) - d(t)$.
-\item \verb+imeas = 3+, cross-correlation traveltime difference for a (banana-doughtnut) sensitivity kernel. The measurement between data and synthetics {\em is not used} in constructing the adjoint source.
-\item \verb+imeas = 4+, amplitude difference for a (banana-doughtnut) sensitivity kernel. The measurement between data and synthetics {\em is not used} in constructing the adjoint source.
-\item \verb+imeas = 5+, cross-correlation traveltime difference for an event kernel. The measurement between data and synthetics {\em is used} in constructing the adjoint source.
-\item \verb+imeas = 6+, amplitude difference for an event kernel. The measurement between data and synthetics {\em is used} in constructing the adjoint source.
+\eqa
+\phi &=& \frac{1}{2}\int [d(t)-s(t)]^2\,dt \nn\\
+f^\dagger(t) &=& s(t)-d(t) \nn
+\ena
+
+\item \verb+imeas = 3+, cross-correlation traveltime for a (banana-doughtnut) sensitivity kernel. The measurement between data and synthetics {\em is not used} in constructing the adjoint source.
+\eqa
+\phi &=& T_{syn} \nn\\
+\delta \phi &=& \delta{T_{syn}} = - \frac{\int w(t) \dot{s}(t) \delta s(t)\,dt}{\int w(t) \dot{s}^2(t)\,dt}
+\quad\quad\quad\text{TTL (43)}\nn \\
+f^\dagger(t) &=& -\frac{w(t)\dot{s}(t)}{\int w(t) \dot{s}^2(t)\,dt} \nn
+\ena
+where $w(t)$ is the time window over which the measurement is made.
+\red{Note the $-$ sign!.}
+
+\item \verb+imeas = 4+, amplitude difference for a (banana-doughtnut) sensitivity kernel. The measurement between data and synthetics {\em is not used} in constructing the adjoint source.
+\eqa
+\phi &=& \ln A_{syn} \nn\\
+\delta \phi &=& \delta ln A_{syn} = \frac{\int w(t)s(t)\delta s(t)\,dt}{\int w(t) s^2(t)\,dt}
+\quad\quad\quad\text{TTL (64)}
+\nn \\
+f^\dagger(t) &=& \frac{w(t) s(t)}{\int w(t) s^2(t)\,dt} \nn
+\ena
+
+\item \verb+imeas = 5+, cross-correlation traveltime difference for an event kernel.
+\eqa
+\phi &=& \frac{1}{2}\sum_{rp}\left[ \frac{T_{obs}-T_{syn}}{\sigma_{\Delta T}} \right]^2 \nn\\
+\delta \phi &=& - \sum_{rp} \frac{T_{obs}-T_{syn}}{\sigma^2_{\Delta T}}\, \delta{T_{syn}} = \sum_{rp} \frac{ \Delta T_{syn}}{\sigma^2_{\Delta T}}
+ \frac{\int w(t) \dot{s}(t) \delta s(t)\,dt}{\int w(t) \dot{s}^2(t)\,dt} \nn \\
+f^\dagger(t) &=& \sum_{rp} \frac{w(t)\dot{s}(t)}{\int w(t) \dot{s}^2(t)\,dt} \frac{ \Delta T_{syn}}{\sigma^2_{\Delta T}} \nn
+\ena
+where the traveltime \textit{delay} of observed data w.r.t. synthetics $\Delta T_{syn}= T_{obs}-T_{syn}$ {\em is used} in constructing the adjoint source. The summation is performed over receivers ($r$) and phases ($p$).
+
+Note the cross-correlation related measurements, including \verb+tshift, dlnA, cc_max+
+given by \verb+compute_syn_cc()+ correspond to
+\eq
+cc_i = \frac{\sum_j s_j \red{d_{j+i}}}{\sum_j s_j s_j},\quad \Delta T = i_{max}*dt,
+\quad
+\Delta \ln A = \frac{1}{2} \ln \left[\frac{\int d^2(t) dt}{\int s^2(t) dt}\right],
+\en
+and average error estimates are given by
+\eq
+\sigma_{\Delta T} = \sqrt{\frac{\int[d(t)-s_c(t)]^2\,dt}{\int \dot{s}_c(t)\,dt}},\quad
+\sigma_{\Delta lnA} = \sqrt{\frac{\int[d(t)-s_c(t)]^2\,dt}{\dot{s}_c(t)\,dt}}
+\en
+where $s_c(t)$ is the reconstructed synthetics after applying $\Delta T$ and $\Delta \ln A$ corrections.
+
+\item \verb+imeas = 6+, amplitude difference for an event kernel. As
+\eq
+d(t) \sim \exp(\Delta \ln A)\, s(t) \sim (1+\Delta \ln A_{syn})\, s(t) \nn,
+\en
+the misfit of amplitude anomaly may be defined by
+\eqa
+\phi &=& \frac{1}{2}\sum_{rp}\left[\frac{\ln A_{obs} - \ln A_{syn}}{\sigma_{\Delta lnA}}\right]^2
+= \sum_{rp}\frac{1}{2}\left[\frac{\Delta \ln A}{\sigma_{\Delta lnA}}\right]^2 \nn\\
+\delta \phi &=& -\sum_{rp}\frac{\Delta \ln A}{\sigma^2_{\Delta lnA}}\,\delta ln A = -\sum_{rp} \frac{\Delta \ln A}{\sigma^2_{\Delta lnA}} \frac{\int w(t)s(t)\delta s(t)\,dt}{\int w(t) s^2(t)\,dt} \nn \\
+f^\dagger(t) &=& -\sum_{rp}\frac{w(t) s(t)}{\int w(t) s^2(t)\,dt} \frac{\Delta \ln A}{\sigma^2_{\Delta lnA}}\nn
+\ena
+where the measurement amplitude anomaly between data and synthetics $\Delta \ln A = \ln A_{obs}-\ln A_{syn} \sim A_{obs}/A_{syn} -1$ {\em is used} in constructing the adjoint source.
+
\item \verb+imeas = 7+, multitaper traveltime difference for an event kernel. The measurement between data and synthetics {\em is used} in constructing the adjoint source. See \verb+multitaper_notes.pdf+.
-\item \verb+imeas = 8+, multitaper ampltidue difference for an event kernel. The measurement between data and synthetics {\em is used} in constructing the adjoint source. See \verb+multitaper_notes.pdf+.
+\eqa
+&&\phi_P(\bbm)= \frac{1}{2}\sum_{rp} \int W_{P_{rp}}(\om) \left[\tau_{rp}^{obs}(\om) - \tau_{rp}(\om,\bbm)\right]^2\, d\om \nn \\
+&&\delta \phi_P = -\sum_{rp}\int W_{P_{rp}}(\om) \Delta \tau_{rp}(\om, \bbm) \delta \tau_{rp}(\om,\bbm)\,d\om \nn \\
+&&f^\dagger_P(t) = \sum_{rp}\sum_j h_j(t) P_j(t), \quad P_j(\om)=W_{P_{rp}}(\om) \Delta \tau_{rp}(\om) p_j(\om), \nn \\
+&& \quad\quad\quad p_j(\om)=\frac{i\om s_j(\om)}{\sum_k |(i\om) s_k(\om)|^2}, \quad s_j(t)=s(t,\bbm)h_j(t)\nn
+\ena
+where $\Delta \tau_{rp}(\om,\bbm)=\tau^{obs}_{rp}(\om)-\tau_{rp}(\om,\bbm)$ is the frequency-dependent traveltime (\textit{delay}) measurements,
+\eq
+W_{P_{rp}}(\om)=\frac{W_{rp}(\om)}{\sigma_{P_{rp}}^2(\om)\int\,W_{rp}(\om)d\om}\nn
+\en
+is the frequency domain taper scaled by error estimates, and $h_j(t)$ is the $j$'th time-domain single taper applied to synthetics.
+The frequency domain taper between $[f_1,f_2]$ may be defined by
+\eq
+W_{rp}(f)=1-\left[\cos\frac{\pi(f-f_1)}{f_2-f_1}\right]^\gamma, \quad \gamma = 10
+\en
+
+\item \verb+imeas = 8+, multitaper amplitude difference for an event kernel. The measurement between data and synthetics {\em is used} in constructing the adjoint source. See \verb+multitaper_notes.pdf+.
+\eqa
+\phi_Q(\bbm)&=& \frac{1}{2}\sum_{rp} \int W_{Q_{rp}}(\om) \left[ \ln A_{rp}^{obs}(\om) - \ln A_{rp}(\om,\bbm)\right]^2\, d\om \nn \\
+\delta \phi_Q &=& -\sum_{rp}\int W_{Q_{rp}}(\om) \Delta \ln A_{rp}(\om, \bbm) \delta \ln A_{rp}(\om,\bbm)\,d\om \nn \\
+f^\dagger_Q(t) &=& \sum_{rp}\sum_j h_j(t) Q_j(t), \quad Q_j(\om)=W_{Q_{rp}}(\om) \Delta \ln A_{Q_{rp}}(\om) q_j(\om), \nn\\
+&&\quad\quad q_j(\om)= -\frac{s_j}{\sum_k |s_k|^2}= i\om p_j(\om)\nn
+\ena
+where $\Delta \ln A_{rp}(\om,\bbm)=\ln A^{obs}_{rp}(\om)-\ln A_{rp}(\om,\bbm)$ is the frequency-dependent amplitude anomaly measurements, and
+\eq
+W_{Q_{rp}}(\om)=\frac{W_{rp}(\om)}{\sigma_{Q_{rp}}^2(\om)\int\,W_{rp}(\om)d\om}\nn
+\en
+Both $\phi_P$ and $\phi_Q$ are dimensionaless numbers, and the adjoint sources $f^\dagger_P(t)$ and $f^\dagger_Q(t)$ have the units of $1/[\text{syn} \cdot \text{time}]$.
+
\end{enumerate}
+\section{More Notes}
+\subsection{NO\_WAVEFORM\_DIFFERENCE}
+Carl: please add more here.
+
+\subsection{DO\_WEIGHTING}
+It is also possible to weigh windowed measurements of different categories if weighting is set in \verb+ma_weighting+ module (\verb+DO_WEIGHTING+):
+\eq
+\Phi = \sum_\alpha \, W_\alpha \, \sum_i \phi_{\alpha,i}
+\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.
+
+\subsection{Time-domain taper}
+Time domain tapers are first applied to windowed data and synthetics (no matter what other options are)
+\eq
+w(t) = 1- \cos^a (2 \pi t/T), \quad 0\le t \le T
+\en
+where $a$ is an even integer (i.e., $a=10$) or a Welch window
+\eq
+w(t) = 1-\frac{(t-T/2)^2}{(T/2)^2}
+\en
+
+\subsection{Reonstructed synthetics}
+\eq
+\tilde{s}(\om)=s(\om) e^{\Delta \ln A(\om)- i \om \Delta\tau(\om)}
+\en
+with the FFT convention in D\&T (p109)
+
+\subsection{Post-processing of adjoint sources}
+Note since data and synthetics have actually been pre-filtered, therefore the adjoint source that satisfies
+\eq
+\delta \phi = \int f^\dagger (t) s^f(t) dt = \frac{1}{2\pi} \int f^\dagger(\om) B^*(\om)s^*(\om)\,d\om
+=\int F^{-1}[f^\dagger(\om)B^*(\om)](t) s(t)\,dt
+\en
+i.e., the adjoint source also needs to be filtered by the same band-pass filter to produce the exact Fr\'{e}chet
+derivatives.
+
+\subsection{DO\_RAY\_DENSITY\_SOURCE}
+Both \verb+sigma_tau+ and \verb+sigma_dlnA+, as well as measurements \verb+Delta_dtau/DlnA+ are set to be 1,
+and carried all the way to the adjoint source \verb+tr/amp_adj_src(:)+. \verb+DO_RAY_DENSITY_SOURCE+ automatically
+sets \verb+Error_type=0+.
+
+\subsection{MTM}
+\verb+dtau_mtm(:)+ starts from $i > 1$, i.e., ignore the time shift at very long periods, and equivalently,
+the first point of \verb+err_dt(1) = LARGE+ so that it does not contribute to adjoint sources and chi values.
%-------------------------------------------------------------------
\pagebreak
@@ -365,7 +659,7 @@
The following individuals have also contributed to the development of the source code or related scripts: Vala \vala, Min Chen.
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.
+The \verb+measure_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.
@@ -380,8 +674,9 @@
%\small
\begin{spacing}{1}
\addcontentsline{toc}{section}{References}
-\bibliographystyle{agu08}
-\bibliography{preamble,REFERENCES,refs_carl,refs_socal}
+\bibliographystyle{gji}
+%\bibliography{preamble,REFERENCES,refs_carl,refs_socal}
+\bibliography{measure}
\end{spacing}
%\normalsize
%-------------------------------------------------------------------
@@ -441,7 +736,7 @@
\caption[]
{{
{\em Left:} Data (black), synthetics (red), and measurement windows.
-{\em Right:} Adjoint sources for a cross-correlation traveltime difference, $\Delta T$ ({\tt imeas = 5}), .
+{\em Right:} Adjoint sources for a cross-correlation traveltime difference, $\Delta T$ ({\tt imeas = 5}).
\label{fig:iker05}
}}
\end{figure}
@@ -451,7 +746,7 @@
\caption[]
{{
{\em Left:} Data (black), synthetics (red), and measurement windows.
-{\em Right:} Adjoint sources for an amplitude difference, $\Delta \ln A$ ({\tt imeas = 6}), .
+{\em Right:} Adjoint sources for an amplitude difference, $\Delta \ln A$ ({\tt imeas = 6}).
\label{fig:iker06}
}}
\end{figure}
@@ -461,7 +756,7 @@
\caption[]
{{
{\em Left:} Data (black), synthetics (red), and measurement windows.
-{\em Right:} Adjoint sources for a multitaper traveltime difference ({\tt imeas = 7}), .
+{\em Right:} Adjoint sources for a multitaper traveltime difference ({\tt imeas = 7}).
\label{fig:iker07}
}}
\end{figure}
@@ -471,7 +766,7 @@
\caption[]
{{
{\em Left:} Data (black), synthetics (red), and measurement windows.
-{\em Right:} Adjoint sources for a multitaper amplitude difference ({\tt imeas = 8}), .
+{\em Right:} Adjoint sources for a multitaper amplitude difference ({\tt imeas = 8}).
\label{fig:iker08}
}}
\end{figure}
Added: seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/run.latex
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/run.latex (rev 0)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/run.latex 2012-11-27 17:05:45 UTC (rev 21075)
@@ -0,0 +1,5 @@
+latex measure_adj_manual
+bibtex measure_adj_manual
+latex measure_adj_manual
+dvipdf measure_adj_manual.dvi measure_adj_manual.pdf
+evince measure_adj_manual.pdf
Property changes on: seismo/3D/ADJOINT_TOMO/measure_adj/USER_MANUAL/run.latex
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/ma_constants.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/ma_constants.f90 2012-11-27 02:07:13 UTC (rev 21074)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/ma_constants.f90 2012-11-27 17:05:45 UTC (rev 21075)
@@ -28,9 +28,9 @@
! takes waveform of first trace dat_dtw, without taking the difference waveform to the second trace syn_dtw
! this is useful to cissor out later reflections which appear in data (no synthetics needed)
- logical:: NO_WAVEFORM_DIFFERENCE = .false.
+ logical, parameter :: NO_WAVEFORM_DIFFERENCE = .false.
! constructs adjoint sources for a "ray density" kernel, where all misfits are equal to one
- logical:: DO_RAY_DENSITY_SOURCE = .false.
-
+ logical, parameter :: DO_RAY_DENSITY_SOURCE = .false.
+
end module ma_constants
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/ma_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/ma_sub.f90 2012-11-27 02:07:13 UTC (rev 21074)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/ma_sub.f90 2012-11-27 17:05:45 UTC (rev 21075)
@@ -9,103 +9,97 @@
contains
- ! =================================================================================================
- ! subroutine mt_measure()
+ 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)
+
+ ! ======================================================================
+ ! subroutine mt_measure(): making measurements on windowed data and synthetics
! Boxcar/Cosine/Multitaper estimates of the transfer function between data and synthetics
!
! Input:
- ! is_mtm -- taper type: 1 for multitaper, 2 for boxcar taper, and 3 for cosine taper
- ! datafile -- original data file in SAC format
- ! file_prefix -- the output file prefix (usually in STA.NT.CMP.N format)
- ! dat_dt(:), syn_dt(:) t0, dt, npts -- original data and synthetics array
- ! tstart, tend -- start and end of the measurement window (can be from Alessia's code)
+ ! datafile --- note it is dummy right now (LQY)
+ ! filename --- output file prefix (e.g., OUT_DIR/PAS.CI.BHZ.02.mtm)
+ ! dat_dt(:), syn_dt(:), t0, dt, npts -- original data and synthetics array
+ ! tstart, tend -- start and end of the measurement window (from FLEXWIN)
! Output:
- ! istart -- starting index of the windowed portion of original trace
+ ! istart -- starting index of the windowed portion in the original trace
! dat_dtw(:), syn_dtw(:), nlen -- windowed and shifted data, windowed synthetics
! tshift, dlnA, cc_max -- time shift and amplitude cross-correlation measurements
+ ! sigma_dt_cc, sigma_dlnA_cc -- corresponding cc error estimates
+ ! syn_dtw_cc(:) -- reconstructed synthetics after Dt/DlnA corrections (LQY)
+ ! i_pmax_dat, i_pmax_syn: the indices of max amplitude of the (freq.) power for data/syn
! i_right -- the maximum reliable frequency estimate index
- ! trans_w(:) -- estimates of transfer function
+ ! trans_w(:) -- estimates of transfer function (MT only ???)
! dtau_w(:), dlnA_w(:) -- estimates of travel-time and amplitude anomaly
- ! err_dt(:), err_dlnA(:) -- error bar of the travel-time and amplitude estimates (MT only)
+ ! sigma_dt, sigma_dlnA -- average travel time and amplitude error for mt (same as sigma_dt_cc)
+ ! err_dt(:), err_dlnA(:) -- error bar of the travel-time and amplitude estimates (OPTIONAL)
!
- ! original coding in Fortran77 by Ying Zhou
- ! upgraded to Fortran90 by Alessia Maggi
+ ! other module variables used:
+ ! is_mtm -- taper type: 1 for multitaper, 2 for cosine taper, and 3 or 0 for boxcar taper
+ ! original coding in Fortran 77 by Ying Zhou
+ ! upgraded to Fortran 90 by Alessia Maggi
! organized into package form by Qinya Liu
! modifications by Carl Tape and Vala Hjorleifsdottir
!
- ! =================================================================================================
+ ! ====================================================================
- 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)
+ implicit none
- implicit none
- integer, intent(in) :: npts
+ character(len=150), intent(in) :: filename,datafile
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
+ double precision, intent(in) :: t0,dt,tstart,tend
+ integer, intent(in) :: npts ! rarely used
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
complex*16, dimension(:), intent(out) :: trans_w
- double precision, dimension(:), intent(out) :: dtau_w,dlnA_w,syn_dtw,dat_dtw,syn_dtw_cc
- integer, intent(out) :: nlen,i_right,istart,i_pmax_dat,i_pmax_syn
double precision, dimension(:), intent(out), optional :: err_dt,err_dlnA
- !double precision, intent(out), optional :: sigma_dt,sigma_dlnA
+
- double precision, dimension(NPT) :: syn_vtw, syn_dtw_mt, syn_dtw_mt_dt, &
+ ! local variables
+ double precision, dimension(NPT) :: syn_dtw_mt, 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
- integer :: ishift,i,ictaper,j,fnum,i_amp_max_unw,i_amp_max,i_right_stop,idf_new,iom
+ integer :: ishift,i,ictaper,j,fnum,i_amp_max_unw,i_amp_max,i_right_stop,idf_new,iom,is_mtm_av
complex*16, dimension(NPT) :: syn_dtwo, dat_dtwo, syn_dtw_ho, dat_dtw_ho, &
- top_mtm, bot_mtm, trans_mtm, wseis_recon
+ top_mtm, bot_mtm, trans_mtm
double precision, dimension(NPT) :: wvec, ey1, ey2, dtau_mtm, dlnA_mtm, &
phi_w, abs_w, err_phi, err_abs, phi_mtm, abs_mtm
double precision :: eph_ave,edt_ave,eabs_ave,eabs2_ave,eph_iom,edt_iom,eabs_iom,eabs2_iom
double precision, dimension(:,:),allocatable :: tas,phi_mul,abs_mul,dtau_mul,dlnA_mul
- character(len=150) :: filename
- logical :: output_logical,display_logical
!-------------------------------------------------------------
-
+ ! check the window [tstart, tend] (the first two unnecessary)
if ( tstart < t0 .or. tend > t0+(npts-1)*dt .or. tstart >= tend) then
print *, 'tstart, t0, tend, t0+(npts-1)*dt:'
print *, tstart, t0, tend, t0+(npts-1)*dt
stop 'Check tstart and tend'
+ else
+ print *, ' start and end time of window: ' ,sngl(tstart),sngl(tend)
endif
! initializes i_right
i_right = 0
- ! LQY -- is this too small ???
- wtr_mtm = 1.e-10
+ ! water level in spectral divisions for mtm measurements
+ wtr_mtm = 1.e-10 ! LQY -- move to par file? too small?
- filename = trim(file_prefix)
-
if (DISPLAY_DETAILS) then
- call dwascii(trim(file_prefix)//'_data',dat_dt,npts,t0,dt)
- call dwascii(trim(file_prefix)//'_syn',syn_dt,npts,t0,dt)
- !call dwsac1(trim(file_prefix)//'_data.sac',dat_dt,npts,t0,dt)
- !call dwsac1(trim(file_prefix)//'_syn.sac',syn_dt,npts,t0,dt)
-
-!!$ call dwrite_sacfile_f(datafile,trim(file_prefix)//'_data.sac',t0,npts,dat_dt)
-!!$ call dwrite_sacfile_f(datafile,trim(file_prefix)//'_syn.sac',t0,npts,syn_dt)
-!!$ !call dwrite_sacfile_f(datafile,trim(file_prefix)//'_diff_dms.sac',t0,npts,dat_dt-syn_dt)
-!!$ !call dwrite_ascfile_f(trim(file_prefix)//'_data.txt',t0,dt,npts,dat_dt)
-!!$ !call dwrite_ascfile_f(trim(file_prefix)//'_syn.txt',t0,dt,npts,syn_dt)
+ print *, ' original data/syn length', npts
+ call dwascii(trim(filename)//'.obs',dat_dt,npts,t0,dt)
+ call dwascii(trim(filename)//'.syn',syn_dt,npts,t0,dt)
endif
-
!--------------------------------------------------------------------------
! window and interpolate data and synthetics
!--------------------------------------------------------------------------
- ! interpolate data and synthetics, and also extract time-windowed records
call interpolate_dat_and_syn(dat_dt,syn_dt,tstart,tend,t0,dt,NPT,dat_dtw,syn_dtw,nlen,istart)
- if (nlen <= 1) stop 'Check the length of the data and syn arrays'
- if (nlen > NPT) stop 'Check the dimension of data and syn arrays'
-
! some constants
sfac1 = (2./dble(nlen))**2 ! for Welch window
ipwr_t = 10 ! for time-domain cosine taper: 1 - [cos(t)]^(ipwr)
@@ -121,53 +115,34 @@
enddo
if (DISPLAY_DETAILS) then
- print *, ' NPTs = ', NPT, ' new nlen = ', nlen
- call dwsac1(trim(file_prefix)//'.obs.sac',dat_dtw,nlen,tstart,dt)
- call dwsac1(trim(file_prefix)//'.syn.sac',syn_dtw,nlen,tstart,dt)
-!!$ call dwrite_sacfile_f(datafile,trim(file_prefix)//'.obs.sac',tstart,nlen,dat_dtw)
-!!$ call dwrite_sacfile_f(datafile,trim(file_prefix)//'.syn.sac',tstart,nlen,syn_dtw)
+ print *, ' new nlen of dat/syn_dtw = ', nlen, '(',sngl(tstart),sngl(tend),')'
+ call dwsac1(trim(filename)//'.obs.sac',dat_dtw,nlen,tstart,dt)
+ call dwsac1(trim(filename)//'.syn.sac',syn_dtw,nlen,tstart,dt)
endif
- ! save a copy of the windowed data
- !dat_dtwc(:) = dat_dtw(:)
-
!------------------------------------------------------------------
! cross-correlation traveltime and amplitude measurements
!------------------------------------------------------------------
! 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
+ ! LQY: Ying suggested to align them at relatively long periods first
call compute_cc(syn_dtw, dat_dtw, nlen, dt, ishift, tshift, dlnA, cc_max)
- ! compute velocity of synthetics
- do i = 2, nlen-1
- syn_vtw(i) = (syn_dtw(i+1) - syn_dtw(i-1)) / (2.0*dt)
- enddo
- syn_vtw(1) = (syn_dtw(2) - syn_dtw(1)) / dt
- syn_vtw(nlen) = (syn_dtw(nlen) - syn_dtw(nlen-1)) /dt
+ ! no need to compute velocity and acceleration of syn in mt_measure()
- ! acceleration
- !do i = 2, nlen-1
- ! syn_atw(i) = (syn_vtw(i+1) - syn_vtw(i-1)) / (2.0*dt)
- !enddo
- !syn_atw(1) = (syn_vtw(2) - syn_vtw(1)) / dt
- !syn_atw(nlen) = (syn_vtw(nlen) - syn_vtw(nlen-1)) / dt
-
! deconstruct data using (negative) cross-correlation measurments
call deconstruct_dat_cc(filename,dat_dtw,tstart,dt,nlen, &
ishift,tshift,dlnA,dat_dtw_cc)
- ! 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)
+ ! reconstruct synthetics using cross-correlation measurements (plotting purposes only)
+ ! syn_dtw_cc (reconstructed syn with both dt/dlnA), syn_dtw_cc_dt (with only dt)
+ call reconstruct_syn_cc(syn_dtw,tstart,dt,nlen, &
+ ishift,tshift,dlnA,syn_dtw_cc,syn_dtw_cc_dt)
if (OUTPUT_MEASUREMENT_FILES) then
call dwsac1(trim(filename)//'.recon_syn_cc.sac',syn_dtw_cc,nlen,tstart,dt)
call dwsac1(trim(filename)//'.recon_syn_cc_dt.sac',syn_dtw_cc_dt,nlen,tstart,dt)
-!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_cc.sac',tstart,nlen,syn_dtw_cc)
-!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_cc_dt.sac',tstart,nlen,syn_dtw_cc_dt)
-!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn',tstart,dt,nlen,syn_dtw_cc)
-!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn_dt',tstart,dt,nlen,syn_dtw_dt)
endif
! compute the estimated uncertainty for the cross-correlation measurment
@@ -176,7 +151,8 @@
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,imeas,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc)
+ is_mtm_av=2
+ call write_average_meas(filename,is_mtm_av,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc)
!========================================
@@ -187,13 +163,13 @@
! set up FFT for the frequency domain
!-----------------------------------------------------------------------------
- ! calculate frequency step and number of frequencies
+ ! calculate frequency step and number of frequencies used for FFT
f0 = 0.
df = 1./(NPT*dt)
dw = TWOPI * df
fnum = NPT/2 + 1
- ! calculate frequency spacing of sampling points
+ ! calculate frequency spacing for the actual time window
df_new = 1.0 / (tend-tstart)
idf_new = df_new / df
@@ -210,15 +186,13 @@
! create complex synthetic seismogram and CC-deconstructed data seismogram
syn_dtwo = cmplx(0.,0.)
dat_dtwo = cmplx(0.,0.)
- !syn_dtwo(1:nlen) = syn_dtw(1:nlen)
- !dat_dtwo(1:nlen) = dat_dtw_cc(1:nlen)
syn_dtwo(1:nlen) = cmplx(syn_dtw(1:nlen))
dat_dtwo(1:nlen) = cmplx(dat_dtw_cc(1:nlen))
call fft(LNPT,syn_dtwo,FORWARD_FFT,dt)
call fft(LNPT,dat_dtwo,FORWARD_FFT,dt)
- ! index of the freq of the max power in the windowed data
+ ! index of the freq of the max power in the windowed data/syn
ampmax_unw = 0.
i_pmax_dat = 1
do i = 1, fnum ! loop over frequencies
@@ -228,55 +202,51 @@
endif
enddo
- ! water level based untapered synthetics
- ! used to determine the i_right values (maximum frequency for measurement)
ampmax_unw = 0.
- do i = 1, fnum ! loop over frequencies
+ i_amp_max_unw = 1
+ do i = 1, fnum
if( abs(syn_dtwo(i)) > ampmax_unw) then
ampmax_unw = abs(syn_dtwo(i))
i_amp_max_unw = i
endif
enddo
+ i_pmax_syn = i_amp_max_unw
+
+ ! water level based untapered synthetics
+ ! used to determine the i_right values (maximum frequency for measurement)
wtr_use_unw = cmplx(ampmax_unw * WTR, 0.)
- ! index of the freq of the max power in the windowed synthetics
- i_pmax_syn = i_amp_max_unw
-
i_right = fnum
i_right_stop = 0
do i = 1,fnum
if( abs(syn_dtwo(i)) <= abs(wtr_use_unw) .and. i_right_stop==0 .and. i > i_amp_max_unw ) then
- i_right_stop = 1
+ i_right_stop = 1 ! power dips below water-level
i_right = i
endif
if( abs(syn_dtwo(i)) >= 10.*abs(wtr_use_unw) .and. i_right_stop==1 .and. i > i_amp_max_unw) then
- i_right_stop = 0
+ i_right_stop = 0 ! power goes above 10*water-level
i_right = i
endif
enddo
if (DISPLAY_DETAILS) then
- print *, 'Frequency of max power in windowed synthetic (Hz):'
- print *, ' i_pmax_syn = ', i_pmax_syn, ', f_pmax = ', sngl(i_pmax_syn * df), ', T_pmax = ', sngl(1./(i_pmax_syn*df))
- print *, 'FFT freq spacing df = ', sngl(df)
- print *, 'measurement spacing df_new = ', sngl(df_new)
- print *, ' i_right = ', i_right, ', stopping freq = ', sngl(i_right * df)
+ print *, ' frequency of max power in windowed synthetic (Hz):'
+ print *, ' i_pmax_syn = ', i_pmax_syn, ', f_pmax = ', sngl(i_pmax_syn * df), ', T_pmax = ', sngl(1./(i_pmax_syn*df))
+ print *, ' FFT freq spacing df = ', sngl(df)
+ print *, ' measurement freq spacing df_new = ', sngl(df_new)
+ print *, ' i_right = ', i_right, ', stopping freq/period = ', sngl(i_right * df),'/', &
+ sngl(1./(df*i_right))
! write out power for each signal
- call dwascii(trim(file_prefix)//'.obs.power',abs(dat_dtwo(1:i_right)),i_right,df,df)
- call dwascii(trim(file_prefix)//'.syn.power',abs(syn_dtwo(1:i_right)),i_right,df,df)
- !call dwsac1(trim(file_prefix)//'.obs.power.sac',abs(dat_dtwo(1:i_right)),i_right,df,df)
- !call dwsac1(trim(file_prefix)//'.syn.power.sac',abs(syn_dtwo(1:i_right)),i_right,df,df)
-!!$ call dwrite_ascfile_f(trim(file_prefix)//'.obs.power',df,df,i_right,abs(dat_dtwo(1:i_right)) )
-!!$ call dwrite_ascfile_f(trim(file_prefix)//'.syn.power',df,df,i_right,abs(syn_dtwo(1:i_right)) )
-
+ call dwascii(trim(filename)//'.obs.power',abs(dat_dtwo(1:i_right)),i_right,df,df)
+ call dwascii(trim(filename)//'.syn.power',abs(syn_dtwo(1:i_right)),i_right,df,df)
endif
!-------------------------------------------------------------------------------
! single-taper estimation of transfer function
!-------------------------------------------------------------------------------
- ! assign number of tapers
+ ! assign number of tapers (is_mtm can only be 1, 2 or 3 at this point)
if (is_mtm == 1) then
ntaper = int(NPI * 2.0)
else
@@ -292,13 +262,6 @@
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.)
@@ -324,16 +287,16 @@
call fft(LNPT,dat_dtw_ho,FORWARD_FFT,dt)
! compute water level for single taper measurement by finding max spectral power
- ! in the tapered synthetics record
+ ! in the single-tapered synthetics record
ampmax = 0.
- do i = 1, fnum ! loop over frequencies
- if( abs(syn_dtw_ho(i)) > ampmax) then ! syn, single_tapered
+ i_amp_max = 1
+ do i = 1, fnum
+ if( abs(syn_dtw_ho(i)) > ampmax) then
ampmax = abs(syn_dtw_ho(i))
i_amp_max = i
endif
enddo
wtr_use = cmplx(ampmax * WTR, 0.)
- !print *, ' wtr_use :', wtr_use
! calculate top and bottom of MT transfer function
do i = 1, fnum
@@ -342,31 +305,28 @@
! calculate transfer function for single taper measurement using water level
if (is_mtm /= 1) then
- if(abs(syn_dtw_ho(i)) > abs(wtr_use)) trans_w(i) = dat_dtw_ho(i) / syn_dtw_ho(i)
- if(abs(syn_dtw_ho(i)) <= abs(wtr_use)) trans_w(i) = dat_dtw_ho(i) / (syn_dtw_ho(i)+wtr_use)
+ if (abs(syn_dtw_ho(i)) > abs(wtr_use)) trans_w(i) = dat_dtw_ho(i) / syn_dtw_ho(i)
+ if (abs(syn_dtw_ho(i)) <= abs(wtr_use)) trans_w(i) = dat_dtw_ho(i) / (syn_dtw_ho(i)+wtr_use)
endif
enddo
- ! 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,
- ! 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
- call write_trans(filename,trans_w,wvec,fnum,i_right,idf_new,df,tshift,dlnA, &
- phi_w,abs_w,dtau_w,dlnA_w,dtau_wa,dlnA_wa)
- call reconstruct_syn(filename,syn_dtwo,wvec,dtau_w,dlnA_w, &
- 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, imeas, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
- endif
-
enddo ! ictapers
- ! for single taper, pass back the transfer function
- if (is_mtm /= 1) return
+ ! 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,
+ ! 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
+ call write_trans(filename,trans_w,wvec,fnum,i_right,idf_new,df,tshift,dlnA, &
+ phi_w,abs_w,dtau_w,dlnA_w,dtau_wa,dlnA_wa)
+ call reconstruct_syn(filename,syn_dtwo,wvec,dtau_w,dlnA_w, &
+ 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)
+ return ! over for single taper (no further error estimates)
+ endif
+
!-------------------------------------------------------------------------------
! multitaper estimation of transfer function
!-------------------------------------------------------------------------------
@@ -374,20 +334,20 @@
! water level for multitaper measurements
ampmax = 0.
do i = 1, fnum
- if( abs(bot_mtm(i)) > ampmax) then
- ampmax = abs(bot_mtm(i))
- i_amp_max = i
- endif
+ if( abs(bot_mtm(i)) > ampmax) then
+ ampmax = abs(bot_mtm(i))
+ i_amp_max = i
+ endif
enddo
wtr_use = cmplx(ampmax * wtr_mtm**2, 0.)
!wtr_use = cmplx(ampmax * WTR, 0.)
-
+
! calculate MT transfer function using water level
do i = 1, fnum
- if(abs(bot_mtm(i)) > abs(wtr_use)) trans_mtm(i) = top_mtm(i) / bot_mtm(i)
- if(abs(bot_mtm(i)) < abs(wtr_use)) trans_mtm(i) = top_mtm(i) / (bot_mtm(i)+wtr_use)
+ if(abs(bot_mtm(i)) > abs(wtr_use)) trans_mtm(i) = top_mtm(i) / bot_mtm(i)
+ if(abs(bot_mtm(i)) < abs(wtr_use)) trans_mtm(i) = top_mtm(i) / (bot_mtm(i)+wtr_use)
enddo
-
+
! multitaper phase, abs, tt, and amp (freq)
call write_trans(filename,trans_mtm,wvec,fnum,i_right,idf_new,df,tshift,dlnA, &
phi_mtm,abs_mtm,dtau_mtm,dlnA_mtm,dtau_wa,dlnA_wa)
@@ -400,25 +360,19 @@
!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)
! CHT: estimate error using the same procedure as for the cross-correlation error estimate
- !sigma_dt = 1. ; sigma_dlnA = 1.
- !call compute_average_error(dat_dtw,syn_dtw_mt,syn_dtw_mt_dt,nlen,dt,sigma_dt,sigma_dlnA)
sigma_dt = sigma_dt_cc ; sigma_dlnA = sigma_dlnA_cc
- ! write average multitaper measurement to file
- call write_average_meas(file_prefix, imeas, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
+ ! write average multitaper measurement to file
+ is_mtm_av=1
+ call write_average_meas(filename, is_mtm_av, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
!-------------------------------------------------------------------------------
! multitaper error estimation
!-------------------------------------------------------------------------------
- if (ntaper > 1) then
+ if (ntaper > 1) then ! this should always be true
- ! save a copy of the control logicals
- output_logical = OUTPUT_MEASUREMENT_FILES
- display_logical = DISPLAY_DETAILS
- ! avoid I/O output for MT error estimates
- OUTPUT_MEASUREMENT_FILES = .false.
- DISPLAY_DETAILS = .false.
+ ! no need to save a copy of the control logicals any more (output/display)
! allocate Jacknife MT estimates
allocate(phi_mul(NPT,ntaper))
@@ -463,7 +417,7 @@
i_amp_max = i
endif
enddo
- wtr_use = cmplx(ampmax * wtr_mtm ** 2, 0.)
+ wtr_use = cmplx(ampmax * wtr_mtm ** 2, 0.) ! again very small wtr_mtm is used
! calculate transfer function using water level
do i = 1, fnum
@@ -474,21 +428,23 @@
call write_trans(filename,trans_mtm,wvec,fnum,i_right,idf_new,df,tshift,dlnA, &
phi_mul(:,iom),abs_mul(:,iom),dtau_mul(:,iom),dlnA_mul(:,iom))
- enddo ! iom
+ enddo ! ntaper
!----------------------
- open(10,file=trim(filename)//'.err_ph')
- open(20,file=trim(filename)//'.err_dt')
- open(30,file=trim(filename)//'.err_abs')
- open(40,file=trim(filename)//'.err_dlnA')
+ if (OUTPUT_MEASUREMENT_FILES) then
+ open(10,file=trim(filename)//'.err_ph')
+ open(20,file=trim(filename)//'.err_dt')
+ open(30,file=trim(filename)//'.err_abs')
+ open(40,file=trim(filename)//'.err_dlnA')
! CHT: Since all freq. domain points are used in constructing the
! adjoint source, we also want to show the entire sigma(f) functions,
! not just the sub-sampled version.
- open(50,file=trim(filename)//'.err_dt_full')
- open(60,file=trim(filename)//'.err_dlnA_full')
-
+ open(50,file=trim(filename)//'.err_dt_full')
+ open(60,file=trim(filename)//'.err_dlnA_full')
+ endif
+
err_phi = 0.
err_dt = 0.
err_abs = 0.
@@ -537,7 +493,7 @@
err_dlnA(i) = sqrt( err_dlnA(i) / (ntaper * (ntaper-1) ) )
! only write out the errors for the 'independent' freq-domain sampling points
- if (mod(i,idf_new) == 0) then
+ if (mod(i,idf_new) == 0 .and. OUTPUT_MEASUREMENT_FILES) then
write(10,*) df*(i-1), phi_mtm(i), err_phi(i)
if (i > 1) write(20,*) df*(i-1), dtau_mtm(i), err_dt(i)
write(30,*) df*(i-1), abs_mtm(i), err_abs(i)
@@ -545,8 +501,10 @@
endif
! CHT: write out the entire dt(f) and dlnA(f) for adjoint sources
- write(50,*) df*(i-1), dtau_mtm(i), err_dt(i)
- write(60,*) df*(i-1), dlnA_mtm(i), err_dlnA(i)
+ if (OUTPUT_MEASUREMENT_FILES) then
+ write(50,*) df*(i-1), dtau_mtm(i), err_dt(i)
+ write(60,*) df*(i-1), dlnA_mtm(i), err_dlnA(i)
+ endif
enddo ! i_right
@@ -562,14 +520,10 @@
dtau_w = dtau_mtm
dlnA_w = dlnA_mtm
- ! reset the control parameters
- OUTPUT_MEASUREMENT_FILES = output_logical
- DISPLAY_DETAILS = display_logical
-
endif
! ------------------------------------------------------------------
- ! End error calculation loop
+ ! End error calculation
! ------------------------------------------------------------------
end subroutine mt_measure
@@ -581,21 +535,24 @@
! adjoint sources by assimulate the measurements passed from mt_measure()
!
! Input:
- ! 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
+ ! istart -- starting index of the windowed portion of original trace, used to assign
+ ! tr/amp_adj_src(:) corresponding to the original synthetics
! dat_dtw(:), syn_dtw(:), nlen, dt -- windowed data and synthetics
! with length nlen and sampling rate dt
- ! tshift, dlnA -- cross-correlation traveltime and amplitude measurements
- ! dtau_w(:), dlnA_w(:), err_dtau(:), err_dlnA(:), i_right -- traveltime and amplitude measurements
- ! and corresponding error bars as a function of frequency (1: i_right)
+ ! tshift, dlnA, sigma_dt_cc,sigma_dlnA_cc -- CC traveltime/amplitude measurements/error estimates
+ ! dtau_w(:), dlnA_w(:), err_dtau(:), err_dlnA(:) -- frequency-dependent tt/amp measurements/error
+ ! i_left, i_right -- range of indices of valid freqs for frequency-dependent tt/amp measurements
!
! Output:
+ ! window_chi(:) -- all available scalar measurement values and chi values (20 in total)
! tr_adj_src(:), tr_chi -- travel-time adjoint source and chi value
- ! am_adj_src(:), am_chi -- amplitude adjoint source and chi value
- ! window_chi(:) -- all available scalar measurement values and chi values
+ ! am_adj_src(:), am_chi -- amplitude adjoint source and chi value (optional)
+
+ ! other global vars used:
+ ! imeas -- adjoint source type: 1/2 for waveform, 3/4 for cc/banana-doughnut, 5/6 for cc, 7/8 for mtm
+ ! this subroutine is only called when COMPUTE_ADJOINT_SOURCE = true.
!
- ! original coding by Carl Tape, finalized by Qinya Liu
+ ! original coding by Carl Tape and finalized by Qinya Liu
! ======================================================================================================
subroutine mt_adj(istart,dat_dtw,syn_dtw,nlen,dt,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc, &
@@ -608,9 +565,9 @@
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(NCHI), intent(inout) :: window_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
@@ -635,19 +592,62 @@
! waveform
if(imeas==1 .or. imeas==2) then
- print *, 'computing waveform adjoint source'
+ print *, ' computing waveform adjoint source'
elseif(imeas==3 .or. imeas==4) then
- print *, 'computing banana-doughtnut adjoint source'
+ print *, ' computing banana-doughtnut adjoint source'
elseif(imeas==5 .or. imeas==6) then
- print *, 'computing cross-correlation adjoint source'
+ print *, ' computing cross-correlation adjoint source'
elseif(imeas==7 .or. imeas==8) then
- print *, 'computing multitaper adjoint source'
+ print *, ' computing multitaper adjoint source'
endif
- ! ----------------------
- ! TAPERS
- ! ----------------------
+ ! define post-processing time-domain taper
+ ! NOTE: If the adjoint sources will be band-pass filtered at the end,
+ ! then perhaps time_window is not necessary (i.e., use boxcar).
+ ! However, if you are using a waveform difference, then you want
+ ! to make sure that the endpoints of the windows are at zero, since
+ ! you would NOT apply the post-processing band-pass filter.
+ ! Note also that the dat/syn_dtw(:) has already been time-windowed
+ time_window(:) = 0.
+ ipwr_t = 10
+ do i = 1,nlen
+ fac = 1. ! boxcar window
+ !fac = 1 - sfac2*((i-1) - dble(nlen1)/2.0)**2 ! welch window
+ !fac = 1. - cos(PI*(i-1)/(nlen-1))**ipwr_t ! cosine window
+ time_window(i) = fac
+ enddo
+
+ ! ----------------------------------
+ ! CROSS CORRELATION ADJOINT SOURCES LQY: does time_window needs to be applied here?
+ ! ----------------------------------
+ if( (imeas >= 3).and.(imeas <= 6) ) then
+
+ ! compute synthetic velocity
+ do i = 2, nlen-1
+ syn_vtw(i) = (syn_dtw(i+1) - syn_dtw(i-1)) / (2.0*dt)
+ enddo
+ syn_vtw(1) = (syn_dtw(2) - syn_dtw(1)) / dt
+ syn_vtw(nlen) = (syn_dtw(nlen) - syn_dtw(nlen-1)) / dt
+
+ ! compute CC traveltime and amplitude banana-dougnut kernels
+ ft_bar_t = 0.
+ Nnorm = dt * sum( syn_vtw(1:nlen) * syn_vtw(1:nlen) )
+ ft_bar_t(1:nlen) = -syn_vtw(1:nlen) / Nnorm
+
+ fa_bar_t = 0.
+ Mnorm = dt * sum( syn_dtw(1:nlen) * syn_dtw(1:nlen) )
+ fa_bar_t(1:nlen) = syn_dtw(1:nlen) / Mnorm
+ endif
+
+ ! ----------------------------------------------
+ ! FREQUENCY-DOMAIN TAPERS FOR MT ADJOINT SOURCES
+ ! ----------------------------------------------
if( is_mtm == 1 ) then
+
+ ! initialize water levels for err_dtau/dlnA division
+ dtau_wtr = WTR * sum(abs(dtau_w(i_left:i_right)))/(i_right-i_left) ! CHT i_left
+ dlnA_wtr = WTR * sum(abs(dlnA_w(i_left:i_right)))/(i_right-i_left) ! CHT i_left
+
! frequency-domain tapers
! THIS CHOICE WILL HAVE AN EFFECT ON THE ADJOINT SOURCES
ipwr_w = 10
@@ -663,16 +663,14 @@
! 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 *, 'Taper normalization factor, ffac = ', ffac
+ if (DISPLAY_DETAILS) print *, ' f-dom taper normalization factor, ffac = ', ffac
! wp_taper and wq_taper are modified frequency-domain tapers
! Notice the option to include the frequency-dependent error.
wp_taper(:) = 0.
wq_taper(:) = 0.
- dtau_wtr = WTR * sum(abs(dtau_w(i_left:i_right)))/(i_right-i_left) ! CHT i_left
- dlnA_wtr = WTR * sum(abs(dlnA_w(i_left:i_right)))/(i_right-i_left) ! CHT i_left
- do i = i_left, i_right ! CHT: 1 --> i_left
+ do i = i_left, i_right ! CHT: i_left = 1
if (ERROR_TYPE == 0 .or. DO_RAY_DENSITY_SOURCE ) then
! no error estimate
@@ -680,13 +678,12 @@
wp_taper(i) = w_taper(i) / ffac
wq_taper(i) = w_taper(i) / ffac
- elseif (ERROR_TYPE == 1) then
+ elseif (ERROR_TYPE == 1) then ! CC error as a constant for all freqs.
! 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 (ERROR_TYPE == 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)
@@ -702,52 +699,6 @@
!!$ enddo
!!$ close(88)
- endif ! is_mtm == 1
-
-
- ! post-processing time-domain taper
- ! NOTE: If the adjoint sources will be band-pass filtered at the end,
- ! then perhaps time_window is not necessary (i.e., use boxcar).
- ! However, if you are using a waveform difference, then you want
- ! to make sure that the endpoints of the windows are at zero, since
- ! you would NOT apply the post-processing band-pass filter.
- time_window(:) = 0.
- ipwr_t = 10
- do i = 1,nlen
- fac = 1. ! boxcar window
- !fac = 1 - sfac2*((i-1) - dble(nlen1)/2.0)**2 ! welch window
- !fac = 1. - cos(PI*(i-1)/(nlen-1))**ipwr_t ! cosine window
- time_window(i) = fac
- enddo
-
- ! ----------------------------------
- ! CROSS CORRELATION ADJOINT SOURCES
- ! ----------------------------------
- if( (imeas >= 3).and.(imeas <= 6) ) then
-
- ! compute synthetic velocity
- do i = 2, nlen-1
- syn_vtw(i) = (syn_dtw(i+1) - syn_dtw(i-1)) / (2.0*dt)
- enddo
- syn_vtw(1) = (syn_dtw(2) - syn_dtw(1)) / dt
- syn_vtw(nlen) = (syn_dtw(nlen) - syn_dtw(nlen-1)) / dt
-
- ! compute CC traveltime and amplitude banana-dougnut kernels
- ft_bar_t = 0.
- Nnorm = dt * sum( syn_vtw(1:nlen) * syn_vtw(1:nlen) )
- ft_bar_t(1:nlen) = -syn_vtw(1:nlen) / Nnorm
-
- fa_bar_t = 0.
- Mnorm = dt * sum( syn_dtw(1:nlen) * syn_dtw(1:nlen) )
- fa_bar_t(1:nlen) = syn_dtw(1:nlen) / Mnorm
- endif
-
- ! -------------------------------
- ! MULTITAPER ADJOINT SOURCES
- ! -------------------------------
-
- if ( is_mtm == 1 .and. COMPUTE_ADJOINT_SOURCE ) then
-
! allocate MT variables
ntaper = int(NPI * 2.0)
allocate(tas(NPT,ntaper))
@@ -779,7 +730,7 @@
syn_dtw_ho_all(1:nlen,ictaper) = cmplx(syn_dtw_h(1:nlen),0.)
syn_vtw_ho_all(1:nlen,ictaper) = cmplx(syn_vtw_h(1:nlen),0.)
- ! apply FFT get complex spectra
+ ! apply FFT get complex spectra LQY: can't we use disp->velo to save one fft operation?
call fft(LNPT,syn_dtw_ho_all(:,ictaper),FORWARD_FFT,DT)
call fft(LNPT,syn_vtw_ho_all(:,ictaper),FORWARD_FFT,DT)
@@ -797,14 +748,14 @@
pwc_adj(:) = cmplx(0.,0.)
qwc_adj(:) = cmplx(0.,0.)
- do i = 1, i_right
+ do i = 1, i_right ! LQY: i_left = 1
pwc_adj(i) = syn_vtw_ho_all(i,ictaper) / v_bot_mtm(i)
qwc_adj(i) = -syn_dtw_ho_all(i,ictaper) / d_bot_mtm(i)
enddo
! compute P_j(w) and Q_j(w)
! NOTE: the MT measurement is incorporated here
- ! also note that wp_taper and wq_taper can contain uncertainty estimations
+ ! also note that wp_taper and wq_taper can contain uncertainty estimations
if( DO_RAY_DENSITY_SOURCE ) then
! uses a misfit measurement dtau, dlnA = 1
pwc_adj(:) = pwc_adj(:) * cmplx(1.0,0.) * cmplx(wp_taper(:),0.)
@@ -819,8 +770,7 @@
call fftinv(LNPT,pwc_adj,REVERSE_FFT,dt,dtau_pj_t)
call fftinv(LNPT,qwc_adj,REVERSE_FFT,dt,dlnA_qj_t)
- ! create adjoint source
- ! applies taper to time signal
+ ! create adjoint source: applies taper to time signal
fp(:) = fp(:) + tas(:,ictaper) * dtau_pj_t(:)
fq(:) = fq(:) + tas(:,ictaper) * dlnA_qj_t(:)
@@ -829,7 +779,7 @@
endif ! MT adjoint source
! -------------------------------------
- ! COMPUTE ADJOINT SOURCE
+ ! ASSEMBLE VARIOUS ADJOINT SOURCES
! -------------------------------------
tr_adj_src = 0.
@@ -847,49 +797,45 @@
waveform_s2 = waveform_temp2
waveform_chi = waveform_temp3
- ! compute traveltime and amplitude adjoint sources
- if (COMPUTE_ADJOINT_SOURCE) then
+ ! compute traveltime and amplitude adjoint sources for imeas
+ do i = 1,nlen
+ i1 = istart + i ! start index in the full adjoint source array(1:npts)
- do i = 1,nlen
- i1 = istart + 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)
+ ! waveform
+ if(imeas==1 .or. imeas==2) then
+ tr_adj_src(i1) = -dat_dtw(i)/waveform_d2 * time_window(i) ! imeas=1
! consider normalizing this by waveform_d2
-
+ am_adj_src(i1) = ( syn_dtw(i) - dat_dtw(i) ) * time_window(i) ! imeas=2
+
! use pure data waveform in time window
if( NO_WAVEFORM_DIFFERENCE ) then
- tr_adj_src(i1) = dat_dtw(i) * time_window(i) ! waveform misfit
+ tr_adj_src(i1) = dat_dtw(i) * time_window(i) ! waveform misfit (imeas=1)
endif
-
- ! 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)
-
- ! cross-correlation
- elseif(imeas==5 .or. imeas==6) then
+
+ ! 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) ! imeas=3
+ am_adj_src(i1) = fa_bar_t(i) * time_window(i) ! imreas=4
+
+ ! 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)
-
+
! ray density
if( DO_RAY_DENSITY_SOURCE ) then
- ! uses a misfit measurement of 1
- tr_adj_src(i1) = - (1.0) * ft_bar_t(i) * time_window(i)
- am_adj_src(i1) = - (1.0) * fa_bar_t(i) * time_window(i)
+ ! uses a misfit measurement of 1
+ tr_adj_src(i1) = - (1.0) * ft_bar_t(i) * time_window(i)
+ am_adj_src(i1) = - (1.0) * fa_bar_t(i) * time_window(i)
endif
-
- ! multitaper
- elseif(imeas==7 .or. imeas==8) then
+
+ ! 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
+ enddo
- endif
-
! -------------------------------------
! COMPUTE MISFIT FUNCTION VALUE
! -------------------------------------
@@ -902,14 +848,13 @@
! 4: cross-correlation, dlnA
!window_chi(:) = 0.
-
! misfit function value
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)
+ ! cc/averaged mt measurement (no uncertainty estimates)
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
@@ -918,20 +863,18 @@
! replaces misfit function values
if( DO_RAY_DENSITY_SOURCE ) then
! uses misfit measurements equal to 1
- ! misfit function value
if(is_mtm==1) window_chi(1) = 0.5 * 2.0 * df * sum( (1.0)**2 * wp_taper(1:i_right) )
if(is_mtm==1) window_chi(2) = 0.5 * 2.0 * df * sum( (1.0)**2 * wq_taper(1:i_right) )
window_chi(3) = 0.5 * (1.0)**2
window_chi(4) = 0.5 * (1.0)**2
- ! measurement (no uncertainty estimates)
if(is_mtm==1) window_chi(5) = sum( 1.0 * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
if(is_mtm==1) window_chi(6) = sum( 1.0 * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
window_chi(7) = 1.0
window_chi(8) = 1.0
endif
- ! estimated mesurement uncertainties
+ ! estimated measurement uncertainties
if(is_mtm==1) window_chi(9) = sigma_dt
if(is_mtm==1) window_chi(10) = sigma_dlnA
window_chi(11) = sigma_dt_cc
@@ -943,12 +886,6 @@
window_chi(15) = 0.5 * waveform_chi
window_chi(16) = nlen*dt
-!!$ open(88,file='testing.dat')
-!!$ do i = 1,i_right
-!!$ write(88,'(5e18.6)') df*i, dtau_w(i), dlnA_w(i), wp_taper(i), wq_taper(i)
-!!$ enddo
-!!$ close(88)
-
if(imeas <= 2) then ! waveform
tr_chi = 0.5 * waveform_chi
am_chi = 0.5 * waveform_chi
@@ -966,11 +903,11 @@
end subroutine mt_adj
!==============================================================================
+ ! OTHER MAJOR SURBROUTINES
!==============================================================================
- !----------------------------------------------------------------------
-
subroutine bandpass(x,n,delta_t,f1,f2)
+ ! a double-precision wrapper around sac xapiir()
! modified from FLEXWIN subroutines on 26-July-2009
implicit none
@@ -992,6 +929,7 @@
! does band-pass filter
! BU - butterworth
! BP - bandpass
+ ! LQY: Shouldn't delta_t_sngl = sngl(delta_t) still be done? same for f1,f2?
call xapiir(x_sngl,n,'BU',TRBDNDW,APARM,IORD,'BP',f1,f2,delta_t,PASSES)
x(1:n) = dble(x_sngl(1:n))
@@ -1100,7 +1038,7 @@
.or. (dlnA < DLNA_MIN) .or. (dlnA > DLNA_MAX) ) then
! zero the CC measurments
if (DISPLAY_DETAILS) then
- print *, 'Fail if ANY of these is true :'
+ print *, 'CC measurements: failed because ONE 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
@@ -1143,11 +1081,12 @@
Wlen = dt*nlen
if( NCYCLE_IN_WINDOW * T_pmax > Wlen ) then
- print *, 'rejecting trace for too few cycles within time window:'
- print *, ' T_pmax : ', T_pmax
- print *, ' Wlen : ', Wlen
- print *, ' NCYCLE_IN_WINDOW : ', NCYCLE_IN_WINDOW
- print *, ' REJECTION: ', NCYCLE_IN_WINDOW*T_pmax, Wlen, NCYCLE_IN_WINDOW * T_pmax < Wlen
+ print *, ' MTM: rejecting for too few cycles within time window:'
+ print *, ' T_pmax : ', sngl(T_pmax)
+ print *, ' Wlen : ', sngl(Wlen)
+ print *, ' NCYCLE_IN_WINDOW : ', NCYCLE_IN_WINDOW
+ print *, ' REJECTION: ', sngl(NCYCLE_IN_WINDOW*T_pmax), &
+ sngl(Wlen), NCYCLE_IN_WINDOW * T_pmax < Wlen
use_trace = .false.
endif
@@ -1156,7 +1095,8 @@
! DECREASE the frequency range of the measurement (and adjoint source)
! --> note NCYCLE_IN_WINDOW and window length
- ! We subjectively state that we want at least 10 frequency points for the multitaper measurement.
+ ! We subjectively state that we want at least ntaper number of
+ ! frequency points between [fstart, fend] for the multitaper measurement.
fstart = max(fstart, NCYCLE_IN_WINDOW/Wlen)
fend = min(fend, 1./(2.0*dt))
@@ -1165,12 +1105,13 @@
if( ntaper > 10 ) ntaper = 10
if( ntaper < 1 ) ntaper = 10
if( use_trace .and. fstart >= fend - ntaper*df ) then
- print *, 'rejecting trace for frequency range (NCYCLE_IN_WINDOW/Wlen):'
- print *, ' fstart, fend, df, ntaper : ', fstart,fend,df,ntaper
- print *, ' NCYCLE_IN_WINDOW, Wlen : ', NCYCLE_IN_WINDOW,Wlen,NCYCLE_IN_WINDOW/Wlen
- print *, ' REJECTION fstart >= fend - ntaper*df : ', fstart, fend - ntaper*df, fstart >= fend - ntaper*df
+ print *, ' MTM: rejecting for frequency range (NCYCLE_IN_WINDOW/Wlen):'
+ print *, ' fstart, fend, df, ntaper : ', sngl(fstart),sngl(fend),sngl(df),ntaper
+ print *, ' NCYCLE_IN_WINDOW, Wlen : ', NCYCLE_IN_WINDOW,sngl(Wlen), &
+ sngl(NCYCLE_IN_WINDOW/Wlen)
+ print *, ' REJECTION fstart >= fend - ntaper*df : ', sngl(fstart), &
+ sngl(fend - ntaper*df), fstart >= fend - ntaper*df
use_trace = .false.
- !stop 'testing rejection criteria'
endif
! assemble frequency vector (NPT is the FFT length)
@@ -1197,27 +1138,31 @@
! determine the indices that denote the new frequency range (CHT)
! IT SEEMS LIKE THERE SHOULD BE NO NEED FOR THIS, SINCE THE SIGNAL HAS ALREADY
! BEEN BAND-PASSED PRIOR TO MAKING THE MULTITAPER MEASUREMENT.
- if (1==1) then
- i_left_old = i_left
- i_right_old = i_right
- do j = i_left_old, i_right_old
- if (fvec(j) > fstart) then
- i_left = j-1
- exit
- endif
- enddo
- do j = i_left_old, i_right_old
- if (fvec(j) > fend) then
- i_right = j-1
- exit
- endif
- enddo
- if (DISPLAY_DETAILS) then
- write(*,'(a24,2i6,2f14.8)') 'Old frequency bounds :', i_left_old, i_right_old, df*i_left_old, df*i_right_old
- write(*,'(a24,2i6,2f14.8)') 'New frequency bounds :', i_left, i_right, df*i_left, df*i_right
+ ! LQY: this should be useful to constraint the range of which use_trace is checked
+
+ i_left_old = i_left
+ i_right_old = i_right
+ do j = i_left_old, i_right_old
+ if (fvec(j) > fstart) then
+ i_left = j-1
+ exit
endif
+ enddo
+ do j = i_left_old, i_right_old
+ if (fvec(j) > fend) then
+ i_right = j-1
+ exit
+ endif
+ enddo
+
+ if (DISPLAY_DETAILS) then
+ print *, ' determine the validity of MTM over new f bounds:'
+ write(*,'(a24,2i6,2f14.8)') ' Old frequency bounds :', i_left_old, i_right_old, &
+ sngl(df*(i_left_old-1)), sngl(df*(i_right_old-1))
+ write(*,'(a24,2i6,2f14.8)') ' New frequency bounds :', i_left, i_right, &
+ sngl(df*(i_left-1)), sngl(df*(i_right-1))
endif
-
+
! update the frequency limits
fstart = (i_left-1)*df
fend = (i_right-1)*df
@@ -1228,9 +1173,9 @@
dtau_w(:) = 0.
use_trace = .false.
if (DISPLAY_DETAILS) then
- print *, 'rejecting trace for too small a time shift:'
- print *, ' dt = ', dt
- print *, ' tshift = ', tshift
+ print *, 'MTM: rejecting for too small a time shift:'
+ print *, ' dt = ', sngl(dt)
+ print *, ' tshift = ', sngl(tshift)
endif
endif
@@ -1241,13 +1186,14 @@
.or. abs(dtau_w(j)) > DT_MAX_SCALE*abs(tshift))) then
use_trace = .false.
if (DISPLAY_DETAILS) then
- print *, 'rejecting trace (T leads to rejection):'
- print *, ' f = ', fvec(j), j
- print *, 'DT_FAC (lower) : ', abs(dtau_w(j)), 1./(DT_FAC * fvec(j)), abs(dtau_w(j)) > 1./(DT_FAC * fvec(j))
- print *, 'ERR_FAC (lower) : ', err_dt(j), 1./(ERR_FAC * fvec(j)), err_dt(j) > 1./(ERR_FAC * fvec(j))
- print *, 'DT_MAX_SCALE (lower) : ', abs(dtau_w(j)), DT_MAX_SCALE*abs(tshift), &
- abs(dtau_w(j)) > DT_MAX_SCALE*abs(tshift)
- !stop 'testing MT trace rejection'
+ print *, ' MTM: rejecting trace on dtau/err at specific frequency:'
+ print *, ' f = ', sngl(fvec(j)), j, sngl(1/fvec(j))
+ print *, ' DT_FAC (lower) : ', sngl(abs(dtau_w(j))), &
+ sngl(1./(DT_FAC * fvec(j))), abs(dtau_w(j)) > 1./(DT_FAC * fvec(j))
+ print *, ' ERR_FAC (lower) : ', sngl(err_dt(j)), &
+ sngl(1./(ERR_FAC * fvec(j))), err_dt(j) > 1./(ERR_FAC * fvec(j))
+ print *, ' DT_MAX_SCALE (lower) : ', sngl(abs(dtau_w(j))), &
+ sngl(DT_MAX_SCALE*abs(tshift)), abs(dtau_w(j)) > DT_MAX_SCALE*abs(tshift)
endif
endif
enddo
@@ -1257,14 +1203,19 @@
!==============================================================================
! subroutines used in mtm_measure() and mtm_adj()
!==============================================================================
+
+ subroutine interpolate_dat_and_syn(data, syn, tstart, tend, t0, dt, NPT, &
+ dat_win, syn_win, nlen, istart)
- subroutine interpolate_dat_and_syn(data, syn, tstart, tend, t0, dt, NPT, dat_win, syn_win, nlen, istart)
+ ! extract windowed portion from original data/syn traces
implicit none
+
double precision, dimension(NPT), intent(in) :: data, syn
- double precision, dimension(NPT), intent(out) :: dat_win, syn_win
double precision, intent(in) :: tstart, tend, t0, dt
integer, intent(in) :: NPT
+
+ double precision, dimension(NPT), intent(out) :: dat_win, syn_win
integer, intent(out) :: nlen, istart
integer :: ii, i
@@ -1276,17 +1227,20 @@
! limits array bounds
if( nlen > NPT ) nlen = NPT
+ ! move checking inside subroutine
+ if (nlen <= 1) stop 'Check the length of the data and syn arrays'
+
do i = 1, nlen
time = tstart + (i-1) * dt
ii = floor((time-t0)/dt) + 1
- ! checks out-of-bounds
+ ! checks out-of-bounds (very unlikely event!)
if( ii >= NPT ) cycle
t1 = floor((time-t0)/dt) * dt + t0
- dat_win(i) = data(ii) + (data(ii+1)-data(ii)) * (time-t1) / dt ! data
- syn_win(i) = syn(ii) + (syn(ii+1)-syn(ii)) * (time-t1) /dt ! syn
+ dat_win(i) = data(ii) + (data(ii+1)-data(ii)) * (time-t1) / dt
+ syn_win(i) = syn(ii) + (syn(ii+1)-syn(ii)) * (time-t1) /dt
enddo
@@ -1417,7 +1371,7 @@
syn_vtw_cc(1) = (syn_dtw_cc(2) - syn_dtw_cc(1)) / dt
syn_vtw_cc(nlen) = (syn_dtw_cc(nlen) - syn_dtw_cc(nlen-1)) / dt
- ! estimated uncertainty in cross-correlation travltime and amplitude
+ ! estimated uncertainty in cross-correlation traveltime and amplitude
sigma_dt_top = sum( (data_dtw(1:nlen) - syn_dtw_cc(1:nlen) )**2 )
sigma_dlnA_top = sigma_dt_top
sigma_dt_bot = sum( syn_vtw_cc(1:nlen)**2 )
@@ -1425,17 +1379,15 @@
sigma_dt = sqrt( sigma_dt_top / sigma_dt_bot )
sigma_dlnA = sqrt( sigma_dlnA_top / sigma_dlnA_bot )
- if(0==1) then
- print *, ' sigma_dt : ', sigma_dt
- print *, ' sigma_dlnA : ', sigma_dlnA
- open(88,file='tshift.dat')
- do i = 1,nlen
- write(88,'(5e18.6)') (i-1)*dt, data_dtw(i), syn_dtw_cc(i), syn_dtw_cc_dt(i), syn_vtw_cc(i)
- enddo
- close(88)
- stop 'testing'
- endif
-
+ ! testing
+!!$ print *, ' sigma_dt : ', sigma_dt
+!!$ print *, ' sigma_dlnA : ', sigma_dlnA
+!!$ open(88,file='tshift.dat')
+!!$ do i = 1,nlen
+!!$ write(88,'(5e18.6)') (i-1)*dt, data_dtw(i), syn_dtw_cc(i), syn_dtw_cc_dt(i), syn_vtw_cc(i)
+!!$ enddo
+!!$ close(88)
+
! make final adjustments to uncertainty estimate
if (ERROR_TYPE == 0) then
! set uncertainty factors to 1 if you do not want to incorporate them
@@ -1455,25 +1407,26 @@
! ---------------------------------------------------------------------------
- subroutine write_average_meas(filename, imeas, dtau_meas, dlnA_meas, dtau_sigma, dlnA_sigma)
+ subroutine write_average_meas(filename, is_mtm, dtau_meas, dlnA_meas, dtau_sigma, dlnA_sigma)
implicit none
character(len=*), intent(in) :: filename
double precision, intent(in) :: dtau_meas, dlnA_meas, dtau_sigma, dlnA_sigma
- integer, intent(in) :: imeas
+ integer, intent(in) :: is_mtm
+
character(len=40) :: stlab, suffix
- if ( imeas == 7 .or. imeas == 8 ) then
- stlab = 'Multitaper' ; suffix = 'average'
+ if ( is_mtm == 1 ) then
+ stlab = 'multitaper averaged' ; suffix = 'average'
else
- stlab = 'Cross-correlation' ; suffix = 'cc'
- endif
+ stlab = 'cross-correlation' ; suffix = 'cc'
+ endif
- if ( imeas .ge. 3 ) then
+ if ( is_mtm > 0) then ! CC or MT
if (DISPLAY_DETAILS) then
- print *, trim(stlab)//' average measurements:'
- print *, ' traveltime :', sngl(dtau_meas), ' +/- ', sngl(dtau_sigma)
- print *, ' amplitude :', sngl(dlnA_meas), ' +/- ', sngl(dlnA_sigma)
+ print *, ' write ', trim(stlab)//' measurements:'
+ print *, ' traveltime :', sngl(dtau_meas), ' +/- ', sngl(dtau_sigma)
+ print *, ' amplitude :', sngl(dlnA_meas), ' +/- ', sngl(dlnA_sigma)
endif
! write average error estimates to file
@@ -1496,24 +1449,30 @@
! The transfer function maps the synthetics to the CC-deconstructed data;
! the CC measurements then need to be applied to match the original data.
+ ! dummy: fnum
implicit none
character(len=*), intent(in) :: filename
complex*16, intent(in) :: trans(:)
double precision, intent(in) :: wvec(:), df, tshift, dlnA
integer, intent(in) :: fnum, i_right, idf_new
+
double precision, dimension(:), intent(out) :: phi_wt, abs_wt, dtau_wt, dlnA_wt
double precision, intent(out), optional :: dtau_wa, dlnA_wa
integer :: i, j
double precision, dimension(NPT) :: fr
double precision :: smth, smth1, smth2 ! f0
+ logical :: ioactive
+ ioactive = .false. ! IO off
+ if (present(dtau_wa) .and. present(dlnA_wa)) ioactive = .true.
+
abs_wt(:) = 0.
phi_wt(:) = 0.
! note that with the idf_new value, these files are SUB-SAMPLED
- if (OUTPUT_MEASUREMENT_FILES) then
+ if (OUTPUT_MEASUREMENT_FILES .and. ioactive) then
open(10,file=trim(filename)//'.ph')
open(20,file=trim(filename)//'.abs')
open(30,file=trim(filename)//'.dlnA')
@@ -1526,10 +1485,11 @@
phi_wt(i) = atan2( aimag(trans(i)) , real(trans(i)) )
abs_wt(i) = abs(trans(i))
fr(i) = df*(i-1)
- if (mod(i,idf_new).eq.0 .and. OUTPUT_MEASUREMENT_FILES) then
+ if (mod(i,idf_new).eq.0 .and. OUTPUT_MEASUREMENT_FILES .and. ioactive) then
write(10,*) fr(i), phi_wt(i)
write(20,*) fr(i), abs_wt(i)
- write(30,*) fr(i), log(abs_wt(i))
+! LQY: should not we apply dlnA from cc before this?
+! write(30,*) fr(i), log(abs_wt(i))
endif
enddo
@@ -1543,13 +1503,15 @@
smth1 = (phi_wt(i+1) + TWOPI) + phi_wt(i-1) - 2.0 * phi_wt(i)
smth2 = (phi_wt(i+1) - TWOPI) + phi_wt(i-1) - 2.0 * phi_wt(i)
if(abs(smth1).lt.abs(smth).and.abs(smth1).lt.abs(smth2).and. abs(phi_wt(i) - phi_wt(i+1)) > PHASE_STEP)then
- if (DISPLAY_DETAILS) print *, 'phase correction : 2 pi', sngl(fr(i)), sngl(phi_wt(i) - phi_wt(i+1))
+ if (DISPLAY_DETAILS .and. ioactive) &
+ print *, ' phase correction : 2 pi', sngl(fr(i)), sngl(phi_wt(i) - phi_wt(i+1))
do j = i+1, i_right
phi_wt(j) = phi_wt(j) + TWOPI
enddo
endif
if(abs(smth2).lt.abs(smth).and.abs(smth2).lt.abs(smth1).and. abs(phi_wt(i) - phi_wt(i+1)) > PHASE_STEP)then
- if (DISPLAY_DETAILS) print *, 'phase correction : - 2 pi', sngl(fr(i)), sngl(phi_wt(i) - phi_wt(i+1))
+ if (DISPLAY_DETAILS .and. ioactive) &
+ print *, ' phase correction : - 2 pi', sngl(fr(i)), sngl(phi_wt(i) - phi_wt(i+1))
do j = i+1, i_right
phi_wt(j) = phi_wt(j) - TWOPI
enddo
@@ -1559,17 +1521,17 @@
! add the CC measurements to the transfer function
if (i > 1) dtau_wt(i) = (-1./wvec(i)) * phi_wt(i) + tshift
dlnA_wt(i) = log(abs_wt(i)) + dlnA
- !dlnA_wt(i) = log(abs_wt(i))
- !!dlnA_wt(i) = abs_wt(i) - 1.
+ !!dlnA_wt(i) = abs_wt(i) - 1. + dlnA
- if(mod(i,idf_new).eq.0 .and. OUTPUT_MEASUREMENT_FILES) then
+ if(mod(i,idf_new).eq.0 .and. OUTPUT_MEASUREMENT_FILES .and. ioactive) then
+ write(30,*) fr(i), dlnA_wt(i)
write(40,*) fr(i), phi_wt(i)
write(50,*) fr(i), dtau_wt(i)
endif
enddo
- if (OUTPUT_MEASUREMENT_FILES) then
+ if (OUTPUT_MEASUREMENT_FILES .and. ioactive) then
close(10)
close(20)
close(30)
@@ -1599,11 +1561,13 @@
subroutine deconstruct_dat_cc(filename,dat_dtw,tstart,dt,nlen,&
ishift,tshift,dlnA,dat_dtw_cc)
-
+
! Using CC measurements, map the data to the synthetics;
! because the windows are picked based on the synthetics,
! we apply the transfer function from the synthetics to the
! CC-deconstructed data.
+ ! dummy inputs: tshift, filename, dt
+
implicit none
character(len=*), intent(in) :: filename
double precision, dimension(NPT), intent(in) :: dat_dtw
@@ -1637,6 +1601,7 @@
! reconstruct the synthetics using cross-correlation measurements:
! (1) apply dT to get syn_dtw_cc_dt
! (2) apply dT and dlnA to get syn_dtw_cc
+ ! dt, tshift are dummy
implicit none
double precision, dimension(NPT), intent(in) :: syn_dtw
integer, intent(in) :: ishift, nlen
@@ -1674,6 +1639,7 @@
double precision, dimension(:), intent(in) :: wvec, dtau_wt, dlnA_wt
integer, intent(in) :: i_right, nlen
double precision, intent(in) :: tstart, dt
+
double precision, dimension(:), intent(out) :: syn_dtw_mt, syn_dtw_mt_dt
complex*16, dimension(NPT) :: wseis_recon
@@ -1703,10 +1669,6 @@
if (OUTPUT_MEASUREMENT_FILES) then
call dwsac1(trim(filename)//'.recon_syn.sac',syn_dtw_mt,nlen,tstart,dt)
call dwsac1(trim(filename)//'.recon_syn_dt.sac',syn_dtw_mt_dt,nlen,tstart,dt)
-!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn.sac',tstart,nlen,syn_dtw_mt)
-!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_dt.sac',tstart,nlen,syn_dtw_mt_dt)
-!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn',tstart,dt,nlen,syn_dtw_mt)
-!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn_dt',tstart,dt,nlen,syn_dtw_mt_dt)
endif
end subroutine reconstruct_syn
@@ -1792,6 +1754,8 @@
! saves interpolated values to output trace
syn(1:npt2) = syn1(1:npt2)
+ ! LQY: zero out any thing beyond npts
+ if (npt1 > npt2) syn(npt2+1:npt1)=0.
end subroutine interpolate_syn
@@ -1805,11 +1769,14 @@
double precision :: Wt
integer :: i !,imax
+ if (2*itmax > npt) stop 'Check taper_start of adjoint source'
+
!imax = maxloc(abs(syn),dim=1) ! index of the max value
!Wt = TWOPI / (2.0*(imax-1)) ! period of the taper
+
Wt = TWOPI / (2.0*(itmax-1)) ! period of the taper
- if(DISPLAY_DETAILS) print *, 'tapering start of record from index 1 to index ', itmax
+ if(DISPLAY_DETAILS) print *, 'tapering start of adjoint source from index 1 to index ', itmax
! apply a cosine taper from the start to the max value,
! such that the starting point is exactly zero
@@ -1835,6 +1802,7 @@
OUT_DIR = 'OUTPUT_FILES' ! default
open(10,file='MEASUREMENT.PAR',status='old',iostat=ios)
+ if ( ios /= 0) stop 'Error opening MEASUREMENT.PAR file'
read(10,*) tt,dtt,nn
read(10,*) imeas0
read(10,*) chan
@@ -1861,27 +1829,27 @@
! check the read-in values
print *, 'INPUTS FROM MEASUREMENT.PAR :'
- print *, ' tt, dtt, nn : ',tt,dtt,nn
+ print *, ' tt, dtt, nn : ',sngl(tt),sngl(dtt),nn
print *, ' imeas : ',imeas
print *, ' chan : ',chan
- print *, ' TLONG, TSHORT : ',TLONG, TSHORT
+ print *, ' TLONG, TSHORT : ',sngl(TLONG), sngl(TSHORT)
fstart0 = 1./TLONG ; fend0 = 1./TSHORT
- print *, ' fstart, fend : ', fstart0, fend0
+ print *, ' fstart, fend : ', sngl(fstart0), sngl(fend0)
print *, ' RUN_BANDPASS : ',RUN_BANDPASS
print *, ' DISPLAY_DETAILS : ',DISPLAY_DETAILS
print *, ' OUTPUT_MEASUREMENT_FILES : ',OUTPUT_MEASUREMENT_FILES
print *, ' COMPUTE_ADJOINT_SOURCE : ',COMPUTE_ADJOINT_SOURCE
- print *, ' TSHIFT_MIN, TSHIFT_MAX : ',TSHIFT_MIN, TSHIFT_MAX
- print *, ' DLNA_MIN, DLNA_MAX : ',DLNA_MIN, DLNA_MAX
- print *, ' CC_MIN : ',CC_MIN
+ print *, ' TSHIFT_MIN, TSHIFT_MAX : ',sngl(TSHIFT_MIN), sngl(TSHIFT_MAX)
+ print *, ' DLNA_MIN, DLNA_MAX : ',sngl(DLNA_MIN), sngl(DLNA_MAX)
+ print *, ' CC_MIN : ',sngl(CC_MIN)
print *, ' ERROR_TYPE : ',ERROR_TYPE
- print *, ' DT_SIGMA_MIN : ',DT_SIGMA_MIN
- print *, ' DLNA_SIGMA_MIN : ',DLNA_SIGMA_MIN
+ print *, ' DT_SIGMA_MIN : ',sngl(DT_SIGMA_MIN)
+ print *, ' DLNA_SIGMA_MIN : ',sngl(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 *, ' WTR, NPI : ',sngl(WTR),NPI
+ print *, ' DT_FAC : ',sngl(DT_FAC)
+ print *, ' ERR_FAC : ',sngl(ERR_FAC)
+ print *, ' DT_MAX_SCALE : ',sngl(DT_MAX_SCALE)
print *, ' NCYCLE_IN_WINDOW : ',NCYCLE_IN_WINDOW
!stop 'checking PAR file input'
@@ -1942,35 +1910,30 @@
if( DO_RAY_DENSITY_SOURCE ) ERROR_TYPE = 0
! assign additional parameters and stop for certain inconsistencies
- if (fstart0.ge.fend0) then
- print *, 'Check input frequency range of the signal'
- stop
- endif
+ if (fstart0 >= fend0) &
+ stop 'Check input frequency range of the signal'
+
+ if (nn > NDIM) &
+ stop 'Error: Change interpolation nn or NDIM'
- if (nn > NDIM) then
- print *, 'Error: Change interpolation nn or NDIM'
- stop
- endif
+ ! LQY: what happens if imeas = 7/8, and itaper = 2,3, is that permitted?
+ ! LQY: how about imeas = 1/2, itaper = 1,2,3 matters?
- ! for CC kernels, ITAPER must be a single taper (2 or 3)
- if ( (ITAPER==1) .and. ((imeas.ge.3).and.(imeas.le.6)) ) then
- print *, 'Error: Change ITAPER to 2 or 3'
- stop
- endif
-
- if ( (imeas==1).or.(imeas==2) ) then
+ if ( imeas == 1 .or. imeas == 2 ) then ! waveforms
is_mtm0 = 0
- elseif ( (imeas.ge.3).and.(imeas.le.6) ) then
- is_mtm0 = ITAPER ! 2 or 3
- elseif ( (imeas==7).or.(imeas==8) ) then
+ elseif ( imeas >= 3 .and. imeas <= 6 ) then ! CC Dt/DlnA
+ ! for CC kernels, ITAPER must be a single taper (2 or 3)
+ if ( ITAPER == 1 ) stop 'Error: Change ITAPER to 2/3 for CC measurements'
+ is_mtm0 = ITAPER ! 2 or 3 for CC adjoint sources
+ elseif ( imeas==7 .or. imeas==8 ) then
is_mtm0 = 1 ! multitaper required for MT adjoint source
else
- print *, 'Error: imeas must by 1-8'
- stop
+ stop 'Error: imeas must by 1-8'
endif
is_mtm = is_mtm0
print *, ' is_mtm :',is_mtm
+ print *, ' '
end subroutine read_par_file
@@ -1984,14 +1947,17 @@
integer,intent(out):: yr,jda,ho,mi
double precision,intent(out):: sec,dist,az,baz,slat,slon
- !real*8,intent(out):: sec,dist,az,baz,slat,slon
character(len=*),intent(out) :: ntw,sta,comp
- !real*8 :: tmp
real :: tmp
- integer :: nsec,msec !,i,klen
+ integer :: nsec,msec
integer :: nerr
+ ! note here data_file is a dummy argument, and we rely on the earlier
+ ! call to rsac1() to retain header and waveform info in computer memory
+ ! for data_file trace.
+
+
! integer header variables
call getnhv('nzyear',yr,nerr)
call getnhv('nzjday',jda,nerr)
@@ -2079,4 +2045,214 @@
end subroutine get_sacfile_header
+!!==================================================
+
+subroutine setup_weighting(chan_syn)
+ !
+ ! determines weights based on number of window picks on Z/R/T components
+ !
+ use ma_weighting
+
+ use ma_constants,only: NDIM
+! use ma_sub,only: get_sacfile_header,drsac1
+ use ma_sub2,only: TOL
+
+ implicit none
+ character(len=10),intent(in) :: chan_syn
+
+ ! local parameters
+ integer :: npairs,ios,ipair,iposition,ipicks
+ character(len=150) :: datafile,synfile !,dummy
+ character(len=4) :: comp_T,comp_Z,comp_R
+ integer :: picks_T, picks_Z, picks_R,npicks
+ ! sac header information
+ integer :: yr,jda,ho,mi
+ double precision :: sec,dist,az,baz,slat,slon,T_surfacewaves
+ character(len=10) :: net,sta,chan_dat,chan,cmp
+ double precision :: t01, dt1, t02, dt2, t0, dt, tstart, tend
+ integer :: npt1, npt2, npts
+ double precision, dimension(NDIM) :: data, syn
+
+ ! initializes
+ picks_R = 0
+ picks_Z = 0
+ picks_T = 0
+
+ num_P_SV_V = 0.d0
+ num_P_SV_R = 0.d0
+ num_SH_T = 0.d0
+
+ num_Rayleigh_V = 0.d0
+ num_Rayleigh_R = 0.d0
+ num_Love_T = 0.d0
+
+ ! substrings (synthetics components)
+ comp_T = trim(chan_syn)//"T."
+ comp_R = trim(chan_syn)//"R."
+ comp_Z = trim(chan_syn)//"Z."
+
+ ! opens measurement windows
+ open(21,file='MEASUREMENT.WINDOWS',status='old',iostat=ios)
+ if (ios /= 0) stop 'Error opening input file: MEASUREMENT WINDOWS'
+ read(21,*,iostat=ios) npairs
+ if (ios /= 0) stop 'Error reading number of pairs of data/syn'
+ ! loops through windows
+ do ipair=1,npairs
+
+ ! reads in file names
+ read(21,'(a)',iostat=ios) datafile
+ if (ios /= 0) stop 'Error reading windows datafile'
+ read(21,'(a)',iostat=ios) synfile
+ if (ios /= 0) stop 'Error reading windows synfile'
+
+ ! read data and syn (read datafile last to take its header later)
+ call drsac1(synfile,syn,npt2,t02,dt2)
+ call drsac1(datafile,data,npt1,t01,dt1)
+
+ if (max(npt1,npt2) > NDIM) &
+ stop 'Error: Too many npts in data or syn'
+
+ ! check if t0 and dt match
+ if (abs(dt1-dt2) > TOL) stop 'Error: check if dt match'
+ dt = dt1
+ npts = min(npt1,npt2)
+ if (abs(t01-t02) > dt) then
+ print *,'data t0: ',t01
+ print *,'syn t0: ',t02
+ stop 'Check if t0 match'
+ endif
+ t0 = t01
+
+ ! figure out station/network/comp names, etc
+ call get_sacfile_header(trim(datafile),yr,jda,ho,mi,sec,net,sta, &
+ chan_dat,dist,az,baz,slat,slon)
+ chan = chan_dat
+ cmp = chan_dat(3:3)
+
+ ! theoretical surface wave arrival time
+ T_surfacewaves = dist / surface_vel
+
+ ! debug output
+ !if (DISPLAY_DETAILS) then
+ !print*,'debug: '
+ !print*,' yr,jda,ho,mi,sec : ',yr,jda,ho,mi,sec
+ !print*,' net,sta,chan_dat : ',net,sta,chan_dat
+ !print*,' dist,az,baz,slat,slon : ',dist,az,baz,slat,slon
+ !print*,' cmp = ',cmp
+ !print*,' dist = ',dist
+ !print*,' T_surfacewaves = ',T_surfacewaves
+ !print*
+ !endif
+
+ ! reads in window picks
+ read(21,*,iostat=ios) npicks
+ if (ios /= 0) stop 'Error reading windows npicks'
+
+ ! loops/skips over picks (start/end times)
+ do ipicks=1,npicks
+
+ read(21,*,iostat=ios) tstart, tend
+ if (ios /= 0) stop 'Error reading window pick: tstart and tend'
+
+ tstart = max(tstart,t0)
+ tend = min(tend, t0+(npts-1)*dt)
+
+ ! body wave picks
+ if( tend <= T_surfacewaves ) then
+ if( cmp(1:1) == "Z" ) num_P_SV_V = num_P_SV_V + 1.d0
+ if( cmp(1:1) == "R" ) num_P_SV_R = num_P_SV_R + 1.d0
+ if( cmp(1:1) == "T" ) num_SH_T = num_SH_T + 1.d0
+ else
+ ! surface wave picks
+ if( cmp(1:1) == "Z" ) num_Rayleigh_V = num_Rayleigh_V + 1.d0
+ if( cmp(1:1) == "R" ) num_Rayleigh_R = num_Rayleigh_R + 1.d0
+ if( cmp(1:1) == "T" ) num_Love_T = num_Love_T + 1.d0
+ endif
+
+ enddo
+
+ ! determines all picks on a trace component
+ ! (also cross-check comp name in filename)
+ ! transverse
+ iposition = INDEX( trim(synfile), comp_T, .false. )
+ if( iposition > 3 .and. iposition < len_trim( synfile) ) then
+ if( cmp(1:1) /= "T" ) stop 'error T component pick'
+ picks_T = picks_T + npicks
+ else
+ ! radial
+ iposition = INDEX( trim(synfile), comp_R, .false. )
+ if( iposition > 3 .and. iposition < len_trim( synfile) ) then
+ if( cmp(1:1) /= "R" ) stop 'error R component pick'
+ picks_R = picks_R + npicks
+ else
+ ! vertical
+ iposition = INDEX( trim(synfile), comp_Z, .false. )
+ if( iposition > 3 .and. iposition < len_trim( synfile) ) then
+ if( cmp(1:1) /= "Z" ) stop 'error Z component pick'
+ picks_Z = picks_Z + npicks
+ endif
+ endif
+ endif
+
+ enddo
+ close(21)
+
+
+ ! check with total number of picks per component
+ if( nint( num_P_SV_R + num_Rayleigh_R ) /= picks_R ) stop 'error R picks'
+ if( nint( num_P_SV_V + num_Rayleigh_V ) /= picks_Z ) stop 'error Z picks'
+ if( nint( num_SH_T + num_Love_T ) /= picks_T ) stop 'error T picks'
+
+ if( DISPLAY_DETAILS ) then
+ print*
+ print*,'weighting measurements: '
+ print*,' picks T:',picks_T
+ print*,' picks R:',picks_R
+ print*,' picks Z:',picks_Z
+ print*
+ print*,' picks P_SV_R: ',nint(num_P_SV_R)
+ print*,' picks P_SV_V: ',nint(num_P_SV_V)
+ print*,' picks SH_T : ',nint(num_SH_T)
+ print*,' picks Rayleigh_R: ',nint(num_Rayleigh_R)
+ print*,' picks Rayleigh_V: ',nint(num_Rayleigh_V)
+ print*,' picks Love_T : ',nint(num_Love_T)
+ print*
+ endif
+
+
+ ! sets up weights based on picks
+ weight_T = 1.0d0
+ weight_R = 1.0d0
+ weight_Z = 1.0d0
+
+ ! weighting tries to balance love waves (tranverse) versus rayleigh waves (radial + vertical)
+ !if( picks_T > 0 ) then
+ ! if( picks_R + picks_Z > 0 ) weight_T = dble(picks_R + picks_Z)/dble(picks_T)
+ !endif
+
+ ! use normalization as weights
+ if( picks_T > 0 ) weight_T = 1.d0 / picks_T
+ if( picks_R > 0 ) weight_R = 1.d0 / picks_R
+ if( picks_Z > 0 ) weight_Z = 1.d0 / picks_Z
+
+ ! use normalization (no traces means zero weights)
+ if( num_P_SV_R > 0. ) num_P_SV_R = 1.d0 / num_P_SV_R
+ if( num_P_SV_V > 0. ) num_P_SV_V = 1.d0 / num_P_SV_V
+ if( num_SH_T > 0. ) num_SH_T = 1.d0 / num_SH_T
+ if( num_Rayleigh_R > 0. ) num_Rayleigh_R = 1.d0 / num_Rayleigh_R
+ if( num_Rayleigh_V > 0. ) num_Rayleigh_V = 1.d0 / num_Rayleigh_V
+ if( num_Love_T > 0. ) num_Love_T = 1.d0 / num_Love_T
+
+ print*,' weight of P_SV_R:',num_P_SV_R
+ print*,' weight of P_SV_V:',num_P_SV_V
+ print*,' weight of SH_T :',num_SH_T
+ print*,' weight of Rayleigh_R:',num_Rayleigh_R
+ print*,' weight of Rayleigh_V:',num_Rayleigh_V
+ print*,' weight of Love_T :',num_Love_T
+ print*
+
+end subroutine setup_weighting
+
+
+
end module ma_sub
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/measure_adj.f90 2012-11-27 02:07:13 UTC (rev 21074)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/measure_adj.f90 2012-11-27 17:05:45 UTC (rev 21075)
@@ -18,7 +18,7 @@
! synthetics is used in constructing the adjoint source.
! 7. imeas = 7, multitaper traveltime difference for an event kernel. The measurement between
! data and synthetics is used in constructing the adjoint source. See multitaper_notes.pdf.
- ! 8. imeas = 8, multitaper ampltidue difference for an event kernel. The measurement between
+ ! 8. imeas = 8, multitaper amplitude difference for an event kernel. The measurement between
! data and synthetics is used in constructing the adjoint source. See multitaper_notes.pdf.
use ma_variables
@@ -42,7 +42,7 @@
double precision :: sec,dist,az,baz,slat,slon
character(len=10) :: net,sta,chan_dat,chan,cmp,chan_syn
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 :: all_chi, 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
complex*16, dimension(NPT) :: trans_mtm
@@ -51,27 +51,24 @@
logical :: use_trace
!double precision :: trbdndw, a
!integer :: iord, passes
-
integer :: ipick_type
double precision :: T_surfacewaves
!********* PROGRAM STARTS HERE *********************
-
! read in MEASUREMENT.PAR (see ma_sub.f90 and write_par_file.pl)
- ! most variables are global (see ma_variables.f90)
+ ! most parameters are global (see ma_variables.f90)
call read_par_file(fstart0,fend0,tt,dtt,nn,chan)
! uses weights to balance love and rayleigh measurements
! we do a normalization of P_SV, P_SH, Love, Rayleigh with the number of measurement picks
- if( DO_WEIGHTING ) then
- call setup_weighting(chan)
- endif
+ if( DO_WEIGHTING ) call setup_weighting(chan)
! 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
+ if (ios /= 0) stop 'Error opening input file: MEASUREMENT WINDOWS'
+
read(11,*,iostat=ios) npairs
- if (ios /= 0) then ; print *, 'Error reading number of pairs of data/syn' ; stop ; endif
+ if (ios /= 0) stop 'Error reading number of pairs of data/syn'
print *, 'reading in the data and synthetics...'
@@ -79,7 +76,7 @@
open(12,file='window_index',status='unknown',iostat=ios)
open(13,file='window_chi',status='unknown',iostat=ios)
- nwin = 0
+ nwin = 0; all_chi=0.
do ipair = 1, npairs
data(:) = 0.0
@@ -90,61 +87,50 @@
! reads in file names for data and synthetics
read(11,'(a)',iostat=ios) datafile
- if (ios /= 0) then ; print *, 'Error reading windows file' ; stop ; endif
+ if (ios /= 0) stop 'Error reading windows file'
read(11,'(a)',iostat=ios) synfile
- if (ios /= 0) then ; print *, 'Error reading windows file' ; stop ; endif
+ if (ios /= 0) stop 'Error reading windows file'
- ! read data and syn
- call drsac1(datafile,data,npt1,t01,dt1)
+ ! read data and syn (in double precision)
+ ! LQY: data is read last to be stored in memory for header retrieval later
call drsac1(synfile,syn,npt2,t02,dt2)
- !print *, npt1,t01,dt1,NDIM ! check: double precision
+ call drsac1(datafile,data,npt1,t01,dt1)
- ! user output
- print *
- print *, 'data: ',trim(datafile)
- print *, ' min/max: ',minval(data(:)),maxval(data(:))
- !print *, ' ',data(1:5)
-
- print *, 'syn: ',trim(synfile)
- print *, ' min/max: ',minval(syn(:)),maxval(syn(:))
- !print *, ' ',syn(1:5)
-
- if (max(npt1,npt2) > NDIM) then
- print *, 'Error: Too many number of points in data or syn'
- stop
+ if (DISPLAY_DETAILS) then
+ print *
+ print *, 'data: ',trim(datafile)
+ print *, ' min/max: ',sngl(minval(data(:))),sngl(maxval(data(:)))
+
+ print *, 'syn: ',trim(synfile)
+ print *, ' min/max: ',sngl(minval(syn(:))),sngl(maxval(syn(:)))
endif
- ! check if t0 and dt match
- if (abs(dt1-dt2) > TOL) then ; print *, 'Error: check if dt match' ; stop ; endif
- dt = dt1
+ ! check if npts, dt, t0 matches for data and synthetics
+ ! no interpolation is necessary at any point
+ if (max(npt1,npt2) > NDIM) &
+ stop 'Error: Too many number of points in data or syn'
npts = min(npt1,npt2)
- ! check
- if (abs(t01-t02) > dt) then
- print*,'data t0: ',t01
- print*,'syn t0: ',t02
- stop 'Check if t0 match'
- else
- print*,' time:',t01
- print*,' dt:',dt
- print*,' npts: ',npts
- endif
+ if (abs(dt1-dt2) > TOL) stop 'Error: check if dt match'
+ dt = dt1
- ! apply bandpass filter to data and synthetics, if desired
+ if (abs(t01-t02) > dt) stop 'Check if t0 match'
+ t0 = t01
+
+ if (DISPLAY_DETAILS) print *,' time, dt, npts :',sngl(t01), sngl(dt), npts
+
+
+ ! apply bandpass filter to data and synthetics with saclib, if desired
! http://www.iris.washington.edu/pipermail/sac-help/2008-March/000376.html
! Access to the kidate, xapiir, and getfil is not simple and not
! supported under the current state of the SAC code base.
- t0 = t01
- if(.not. RUN_BANDPASS) then
- t0 = t01
- else
- !call interpolate_syn(syn,t02,dt,npt2,t01,dt,npts)
- t0 = t01
+ if(RUN_BANDPASS) then
call bandpass(data,npts,dt,fstart0,fend0)
+ call bandpass(syn,npts,dt,fstart0,fend0)
endif
- ! figure out station name, network name, component name, etc
+ ! find out station/network/comp names,etc
call get_sacfile_header(trim(datafile),yr,jda,ho,mi,sec,net,sta, &
chan_dat,dist,az,baz,slat,slon)
@@ -155,60 +141,36 @@
cmp = chan_dat(3:3)
chan_syn = trim(chan)//trim(cmp)
+ ! example: OUT/PAS.CI.BHZ
file_prefix0 = trim(sta)//'.'//trim(net)//'.'//trim(chan_syn)
file_prefix2 = trim(OUT_DIR)//'/'//trim(file_prefix0)
print *
- print *, trim(file_prefix2), ' --- '
+ 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', imeas0
-!!$ if (imeas == 0) then
-!!$ adj_file_prefix = trim(file_prefix2) // '.ik'
-!!$ else if (imeas == 1) then
-!!$ adj_file_prefix = trim(file_prefix2) // '.mtm'
-!!$ else if (imeas == 2) then
-!!$ adj_file_prefix = trim(file_prefix2) // '.cc'
-!!$ else if (imeas == 3) then
-!!$ adj_file_prefix = trim(file_prefix2) // '.bdcc'
-!!$ else
-!!$ print *, 'imeas = ', imeas
-!!$ print *, 'Error: imeas must be 0, 1, 2, or 3'
-!!$ stop
-!!$ endif
! reads number of measurement windows
read(11,*,iostat=ios) num_meas
- if (ios /= 0) then ; print *, 'Error reading num_meas' ; stop ; endif
+ if (ios /= 0) stop 'Error reading num_meas'
do j = 1, num_meas
- ! reads in start and end time of measurement window
+ ! reads in start and end time of the measurement window
read(11,*,iostat=ios) tstart, tend
- if (ios /= 0) then ; print *, 'Error reading tstart and tend' ; stop ; endif
+ if (ios /= 0) stop 'Error reading tstart and tend'
! checks start and end times of window compared to trace lengths
tstart = max(tstart,t0)
tend = min(tend, t0+(npts-1)*dt)
- nlen = floor((tend-tstart)/dt) + 1 ! see subroutine interpolate_data_and_syn
+ nlen = floor((tend-tstart)/dt) + 1 ! dummy, replaced later in mt_measure()
- ! body wave picks
- ipick_type = 0
- if( tend <= T_surfacewaves ) then
- if( cmp(1:1) == "Z" ) ipick_type = P_SV_V
- if( cmp(1:1) == "R" ) ipick_type = P_SV_R
- if( cmp(1:1) == "T" ) ipick_type = SH_T
- else
- ! surface wave picks
- if( cmp(1:1) == "Z" ) ipick_type = Rayleigh_V
- if( cmp(1:1) == "R" ) ipick_type = Rayleigh_R
- if( cmp(1:1) == "T" ) ipick_type = Love_T
- endif
-
! write values to output file
nwin = nwin + 1 ! overall window counter
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
+ ! add taper type to file prefix: OUT/PAS.CI.BHZ.01.mtm
write(file_prefix,'(a,i2.2)') trim(file_prefix2)//'.', j
+
if (is_mtm == 1) then
measure_file_prefix = trim(file_prefix) // '.mtm' ! multitaper taper
elseif (is_mtm == 2) then
@@ -218,7 +180,7 @@
endif
print *
- print *, ' Measurements No.', j, ' ... '
+ print *, ' Measurement window No.', j, ' ... '
! initialize the measurements
window_chi(:) = 0.
@@ -233,45 +195,41 @@
window_chi(19) = 0.5 * sum( (data-syn)**2 )
window_chi(20) = npts*dt
- ! get the starting frequency to avoid measuring long periods not present
- !fstart = fstart0 ; fend = fend0
- !call mt_measure_select_0(nlen,dt,i_left,i_right,fstart,fend,use_trace)
-
! make measurements
! 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_right = i_right0
- i_left = 1
+ i_left = 1 ! LQY: is it feasible that i_left is not 1? mt_adj() inherently assumes it.
! period of the max power of the synthetic record
T_pmax_dat = (dt*NPT) / dble(i_pmax_dat)
T_pmax_syn = (dt*NPT) / dble(i_pmax_syn)
- ! adjust measurements for MT adjoint source
+ ! adjust frequency ranges for MT measurements
+ ! fstart is constrained by NCYCLE_IN_WINDOW/tlen, fend constrained by i_right
if (is_mtm == 1) then
fstart = fstart0 ; fend = fend0
call mt_measure_select(nlen,tshift,i_pmax_syn,dtau_w,err_dt, &
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
+ print *, ' Tlong/Tshort (input) :', sngl(1/fstart0), sngl(1/fend0)
+ print *, ' Tlong/Tshort (adjusted) :', sngl(1/fstart), sngl(1/fend)
+ print *, ' period of max data/syn power :', sngl(T_pmax_dat), sngl(T_pmax_syn)
! if MT measurement window is rejected by mt_measure_select, then use a CC measurement
if(.not. use_trace) then
!stop 'Check why this MT measurement was rejected'
- print *, 'Reverting from multitaper measurement to cross-correlation measurement'
+ print *, ' reverting from MT measurement to CC measurement...'
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,&
+ is_mtm = 3 ! LQY: WHY not is_mtm = 2?
+ 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.
+ else
+ print *, ' using this MTM. '
endif
-
- else
- use_trace = .true.
endif
! check that the CC measurements are within the specified input range
@@ -286,17 +244,19 @@
endif
! compute adjoint sources and misfit function values and also the CC-reconstructed records
- if (use_trace) then
+ if (COMPUTE_ADJOINT_SOURCE) then
+ print *, ' Generating adjoint source and chi value for imeas = ', imeas
- ! banana-doughnut kernel (needs only synthetic trace)
+ ! banana-doughnut kernel (needs only synthetic trace)
+ ! LQY: what is this section intended to do?
+ ! reset imeas == 3 for adjoint sources without time shift and uncertainty scaling
+ ! (pure cross-correlation adjoint source for banana-doughnuts)
if(imeas == 5 .and. trim(datafile) == trim(synfile) ) then
- print*,'cross-correlation measurement:'
- print*,' only synthetic file: ',trim(synfile)
- print*,' without traveltime difference/uncertainty'
- print*
- ! uses imeas == 3 for adjoint sources without time shift and uncertainty scaling
- ! (pure cross-correlation adjoint source for banana-doughnuts)
- imeas = 3
+ print*,'cross-correlation measurement:'
+ print*,' only synthetic file: ',trim(synfile)
+ print*,' without traveltime difference/uncertainty'
+ print*
+ imeas = 3
endif
tr_chi = 0.0 ; am_chi = 0.0 ! must be initialized
@@ -315,54 +275,64 @@
write(13,'(a14,a8,a3,a5,i4,i4,2e14.6,20e14.6,2e14.6,2f14.6)') &
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)
+ print *, ' tr_chi = ', sngl(tr_chi), ' am_chi = ', sngl(am_chi)
! uses weighting to balance love / rayleigh measurements
if( DO_WEIGHTING ) then
- ! weights by transverse/radial/vertical
- !if( cmp(1:1) == "T") then
- ! tr_adj_src(:) = tr_adj_src(:) * weight_T
- !else
- ! if( cmp(1:1) == "R") then
- ! tr_adj_src(:) = tr_adj_src(:) * weight_R
- ! else
- ! if( cmp(1:1) == "Z") then
- ! tr_adj_src(:) = tr_adj_src(:) * weight_Z
- ! endif
- ! endif
- !endif
+ ipick_type = 0
+ if( tend <= T_surfacewaves ) then
+ ! body wave picks
+ if( cmp(1:1) == "Z" ) ipick_type = P_SV_V
+ if( cmp(1:1) == "R" ) ipick_type = P_SV_R
+ if( cmp(1:1) == "T" ) ipick_type = SH_T
+ else
+ ! surface wave picks
+ if( cmp(1:1) == "Z" ) ipick_type = Rayleigh_V
+ if( cmp(1:1) == "R" ) ipick_type = Rayleigh_R
+ if( cmp(1:1) == "T" ) ipick_type = Love_T
+ endif
+
+ ! LQY: shouldn't chi values be changed accordingly?????
+ ! No total chi value is calculated ...
+
! weights by phase types
select case(ipick_type)
case( P_SV_V )
tr_adj_src(:) = tr_adj_src(:) * num_P_SV_V
+ tr_chi = tr_chi * num_P_SV_V
case( P_SV_R )
tr_adj_src(:) = tr_adj_src(:) * num_P_SV_R
+ tr_chi = tr_chi * num_P_SV_R
case( SH_T )
tr_adj_src(:) = tr_adj_src(:) * num_SH_T
+ tr_chi = tr_chi * num_SH_T
case( Rayleigh_V )
tr_adj_src(:) = tr_adj_src(:) * num_Rayleigh_V
+ tr_chi = tr_chi * num_Rayleigh_V
case( Rayleigh_R )
tr_adj_src(:) = tr_adj_src(:) * num_Rayleigh_R
+ tr_chi = tr_chi * num_Rayleigh_R
case( Love_T )
tr_adj_src(:) = tr_adj_src(:) * num_Love_T
+ tr_chi = tr_chi * num_Love_T
case default
stop 'error ipick_type unknown'
end select
- endif
+ endif
! combine adjoint sources from different measurement windows
- if (COMPUTE_ADJOINT_SOURCE) then
- 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(:) ! imeas = 2,4,6,8
- endif
- endif
-
+ if (mod(imeas,2)==1) then
+ adj_syn_all(:) = adj_syn_all(:) + tr_adj_src(:) ! imeas = 1,3,5,7
+ all_chi = all_chi + tr_chi
+ else
+ adj_syn_all(:) = adj_syn_all(:) + am_adj_src(:) ! imeas = 2,4,6,8
+ all_chi = all_chi + am_chi
+ endif
+
! combine CC-reconstructed records
recon_cc_all(istart:istart+nlen-1) = recon_cc_all(istart:istart+nlen-1) + syn_dtw_cc(1:nlen)
- endif
+ endif ! COMPUTE_ADJOINT_SOURCE
! CHT: (re-)set to multitaper parameters, if originally specified
if (is_mtm0 == 1) then
@@ -370,13 +340,17 @@
is_mtm = is_mtm0
endif
- enddo ! nmeas
+ enddo ! nmeas
!----------------------------
- ! write out the adjoint source
+ ! write out the adjoint source for the trace (STA.NI.CMP) by combining contribution from all wins
if (COMPUTE_ADJOINT_SOURCE) then
+ ! 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)
+
! 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.
@@ -387,7 +361,7 @@
! then kernels,
! i.e. for a traveltime measurement: DeltaT = 1/N * int F(d/dt s) F(ds)
! should contain this filter as well.
- !
+ !
! when we construct the adjoint source here,it is initially a filtered version
! as well F(s_adj) since we use/depend on filtered synthetics F(s).
! however, for kernel simulations, we do run with a reconstructed forward wavefield,
@@ -404,7 +378,7 @@
! we do use a bandpass filter here again on the adjoint source. this is slightly different
! to the transfer function filter in SAC used initially to filter data & synthetics.
! but this seems to be the best and fairly easy what we can do here...
- call bandpass(adj_syn_all,npts,dt,fstart0,fend0)
+ call bandpass(adj_syn_all,npts,dt,fstart0,fend0) ! sac butterworth filter
! cut and interpolate to match time-stepping for SEM
! NOTE: This can leave a non-zero value to start the record,
@@ -423,26 +397,14 @@
call dwascii(trim(adj_file_prefix)//'.density.adj',adj_syn_all,nn,tt,dtt)
else
call dwascii(trim(adj_file_prefix)//'.adj',adj_syn_all,nn,tt,dtt)
+ ! LQY add the sum of chi values (total misfit), weights included if DO_WEIGHTING on
+ open(14,file='window_chi_sum',status='unknown')
+ write(14,*) all_chi
+ close(14)
endif
- !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
-
- ! cut and interpolate, then write ASCII file
- !call interpolate_syn(recon_cc_all,t0,dt,npts,tt,dtt,nn)
- !call dwrite_ascfile_f(trim(file_prefix2)//'.recon.cc',tt,dtt,nn,recon_cc_all)
-
- ! write SAC file
- call dwsac1(trim(file_prefix2)//'.recon.cc.sac',recon_cc_all,npts,t0,dt)
-
- !if (nerr > 0) then ; print *, 'Error writing reconstructed CC file' ; stop ; endif
- !call dwrite_sacfile_f(trim(datafile),trim(file_prefix2)//'.recon.cc',t0,npts,recon_cc_all)
-
enddo ! npairs
close(11) ! read: MEASUREMENT.WINDOWS
@@ -451,214 +413,3 @@
end program measure_adj
-subroutine setup_weighting(chan_syn)
- !
- ! determines weights based on number of window picks on radial, transverse and vertical components
- !
- use ma_weighting
-
- use ma_constants,only: NDIM
- use ma_sub,only: get_sacfile_header,drsac1
- use ma_sub2,only: TOL
-
- implicit none
- character(len=10) :: chan_syn
-
- ! local parameters
- integer :: npairs,ios,ipair,iposition,ipicks
- character(len=150) :: datafile,synfile !,dummy
- character(len=4) :: comp_T,comp_Z,comp_R
- integer :: picks_T, picks_Z, picks_R,npicks
- ! sac header information
- integer :: yr,jda,ho,mi
- double precision :: sec,dist,az,baz,slat,slon,T_surfacewaves
- character(len=10) :: net,sta,chan_dat,chan,cmp
- double precision :: t01, dt1, t02, dt2, t0, dt, tstart, tend
- integer :: npt1, npt2, npts
- double precision, dimension(NDIM) :: data, syn
-
- ! initializes
- picks_R = 0
- picks_Z = 0
- picks_T = 0
-
- num_P_SV_V = 0.d0
- num_P_SV_R = 0.d0
- num_SH_T = 0.d0
-
- num_Rayleigh_V = 0.d0
- num_Rayleigh_R = 0.d0
- num_Love_T = 0.d0
-
- ! substrings (synthetics components)
- comp_T = trim(chan_syn)//"T."
- comp_R = trim(chan_syn)//"R."
- comp_Z = trim(chan_syn)//"Z."
-
- ! opens measurement windows
- open(21,file='MEASUREMENT.WINDOWS',status='old',iostat=ios)
- if (ios /= 0) then ; print *, 'Error opening input file: MEASUREMENT WINDOWS' ; stop ; endif
- read(21,*,iostat=ios) npairs
- if (ios /= 0) then ; print *, 'Error reading number of pairs of data/syn' ; stop ; endif
-
- ! loops through windows
- do ipair=1,npairs
-
- ! reads in file names
- read(21,'(a)',iostat=ios) datafile
- if (ios /= 0) then ; print *, 'Error reading windows datafile' ; stop ; endif
- read(21,'(a)',iostat=ios) synfile
- if (ios /= 0) then ; print *, 'Error reading windows synfile' ; stop ; endif
-
- ! read data and syn
- call drsac1(datafile,data,npt1,t01,dt1)
- call drsac1(synfile,syn,npt2,t02,dt2)
-
- if (max(npt1,npt2) > NDIM) then
- print *, 'Error: Too many number of points in data or syn'
- stop
- endif
- ! check if t0 and dt match
- if (abs(dt1-dt2) > TOL) then ; print *, 'Error: check if dt match' ; stop ; endif
- dt = dt1
- npts = min(npt1,npt2)
- if (abs(t01-t02) > dt) then
- print*,'data t0: ',t01
- print*,'syn t0: ',t02
- stop 'Check if t0 match'
- endif
- t0 = t01
-
- ! 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
- cmp = chan_dat(3:3)
-
-
- ! theoretical surface wave arrival time
- T_surfacewaves = dist / surface_vel
-
- ! debug output
- !print*
- !print*,'debug: '
- !print*,' yr,jda,ho,mi,sec : ',yr,jda,ho,mi,sec
- !print*,' net,sta,chan_dat : ',net,sta,chan_dat
- !print*,' dist,az,baz,slat,slon : ',dist,az,baz,slat,slon
- !print*,' cmp = ',cmp
- !print*,' dist = ',dist
- !print*,' T_surfacewaves = ',T_surfacewaves
- !print*
-
- ! reads in window picks
- read(21,*,iostat=ios) npicks
- if (ios /= 0) then ; print *, 'Error reading windows npicks' ; stop ; endif
-
- ! loops/skips over picks (start/end times)
- do ipicks=1,npicks
- !read(21,'(a)',iostat=ios) dummy
- !if (ios /= 0) then ; print *, 'Error reading window pick' ; stop ; endif
-
- read(21,*,iostat=ios) tstart, tend
- if (ios /= 0) then ; print *, 'Error reading window pick: tstart and tend' ; stop ; endif
-
- tstart = max(tstart,t0)
- tend = min(tend, t0+(npts-1)*dt)
- !nlen = floor((tend-tstart)/dt) + 1 ! see subroutine interpolate_data_and_syn
-
- ! body wave picks
- if( tend <= T_surfacewaves ) then
- if( cmp(1:1) == "Z" ) num_P_SV_V = num_P_SV_V + 1.d0
- if( cmp(1:1) == "R" ) num_P_SV_R = num_P_SV_R + 1.d0
- if( cmp(1:1) == "T" ) num_SH_T = num_SH_T + 1.d0
- else
- ! surface wave picks
- if( cmp(1:1) == "Z" ) num_Rayleigh_V = num_Rayleigh_V + 1.d0
- if( cmp(1:1) == "R" ) num_Rayleigh_R = num_Rayleigh_R + 1.d0
- if( cmp(1:1) == "T" ) num_Love_T = num_Love_T + 1.d0
- endif
-
- enddo
-
- ! determines all picks on a trace component
- ! transverse
- iposition = INDEX( trim(synfile), comp_T, .false. )
- if( iposition > 3 .and. iposition < len_trim( synfile) ) then
- if( cmp(1:1) /= "T" ) stop 'error T component pick'
- picks_T = picks_T + npicks
- else
- ! radial
- iposition = INDEX( trim(synfile), comp_R, .false. )
- if( iposition > 3 .and. iposition < len_trim( synfile) ) then
- if( cmp(1:1) /= "R" ) stop 'error R component pick'
- picks_R = picks_R + npicks
- else
- ! vertical
- iposition = INDEX( trim(synfile), comp_Z, .false. )
- if( iposition > 3 .and. iposition < len_trim( synfile) ) then
- if( cmp(1:1) /= "Z" ) stop 'error Z component pick'
- picks_Z = picks_Z + npicks
- endif
- endif
- endif
-
- enddo
- close(21)
-
-
- ! check with total number of picks per component
- if( nint( num_P_SV_R + num_Rayleigh_R ) /= picks_R ) stop 'error R picks'
- if( nint( num_P_SV_V + num_Rayleigh_V ) /= picks_Z ) stop 'error Z picks'
- if( nint( num_SH_T + num_Love_T ) /= picks_T ) stop 'error T picks'
-
- if( DO_WEIGHTING ) then
- print*
- print*,'weighting measurements: '
- print*,' picks T:',picks_T
- print*,' picks R:',picks_R
- print*,' picks Z:',picks_Z
- print*
- print*,' picks P_SV_R: ',nint(num_P_SV_R)
- print*,' picks P_SV_V: ',nint(num_P_SV_V)
- print*,' picks SH_T : ',nint(num_SH_T)
- print*,' picks Rayleigh_R: ',nint(num_Rayleigh_R)
- print*,' picks Rayleigh_V: ',nint(num_Rayleigh_V)
- print*,' picks Love_T : ',nint(num_Love_T)
- print*
- endif
-
-
- ! sets up weights based on picks
- weight_T = 1.0d0
- weight_R = 1.0d0
- weight_Z = 1.0d0
-
- ! weighting tries to balance love waves (tranverse) versus rayleigh waves (radial + vertical)
- !if( picks_T > 0 ) then
- ! if( picks_R + picks_Z > 0 ) weight_T = dble(picks_R + picks_Z)/dble(picks_T)
- !endif
-
- ! use normalization as weights
- if( picks_T > 0 ) weight_T = 1.d0 / picks_T
- if( picks_R > 0 ) weight_R = 1.d0 / picks_R
- if( picks_Z > 0 ) weight_Z = 1.d0 / picks_Z
-
- ! use normalization
- if( num_P_SV_R > 0. ) num_P_SV_R = 1.d0 / num_P_SV_R
- if( num_P_SV_V > 0. ) num_P_SV_V = 1.d0 / num_P_SV_V
- if( num_SH_T > 0. ) num_SH_T = 1.d0 / num_SH_T
- if( num_Rayleigh_R > 0. ) num_Rayleigh_R = 1.d0 / num_Rayleigh_R
- if( num_Rayleigh_V > 0. ) num_Rayleigh_V = 1.d0 / num_Rayleigh_V
- if( num_Love_T > 0. ) num_Love_T = 1.d0 / num_Love_T
-
- if( DO_WEIGHTING ) then
- print*,' weight num_P_SV_R:',num_P_SV_R
- print*,' weight num_P_SV_V:',num_P_SV_V
- print*,' weight num_SH_T :',num_SH_T
- print*,' weight num_Rayleigh_R:',num_Rayleigh_R
- print*,' weight num_Rayleigh_V:',num_Rayleigh_V
- print*,' weight num_Love_T :',num_Love_T
- print*
- endif
-
-end subroutine setup_weighting
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/prepare_adj_src.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/prepare_adj_src.pl 2012-11-27 02:07:13 UTC (rev 21074)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/prepare_adj_src.pl 2012-11-27 17:05:45 UTC (rev 21075)
@@ -1,4 +1,4 @@
-#!/usr/bin/perl -w
+#!/usr/bin/perl
#
#--------------------------------------------------------------------
# prepare_adj_src.pl
@@ -22,7 +22,7 @@
use DELAZ5;
#Reading input arguments:
- at ARGV>0 || die("prepare_adj_src.pl -m CMT -z BH -s STATION -o OUTDIR -i iadj all_adj_files\n");
+ at 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");}
@@ -117,4 +117,4 @@
}
# write out STATIONS_ADJOINT file
-system("echo $nstation > STATIONS_ADJOINT; cat ${stafile_temp} >> STATIONS_ADJOINT");
+system("echo $nstation > STATIONS_ADJOINT; cat ${stafile_temp} >> STATIONS_ADJOINT; rm -f ${stafile_temp}");
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/rotate_adj_src.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/rotate_adj_src.f90 2012-11-27 02:07:13 UTC (rev 21074)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/rotate_adj_src.f90 2012-11-27 17:05:45 UTC (rev 21075)
@@ -22,7 +22,7 @@
if (trim(bazch)=='' .or. trim(zfile) =='' .or. trim(tfile)=='' .or. &
trim(rfile) == '' .or. trim(efile) =='' .or. trim(nfile) =='') then
- stop 'rotate_adj_src baz zfile tfile rfile efile nfile'
+ stop 'rotate_adj_src baz(radian!) zfile tfile rfile efile nfile'
endif
read(bazch,*) baz
@@ -31,30 +31,30 @@
inquire(file=trim(rfile),exist=r_exist)
inquire(file=trim(zfile),exist=z_exist)
+ ! initialization
+ rdata = 0; tdata = 0
+
+ ! 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'
-
- call drascii(zfile,tdata,npts,t0,dt)
- !call dread_ascfile_f(zfile,t0,dt,npts,tdata)
- tdata = 0.
- rdata = 0.
+ ! no need to read Z comp adjoint source (keep the same)
+ ! call drascii(zfile,tdata,npts,t0,dt)
endif
+ ! read in T file
if (t_exist) then
call drascii(tfile,tdata,nptst,t0t,dtt)
- !call dread_ascfile_f(tfile,t0t,dtt,nptst,tdata)
- if (.not. r_exist) rdata = 0.
endif
+ ! read in R file
if (r_exist) then
call drascii(rfile,rdata,nptsr,t0r,dtr)
- !call dread_ascfile_f(rfile,t0r,dtr,nptsr,rdata)
- if (.not. t_exist) tdata = 0.
endif
+
+ ! check consistency of t0,dt,npts
if (t_exist .and. r_exist) then
if (abs(t0t-t0r)>1.0e-2 .or. abs(dtt-dtr)>1.0e-2 .or. nptst /= nptsr) &
stop 'check t0 and npts'
endif
-
if (t_exist) then
t0 =t0t; dt = dtt; npts = nptst
else if (r_exist) then
@@ -62,20 +62,20 @@
endif
if (NDIM < npts) stop 'Increase NDIM'
+ ! rotate T,R to E,N based on baz (note in radian!)
costh = cos(baz)
sinth = sin(baz)
edata(1:npts) = -costh * tdata(1:npts) - sinth * rdata(1:npts)
ndata(1:npts) = sinth * tdata(1:npts) - costh * rdata(1:npts)
+ ! write E,N files
call dwascii(efile,edata,npts,t0,dt)
call dwascii(nfile,ndata,npts,t0,dt)
- !call dwrite_ascfile_f(efile,t0,dt,npts,edata)
- !call dwrite_ascfile_f(nfile,t0,dt,npts,ndata)
+ ! write Z file if did not exist
if (.not. z_exist) then
tdata = 0.
call dwascii(zfile,tdata,npts,t0,dt)
- !call dwrite_ascfile_f(zfile,t0,dt,npts,tdata)
endif
end program rotate_adj_src
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/write_par_file.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/write_par_file.pl 2012-11-27 02:07:13 UTC (rev 21074)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/write_par_file.pl 2012-11-27 17:05:45 UTC (rev 21075)
@@ -22,7 +22,11 @@
#
#==========================================================
-if (@ARGV < 8) {die("Usage: write_par_file.pl xxx\n");}
+if (@ARGV < 8) {die("Usage: write_par_file.pl tstart/dt/npts imeas chan Tshort/Tlong
+ RUN_BANDPASS/DISPLAY_DETAILS/OUTPUT_MEASUREMENT_FILES/COMPUTE_ADJONIT_SOURCE
+ tshift_min/tshift_max/dlnA_min/dlnA_max/cc_min error_type/dt_sigma_min/dlnA_sigma_min
+ itaper/wtr/npi/dt_fac/err_fac/dt_max_scale/ncycle_in_window
+ **** writes MEASUREMENT.PAR with given parameters *** \n\n");}
($tvec,$imeas,$chan,$Ts,$iparbools,$par1,$par2,$par3) = @ARGV;
# extract variables
More information about the CIG-COMMITS
mailing list