[cig-commits] r17886 - in seismo/3D/ADJOINT_TOMO/measure_adj: . PLOTS scripts_meas
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Wed Feb 16 17:13:56 PST 2011
Author: danielpeter
Date: 2011-02-16 17:13:56 -0800 (Wed, 16 Feb 2011)
New Revision: 17886
Modified:
seismo/3D/ADJOINT_TOMO/measure_adj/Makefile
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl
seismo/3D/ADJOINT_TOMO/measure_adj/README
seismo/3D/ADJOINT_TOMO/measure_adj/mt_constants.f90
seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub2.f90
seismo/3D/ADJOINT_TOMO/measure_adj/mt_variables.f90
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.pl
Log:
adds possibility to normalize misfits for different wave types and to calculate ray density adjoint sources
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/Makefile
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/Makefile 2011-02-16 21:49:17 UTC (rev 17885)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/Makefile 2011-02-17 01:13:56 UTC (rev 17886)
@@ -1,13 +1,20 @@
F90 = gfortran
-F90_FLAGS = -O2 -fbounds-check
+#64-bit
+F90_FLAGS = -O2 -fbounds-check #-Wall
+
+#32-bit
+#F90_FLAGS = -O2 -fbounds-check -m32 #-Wall
+
# new version with default sac libraries
SACLIBDIR = ${SACHOME}/lib
LIB = -lsacio -lsac
#LIB = -L/opt/seismo/lib -lDRWFiles -lf90recipes -lDSacio -lDSacLib -lSacTools -lm
+
+
# NOTE: order matters if modules depend on other modules
-MOD = mt_constants mt_variables ascii_rw mt_sub2 mt_sub
+MOD = mt_constants mt_variables ascii_rw mt_sub2 mt_sub
SRC_DIR = .
MOD_DIR = mod
@@ -26,10 +33,10 @@
$(F90) -o $(BIN_DIR)/$* $(F90_FLAGS) $(SRC_DIR)/$*.f90 -I$(MOD_DIR) -J$(MOD_DIR) $(OBJ) -L${SACLIBDIR} $(LIB)
$(F90_OBJ): $(OBJ_DIR)/%.o : $(SRC_DIR)/%.f90
- $(F90) -o $@ $(F90_FLAGS) -c $(SRC_DIR)/$*.f90
+ $(F90) -o $@ $(F90_FLAGS) -c $(SRC_DIR)/$*.f90
$(MOD_OBJ): $(OBJ_DIR)/%.o : $(SRC_DIR)/%.f90
- $(F90) -o $@ $(F90_FLAGS) -c $(SRC_DIR)/$*.f90 -I$(MOD_DIR) -J$(MOD_DIR)
+ $(F90) -o $@ $(F90_FLAGS) -c $(SRC_DIR)/$*.f90 -I$(MOD_DIR) -J$(MOD_DIR)
#rotate_adj_src: rotate_adj_src.f90
# $(F90) -o rotate_adj_src rotate_adj_src.f90 -L/opt/seismo/lib -lDRWFiles
@@ -37,6 +44,6 @@
.PHONY : clean
clean:
- \rm -f *.o *.mod *~ $(OBJ_DIR)/*.o $(MOD_DIR)/*.mod *.txt* STA.* OUTPUT_FILES/* *.sac rotate_adj_src mt_measure_adj
+ \rm -f *.o *.mod $(OBJ_DIR)/*.o $(MOD_DIR)/*.mod *.txt* STA.* OUTPUT_FILES/* *.sac rotate_adj_src mt_measure_adj
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 2011-02-16 21:49:17 UTC (rev 17885)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl 2011-02-17 01:13:56 UTC (rev 17886)
@@ -36,6 +36,8 @@
if($opt_i) {$smodel=$opt_i;}
if($opt_j) {$Ts=$opt_j;}
+print "\nplot_win_adj_all.pl:\n";
+
# loop over STATIONS -- first line is the number of stations
$ns = `sed -n 1p $opt_a`;
$nline = 1;
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/README 2011-02-16 21:49:17 UTC (rev 17885)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/README 2011-02-17 01:13:56 UTC (rev 17886)
@@ -19,7 +19,9 @@
Qinya Liu, 12/2/06
Carl Tape, 2/7/08
-This package measures the transfer function for given data and synthetics, and then assimilates the measurements to create the adjoint source and misfit function values necessary for a conjugate gradient inversion.
+This package measures the transfer function for given data and synthetics,
+and then assimilates the measurements to create the adjoint source and
+misfit function values necessary for a conjugate gradient inversion.
Input: MEASUREMENT.WINDOWS which has the format:
@@ -39,9 +41,12 @@
Output: *.txt *.sac and OUTPUT_FILES/*
-For examples of different dataset, try other INPUT_FILES/MEASUREMENT_WINDOWS* files by copying them to MEASUREMENT.WINDOWS
+For examples of different dataset, try other INPUT_FILES/MEASUREMENT_WINDOWS* files
+by copying them to MEASUREMENT.WINDOWS
-One thing tricky: change the measurement window length, while keeping the same number of tapers is going to change the frequency domain averaging width W, therefore, affect the smoothness (sometimes amplitude as well) of the adjoint source.
+One thing tricky: change the measurement window length, while keeping
+the same number of tapers is going to change the frequency domain averaging width W,
+therefore, affect the smoothness (sometimes amplitude as well) of the adjoint source.
General scripts are found in the main directory:
@@ -57,7 +62,9 @@
==============================
Carl Tape, 2/7/08
-Several scripts have been added to facilitate running multiple events and generating multiple plots for the purpose of a tomography inversion. These are by no means in the "final" user-friendly state.
+Several scripts have been added to facilitate running multiple events
+and generating multiple plots for the purpose of a tomography inversion.
+These are by no means in the "final" user-friendly state.
These scripts have been placed in the directory scripts_tomo:
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_constants.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_constants.f90 2011-02-16 21:49:17 UTC (rev 17885)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_constants.f90 2011-02-17 01:13:56 UTC (rev 17886)
@@ -11,7 +11,7 @@
double precision, parameter :: LARGE_VAL = 1.0d8
! FFT parameters
- integer, parameter :: LNPT = 13, NPT = 2**LNPT, NDIM = 40000
+ integer, parameter :: LNPT = 15, NPT = 2**LNPT, NDIM = 80000
double precision, parameter :: FORWARD_FFT = 1.0
double precision, parameter :: REVERSE_FFT = -1.0
@@ -26,4 +26,11 @@
integer, parameter :: IORD = 4
integer, parameter :: PASSES = 2
+ ! 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.
+
+ ! constructs adjoint sources for a "ray density" kernel, where all misfits are equal to one
+ logical:: DO_RAY_DENSITY_SOURCE = .false.
+
end module mt_constants
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2011-02-16 21:49:17 UTC (rev 17885)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2011-02-17 01:13:56 UTC (rev 17886)
@@ -1,20 +1,40 @@
program mt_measure_adj
- ! main program that calls the subroutines to make the MT measurements
+ ! main program that calls the subroutines to make the MT measurements
! and compute the corresponding adjoint sources
+ ! input parameter:
+ ! 1. imeas = 1, normalized waveform difference. Adjoint source is constructed from the data
+ ! only, with the form −d(t)/ || d(t) || 2
+ ! 2. imeas = 2, waveform difference, s(t) − d(t).
+ ! 3. imeas = 3, cross-correlation traveltime difference for a (banana-doughtnut) sensitivity ker-
+ ! nel. The measurement between data and synthetics is not used in constructing the adjoint
+ ! source.
+ ! 4. imeas = 4, amplitude difference for a (banana-doughtnut) sensitivity kernel. The measure-
+ ! ment between data and synthetics is not used in constructing the adjoint source.
+ ! 5. imeas = 5, cross-correlation traveltime difference for an event kernel. The measurement
+ ! between data and synthetics is used in constructing the adjoint source.
+ ! 6. imeas = 6, amplitude difference for an event kernel. The measurement between data and
+ ! 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
+ ! data and synthetics is used in constructing the adjoint source. See multitaper_notes.pdf.
+
+
use mt_variables
use mt_constants
use ascii_rw ! dwascii()
use mt_sub2 ! fft(), fftinv()
use mt_sub ! mt_measure(), mt_adj()
+ use mt_weighting
implicit none
character(len=150) :: datafile,synfile,file_prefix,file_prefix0,file_prefix2,measure_file_prefix,adj_file_prefix
- integer :: num_files, num_meas, i, j, k, ios, npt1, npt2, npts, nn, nerr
+ integer :: num_meas, j, ios, npt1, npt2, npts, nn
double precision, dimension(NDIM) :: data, syn, adj_syn_all, tr_adj_src, am_adj_src, recon_cc_all, syn_dtw_cc
- double precision :: t01, dt1, t02, dt2, t0, dt, tstart, tend, chi, tt, dtt, df
+ double precision :: t01, dt1, t02, dt2, t0, dt, tstart, tend, tt, dtt, df
double precision, dimension(NCHI) :: window_chi
double precision :: fend0, fstart0, fend, fstart
!double precision :: TSHORT, TLONG ! mtm_variables.f90
@@ -28,23 +48,33 @@
!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
- integer :: nlen, nlen0, i_left, i_pmax_dat, i_pmax_syn, i_right, i_right0, istart, &
- sta_len,net_len,chan_len, ipair, npairs, nwin, itmax
- logical :: use_trace, flag_short_window
+ integer :: nlen, i_left, i_pmax_dat, i_pmax_syn, i_right, i_right0, istart, &
+ ipair, npairs, nwin, itmax
+ 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 mt_sub.f90 and write_par_file.pl)
! most variables are global (see mt_variables.f90)
call read_par_file(fstart0,fend0,tt,dtt,nn,chan)
+ ! 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
+
! 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
read(11,*,iostat=ios) npairs
if (ios /= 0) then ; print *, 'Error reading number of pairs of data/syn' ; stop ; endif
+
print *, 'reading in the data and synthetics...'
! output files
@@ -60,6 +90,7 @@
adj_syn_all(:) = 0.0
recon_cc_all(:) = 0.0
+ ! 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
read(11,'(a)',iostat=ios) synfile
@@ -70,6 +101,17 @@
call drsac1(synfile,syn,npt2,t02,dt2)
!print *, npt1,t01,dt1,NDIM ! check: double precision
+ ! 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
@@ -78,11 +120,24 @@
! 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)
+ 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
+
+
! apply bandpass filter to data and synthetics, 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
+ ! Access to the kidate, xapiir, and getfil is not simple and not
! supported under the current state of the SAC code base.
t0 = t01
@@ -95,8 +150,12 @@
endif
! 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)
+ call get_sacfile_header(trim(datafile),yr,jda,ho,mi,sec,net,sta, &
+ chan_dat,dist,az,baz,slat,slon)
+ ! theoretical surface wave arrival time
+ T_surfacewaves = dist / surface_vel
+
! synthetics always have the form BH_ or LH_, but the data may not (HH_, LH_, BL_, etc).
cmp = chan_dat(3:3)
chan_syn = trim(chan)//trim(cmp)
@@ -122,18 +181,34 @@
!!$ stop
!!$ endif
+ ! reads number of measurement windows
read(11,*,iostat=ios) num_meas
if (ios /= 0) then ; print *, 'Error reading num_meas' ; stop ; endif
- do j = 1, num_meas
-
+ do j = 1, num_meas
+ ! reads in start and end time of measurement window
read(11,*,iostat=ios) tstart, tend
if (ios /= 0) then ; print *, 'Error reading tstart and tend' ; stop ; endif
+ ! 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
+ ! 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
@@ -164,7 +239,7 @@
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
+ ! 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)
@@ -183,7 +258,8 @@
! adjust measurements for MT adjoint source
if (is_mtm == 1) then
fstart = fstart0 ; fend = fend0
- call mt_measure_select(nlen,tshift,i_pmax_syn,dtau_w,dlnA_w,err_dt,err_dlnA,dt,i_left,i_right,fstart,fend,use_trace)
+ 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
@@ -207,17 +283,28 @@
! check that the CC measurements are within the specified input range
if (imeas >= 5) call cc_measure_select(tshift,dlnA,cc_max)
- ! write frequency limits to file
- if (OUTPUT_MEASUREMENT_FILES) then
- df = 1./(dt*NPT)
- open(71,file=trim(measure_file_prefix)//'.freq_limits')
- write(71,'(6f18.8)') fstart0, fend0, df, i_right0*df, fstart, fend
- close(71)
- endif
+ ! write frequency limits to file
+ if (OUTPUT_MEASUREMENT_FILES) then
+ df = 1./(dt*NPT)
+ open(71,file=trim(measure_file_prefix)//'.freq_limits')
+ write(71,'(6f18.8)') fstart0, fend0, df, i_right0*df, fstart, fend
+ close(71)
+ endif
! compute adjoint sources and misfit function values and also the CC-reconstructed records
if (use_trace) then
+ ! banana-doughnut kernel (needs only synthetic trace)
+ 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
+ endif
+
tr_chi = 0.0 ; am_chi = 0.0 ! must be initialized
call mt_adj(istart,data_dtw,syn_dtw,nlen,dt,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc,&
dtau_w,dlnA_w,err_dt,err_dlnA,sigma_dt,sigma_dlnA,i_left,i_right,&
@@ -236,6 +323,39 @@
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)
+ ! 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
+ ! weights by phase types
+ select case(ipick_type)
+ case( P_SV_V )
+ tr_adj_src(:) = tr_adj_src(:) * num_P_SV_V
+ case( P_SV_R )
+ tr_adj_src(:) = tr_adj_src(:) * num_P_SV_R
+ case( SH_T )
+ tr_adj_src(:) = tr_adj_src(:) * num_SH_T
+ case( Rayleigh_V )
+ tr_adj_src(:) = tr_adj_src(:) * num_Rayleigh_V
+ case( Rayleigh_R )
+ tr_adj_src(:) = tr_adj_src(:) * num_Rayleigh_R
+ case( Love_T )
+ tr_adj_src(:) = tr_adj_src(:) * num_Love_T
+ case default
+ stop 'error ipick_type unknown'
+ end select
+ endif
+
! combine adjoint sources from different measurement windows
if (COMPUTE_ADJOINT_SOURCE) then
if (mod(imeas,2)==1) then
@@ -263,25 +383,54 @@
if (COMPUTE_ADJOINT_SOURCE) then
- ! OPTIONAL: A conservative choice is to filter the adjoint source,
- ! since higher frequencies could enter from the tapering operations.
- ! Note: time_window in mt_adj.f90 tapers the windows.
- !call bandpass(adj_syn_all,npts,dt,fstart0,fend0)
+ ! 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.
- ! cut and interpolate to match time-stepping for SEM
- ! NOTE: This can leave a non-zero value to start the record,
- ! which is NOT GOOD for the SEM simulation.
- call interpolate_syn(adj_syn_all,t0,dt,npts,tt,dtt,nn)
+ ! note also:
+ ! measurements are done on filtered synthetics F(s) and filtered data F(d), such that DeltaT
+ ! is given for filtered data & synthetics.
+ ! 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,
+ ! which is unfiltered (only filter there is by source half-time), but we want to convolve
+ ! K = int F*(s_adj) F(s)
+ ! using the same (bandpass) filter F() as used for filtereing data & synthetics in the meausurements
+ ! We can write the kernel expression as K = int F*{F* (s_adj)} s
+ ! thus we should apply the filter F() twice on the adjoint source
+ !
+ ! why is this important? the filter, like bandpassing, is usually acausal, that is, it can
+ ! introduce a slight phase-shift to the data. but, phase-shifts is what we are interested in
+ ! and invert for. so, filtering might affect our inversions...
- ! Taper the start of the adjoint source, since cutting the record
- ! may have left a non-zero value to start the record,
- ! which is not good for the SEM simulation.
- itmax = int(TSHORT/dtt)
- call taper_start(adj_syn_all,nn,itmax)
+ ! we 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)
- ! we can have choices of outputing the adjoint source as ASCII or SAC format
- call dwascii(trim(adj_file_prefix)//'.adj',adj_syn_all,nn,tt,dtt)
- !call dwsac1(trim(adj_file_prefix)//'.adj.sac',adj_syn_all,nn,tt,dtt)
+
+ ! cut and interpolate to match time-stepping for SEM
+ ! NOTE: This can leave a non-zero value to start the record,
+ ! which is NOT GOOD for the SEM simulation.
+ call interpolate_syn(adj_syn_all,t0,dt,npts,tt,dtt,nn)
+
+ ! Taper the start of the adjoint source, since cutting the record
+ ! may have left a non-zero value to start the record,
+ ! which is not good for the SEM simulation.
+ itmax = int(TSHORT/dtt)
+ call taper_start(adj_syn_all,nn,itmax)
+
+ ! we can have choices of outputing the adjoint source as ASCII or SAC format
+ if( DO_RAY_DENSITY_SOURCE ) then
+ 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)
+ 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)
@@ -296,7 +445,8 @@
! 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
+
+ !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
@@ -306,4 +456,215 @@
close(13) ! write: window_chi
end program mt_measure_adj
-
+
+subroutine setup_weighting(chan_syn)
+ !
+ ! determines weights based on number of window picks on radial, transverse and vertical components
+ !
+ use mt_weighting
+
+ use mt_constants,only: NDIM
+ use mt_sub,only: get_sacfile_header,drsac1
+ use mt_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
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2011-02-16 21:49:17 UTC (rev 17885)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2011-02-17 01:13:56 UTC (rev 17886)
@@ -14,7 +14,7 @@
! 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
+ ! 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
@@ -28,16 +28,18 @@
! 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)
!
- ! original coding in Fortran77 by Ying Zhou
+ ! original coding in Fortran77 by Ying Zhou
! upgraded to Fortran90 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
integer, intent(in) :: npts
double precision, dimension(:), intent(in) :: dat_dt,syn_dt
double precision, intent(in) :: tstart,tend,t0,dt
@@ -52,9 +54,9 @@
double precision, dimension(NPT) :: syn_vtw, syn_dtw_mt, syn_dtw_mt_dt, &
syn_dtw_cc_dt, dat_dtw_cc, syn_dtw_h, dat_dtw_h
- double precision :: sfac1,fac,omega,f0,df,df_new,dw,&
- ampmax_unw,wtr_use_unw,ampmax,wtr_use,wtr_mtm,dtau_wa,dlnA_wa
- integer :: ishift,i,ictaper,j,k,fnum,i_amp_max_unw,i_amp_max,i_right_stop,idf_new,iom
+ 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
complex*16, dimension(NPT) :: syn_dtwo, dat_dtwo, syn_dtw_ho, dat_dtw_ho, &
top_mtm, bot_mtm, trans_mtm, wseis_recon
@@ -77,7 +79,7 @@
i_right = 0
! LQY -- is this too small ???
- wtr_mtm = 1.e-10
+ wtr_mtm = 1.e-10
filename = trim(file_prefix)
@@ -90,8 +92,8 @@
!!$ 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)
+!!$ !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)
endif
!--------------------------------------------------------------------------
@@ -154,18 +156,18 @@
! deconstruct data using (negative) cross-correlation measurments
call deconstruct_dat_cc(filename,dat_dtw,tstart,dt,nlen, &
- ishift,tshift,dlnA,dat_dtw_cc)
+ 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)
+ 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)
+!!$ !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
@@ -214,41 +216,41 @@
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)
+ call fft(LNPT,dat_dtwo,FORWARD_FFT,dt)
! index of the freq of the max power in the windowed data
ampmax_unw = 0.
i_pmax_dat = 1
do i = 1, fnum ! loop over frequencies
- if( abs(dat_dtwo(i)) > ampmax_unw) then
+ if( abs(dat_dtwo(i)) > ampmax_unw) then
ampmax_unw = abs(dat_dtwo(i))
i_pmax_dat = i
endif
enddo
- ! water level based untapered synthetics
- ! used to determine the i_right values (maximum frequency for measurement)
+ ! 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
- if( abs(syn_dtwo(i)) > ampmax_unw) then
+ if( abs(syn_dtwo(i)) > ampmax_unw) then
ampmax_unw = abs(syn_dtwo(i))
i_amp_max_unw = i
endif
enddo
- wtr_use_unw = cmplx(ampmax_unw * WTR, 0.)
+ 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 = 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 = 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
+ 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 = i
endif
enddo
@@ -256,7 +258,7 @@
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 *, '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)
@@ -277,7 +279,7 @@
! assign number of tapers
if (is_mtm == 1) then
ntaper = int(NPI * 2.0)
- else
+ else
ntaper = 1
endif
allocate(tas(NPT,ntaper))
@@ -303,7 +305,7 @@
bot_mtm(:) = cmplx(0.,0.)
trans_mtm(:) = cmplx(0.,0.)
- do ictaper = 1, ntaper
+ do ictaper = 1, ntaper
syn_dtw_ho(:) = cmplx(0.,0.) ! note: this has to be initialized inside the loop
dat_dtw_ho(:) = cmplx(0.,0.)
@@ -325,12 +327,12 @@
! in the tapered synthetics record
ampmax = 0.
do i = 1, fnum ! loop over frequencies
- if( abs(syn_dtw_ho(i)) > ampmax) then ! syn, single_tapered
+ if( abs(syn_dtw_ho(i)) > ampmax) then ! syn, single_tapered
ampmax = abs(syn_dtw_ho(i))
i_amp_max = i
endif
enddo
- wtr_use = cmplx(ampmax * WTR, 0.)
+ wtr_use = cmplx(ampmax * WTR, 0.)
!print *, ' wtr_use :', wtr_use
! calculate top and bottom of MT transfer function
@@ -340,7 +342,7 @@
! 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)
if(abs(syn_dtw_ho(i)) <= abs(wtr_use)) trans_w(i) = dat_dtw_ho(i) / (syn_dtw_ho(i)+wtr_use)
endif
enddo
@@ -349,12 +351,12 @@
! 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
+ 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)
+ 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)
@@ -372,15 +374,15 @@
! water level for multitaper measurements
ampmax = 0.
do i = 1, fnum
- if( abs(bot_mtm(i)) > ampmax) then
+ 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_mtm**2, 0.)
!wtr_use = cmplx(ampmax * WTR, 0.)
- ! calculate MT transfer function using water level
+ ! 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)
@@ -392,7 +394,7 @@
! apply transfer function to the syn
call reconstruct_syn(filename,syn_dtwo,wvec,dtau_mtm,dlnA_mtm, &
- i_right,tstart,dt,nlen,syn_dtw_mt,syn_dtw_mt_dt)
+ i_right,tstart,dt,nlen,syn_dtw_mt,syn_dtw_mt_dt)
! check quality
!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)
@@ -424,7 +426,7 @@
allocate(dtau_mul(NPT,ntaper))
allocate(dlnA_mul(NPT,ntaper))
- do iom = 1, ntaper
+ do iom = 1, ntaper
top_mtm(:) = cmplx(0.,0.)
bot_mtm(:) = cmplx(0.,0.)
@@ -453,10 +455,10 @@
enddo
enddo ! ictaper
- ! water level
+ ! water level
ampmax = 0.
do i = 1, fnum
- if( abs(bot_mtm(i)).gt.ampmax) then
+ if( abs(bot_mtm(i)).gt.ampmax) then
ampmax = abs(bot_mtm(i))
i_amp_max = i
endif
@@ -472,7 +474,7 @@
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 ! iom
!----------------------
@@ -527,7 +529,7 @@
err_phi(i) = sqrt( err_phi(i) / (ntaper * (ntaper-1) ) )
err_dt(i) = sqrt( err_dt(i) / (ntaper * (ntaper-1) ) )
- ! set the error bar for the first point corresponding to
+ ! set the error bar for the first point corresponding to
! static offset to be large, which makes no contribution to
! the adjoint source
if (i == 1) err_dt(i) = LARGE_VAL
@@ -561,8 +563,8 @@
dlnA_w = dlnA_mtm
! reset the control parameters
- OUTPUT_MEASUREMENT_FILES = output_logical
- DISPLAY_DETAILS = display_logical
+ OUTPUT_MEASUREMENT_FILES = output_logical
+ DISPLAY_DETAILS = display_logical
endif
@@ -585,14 +587,14 @@
! 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
+ ! 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)
!
- ! Output:
- ! tr_adj_src(:), tr_chi -- travel-time adjoint source and chi value
+ ! Output:
+ ! 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
- !
+ !
! original coding by Carl Tape, finalized by Qinya Liu
! ======================================================================================================
@@ -600,6 +602,7 @@
dtau_w,dlnA_w,err_dtau,err_dlnA,sigma_dt,sigma_dlnA,i_left,i_right, &
window_chi,tr_adj_src,tr_chi,am_adj_src,am_chi)
+ implicit none
integer, intent(in) :: istart, nlen, i_left, i_right
double precision, dimension(:), intent(in) :: dat_dtw, syn_dtw
double precision, intent(in) :: dt, tshift, dlnA, sigma_dt_cc, sigma_dlnA_cc, sigma_dt, sigma_dlnA
@@ -614,7 +617,7 @@
double precision, dimension(NPT) :: syn_vtw, syn_vtw_h, syn_dtw_h, ey1, ey2
double precision, dimension(NPT) :: ft_bar_t, fa_bar_t, fp, fq, wp_taper, wq_taper
complex*16, dimension(NPT) :: d_bot_mtm, v_bot_mtm
- integer :: i, i1, ictaper, ntaper, i2, ishift
+ integer :: i, i1, ictaper, ntaper
double precision :: df,Nnorm,Mnorm,fac,ffac,w_taper(NPT), time_window(NPT)
double precision, dimension(:,:), allocatable :: tas
complex*16, dimension(:,:),allocatable :: syn_dtw_ho_all, syn_vtw_ho_all
@@ -671,15 +674,19 @@
do i = i_left, i_right ! CHT: 1 --> i_left
- if (ERROR_TYPE == 0) then ! no error estimate
+ if (ERROR_TYPE == 0 .or. DO_RAY_DENSITY_SOURCE ) then
+ ! no error estimate
+ ! only adds normalization factor
wp_taper(i) = w_taper(i) / ffac
- wq_taper(i) = w_taper(i) / ffac
+ wq_taper(i) = w_taper(i) / ffac
- elseif (ERROR_TYPE == 1) then ! MT error estimate is assigned the CC error estimate
+ elseif (ERROR_TYPE == 1) then
+ ! MT error estimate is assigned the CC error estimate
wp_taper(i) = w_taper(i) / ffac / (sigma_dt ** 2)
wq_taper(i) = w_taper(i) / ffac / (sigma_dlnA ** 2)
- elseif (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)
@@ -742,7 +749,7 @@
if ( (is_mtm == 1).and.COMPUTE_ADJOINT_SOURCE ) then
! allocate MT variables
- ntaper = int(NPI * 2.0)
+ ntaper = int(NPI * 2.0)
allocate(tas(NPT,ntaper))
allocate(syn_dtw_ho_all(NPT,ntaper))
allocate(syn_vtw_ho_all(NPT,ntaper))
@@ -776,35 +783,44 @@
call fft(LNPT,syn_dtw_ho_all(:,ictaper),FORWARD_FFT,DT)
call fft(LNPT,syn_vtw_ho_all(:,ictaper),FORWARD_FFT,DT)
- d_bot_mtm(:) = d_bot_mtm(:) + syn_dtw_ho_all(:,ictaper) * conjg(syn_dtw_ho_all(:,ictaper))
+ d_bot_mtm(:) = d_bot_mtm(:) + syn_dtw_ho_all(:,ictaper) * conjg(syn_dtw_ho_all(:,ictaper))
v_bot_mtm(:) = v_bot_mtm(:) + syn_vtw_ho_all(:,ictaper) * conjg(syn_vtw_ho_all(:,ictaper))
enddo ! ictaper
! compute p_j, q_j, P_j, Q_j and adjoint source fp, fq
- fp = 0.
+ fp = 0.
fq = 0.
do ictaper = 1,ntaper
- ! compute p_j(w) and q_j(w)
+ ! compute p_j(w) and q_j(w)
pwc_adj(:) = cmplx(0.,0.)
qwc_adj(:) = cmplx(0.,0.)
do i = 1, i_right
- pwc_adj(i) = syn_vtw_ho_all(i,ictaper) / v_bot_mtm(i)
+ 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
- pwc_adj(:) = pwc_adj(:) * cmplx(dtau_w(:),0.) * cmplx(wp_taper(:),0.)
- qwc_adj(:) = qwc_adj(:) * cmplx(dlnA_w(:),0.) * cmplx(wq_taper(:),0.)
+ ! 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.)
+ qwc_adj(:) = qwc_adj(:) * cmplx(1.0,0.) * cmplx(wq_taper(:),0.)
+ else
+ ! adds misfit measurement dtau, dlnA
+ pwc_adj(:) = pwc_adj(:) * cmplx(dtau_w(:),0.) * cmplx(wp_taper(:),0.)
+ qwc_adj(:) = qwc_adj(:) * cmplx(dlnA_w(:),0.) * cmplx(wq_taper(:),0.)
+ endif
! IFFT into the time domain
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
fp(:) = fp(:) + tas(:,ictaper) * dtau_pj_t(:)
fq(:) = fq(:) + tas(:,ictaper) * dlnA_qj_t(:)
@@ -821,7 +837,7 @@
! integrated waveform difference squared
waveform_temp1 = 0. ; waveform_temp2 = 0. ; waveform_temp3 = 0.
- do i = 1,nlen
+ do i = 1,nlen
waveform_temp1 = waveform_temp1 + ( dat_dtw(i) * time_window(i) )**2
waveform_temp2 = waveform_temp2 + ( syn_dtw(i) * time_window(i) )**2
waveform_temp3 = waveform_temp3 + (( dat_dtw(i) - syn_dtw(i) ) * time_window(i) )**2
@@ -831,34 +847,47 @@
waveform_s2 = waveform_temp2
waveform_chi = waveform_temp3
+ ! compute traveltime and amplitude adjoint sources
if (COMPUTE_ADJOINT_SOURCE) then
- do i = 1,nlen
- i1 = istart + i
+ 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)
- ! consider normalizing this by waveform_d2
+ ! waveform
+ if(imeas==1 .or. imeas==2) then
+ tr_adj_src(i1) = -dat_dtw(i)/waveform_d2 * time_window(i)
+ am_adj_src(i1) = ( syn_dtw(i) - dat_dtw(i) ) * time_window(i)
+ ! consider normalizing this by waveform_d2
- ! 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)
+ ! use pure data waveform in time window
+ if( NO_WAVEFORM_DIFFERENCE ) then
+ tr_adj_src(i1) = dat_dtw(i) * time_window(i) ! waveform misfit
+ endif
- ! 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)
+ ! banana-doughnut kernel adjoint source (no measurement)
+ elseif(imeas==3 .or. imeas==4) then
+ tr_adj_src(i1) = ft_bar_t(i) * time_window(i)
+ am_adj_src(i1) = fa_bar_t(i) * time_window(i)
- ! multitaper
- elseif(imeas==7 .or. imeas==8) then
- tr_adj_src(i1) = fp(i) * time_window(i)
- am_adj_src(i1) = fq(i) * time_window(i)
+ ! cross-correlation
+ elseif(imeas==5 .or. imeas==6) then
+ tr_adj_src(i1) = -(tshift / sigma_dt_cc**2 ) * ft_bar_t(i) * time_window(i)
+ am_adj_src(i1) = -(dlnA / sigma_dlnA_cc**2 ) * fa_bar_t(i) * time_window(i)
+
+ ! 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)
endif
- enddo
+ ! multitaper
+ elseif(imeas==7 .or. imeas==8) then
+ tr_adj_src(i1) = fp(i) * time_window(i)
+ am_adj_src(i1) = fq(i) * time_window(i)
+ endif
+ enddo
+
endif
! -------------------------------------
@@ -873,6 +902,7 @@
! 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) )
@@ -885,12 +915,28 @@
window_chi(7) = tshift
window_chi(8) = dlnA
+ ! 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
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 * waveform_d2
window_chi(14) = 0.5 * waveform_s2
@@ -921,12 +967,13 @@
!==============================================================================
!==============================================================================
-
+
!----------------------------------------------------------------------
-
+
subroutine bandpass(x,n,delta_t,f1,f2)
! modified from FLEXWIN subroutines on 26-July-2009
+ implicit none
integer, intent(in) :: n
double precision, intent(inout), dimension(*) :: x
double precision, intent(in) :: delta_t,f1,f2
@@ -943,6 +990,8 @@
! new version, uses subroutines in libsac.a
! does band-pass filter
+ ! BU - butterworth
+ ! BP - bandpass
call xapiir(x_sngl,n,'BU',TRBDNDW,APARM,IORD,'BP',f1,f2,delta_t,PASSES)
x(1:n) = dble(x_sngl(1:n))
@@ -956,16 +1005,17 @@
subroutine drsac1(datafile,data,npt1,b1,dt1)
! read sac file and convert to double precision
+ implicit none
character(len=*),intent(in) :: datafile
real, dimension(NDIM) :: dat_sngl
double precision, dimension(NDIM), intent(out) :: data
integer :: npt1, nerr
real :: b1_sngl,dt1_sngl
double precision :: b1,dt1
-
+
! read file as single precision
call rsac1(datafile,dat_sngl,npt1,b1_sngl,dt1_sngl,NDIM,nerr)
- if (nerr > 0) then
+ if (nerr > 0) then
print *, 'Error reading sac file', trim(datafile)
stop
endif
@@ -984,6 +1034,7 @@
! --> includes an option to add minmax values to sac file,
! which are used in the plotting scripts
+ implicit none
character(len=*),intent(in) :: datafile
integer, intent(in) :: npt1
double precision, dimension(npt1), intent(in) :: data
@@ -993,7 +1044,7 @@
real, dimension(npt1) :: dat_sngl,ti_sngl
real :: b1_sngl,dt1_sngl,xmin_sngl,xmax_sngl
integer :: nerr,i
-
+
! convert to single precision
b1_sngl = real(b1)
dt1_sngl = real(dt1)
@@ -1006,6 +1057,8 @@
ti_sngl(i) = b1_sngl + (i-1)*dt1_sngl
enddo
+ !call newhdr() ! create a new header
+
! set minmax values in sac file
xmin_sngl = minval(dat_sngl)
xmax_sngl = maxval(dat_sngl)
@@ -1013,6 +1066,10 @@
call setfhv('depmax',xmax_sngl,nerr)
call setnhv('npts',npt1,nerr) ! sets number of points
+ !call setfhv('b',ti_sngl(1),nerr) ! sets begin
+ !call setfhv('e',ti_sngl(npt1),nerr) ! sets end
+ !call setlhv('leven',.false.,nerr) ! sets un-even sampling
+ !call setihv('iftype','itime',nerr) ! sets file type: time file
! write file with headers
call wsac0(datafile,ti_sngl,dat_sngl,nerr)
@@ -1031,18 +1088,17 @@
subroutine cc_measure_select(tshift,dlnA,cc_max)
- ! CHT: If the CC timeshift is for some reason larger than the allowable max,
+ ! CHT: If the CC timeshift is for some reason larger than the allowable max,
! then effectively eliminate the window by zeroing the
! cross-correlation traveltime and amplitude measurements.
! See subroutine compute_cc in mt_sub.f90.
+ implicit none
double precision, intent(inout) :: tshift, dlnA, cc_max
if( (cc_max < CC_MIN) .or. (tshift < TSHIFT_MIN) .or. (tshift > TSHIFT_MAX) &
.or. (dlnA < DLNA_MIN) .or. (dlnA > DLNA_MAX) ) then
! zero the CC measurments
- tshift = 0.0
- dlnA = 0.0
if (DISPLAY_DETAILS) then
print *, 'Fail if ANY of these is true :'
print *, ' cc_max : ', cc_max, CC_MIN, cc_max < CC_MIN
@@ -1051,39 +1107,44 @@
print *, ' dlnA : ', dlnA, DLNA_MIN, dlnA < DLNA_MIN
print *, ' dlnA : ', dlnA, DLNA_MAX, dlnA > DLNA_MAX
endif
+
+ ! zero the CC measurments
+ tshift = 0.0
+ dlnA = 0.0
endif
end subroutine cc_measure_select
!-----------------------------------------------------------------------------
- subroutine mt_measure_select(nlen,tshift,i_pmax_syn,dtau_w,dlnA_w,err_dt,err_dlnA, &
- dt,i_left,i_right,fstart,fend,use_trace)
+ subroutine mt_measure_select(nlen,tshift,i_pmax_syn,dtau_w,err_dt, &
+ dt,i_left,i_right,fstart,fend,use_trace)
! an important subroutine to determine whether an MT measurement should be rejected,
! in which case a CC measurement is used -- several choices are made here
+ implicit none
integer, intent(in) :: nlen, i_pmax_syn
double precision, intent(in) :: tshift, dt
- double precision, dimension(:), intent(inout) :: dtau_w, dlnA_w, err_dt, err_dlnA
+ double precision, dimension(:), intent(inout) :: dtau_w, err_dt
double precision, intent(inout) :: fstart, fend
integer,intent(inout) :: i_left, i_right
logical,intent(out) :: use_trace
- double precision :: df, fvec(NPT), dt_a, f_pmax, T_pmax, Wlen
+ double precision :: df, fvec(NPT), f_pmax, T_pmax, Wlen
integer :: i_right_old, i_left_old
- integer :: j
- logical :: stop_freq
+ integer :: j,ntaper
+ !logical :: stop_freq
use_trace = .true.
df = 1./(dt*NPT)
f_pmax = df * i_pmax_syn
- T_pmax = 1./ f_pmax
+ T_pmax = 1./ f_pmax
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 *, ' 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
@@ -1098,10 +1159,16 @@
! We subjectively state that we want at least 10 frequency points for the multitaper measurement.
fstart = max(fstart, NCYCLE_IN_WINDOW/Wlen)
fend = min(fend, 1./(2.0*dt))
- if (fstart >= fend - 10.0*df ) then
- print *, 'rejecting trace for frequency range (NCYCLE_IN_WINDOW):'
- print *, ' fstart >= fend - 10*df : '
- print *, ' fstart, fend, df, fend - 10*df : ', fstart, fend - 10.0*df
+
+ ! number of tapers (slepian tapers, type = 1)
+ ntaper = int(NPI * 2.0)
+ 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
use_trace = .false.
!stop 'testing rejection criteria'
endif
@@ -1177,7 +1244,7 @@
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 *, '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'
@@ -1193,6 +1260,7 @@
subroutine interpolate_dat_and_syn(data, syn, tstart, tend, t0, dt, NPT, dat_win, syn_win, nlen, istart)
+ 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
@@ -1205,13 +1273,21 @@
nlen = floor((tend-tstart)/dt) + 1
istart = floor((tstart-t0)/dt)
+ ! limits array bounds
+ if( nlen > NPT ) nlen = NPT
+
do i = 1, nlen
time = tstart + (i-1) * dt
ii = floor((time-t0)/dt) + 1
+
+ ! checks out-of-bounds
+ 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
+
enddo
end subroutine interpolate_dat_and_syn
@@ -1223,13 +1299,14 @@
! time shift MEASUREMENT between data (data) and synthetics (syn)
! CHT: modified the subroutine to resemble the one used in FLEXWIN
+ implicit none
double precision, dimension(*), intent(in) :: syn, data
integer, intent(in) :: nlen
double precision, intent(in) :: dt
double precision, intent(out) :: tshift, dlnA, cc_max
integer, intent(out) :: ishift
- double precision :: cr_shift, cc, norm_s, norm
+ double precision :: cc, norm_s, norm ! cr_shift
integer i1, i2, i, j, i_left, i_right, id_left, id_right
!!$ ! these choices will slide the entire windowed record past the other
@@ -1238,12 +1315,12 @@
!!$ i_right = floor( cr_shift / dt )
!!$
!!$ ! cross-correlation
-!!$ ishift = 0
+!!$ ishift = 0
!!$ do i = i_left, i_right, 1
!!$
!!$ cc = 0.
!!$ do j = 1, nlen
-!!$ if((j+i) > 1 .and. (j+i) < nlen) cc = cc + syn(j) * data(j+i)
+!!$ if((j+i) > 1 .and. (j+i) < nlen) cc = cc + syn(j) * data(j+i)
!!$ enddo
!!$
!!$ !if(cc > cc_max) then
@@ -1254,7 +1331,7 @@
!!$ endif
!!$
!!$ enddo
-!!$ tshift = ishift*dt
+!!$ tshift = ishift*dt
! initialise shift and cross correlation to zero
ishift = 0
@@ -1309,7 +1386,7 @@
! The previously used expression for dlnA of Dahlen and Baig (2002),
! is a first-order perturbation of ln(A1/A2) = (A1-A2)/A2 .
! The new expression is better suited to getting Gaussian-distributed
- ! values between -1 and 1, with dlnA = 0 indicating perfect fit, as before.
+ ! values between -1 and 1, with dlnA = 0 indicating perfect fit, as before.
dlnA = 0.5 * log( sum(data(i1:i2) * data(i1:i2)) / sum(syn(i1:i2) * syn(i1:i2)) )
end subroutine compute_cc
@@ -1323,6 +1400,7 @@
! and the reconstructed synthetics.
! NOTE: We implement the exact equations that are in the Latex notes.
+ implicit none
double precision, dimension(*), intent(in) :: data_dtw, syn_dtw_cc, syn_dtw_cc_dt
integer, intent(in) :: nlen
double precision, intent(in) :: dt
@@ -1342,7 +1420,7 @@
! estimated uncertainty in cross-correlation travltime 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 )
+ sigma_dt_bot = sum( syn_vtw_cc(1:nlen)**2 )
sigma_dlnA_bot = sum( (syn_dtw_cc_dt(1:nlen))**2 )
sigma_dt = sqrt( sigma_dt_top / sigma_dt_bot )
sigma_dlnA = sqrt( sigma_dlnA_top / sigma_dlnA_bot )
@@ -1377,8 +1455,9 @@
! ---------------------------------------------------------------------------
- subroutine write_average_meas(filename, imeas, dtau_meas, dlnA_meas, dtau_sigma, dlnA_sigma)
+ subroutine write_average_meas(filename, imeas, 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
@@ -1413,11 +1492,12 @@
! ---------------------------------------------------------------------------
subroutine write_trans(filename, trans, wvec, fnum, i_right, idf_new, df, tshift, dlnA, &
- phi_wt, abs_wt, dtau_wt, dlnA_wt, dtau_wa, dlnA_wa)
+ phi_wt, abs_wt, dtau_wt, dlnA_wt, dtau_wa, dlnA_wa)
! The transfer function maps the synthetics to the CC-deconstructed data;
! the CC measurements then need to be applied to match the original data.
+ implicit none
character(len=*), intent(in) :: filename
complex*16, intent(in) :: trans(:)
double precision, intent(in) :: wvec(:), df, tshift, dlnA
@@ -1427,7 +1507,7 @@
integer :: i, j
double precision, dimension(NPT) :: fr
- double precision :: f0, smth, smth1, smth2
+ double precision :: smth, smth1, smth2 ! f0
abs_wt(:) = 0.
phi_wt(:) = 0.
@@ -1441,7 +1521,7 @@
open(50,file=trim(filename)//'.dt')
endif
- ! loop to calculate phase and amplitude
+ ! loop to calculate phase and amplitude
do i = 1, i_right
phi_wt(i) = atan2( aimag(trans(i)) , real(trans(i)) )
abs_wt(i) = abs(trans(i))
@@ -1462,7 +1542,7 @@
smth = phi_wt(i+1) + phi_wt(i-1) - 2.0 * phi_wt(i)
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(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))
do j = i+1, i_right
phi_wt(j) = phi_wt(j) + TWOPI
@@ -1518,13 +1598,13 @@
! --------------------------------------------------------------------
subroutine deconstruct_dat_cc(filename,dat_dtw,tstart,dt,nlen,&
- ishift,tshift,dlnA,dat_dtw_cc)
+ 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.
-
+ implicit none
character(len=*), intent(in) :: filename
double precision, dimension(NPT), intent(in) :: dat_dtw
integer, intent(in) :: ishift, nlen
@@ -1546,18 +1626,18 @@
!if (DISPLAY_DETAILS) then
! call dwrite_sacfile_f(datafile,'windowed_shifted_data.sac',tstart,nlen,dat_dtw0)
- !endif
+ !endif
end subroutine deconstruct_dat_cc
! --------------------------------------------------------------------
- subroutine reconstruct_syn_cc(syn_dtw,tstart,dt,nlen,ishift,tshift,dlnA,syn_dtw_cc,syn_dtw_cc_dt)
+ subroutine reconstruct_syn_cc(syn_dtw,tstart,dt,nlen,ishift,tshift,dlnA,syn_dtw_cc,syn_dtw_cc_dt)
! 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
-
+ implicit none
double precision, dimension(NPT), intent(in) :: syn_dtw
integer, intent(in) :: ishift, nlen
double precision, intent(in) :: tshift, dlnA, tstart, dt
@@ -1588,7 +1668,7 @@
! reconstruct the synthetics using multitaper measurements:
! (1) apply dtau(w) and dlnA(w) to get syn_dtw_mt0
-
+ implicit none
character(len=*), intent(in) :: filename
complex*16, dimension(:), intent(in) :: syn_dtwo
double precision, dimension(:), intent(in) :: wvec, dtau_wt, dlnA_wt
@@ -1604,7 +1684,7 @@
syn_dtw_mt(:) = 0.
wseis_recon(:) = cmplx(0.,0.)
do i = 1,i_right
- omega = wvec(i)
+ omega = wvec(i)
wseis_recon(i) = syn_dtwo(i) * exp(dlnA_wt(i)) * exp(-CCI*omega*dtau_wt(i))
!wseis_recon(i) = syn_dtwo(i) * (1.+ dlnA_wt(i)) * exp(-CCI*omega*dtau_wt(i))
!wseis_recon(i) = syn_dtwo(i) * trans_mtm(i) * exp(-CCI*omega*tshift)
@@ -1615,7 +1695,7 @@
syn_dtw_mt_dt(:) = 0.
wseis_recon(:) = cmplx(0.,0.)
do i = 1,i_right
- omega = wvec(i)
+ omega = wvec(i)
wseis_recon(i) = syn_dtwo(i) * exp(-CCI*omega*dtau_wt(i))
enddo
call fftinv(LNPT,wseis_recon,REVERSE_FFT,dt,syn_dtw_mt_dt)
@@ -1625,8 +1705,8 @@
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)
+!!$ !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
@@ -1674,24 +1754,43 @@
subroutine interpolate_syn(syn,t1,dt1,npt1,t2,dt2,npt2)
+ implicit none
double precision, dimension(:),intent(inout) :: syn
integer,intent(in) :: npt1,npt2
double precision,intent(in) :: t1,dt1,t2,dt2
double precision :: syn1(NDIM), time, tt
integer i, ii
-
+
+ ! initializes trace holding interpolated values
syn1(1:npt2) = 0.
+ ! loops over number of time steps in complete trace
do i = 1, npt2
+
+ ! sets time (in s) at this time step:
+ ! t2 : start time of trace
+ ! dt2: delta_t of a single time step
time = t2 + (i-1) * dt2
+
+ ! checks if time is within measurement window
+ ! t1: start time of measurement window
+ ! npt1: number of time steps in measurement window
+ ! dt1: delta_t of a single time step in measurement window
if (time > t1 .and. time < t1 + (npt1-1)*dt1) then
+
+ ! sets index of time steps within this window: is 1 at the beginning of window
ii = floor((time-t1)/dt1) + 1
+
+ ! time increment within this single time step to match the exact value of time
tt = time - ((ii-1)*dt1 + t1)
+
+ ! interpolates value of trace for the exact time
syn1(i) = (syn(ii+1)-syn(ii)) * tt/dt1 + syn(ii)
endif
enddo
+ ! saves interpolated values to output trace
syn(1:npt2) = syn1(1:npt2)
end subroutine interpolate_syn
@@ -1700,11 +1799,12 @@
subroutine taper_start(syn,npt,itmax)
+ implicit none
double precision, dimension(:),intent(inout) :: syn
integer,intent(in) :: npt, itmax
double precision :: Wt
- integer :: imax, i
-
+ integer :: i !,imax
+
!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
@@ -1721,8 +1821,10 @@
!-------------------------------------------------------------------
+
subroutine read_par_file(fstart0,fend0,tt,dtt,nn,chan)
+ implicit none
double precision, intent(out) :: fstart0,fend0,tt,dtt
integer, intent(out) :: nn
character(len=10), intent(out) :: chan
@@ -1739,7 +1841,7 @@
read(10,*) TLONG, TSHORT
read(10,*) RUN_BANDPASS
read(10,*) DISPLAY_DETAILS
- read(10,*) OUTPUT_MEASUREMENT_FILES
+ read(10,*) OUTPUT_MEASUREMENT_FILES
read(10,*) COMPUTE_ADJOINT_SOURCE
read(10,*) TSHIFT_MIN, TSHIFT_MAX
read(10,*) DLNA_MIN, DLNA_MAX
@@ -1754,34 +1856,90 @@
read(10,*) DT_MAX_SCALE
read(10,*) NCYCLE_IN_WINDOW
close(10)
-
+
imeas = imeas0
! check the read-in values
print *, 'INPUTS FROM MEASUREMENT.PAR :'
print *, ' tt, dtt, nn : ',tt,dtt,nn
- print *, ' chan :',chan
+ print *, ' chan : ',chan
print *, ' imeas : ',imeas
print *, ' TLONG, TSHORT : ',TLONG, TSHORT
fstart0 = 1./TLONG ; fend0 = 1./TSHORT
- print *, ' fstart, fend :', fstart0, fend0
- print *, ' RUN_BANDPASS :',RUN_BANDPASS
- print *, ' DISPLAY_DETAILS :',DISPLAY_DETAILS
- print *, ' OUTPUT_MEASUREMENT_FILES :',OUTPUT_MEASUREMENT_FILES
- print *, ' TSHIFT_MIN, TSHIFT_MAX :',TSHIFT_MIN, TSHIFT_MAX
- print *, ' DLNA_MIN, DLNA_MAX :',DLNA_MIN, DLNA_MAX
- print *, ' CC_MIN :'
- print *, ' ERROR_TYPE :',ERROR_TYPE
- print *, ' DT_SIGMA_MIN :',DT_SIGMA_MIN
- print *, ' DLNA_SIGMA_MIN :',DLNA_SIGMA_MIN
+ print *, ' fstart, fend : ', fstart0, fend0
+ print *, ' RUN_BANDPASS : ',RUN_BANDPASS
+ print *, ' DISPLAY_DETAILS : ',DISPLAY_DETAILS
+ print *, ' OUTPUT_MEASUREMENT_FILES : ',OUTPUT_MEASUREMENT_FILES
+ print *, ' TSHIFT_MIN, TSHIFT_MAX : ',TSHIFT_MIN, TSHIFT_MAX
+ print *, ' DLNA_MIN, DLNA_MAX : ',DLNA_MIN, DLNA_MAX
+ print *, ' CC_MIN : ',CC_MIN
+ print *, ' ERROR_TYPE : ',ERROR_TYPE
+ print *, ' DT_SIGMA_MIN : ',DT_SIGMA_MIN
+ print *, ' DLNA_SIGMA_MIN : ',DLNA_SIGMA_MIN
print *, ' ITAPER : ',ITAPER
print *, ' WTR, NPI : ',WTR,NPI
- print *, ' DT_FAC :',DT_FAC
- print *, ' ERR_FAC :',ERR_FAC
- print *, ' DT_MAX_SCALE :',DT_MAX_SCALE
- print *, ' NCYCLE_IN_WINDOW :',NCYCLE_IN_WINDOW
+ print *, ' DT_FAC : ',DT_FAC
+ print *, ' ERR_FAC : ',ERR_FAC
+ print *, ' DT_MAX_SCALE : ',DT_MAX_SCALE
+ print *, ' NCYCLE_IN_WINDOW : ',NCYCLE_IN_WINDOW
!stop 'checking PAR file input'
-
+
+ ! old format way..
+ ! open(10,file='MEASUREMENT.PAR',status='old',iostat=ios)
+ ! read(10,'(a)') out_dir
+ ! read(10,*) is_mtm0
+ ! read(10,*) wtr,npi
+ ! read(10,*) iker0
+ ! read(10,*) RUN_BANDPASS
+ ! read(10,*) TLONG, TSHORT
+ ! read(10,*) tt,dtt,nn
+ ! read(10,*) DISPLAY_DETAILS
+ ! read(10,*) OUTPUT_MEASUREMENT_FILES
+ ! read(10,*) INCLUDE_ERROR
+ ! read(10,*) DT_FAC
+ ! read(10,*) ERR_FAC
+ ! read(10,*) DT_MAX_SCALE
+ ! read(10,*) NCYCLE_IN_WINDOW
+ ! read(10,*) BEFORE_QUALITY, AFTER_QUALITY
+ ! read(10,*) BEFORE_TSHIFT, AFTER_TSHIFT
+ ! read(10,*) DT_SIGMA_MIN, DLNA_SIGMA_MIN
+ ! close(10)
+ !
+ ! out_dir = adjustl(out_dir)
+ ! iker = iker0
+ ! is_mtm = is_mtm0
+ !
+ ! ! check the read-in values
+ ! print *, 'INPUTS FROM MEASUREMENT.PAR :'
+ ! print *, ' is_mtm : ',is_mtm
+ ! print *, ' wtr, npi : ',wtr,npi
+ ! print *, ' iker : ',iker
+ ! print *, ' RUN_BANDPASS :',RUN_BANDPASS
+ ! print *, ' TLONG, TSHORT : ',TLONG, TSHORT
+ ! fstart0 = 1./TLONG ; fend0 = 1./TSHORT
+ ! print *, ' fstart, fend :', fstart0, fend0
+ ! print *, ' tt, dtt, nn : ',tt,dtt,nn
+ ! print *, ' out_dir : ',trim(out_dir)
+ ! print *, ' DISPLAY_DETAILS :',DISPLAY_DETAILS
+ ! print *, ' OUTPUT_MEASUREMENT_FILES :',OUTPUT_MEASUREMENT_FILES
+ ! print *, ' INCLUDE_ERROR :',INCLUDE_ERROR
+ ! print *, ' DT_FAC :',DT_FAC
+ ! print *, ' ERR_FAC :',ERR_FAC
+ ! print *, ' DT_MAX_SCALE :',DT_MAX_SCALE
+ ! print *, ' NCYCLE_IN_WINDOW :',NCYCLE_IN_WINDOW
+ ! print *, ' BEFORE_QUALITY, AFTER_QUALITY :',BEFORE_QUALITY, AFTER_QUALITY
+ ! print *, ' BEFORE_TSHIFT, AFTER_TSHIFT :',BEFORE_TSHIFT, AFTER_TSHIFT
+ ! print *, ' DT_SIGMA_MIN, DLNA_SIGMA_MIN :',DT_SIGMA_MIN, DLNA_SIGMA_MIN
+ ! !stop 'checking PAR file input'
+ ! apply filter (this should EXACTLY match the filter used in the windowing code)
+ !trbdndw = 0.3
+ !a = 30.
+ !iord = 4
+ !passes = 2
+
+ ! ray density
+ 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'
@@ -1817,16 +1975,21 @@
!-------------------------------------------------------------------
- subroutine get_sacfile_header(data_file,yr,jda,ho,mi,sec,ntw,sta,comp,dist,az,baz,slat,slon)
+ subroutine get_sacfile_header(data_file,yr,jda,ho,mi,sec,ntw,sta, &
+ comp,dist,az,baz,slat,slon)
+ implicit none
character(len=*),intent(in) :: data_file
integer,intent(out):: yr,jda,ho,mi
- real*8,intent(out):: sec,dist,az,baz,slat,slon
+ 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
- integer :: klen,nerr
+ integer :: nsec,msec !,i,klen
+ integer :: nerr
! integer header variables
call getnhv('nzyear',yr,nerr)
@@ -1840,30 +2003,73 @@
! string headers
call getkhv('knetwk',ntw,nerr)
+ if(nerr .ne. 0) then
+ write(*,*)'Error reading variable: knetwk'
+ call exit(-1)
+ endif
+
call getkhv('kstnm',sta,nerr)
+ if(nerr .ne. 0) then
+ write(*,*)'Error reading variable: kstnm'
+ call exit(-1)
+ endif
+
call getkhv('kcmpnm',comp,nerr)
+ if(nerr .ne. 0) then
+ write(*,*)'Error reading variable: kcmpnm'
+ call exit(-1)
+ endif
! decimal headers
- call getfhv('dist',dist,nerr)
- call getfhv('az',az,nerr)
- call getfhv('baz',baz,nerr)
- call getfhv('stlo',slon,nerr)
- call getfhv('stla',slat,nerr)
+ call getfhv('dist',tmp,nerr)
+ if(nerr .ne. 0) then
+ write(*,*)'Error reading variable: dist'
+ call exit(-1)
+ endif
+ dist = tmp
+ call getfhv('az',tmp,nerr)
+ if(nerr .ne. 0) then
+ write(*,*)'Error reading variable: az'
+ call exit(-1)
+ endif
+ az = tmp
+
+ call getfhv('baz',tmp,nerr)
+ if(nerr .ne. 0) then
+ write(*,*)'Error reading variable: baz'
+ call exit(-1)
+ endif
+ baz = tmp
+
+ call getfhv('stlo',tmp,nerr)
+ if(nerr .ne. 0) then
+ write(*,*)'Error reading variable: stlo'
+ call exit(-1)
+ endif
+ slon = tmp
+
+ call getfhv('stla',tmp,nerr)
+ if(nerr .ne. 0) then
+ write(*,*)'Error reading variable: stla'
+ call exit(-1)
+ endif
+ slat = tmp
+
!!$ ! integer header variables
-!!$ call saclst_iheader_f(data_file,'nzyear', yr)
-!!$ call saclst_iheader_f(data_file,'nzjday', jda)
-!!$ call saclst_iheader_f(data_file,'nzhour', ho)
-!!$ call saclst_iheader_f(data_file,'nzmin', mi)
-!!$ call saclst_iheader_f(data_file,'nzsec', nsec)
-!!$ call saclst_iheader_f(data_file,'nzmsec', msec)
-!!$
+!!$ call saclst_iheader_f(data_file,'nzyear', yr)
+!!$ call saclst_iheader_f(data_file,'nzjday', jda)
+!!$ call saclst_iheader_f(data_file,'nzhour', ho)
+!!$ call saclst_iheader_f(data_file,'nzmin', mi)
+!!$ call saclst_iheader_f(data_file,'nzsec', nsec)
+!!$ call saclst_iheader_f(data_file,'nzmsec', msec)
+!!$
!!$ sec=nsec+msec/1000.0
!!$
!!$ call saclst_kheader_f(data_file,'knetwk',ntw,klen)
!!$ call saclst_kheader_f(data_file,'kstnm', sta,klen)
!!$ call saclst_kheader_f(data_file,'kcmpnm',comp,klen)
-!!$
+!!$
!!$ call dsaclst_fheader_f(data_file,'dist',dist)
!!$ call dsaclst_fheader_f(data_file,'az', az)
!!$ call dsaclst_fheader_f(data_file,'baz', baz)
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub2.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub2.f90 2011-02-16 21:49:17 UTC (rev 17885)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub2.f90 2011-02-17 01:13:56 UTC (rev 17886)
@@ -1,5 +1,8 @@
module mt_sub2
+
+ use mt_constants
+
implicit none
! TOLERRANCE CONTROL
@@ -24,7 +27,7 @@
double precision :: zzign,zign,flx,v
integer :: lblock,k,fk,jh,ii,istart
integer :: l,iblock,nblock,i,lbhalf,j,lx
-
+
! sign must be +1. or -1.
if(zzign >= 0.) then
zign = 1.
@@ -33,6 +36,12 @@
endif
lx = 2**n
+
+ ! checks bounds
+ if( lx > NPT ) stop 'error fft increase NPT, or decrease n'
+
+
+
do 1 i=1,n
1 m(i) = 2**(n-i)
do 4 l=1,n
@@ -52,6 +61,9 @@
do 2 i=1,lbhalf
j = istart+i
jh = j+lbhalf
+ ! checks bounds
+ if( jh < 1 .or. jh > NPT ) stop 'error fft bounds'
+
q = xi(jh)*wk
xi(jh) = xi(j)-q
xi(j) = xi(j)+q
@@ -66,6 +78,8 @@
do 7 j=1,lx
if(k < j) go to 5
hold = xi(j)
+ ! checks bounds
+ if( k+1 < 1 .or. k+1 > NPT ) stop 'error fft k bounds'
xi(j) = xi(k+1)
xi(k+1) = hold
5 do 6 i=1,n
@@ -274,7 +288,7 @@
double precision, dimension(ndim,*) :: r
integer :: nt,n,ndim,nev,ipar
- double precision, dimension(ndim) :: bb
+ !double precision, dimension(ndim) :: bb
double precision :: q,el,elam,u,umeps,x,ddot,rnorm
integer :: i,j,ik,iag,m,jk,jm1
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_variables.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_variables.f90 2011-02-16 21:49:17 UTC (rev 17885)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_variables.f90 2011-02-17 01:13:56 UTC (rev 17886)
@@ -1,7 +1,7 @@
module mt_variables
use mt_constants
-!
+!
! multi-taper measurements
!
! Ying Zhou: The fit between the recovered data and the data can be improved
@@ -26,3 +26,31 @@
logical :: DISPLAY_DETAILS,OUTPUT_MEASUREMENT_FILES,RUN_BANDPASS,COMPUTE_ADJOINT_SOURCE
end module mt_variables
+
+
+module mt_weighting
+
+! module for weighting/normalizing measurements
+
+ logical,parameter :: DO_WEIGHTING = .false.
+
+ ! transverse, radial and vertical weights
+ double precision :: weight_T, weight_R, weight_Z
+ ! body waves: number of picks on vertical, radial and transverse component
+ double precision :: num_P_SV_V,num_P_SV_R,num_SH_T
+ ! surface waves: number of pick on vertical, radial and transverse
+ double precision :: num_Rayleigh_V,num_Rayleigh_R,num_Love_T
+
+ ! typical surface wave speed in km/s, to calculate surface wave arrival times
+ ! Love waves faster than Rayleigh
+ double precision, parameter :: surface_vel = 4.0
+
+ ! wave type pick
+ integer, parameter :: P_SV_V = 1
+ integer, parameter :: P_SV_R = 2
+ integer, parameter :: SH_T = 3
+ integer, parameter :: Rayleigh_V = 4
+ integer, parameter :: Rayleigh_R = 5
+ integer, parameter :: Love_T = 6
+
+end module mt_weighting
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.pl 2011-02-16 21:49:17 UTC (rev 17885)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.pl 2011-02-17 01:13:56 UTC (rev 17886)
@@ -37,6 +37,8 @@
#open(MATLAB,"> plot1.m");
# NOTE: Type "which matlab" to make sure that matlab is not alised.
open(MATLAB,"|/opt/matlab/bin/matlab -nosplash -nojvm") || die("Error opening matlab program\n");
+#open(MATLAB,"|/usr/local/bin/matlab -nosplash -nojvm") || die("Error opening matlab program\n");
+
print MATLAB "P=path; path(P,'$progdir');\n";
print MATLAB "amps=[];phs=[];dlnAs=[];dts=[];gcarcs=[];azs=[];stats=[];\n";
@@ -48,7 +50,7 @@
($filebase) = split("dlnA",`basename $dlnA_files[$i]`); chop($filebase);
($inputdir) = split(" ",`dirname $dlnA_files[$i]`);
- `sac2ascii.pl $inputdir/$filebase.obs.sac $inputdir/$filebase.syn.sac $inputdir/$filebase.recon_syn.sac $inputdir/$filebase.recon_syn_dt.sac $inputdir/$filebase.recon_syn_cc.sac $inputdir/$filebase.recon_syn_cc_dt.sac`;
+ `sac2ascii.pl $inputdir/$filebase.obs.sac $inputdir/$filebase.syn.sac $inputdir/$filebase.recon_syn.sac $inputdir/$filebase.recon_syn_dt.sac $inputdir/$filebase.recon_syn_cc.sac $inputdir/$filebase.recon_syn_cc_dt.sac`;
$epsfile = "$inputdir/${filebase}_measure.eps";
(undef,$stname,$stnetwk,$comp,$az,$gcarc)=split(" ",`saclst kstnm knetwk kcmpnm az gcarc f $inputdir/${filebase}.obs.sac`);
@@ -59,7 +61,7 @@
print MATLAB "dataf=dlmread('$inputdir/$filebase.obs.power');\n";
print MATLAB "synf=dlmread('$inputdir/$filebase.syn.power');\n";
print MATLAB "new=dlmread('$inputdir/$filebase.recon_syn');\n";
- print MATLAB "new_dt=dlmread('$inputdir/$filebase.recon_syn_dt');\n";
+ print MATLAB "new_dt=dlmread('$inputdir/$filebase.recon_syn_dt');\n";
if ($iall == 1) {
# read sub-sampled error curves
print MATLAB "dlnA=dlmread('$inputdir/$filebase.err_dlnA');\n";
@@ -78,11 +80,11 @@
# read multitaper average measurements
print MATLAB "dlnA_ave=dlmread('$inputdir/$filebase.dlnA_average');\n";
- print MATLAB "dt_ave =dlmread('$inputdir/$filebase.dt_average');\n";
+ print MATLAB "dt_ave =dlmread('$inputdir/$filebase.dt_average');\n";
# read cross-correlation measurements
print MATLAB "dlnA_cc=dlmread('$inputdir/$filebase.dlnA_cc');\n";
- print MATLAB "dt_cc =dlmread('$inputdir/$filebase.dt_cc');\n";
+ print MATLAB "dt_cc =dlmread('$inputdir/$filebase.dt_cc');\n";
# read cross-correlation reconstruction
print MATLAB "new_cc=dlmread('$inputdir/$filebase.recon_syn_cc');\n";
@@ -98,7 +100,11 @@
print MATLAB "print ('-depsc2','$epsfile');\n";
print MATLAB "close;\n";
# print MATLAB "plot_mtm_quality('$inputdir');\n";
+
+ print "\nplotted to: $epsfile \n";
+
}
print MATLAB "\n";
close(MATLAB);
+
More information about the CIG-COMMITS
mailing list