[cig-commits] r15545 - in seismo/3D/ADJOINT_TOMO/measure_adj: . scripts_tomo
carltape at geodynamics.org
carltape at geodynamics.org
Sun Aug 16 03:59:55 PDT 2009
Author: carltape
Date: 2009-08-16 03:59:55 -0700 (Sun, 16 Aug 2009)
New Revision: 15545
Modified:
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/scripts_tomo/run_tomo.pl
Log:
Modified MT code so that the transfer function is applied between the windowed synthetics and the CC-deconstructed windowed data, i.e., shift the data by -DT and amplify by -DlnA. Also ensured consistency between the DlnA conventions for CC and MT, with DlnA = ln(Aobs/Asyn). Simplified several aspects of the code, and also renamed many of the variables in mt_sub.f90. Tested changes with 234 socal events.
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-08-14 22:13:56 UTC (rev 15544)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-08-16 10:59:55 UTC (rev 15545)
@@ -31,7 +31,7 @@
double precision :: tshift_xc, sigma_dt_cc, dlnA, sigma_dlnA_cc, sigma_dt, sigma_dlnA
double precision :: tr_chi, am_chi, cc_max, T_pmax_dat, T_pmax_syn
!double precision :: tshift_f1f2, cc_max_f1f2
- double precision, dimension(NPT) :: dtau_w, dlnA_w, err_dt, err_dlnA, dzr0, dzr20
+ 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
@@ -260,8 +260,8 @@
! make measurements
! CHT: also compute reconstructed synthetics for CC (and MT, if specified) measurements
- call mt_measure(is_mtm,iker,datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,istart,dzr20,dzr0,nlen,&
- tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tseis_recon_cc,&
+ call mt_measure(is_mtm,iker,datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,&
+ istart,data_dtw,syn_dtw,nlen,tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tseis_recon_cc,&
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
@@ -298,7 +298,7 @@
print *, 'Reverting from multitaper measurement to cross-correlation measurement'
iker = 2
is_mtm = 3
- call mt_measure(is_mtm,iker,datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,istart,dzr20,dzr0,nlen,&
+ call mt_measure(is_mtm,iker,datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,istart,data_dtw,syn_dtw,nlen,&
tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tseis_recon_cc,&
i_pmax_dat,i_pmax_syn,i_right,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA)
use_trace = .true.
@@ -335,10 +335,10 @@
tr_chi = 0.0 ; am_chi = 0.0 ! must be initialized
if (iker==0) then
- call mt_adj(iker,istart,dzr20,dzr0,nlen,dt,tshift_xc,dlnA,sigma_dt_cc,sigma_dlnA_cc,&
+ call mt_adj(iker,istart,data_dtw,syn_dtw,nlen,dt,tshift_xc,dlnA,sigma_dt_cc,sigma_dlnA_cc,&
dtau_w,dlnA_w,err_dt,err_dlnA,sigma_dt,sigma_dlnA,i_left,i_right,tr_adj_src,tr_chi,window_chi)
else
- call mt_adj(iker,istart,dzr20,dzr0,nlen,dt,tshift_xc,dlnA,sigma_dt_cc,sigma_dlnA_cc,&
+ call mt_adj(iker,istart,data_dtw,syn_dtw,nlen,dt,tshift_xc,dlnA,sigma_dt_cc,sigma_dlnA_cc,&
dtau_w,dlnA_w,err_dt,err_dlnA,sigma_dt,sigma_dlnA,i_left,i_right,tr_adj_src,tr_chi,window_chi,am_adj_src,am_chi)
endif
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-08-14 22:13:56 UTC (rev 15544)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-08-16 10:59:55 UTC (rev 15545)
@@ -17,56 +17,59 @@
! 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)
- ! data(:), syn(:) t0, dt, npts -- original data and synthetics array
+ ! 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)
! Output:
! istart -- starting index of the windowed portion of original trace
- ! dzr20(:), dzr0(:), nlen -- windowed and shifted data, windowed synthetics
- ! tshift_xc, dlnA, cc_max -- time shift and amplitude cross-correlation measurements
+ ! dat_dtw(:), syn_dtw(:), nlen -- windowed and shifted data, windowed synthetics
+ ! tshift, dlnA, cc_max -- time shift and amplitude cross-correlation measurements
! i_right -- the maximum reliable frequency estimate index
! trans_w(:) -- estimates of transfer function
! 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 Fortran 77 by Ying Zhou
- ! upgraded to Fortran 90 by Alessia Maggi
+ ! 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
!
! Evolution of data and syn arrays:
- ! data -- window + time taper --> dzr2 -- time shift --> dzr20 -- complex FFT --> wseis_dat
- ! dzr20 -- single-tapered --> dzr2 -- complex FFT --> wseis2
- ! syn -- window + time taper --> dzr = dzr20 -- complex FFT --> wseis_syn
- ! dzr0 -- single-tapered --> dzr -- complex FFT --> wseis
+ ! data -- window + time taper --> dataw -- time shift --> dataw0 -- complex FFT --> wseis_dat
+ ! dataw0 -- single-tapered --> dataw -- complex FFT --> wseis2
+ ! syn -- window + time taper --> synw = dataw0 -- complex FFT --> wseis_syn
+ ! synw0 -- single-tapered --> synw -- complex FFT --> wseis
! =================================================================================================
- subroutine mt_measure(is_mtm,iker,datafile,file_prefix,data,syn,t0,dt,npts,tstart,tend, &
- istart,dzr20,dzr0,nlen,tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tseis_recon_cc, &
+ subroutine mt_measure(is_mtm,iker,datafile,file_prefix,dat_dt,syn_dt,t0,dt,npts,tstart,tend, &
+ istart,dat_dtw,syn_dtw,nlen,tshift,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,syn_dtw_cc, &
i_pmax_dat,i_pmax_syn,i_right,trans_w,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,err_dt,err_dlnA)
- integer, intent(in) :: is_mtm, iker, npts
- double precision, dimension(:), intent(in) :: data, syn
- double precision, intent(in) :: tstart, tend, t0, dt
- character(len=150), intent(in) :: file_prefix, datafile
+ integer, intent(in) :: is_mtm,iker,npts
+ double precision, dimension(:), intent(in) :: dat_dt,syn_dt
+ double precision, intent(in) :: tstart,tend,t0,dt
+ character(len=150), intent(in) :: file_prefix,datafile
- double precision, intent(out) :: tshift_xc, sigma_dt_cc, dlnA, sigma_dlnA_cc, cc_max, sigma_dt,sigma_dlnA
+ double precision, intent(out) :: tshift,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,sigma_dt,sigma_dlnA
complex*16, dimension(:), intent(out) :: trans_w
- double precision, dimension(:), intent(out) :: dtau_w, dlnA_w, dzr0, dzr20
- 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, 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) :: dzr, dzr2, dzr3, tseis_recon, tseis_recon_dt, tseis_recon_cc, tseis_recon_cc_dt
- double precision, dimension(NPT) :: syn_displ, dat_displ, syn_veloc, dat_veloc, syn_accel
- double precision :: sfac1, fac, omega, f0, df, df_new, dw, ampmax_unw, wtr_use_unw, ampmax, wtr_use, wtr_mtm, dtau_wa, dlnA_wa
- integer :: ishift, i, ictaper, j, k, fnum, i_amp_max_unw, i_amp_max, i_right_stop, idf_new, iom, iker_temp
+ 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,iker_temp
- complex*16, dimension(NPT) :: wseis_syn, wseis_dat, top_mtm, bot_mtm, wseis, wseis2, trans_mtm, wseis_recon
- 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
+ complex*16, dimension(NPT) :: syn_dtwo, dat_dtwo, syn_dtw_ho, dat_dtw_ho, &
+ top_mtm, bot_mtm, trans_mtm, wseis_recon
+ 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
+ logical :: output_logical,display_logical
!-------------------------------------------------------------
@@ -82,16 +85,16 @@
filename = trim(file_prefix)
if (DISPLAY_DETAILS) then
- call dwascii(trim(file_prefix)//'_data',data,npts,t0,dt)
- call dwascii(trim(file_prefix)//'_syn',syn,npts,t0,dt)
- !call dwsac1(trim(file_prefix)//'_data.sac',data,npts,t0,dt)
- !call dwsac1(trim(file_prefix)//'_syn.sac',syn,npts,t0,dt)
+ 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,data)
-!!$ call dwrite_sacfile_f(datafile,trim(file_prefix)//'_syn.sac',t0,npts,syn)
-!!$ !call dwrite_sacfile_f(datafile,trim(file_prefix)//'_diff_dms.sac',t0,npts,data-syn)
-!!$ !call dwrite_ascfile_f(trim(file_prefix)//'_data.txt',t0,dt,npts,data)
-!!$ !call dwrite_ascfile_f(trim(file_prefix)//'_syn.txt',t0,dt,npts,syn)
+!!$ 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)
endif
!--------------------------------------------------------------------------
@@ -99,7 +102,7 @@
!--------------------------------------------------------------------------
! interpolate data and synthetics, and also extract time-windowed records
- call interpolate_data_and_syn(data, syn, tstart, tend, t0, dt, NPT, dzr2, dzr, nlen, istart)
+ 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'
@@ -114,72 +117,67 @@
!fac = 1 - sfac1*((i-1) - dble(nlen)/2.)**2 ! welch window
fac = 1. - cos(PI*(i-1)/(nlen-1))**ipwr_t ! cosine window
- dzr(i) = dzr(i) * fac ! syn, windowed
- dzr2(i) = dzr2(i) * fac ! dat, windowed
+ syn_dtw(i) = syn_dtw(i) * fac ! syn, windowed
+ dat_dtw(i) = dat_dtw(i) * fac ! dat, windowed
enddo
if (DISPLAY_DETAILS) then
print *, ' NPTs = ', NPT, ' new nlen = ', nlen
- call dwsac1(trim(file_prefix)//'.obs.sac',dzr2,nlen,tstart,dt)
- call dwsac1(trim(file_prefix)//'.syn.sac',dzr,nlen,tstart,dt)
-!!$ call dwrite_sacfile_f(datafile,trim(file_prefix)//'.obs.sac',tstart,nlen,dzr2)
-!!$ call dwrite_sacfile_f(datafile,trim(file_prefix)//'.syn.sac',tstart,nlen,dzr)
+ 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)
endif
! save a copy of the windowed data
- dzr3(:) = dzr2(:)
+ !dat_dtwc(:) = dat_dtw(:)
!------------------------------------------------------------------
! cross-correlation traveltime and amplitude measurements
!------------------------------------------------------------------
! compute cross-correltaion time shift and also amplitude measurmement
+ ! NOTE: records have already been windowed, so no information outside windows is considered
! LQY: Ying suggested to align them at relatively long periods
- call compute_cc(dzr, dzr2, nlen, dt, ishift, tshift_xc, dlnA, cc_max)
+ call compute_cc(syn_dtw, dat_dtw, nlen, dt, ishift, tshift, dlnA, cc_max)
- ! apply time shift to OBSERVED seismogram
- dzr0(:) = dzr(:)
- dzr20(:) = dzr2(:)
- do i = 1, nlen
- if ((ishift+i) > 1 .and. (ishift+i) < nlen ) dzr20(i) = dzr2(i+ishift)
- enddo
- ! fill the missing time window with the endpoint value
- if (ishift < 0) dzr20(1:-ishift+1) = dzr20(-ishift+2)
- if (ishift > 0) dzr20(nlen-ishift:nlen) = dzr20(nlen-ishift-1)
-
- !if (DISPLAY_DETAILS) then
- ! call dwrite_sacfile_f(datafile,'windowed_shifted_data.sac',tstart,nlen,dzr20)
- !endif
-
- ! compute velocity of synthetics and data
- syn_displ(1:nlen) = dzr0(1:nlen)
- dat_displ(1:nlen) = dzr20(1:nlen)
-
+ ! compute velocity of synthetics
do i = 2, nlen-1
- syn_veloc(i) = (syn_displ(i+1) - syn_displ(i-1)) / (2.0*dt)
- dat_veloc(i) = (dat_displ(i+1) - dat_displ(i-1)) / (2.0*dt)
+ syn_vtw(i) = (syn_dtw(i+1) - syn_dtw(i-1)) / (2.0*dt)
enddo
- syn_veloc(1) = (syn_displ(2) - syn_displ(1)) / dt
- syn_veloc(nlen) = (syn_displ(nlen) - syn_displ(nlen-1)) /dt
- dat_veloc(1) = (dat_displ(2) - dat_displ(1)) / dt
- dat_veloc(nlen) = (dat_displ(nlen) - dat_displ(nlen-1)) /dt
+ syn_vtw(1) = (syn_dtw(2) - syn_dtw(1)) / dt
+ syn_vtw(nlen) = (syn_dtw(nlen) - syn_dtw(nlen-1)) /dt
- do i = 2, nlen-1
- syn_accel(i) = (syn_veloc(i+1) - syn_veloc(i-1)) / (2.0*dt)
- enddo
- syn_accel(1) = (syn_veloc(2) - syn_veloc(1)) / dt
- syn_accel(nlen) = (syn_veloc(nlen) - syn_veloc(nlen-1)) / dt
+ ! 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
- call reconstruct_syn_cc(datafile,filename,dzr,dzr2,tstart,dt,nlen,ishift,tshift_xc,dlnA, tseis_recon_cc, tseis_recon_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)
+ endif
+
! compute the estimated uncertainty for the cross-correlation measurment
sigma_dt_cc = 1.
sigma_dlnA_cc = 1.
- call compute_average_error(dzr3, tseis_recon_cc, tseis_recon_cc_dt, nlen, dt, sigma_dt_cc, sigma_dlnA_cc)
+ 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(file_prefix, IKER_CC, tshift_xc,dlnA,sigma_dt_cc,sigma_dlnA_cc)
+ call write_average_meas(filename,IKER_CC,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc)
!========================================
@@ -210,23 +208,23 @@
endif
enddo
- ! create complex synthetic seismogram and shifted-data seismogram
- wseis_syn = cmplx(0.,0.)
- wseis_dat = cmplx(0.,0.)
- !wseis_syn(1:nlen) = dzr0(1:nlen)
- !wseis_dat(1:nlen) = dzr20(1:nlen)
- wseis_syn(1:nlen) = cmplx(dzr0(1:nlen))
- wseis_dat(1:nlen) = cmplx(dzr20(1:nlen))
+ ! 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,wseis_syn,FORWARD_FFT,dt)
- call fft(LNPT,wseis_dat,FORWARD_FFT,dt)
+ 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
ampmax_unw = 0.
i_pmax_dat = 1
do i = 1, fnum ! loop over frequencies
- if( abs(wseis_dat(i)) > ampmax_unw) then
- ampmax_unw = abs(wseis_dat(i))
+ if( abs(dat_dtwo(i)) > ampmax_unw) then
+ ampmax_unw = abs(dat_dtwo(i))
i_pmax_dat = i
endif
enddo
@@ -235,8 +233,8 @@
! used to determine the i_right values (maximum frequency for measurement)
ampmax_unw = 0.
do i = 1, fnum ! loop over frequencies
- if( abs(wseis_syn(i)) > ampmax_unw) then
- ampmax_unw = abs(wseis_syn(i))
+ if( abs(syn_dtwo(i)) > ampmax_unw) then
+ ampmax_unw = abs(syn_dtwo(i))
i_amp_max_unw = i
endif
enddo
@@ -248,11 +246,11 @@
i_right = fnum
i_right_stop = 0
do i = 1,fnum
- if( abs(wseis_syn(i)) <= abs(wtr_use_unw) .and. i_right_stop==0 .and. i > i_amp_max_unw ) then
+ 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(wseis_syn(i)) >= 10.*abs(wtr_use_unw) .and. i_right_stop==1 .and. i > i_amp_max_unw) then
+ 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
@@ -266,12 +264,12 @@
print *, ' i_right = ', i_right, ', stopping freq = ', sngl(i_right * df)
! write out power for each signal
- call dwascii(trim(file_prefix)//'.obs.power',abs(wseis_dat(1:i_right)),i_right,df,df)
- call dwascii(trim(file_prefix)//'.syn.power',abs(wseis_syn(1:i_right)),i_right,df,df)
- !call dwsac1(trim(file_prefix)//'.obs.power.sac',abs(wseis_dat(1:i_right)),i_right,df,df)
- !call dwsac1(trim(file_prefix)//'.syn.power.sac',abs(wseis_syn(1:i_right)),i_right,df,df)
-!!$ call dwrite_ascfile_f(trim(file_prefix)//'.obs.power',df,df,i_right,abs(wseis_dat(1:i_right)) )
-!!$ call dwrite_ascfile_f(trim(file_prefix)//'.syn.power',df,df,i_right,abs(wseis_syn(1:i_right)) )
+ 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)) )
endif
@@ -303,27 +301,28 @@
do ictaper = 1, ntaper
- wseis(:) = cmplx(0.,0.) ! note: this has to be initialized inside the loop
- wseis2(:) = cmplx(0.,0.)
+ syn_dtw_ho(:) = cmplx(0.,0.) ! note: this has to be initialized inside the loop
+ dat_dtw_ho(:) = cmplx(0.,0.)
+ ! apply time-domain taper
do i = 1, nlen
- dzr(i) = dzr0(i) * tas(i,ictaper) ! single-tapered, windowed syn
- dzr2(i) = dzr20(i) * tas(i,ictaper) ! single-tapered, windowed, shifted data
+ syn_dtw_h(i) = syn_dtw(i) * tas(i,ictaper) ! single-tapered, windowed syn
+ dat_dtw_h(i) = dat_dtw_cc(i) * tas(i,ictaper) ! single-tapered, windowed, shifted data
enddo
- wseis(1:nlen) = cmplx(dzr(1:nlen),0.) ! syn
- wseis2(1:nlen) = cmplx(dzr2(1:nlen),0.) ! data
+ syn_dtw_ho(1:nlen) = cmplx(syn_dtw_h(1:nlen),0.)
+ dat_dtw_ho(1:nlen) = cmplx(dat_dtw_h(1:nlen),0.)
! apply FFT to get complex spectra
- call fft(LNPT,wseis, FORWARD_FFT,dt) ! syn
- call fft(LNPT,wseis2,FORWARD_FFT,dt) ! data
+ call fft(LNPT,syn_dtw_ho,FORWARD_FFT,dt)
+ 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
ampmax = 0.
do i = 1, fnum ! loop over frequencies
- if( abs(wseis(i)) > ampmax) then ! syn, single_tapered
- ampmax = abs(wseis(i))
+ if( abs(syn_dtw_ho(i)) > ampmax) then ! syn, single_tapered
+ ampmax = abs(syn_dtw_ho(i))
i_amp_max = i
endif
enddo
@@ -332,13 +331,13 @@
! calculate top and bottom of MT transfer function
do i = 1, fnum
- top_mtm(i) = top_mtm(i) + wseis2(i) * conjg(wseis(i)) ! uses data and syn
- bot_mtm(i) = bot_mtm(i) + wseis(i) * conjg(wseis(i)) ! uses syn only
+ top_mtm(i) = top_mtm(i) + dat_dtw_ho(i) * conjg(syn_dtw_ho(i)) ! uses data and syn
+ bot_mtm(i) = bot_mtm(i) + syn_dtw_ho(i) * conjg(syn_dtw_ho(i)) ! uses syn only
! calculate transfer function for single taper measurement using water level
if (is_mtm /= 1) then
- if(abs(wseis(i)) > abs(wtr_use)) trans_w(i) = wseis2(i) / wseis(i)
- if(abs(wseis(i)) <= abs(wtr_use)) trans_w(i) = wseis2(i) / (wseis(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
@@ -348,10 +347,12 @@
! 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_xc, phi_w,abs_w,dtau_w,dlnA_w,dtau_wa,dlnA_wa)
- call reconstruct_syn(datafile,filename,wseis_syn,wvec,dtau_w,dlnA_w,i_right,tstart,dt,nlen, tseis_recon, tseis_recon_dt)
- !call check_recon_quality(filename,dzr20,dzr0,dzr3,tseis_recon,nlen,dt,tshift_xc,tshift_f1f2,cc_max_f1f2,cc_max)
- !call compute_average_error(dzr3,tseis_recon,tseis_recon_dt,nlen,dt,sigma_dt,sigma_dlnA)
+ 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, IKER_MT, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
endif
@@ -382,17 +383,20 @@
enddo
! multitaper phase, abs, tt, and amp (freq)
- call write_trans(filename,trans_mtm,wvec,fnum,i_right,idf_new,df,tshift_xc,phi_mtm,abs_mtm,dtau_mtm,dlnA_mtm,dtau_wa,dlnA_wa)
+ 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)
! apply transfer function to the syn
- call reconstruct_syn(datafile,filename,wseis_syn,wvec,dtau_mtm,dlnA_mtm,i_right,tstart,dt,nlen,tseis_recon,tseis_recon_dt)
+ call reconstruct_syn(filename,syn_dtwo,wvec,dtau_mtm,dlnA_mtm, &
+ i_right,tstart,dt,nlen,syn_dtw_mt,syn_dtw_mt_dt)
! check quality
- !call check_recon_quality(filename,dzr20,dzr0,dzr3,tseis_recon,nlen,dt,tshift_xc, tshift_f1f2, cc_max_f1f2,cc_max)
+ !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(dzr3,tseis_recon,tseis_recon_dt,nlen,dt,sigma_dt,sigma_dlnA)
+ !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, iker, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
@@ -418,30 +422,30 @@
do iom = 1, ntaper
- top_mtm(:) = cmplx(0,0)
- bot_mtm(:) = cmplx(0,0)
+ top_mtm(:) = cmplx(0.,0.)
+ bot_mtm(:) = cmplx(0.,0.)
do ictaper = 1, ntaper
if(ictaper.eq.iom) cycle
- ! apply ictaper'th taper
- dzr(1:nlen) = dzr0(1:nlen) * tas(1:nlen,ictaper)
- dzr2(1:nlen) = dzr20(1:nlen) * tas(1:nlen,ictaper)
+ ! apply ictaper-th taper
+ syn_dtw_h(1:nlen) = syn_dtw(1:nlen) * tas(1:nlen,ictaper)
+ dat_dtw_h(1:nlen) = dat_dtw_cc(1:nlen) * tas(1:nlen,ictaper)
! complex tapered series
- wseis(:) = cmplx(0,0)
- wseis2(:) = cmplx(0,0)
- wseis(1:nlen) = cmplx(dzr(1:nlen),0)
- wseis2(1:nlen) = cmplx(dzr2(1:nlen),0)
+ syn_dtw_ho(:) = cmplx(0.,0.)
+ dat_dtw_ho(:) = cmplx(0.,0.)
+ syn_dtw_ho(1:nlen) = cmplx(syn_dtw_h(1:nlen),0.)
+ dat_dtw_ho(1:nlen) = cmplx(dat_dtw_h(1:nlen),0.)
! apply f.t. to get complex spectra
- call fft(LNPT,wseis,FORWARD_FFT,dt)
- call fft(LNPT,wseis2,FORWARD_FFT,dt)
+ call fft(LNPT,syn_dtw_ho,FORWARD_FFT,dt)
+ call fft(LNPT,dat_dtw_ho,FORWARD_FFT,dt)
! calculate top and bottom of Jacknife transfer function
do i = 1, fnum
- top_mtm(i) = top_mtm(i) + wseis2(i) * conjg(wseis(i))
- bot_mtm(i) = bot_mtm(i) + wseis(i) * conjg(wseis(i))
+ top_mtm(i) = top_mtm(i) + dat_dtw_ho(i) * conjg(syn_dtw_ho(i))
+ bot_mtm(i) = bot_mtm(i) + syn_dtw_ho(i) * conjg(syn_dtw_ho(i))
enddo
enddo ! ictaper
@@ -461,8 +465,8 @@
if(abs(bot_mtm(i)).le.abs(wtr_use)) trans_mtm(i) = top_mtm(i) /(bot_mtm(i)+wtr_use)
enddo
- call write_trans(filename,trans_mtm,wvec,fnum,i_right,idf_new,df, &
- tshift_xc,phi_mul(:,iom),abs_mul(:,iom),dtau_mul(:,iom),dlnA_mul(:,iom))
+ 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
@@ -492,7 +496,6 @@
eabs2_ave = 0.
do iom = 1, ntaper
-
eph_iom = ntaper*phi_mtm(i) - (ntaper-1)*phi_mul(i,iom)
eph_ave = eph_ave + eph_iom
@@ -504,7 +507,6 @@
eabs2_iom = ntaper*dlnA_mtm(i) - (ntaper-1)*dlnA_mul(i,iom)
eabs2_ave = eabs2_ave + eabs2_iom
-
enddo
eph_ave = eph_ave / (ntaper)
@@ -576,9 +578,9 @@
! iker -- 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
- ! dzr20(:), dzr0(:), nlen, dt -- windowed, shifted data, windowed synthetics with length nlen and
- ! sampling rate dt
- ! tshift_xc, dlnA -- cross-correlation traveltime and amplitude measurements
+ ! 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)
!
@@ -590,13 +592,13 @@
! original coding by Carl Tape, finalized by Qinya Liu
! ======================================================================================================
- subroutine mt_adj(iker,istart,dzr20,dzr0,nlen,dt,tshift_xc,dlnA,sigma_dt_cc,sigma_dlnA_cc,&
+ subroutine mt_adj(iker,istart,dat_dtw,syn_dtw,nlen,dt,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc,&
dtau_w,dlnA_w,err_dtau,err_dlnA,sigma_dt,sigma_dlnA,i_left,i_right, &
tr_adj_src,tr_chi,window_chi,am_adj_src,am_chi)
integer, intent(in) :: iker, istart, nlen, i_left, i_right
- double precision, dimension(:), intent(in) :: dzr20, dzr0
- double precision, intent(in) :: dt, tshift_xc, dlnA, sigma_dt_cc, sigma_dlnA_cc, sigma_dt, sigma_dlnA
+ double precision, dimension(:), intent(in) :: dat_dtw, syn_dtw
+ double precision, intent(in) :: dt, tshift, dlnA, sigma_dt_cc, sigma_dlnA_cc, sigma_dt, sigma_dlnA
double precision, dimension(:), intent(in) :: dtau_w, dlnA_w, err_dtau, err_dlnA
double precision, dimension(:), intent(out) :: tr_adj_src
@@ -605,13 +607,13 @@
double precision, dimension(:), intent(out), optional :: am_adj_src
double precision, intent(out), optional :: am_chi
- double precision, dimension(NPT) :: vzr0, vzr, dzr, dzr2, dzr30, ey1, ey2
+ 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
double precision :: df,Nnorm,Mnorm,fac,ffac,w_taper(NPT), time_window(NPT)
double precision, dimension(:,:), allocatable :: tas
- complex*16, dimension(:,:),allocatable :: wseis, wseis2
+ complex*16, dimension(:,:),allocatable :: syn_dtw_ho_all, syn_vtw_ho_all
complex*16, dimension(NPT) :: pwc_adj,qwc_adj
double precision, dimension(NPT) :: dtau_pj_t, dlnA_qj_t
double precision :: dtau_wtr, dlnA_wtr, err_t, err_A
@@ -629,7 +631,7 @@
! ----------------------
! frequency-domain tapers
- ! THIS CHOICE CAN HAVE A SUBSTANTIAL EFFECT ON THE ADJOINT SOURCES
+ ! THIS CHOICE WILL HAVE AN EFFECT ON THE ADJOINT SOURCES
ipwr_w = 10
w_taper(:) = 0.
do i = i_left, i_right ! CHT: 1 --> i_left
@@ -699,19 +701,19 @@
! compute synthetic velocity
do i = 2, nlen-1
- vzr0(i) = (dzr0(i+1) - dzr0(i-1)) / (2.0*dt)
+ syn_vtw(i) = (syn_dtw(i+1) - syn_dtw(i-1)) / (2.0*dt)
enddo
- vzr0(1) = (dzr0(2) - dzr0(1)) / dt
- vzr0(nlen) = (dzr0(nlen) - dzr0(nlen-1)) / dt
+ 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( vzr0(1:nlen) * vzr0(1:nlen) )
- ft_bar_t(1:nlen) = -vzr0(1:nlen) / Nnorm
+ 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( dzr0(1:nlen) * dzr0(1:nlen) )
- fa_bar_t(1:nlen) = dzr0(1:nlen) / Mnorm
+ Mnorm = dt * sum( syn_dtw(1:nlen) * syn_dtw(1:nlen) )
+ fa_bar_t(1:nlen) = syn_dtw(1:nlen) / Mnorm
! -------------------------------
! MULTITAPER ADJOINT SOURCES
@@ -722,8 +724,8 @@
! allocate MT variables
ntaper = int(NPI * 2.0)
allocate(tas(NPT,ntaper))
- allocate(wseis(NPT,ntaper))
- allocate(wseis2(NPT,ntaper))
+ allocate(syn_dtw_ho_all(NPT,ntaper))
+ allocate(syn_vtw_ho_all(NPT,ntaper))
! get the MT tapers
call staper(nlen, NPI, NTAPER, tas, NPT, ey1, ey2)
@@ -735,27 +737,27 @@
do ictaper = 1,ntaper
! tapered synthetic displacement
- dzr(1:nlen) = dzr0(1:nlen) * tas(1:nlen,ictaper)
+ syn_dtw_h(1:nlen) = syn_dtw(1:nlen) * tas(1:nlen,ictaper)
- ! note: cannot directly apply tapers on vzr0, which is not the same as vzr
+ ! compute velocity of tapered syn
do i = 2, nlen-1
- vzr(i) = (dzr(i+1) - dzr(i-1)) / (2.0*dt)
+ syn_vtw_h(i) = (syn_dtw_h(i+1) - syn_dtw_h(i-1)) / (2.0*dt)
enddo
- vzr(1) = (dzr(2) - dzr(1)) / dt
- vzr(nlen) = (dzr(nlen) - dzr(nlen-1)) /dt
+ syn_vtw_h(1) = (syn_dtw_h(2) - syn_dtw_h(1)) / dt
+ syn_vtw_h(nlen) = (syn_dtw_h(nlen) - syn_dtw_h(nlen-1)) /dt
! single-tapered complex synthetic displacement and velocity
- wseis(:,ictaper) = 0.
- wseis(1:nlen,ictaper) = cmplx(dzr(1:nlen),0.)
- wseis2(:,ictaper) = 0.
- wseis2(1:nlen,ictaper) = cmplx(vzr(1:nlen),0.)
+ syn_dtw_ho_all(:,ictaper) = 0.
+ syn_vtw_ho_all(:,ictaper) = 0.
+ 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
- call fft(LNPT,wseis(:,ictaper),FORWARD_FFT,DT)
- call fft(LNPT,wseis2(:,ictaper),FORWARD_FFT,DT)
+ 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(:) + wseis(:,ictaper) * conjg(wseis(:,ictaper))
- v_bot_mtm(:) = v_bot_mtm(:) + wseis2(:,ictaper) * conjg(wseis2(:,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
@@ -769,11 +771,12 @@
qwc_adj(:) = cmplx(0.,0.)
do i = 1, i_right
- pwc_adj(i) = wseis2(i,ictaper) / v_bot_mtm(i)
- qwc_adj(i) = -wseis(i,ictaper) / d_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.)
@@ -798,22 +801,12 @@
tr_adj_src = 0.
if(iker /= 0) am_adj_src = 0.
- ishift = int(tshift_xc/dt)
-
- ! CHT: Compute the integrated waveform difference squared,
- ! no matter what measurement you make.
- ! To do this, we must UNSHIFT the data.
+ ! CHT: Compute the integrated waveform difference squared
waveform_temp1 = 0. ; waveform_temp2 = 0. ; waveform_temp3 = 0.
do i = 1,nlen
- i1 = istart + i
- if (i-ishift < 1 .or. i-ishift > nlen) then
- i2 = i
- else
- i2 = i-ishift
- endif
- waveform_temp1 = waveform_temp1 + ( dzr20(i2) * time_window(i) )**2
- waveform_temp2 = waveform_temp2 + ( dzr0(i) * time_window(i) )**2
- waveform_temp3 = waveform_temp3 + (( dzr20(i2) - dzr0(i) ) * time_window(i) )**2
+ 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
enddo
! NOTE: does not include DT factor or normalization by duration of window
waveform_d2 = 0.5 * waveform_temp1
@@ -826,21 +819,14 @@
i1 = istart + i
if(iker==0) then ! waveform
+ tr_adj_src(i1) = ( syn_dtw(i) - dat_dtw(i) ) * time_window(i)
- ! shift back the shifted data
- if (i-ishift < 1 .or. i-ishift > nlen) then
- i2 = i
- else
- i2 = i-ishift
- endif
- tr_adj_src(i1) = ( dzr0(i) - dzr20(i2) ) * time_window(i)
-
elseif(iker==1) then ! multitaper
tr_adj_src(i1) = fp(i) * time_window(i)
am_adj_src(i1) = fq(i) * time_window(i)
elseif(iker==2) then ! cross-correlation
- tr_adj_src(i1) = -(tshift_xc / sigma_dt_cc**2 ) * ft_bar_t(i) * time_window(i)
+ 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)
@@ -862,13 +848,13 @@
! misfit function value
if(iker==1) window_chi(1) = 0.5 * 2.0 * df * sum( (dtau_w(1:i_right))**2 * wp_taper(1:i_right) )
if(iker==1) window_chi(2) = 0.5 * 2.0 * df * sum( (dlnA_w(1:i_right))**2 * wq_taper(1:i_right) )
- window_chi(3) = 0.5 * (tshift_xc/sigma_dt_cc)**2
+ window_chi(3) = 0.5 * (tshift/sigma_dt_cc)**2
window_chi(4) = 0.5 * (dlnA/sigma_dlnA_cc)**2
! measurement (no uncertainty estimates)
if(iker==1) window_chi(5) = sum( dtau_w(1:i_right) * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
if(iker==1) window_chi(6) = sum( dlnA_w(1:i_right) * w_taper(1:i_right) ) / sum(w_taper(1:i_right) )
- window_chi(7) = tshift_xc
+ window_chi(7) = tshift
window_chi(8) = dlnA
! estimated mesurement uncertainties
@@ -942,20 +928,20 @@
! read sac file and convert to double precision
character(len=*),intent(in) :: datafile
- real, dimension(NDIM) :: data_sngl
+ 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,data_sngl,npt1,b1_sngl,dt1_sngl,NDIM,nerr)
+ call rsac1(datafile,dat_sngl,npt1,b1_sngl,dt1_sngl,NDIM,nerr)
if (nerr > 0) stop ' Error reading sac file'
! return double precision quantities
b1 = dble(b1_sngl)
dt1 = dble(dt1_sngl)
- data = dble(data_sngl)
+ data = dble(dat_sngl)
end subroutine drsac1
@@ -972,14 +958,14 @@
double precision, intent(in) :: b1,dt1
logical, parameter :: minmax_header = .true.
- real, dimension(npt1) :: data_sngl,ti_sngl
+ 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)
- data_sngl = real(data)
+ dat_sngl = real(data)
if (minmax_header) then
! get time vector
@@ -989,16 +975,16 @@
enddo
! set minmax values in sac file
- xmin_sngl = minval(data_sngl)
- xmax_sngl = maxval(data_sngl)
+ xmin_sngl = minval(dat_sngl)
+ xmax_sngl = maxval(dat_sngl)
call setfhv('depmin',xmin_sngl,nerr)
call setfhv('depmax',xmax_sngl,nerr)
! write file with headers
- call wsac0(datafile,ti_sngl,data_sngl,nerr)
+ call wsac0(datafile,ti_sngl,dat_sngl,nerr)
else
- call wsac1(datafile,data_sngl,npt1,b1_sngl,dt1_sngl,nerr)
+ call wsac1(datafile,dat_sngl,npt1,b1_sngl,dt1_sngl,nerr)
endif
if (nerr > 0) stop ' Error writing sac file'
@@ -1006,10 +992,14 @@
!-----------------------------------------------------------------------------
- subroutine mt_measure_select(nlen,tshift_xc,i_pmax_syn,dtau_w,dlnA_w,err_dt,err_dlnA,dt,i_left,i_right,fstart,fend,use_trace)
+ 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)
+ ! 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
+
integer, intent(in) :: nlen, i_pmax_syn
- double precision, intent(in) :: tshift_xc, dt
+ double precision, intent(in) :: tshift, dt
double precision, dimension(:), intent(inout) :: dtau_w, dlnA_w, err_dt, err_dlnA
double precision, intent(inout) :: fstart, fend
integer,intent(inout) :: i_left, i_right
@@ -1065,7 +1055,7 @@
!!$ do j = 1, i_right
!!$ if (stop_freq) exit
!!$ print *, j, dtau_w(j),stop_freq
-!!$ if (abs(dtau_w(j)) > 3 * abs(tshift_xc)) then
+!!$ if (abs(dtau_w(j)) > 3 * abs(tshift)) then
!!$ dtau_w(j) = 0
!!$ else if (j /= 1) then
!!$ stop_freq = .true.
@@ -1101,13 +1091,14 @@
fend = (i_right-1)*df
! if the cross-correlation time-shift is <= a time-step, set dtau(w) to zero
- if ( abs(tshift_xc) <= 1.01*dt ) then
+ ! NOTE: this should probably be a user parameter
+ if ( abs(tshift) <= 1.01*dt ) then
dtau_w(:) = 0.
use_trace = .false.
if (DISPLAY_DETAILS) then
print *, 'rejecting trace for too small a time shift:'
print *, ' dt = ', dt
- print *, ' tshift_xc = ', tshift_xc
+ print *, ' tshift = ', tshift
endif
endif
@@ -1115,15 +1106,15 @@
! CHT: dtau_w(j) --> abs(dtau_w(j)) for the first criterion
do j = i_left, i_right
if (use_trace .and. (abs(dtau_w(j)) > 1./(DT_FAC*fvec(j)) .or. err_dt(j) > 1./(ERR_FAC*fvec(j)) &
- .or. abs(dtau_w(j)) > DT_MAX_SCALE*abs(tshift_xc))) then
+ .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_xc), &
- abs(dtau_w(j)) > DT_MAX_SCALE*abs(tshift_xc)
+ 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'
endif
endif
@@ -1135,7 +1126,7 @@
! subroutines used in mtm_measure() and mtm_adj()
!==============================================================================
- subroutine interpolate_data_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)
double precision, dimension(NPT), intent(in) :: data, syn
double precision, dimension(NPT), intent(out) :: dat_win, syn_win
@@ -1158,11 +1149,11 @@
syn_win(i) = syn(ii) + (syn(ii+1)-syn(ii)) * (time-t1) /dt ! syn
enddo
- end subroutine interpolate_data_and_syn
+ end subroutine interpolate_dat_and_syn
!-----------------------------------------------------------------------------
- subroutine compute_cc(syn, data, nlen, dt, ishift, tshift_xc, dlnA, cc_max)
+ subroutine compute_cc(syn, data, nlen, dt, ishift, tshift, dlnA, cc_max)
! time shift MEASUREMENT between data (data) and synthetics (syn)
! CHT: modified the subroutine to resemble the one used in FLEXWIN
@@ -1170,7 +1161,7 @@
double precision, dimension(*), intent(in) :: syn, data
integer, intent(in) :: nlen
double precision, intent(in) :: dt
- double precision, intent(out) :: tshift_xc, dlnA, cc_max
+ double precision, intent(out) :: tshift, dlnA, cc_max
integer, intent(out) :: ishift
double precision :: cr_shift, cc, norm_s, norm
@@ -1198,7 +1189,7 @@
!!$ endif
!!$
!!$ enddo
-!!$ tshift_xc = ishift*dt
+!!$ tshift = ishift*dt
! initialise shift and cross correlation to zero
ishift = 0
@@ -1245,12 +1236,12 @@
endif
enddo
- tshift_xc = ishift*dt
+ tshift = ishift*dt
! amplitude MEASUREMENT (definition of Dahlen and Baig (2002), Eq. 3,17,18 : dlnA = Aobs/Asyn - 1)
! note that the records are UN-SHIFTED
- !dlnA = sqrt( (dt * sum( dat_displ(i1:i2) * dat_displ(i1:i2) )) &
- ! / (dt * sum( syn_displ(i1:i2) * syn_displ(i1:i2) )) ) - 1.
+ !dlnA = sqrt( (dt * sum( dat_dtw0(i1:i2) * dat_dtw0(i1:i2) )) &
+ ! / (dt * sum( syn_dtw(i1:i2) * syn_dtw(i1:i2) )) ) - 1.
! CHT revision 15-Oct-2008
! The previously used expression for dlnA is a first-order perturbation of ln(A1/A2) = (A1-A2)/A2
@@ -1262,93 +1253,44 @@
! ---------------------------------------------------------------------------
-!!$ subroutine compute_cc_error(dzr,dzr2,dzr20,nlen,dt,tshift_xc,dlnA,sigma_dt_cc,sigma_dlnA_cc)
-!!$
-!!$ ! CHT: estimate the uncertainty in the measurement
-!!$ ! NOTE: We implement the formula shown in Jeroen's Latex notes, only here
-!!$ ! the data are shifted, not the synthetics.
-!!$
-!!$ double precision, dimension(*), intent(in) :: dzr, dzr2, dzr20
-!!$ integer, intent(in) :: nlen
-!!$ double precision, intent(in) :: dt, tshift_xc, dlnA
-!!$ double precision, intent(inout) :: sigma_dt_cc, sigma_dlnA_cc
-!!$
-!!$ double precision, dimension(nlen) :: vzr
-!!$ double precision :: ntop
-!!$ integer i
-!!$
-!!$ ! compute synthetic velocity (unshifted)
-!!$ do i = 2, nlen-1
-!!$ vzr(i) = (dzr(i+1) - dzr(i-1)) / (2.0*dt)
-!!$ enddo
-!!$ vzr(1) = (dzr(2) - dzr(1)) / dt
-!!$ vzr(nlen) = (dzr(nlen) - dzr(nlen-1)) / dt
-!!$
-!!$ ! estimated uncertainty in cross-correlation travltime and amplitude
-!!$ ntop = sum( (dzr20(1:nlen) - (1.0 + dlnA)*dzr(1:nlen) )**2 )
-!!$ sigma_dt_cc = sqrt( ntop / sum( ((1.0 + dlnA)*vzr(1:nlen))**2 ) )
-!!$ sigma_dlnA_cc = sqrt( ntop / sum( dzr(1:nlen)**2 ) )
-!!$
-!!$ if(0==1) then
-!!$ print *, ' tshift : ', tshift_xc, ' +/- ', sigma_dt_cc
-!!$ print *, ' dlnA : ', dlnA, ' +/- ', sigma_dlnA_cc
-!!$
-!!$ ! time, syn-d, data-d, syn-v, shifted data, stretched syn
-!!$ open(88,file='tshift.dat')
-!!$ do i = 1,nlen
-!!$ write(88,'(6e18.6)') (i-1)*dt, dzr(i), dzr2(i), vzr(i), dzr20(i), (1.0 + dlnA)*dzr(i)
-!!$ enddo
-!!$ close(88)
-!!$ stop 'testing'
-!!$ endif
-!!$
-!!$ ! set uncertainty factors to 1 if you do not want to incorporate them
-!!$ ! into the adjoint sources and the misfit function values
-!!$ if (INCLUDE_ERROR == 0) then
-!!$ sigma_dt_cc = 1.0
-!!$ sigma_dlnA_cc = 1.0
-!!$ endif
-!!$
-!!$ end subroutine compute_cc_error
+ subroutine compute_average_error(data_dtw, syn_dtw_cc, syn_dtw_cc_dt, &
+ nlen, dt, sigma_dt, sigma_dlnA)
- ! ---------------------------------------------------------------------------
-
- subroutine compute_average_error(dat_d, syn_d_dt_dlnA, syn_d_dt, nlen, dt, sigma_dt, sigma_dlnA)
-
- ! CHT: estimate the uncertainty in the measurement (CC or MT),
+ ! CHT: Estimate the uncertainty in the CC measurement
! based on the integrated waveform difference between the data
! and the reconstructed synthetics.
! NOTE: We implement the exact equations that are in the Latex notes.
- double precision, dimension(*), intent(in) :: dat_d, syn_d_dt_dlnA, syn_d_dt
+ double precision, dimension(*), intent(in) :: data_dtw, syn_dtw_cc, syn_dtw_cc_dt
integer, intent(in) :: nlen
double precision, intent(in) :: dt
double precision, intent(inout) :: sigma_dt, sigma_dlnA
- double precision, dimension(nlen) :: syn_v_dt_dlnA
- double precision :: ntop
+ double precision, dimension(nlen) :: syn_vtw_cc
+ double precision :: sigma_dt_top, sigma_dlnA_top, sigma_dt_bot, sigma_dlnA_bot
integer i
- ! compute synthetic velocity (unshifted)
+ ! compute synthetic velocity (shifted and stretched)
do i = 2, nlen-1
- syn_v_dt_dlnA(i) = (syn_d_dt_dlnA(i+1) - syn_d_dt_dlnA(i-1)) / (2.0*dt)
+ syn_vtw_cc(i) = (syn_dtw_cc(i+1) - syn_dtw_cc(i-1)) / (2.0*dt)
enddo
- syn_v_dt_dlnA(1) = (syn_d_dt_dlnA(2) - syn_d_dt_dlnA(1)) / dt
- syn_v_dt_dlnA(nlen) = (syn_d_dt_dlnA(nlen) - syn_d_dt_dlnA(nlen-1)) / dt
+ 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
- ntop = sum( (dat_d(1:nlen) - syn_d_dt_dlnA(1:nlen) )**2 )
- sigma_dt = sqrt( ntop / sum( (syn_v_dt_dlnA(1:nlen))**2 ) )
- sigma_dlnA = sqrt( ntop / sum( (syn_d_dt(1:nlen))**2 ) )
+ 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_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 )
if(0==1) then
print *, ' sigma_dt : ', sigma_dt
print *, ' sigma_dlnA : ', sigma_dlnA
-
- ! time, syn-d, data-d, syn-v, shifted data, stretched syn
open(88,file='tshift.dat')
do i = 1,nlen
- write(88,'(5e18.6)') (i-1)*dt, dat_d(i), syn_d_dt_dlnA(i), syn_d_dt(i), syn_v_dt_dlnA(i)
+ 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'
@@ -1410,12 +1352,15 @@
! ---------------------------------------------------------------------------
- subroutine write_trans(filename,trans, wvec, fnum, i_right, idf_new, df, tshift, &
+ 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)
+ ! The transfer function maps the synthetics to the CC-deconstructed data;
+ ! the CC measurements then need to be applied to match the original data.
+
character(len=*), intent(in) :: filename
complex*16, intent(in) :: trans(:)
- double precision, intent(in) :: wvec(:), df, tshift
+ 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
@@ -1424,7 +1369,6 @@
double precision, dimension(NPT) :: fr
double precision :: f0, smth, smth1, smth2
-
abs_wt(:) = 0.
phi_wt(:) = 0.
@@ -1437,7 +1381,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))
@@ -1445,25 +1389,19 @@
if (mod(i,idf_new).eq.0 .and. OUTPUT_MEASUREMENT_FILES) then
write(10,*) fr(i), phi_wt(i)
write(20,*) fr(i), abs_wt(i)
- write(30,*) fr(i), abs_wt(i) - 1.
+ write(30,*) fr(i), log(abs_wt(i))
endif
enddo
+ ! NOTE: the CC measurements dT (tshift) and dlnA are BOTH included
dtau_wt(1) = tshift
do i = 1, i_right
- if (i > 1) dtau_wt(i) = (-1./wvec(i)) * phi_wt(i) + tshift ! NOTE: tshift is included
- dlnA_wt(i) = abs_wt(i) - 1.
- if(mod(i,idf_new).eq.0 .and. OUTPUT_MEASUREMENT_FILES) then
- write(40,*) fr(i), phi_wt(i)
- write(50,*) fr(i), dtau_wt(i)
- endif
-
if (i > 1 .and. i < i_right) then
! check the smoothness (2nd-order derivative) by 2*pi changes
- 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)
+ 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 (DISPLAY_DETAILS) print *, 'phase correction : 2 pi', sngl(fr(i)), sngl(phi_wt(i) - phi_wt(i+1))
do j = i+1, i_right
@@ -1478,6 +1416,17 @@
endif
endif
+ ! 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.
+
+ if(mod(i,idf_new).eq.0 .and. OUTPUT_MEASUREMENT_FILES) then
+ write(40,*) fr(i), phi_wt(i)
+ write(50,*) fr(i), dtau_wt(i)
+ endif
+
enddo
if (OUTPUT_MEASUREMENT_FILES) then
@@ -1494,118 +1443,141 @@
dlnA_wa = sum( dlnA_wt(1:i_right) ) / i_right
endif
- if (DISPLAY_DETAILS) then
- !print *, ' Taper traveltime measurement average : ', sngl(dtau_wa)
- !print *, ' Taper amplitude measurement average : ', sngl(dlnA_wa)
- !print *, ' i_right : ', i_right
- !f0 = 0.
- !call dwrite_ascfile_f(trim(filename)//'.dt_full',f0,df,i_right,dtau_wt(1:i_right))
- !call dwrite_ascfile_f(trim(filename)//'.dlnA_full',f0,df,i_right,dlnA_wt(1:i_right))
- !call dwrite_ascfile_f(trim(filename)//'.transfer_full',f0,df,i_right,abs(trans(1:i_right)))
- endif
+!!$ if (DISPLAY_DETAILS) then
+!!$ print *, ' Taper traveltime measurement average : ', sngl(dtau_wa)
+!!$ print *, ' Taper amplitude measurement average : ', sngl(dlnA_wa)
+!!$ print *, ' i_right : ', i_right
+!!$ !f0 = 0.
+!!$ !call dwrite_ascfile_f(trim(filename)//'.dt_full',f0,df,i_right,dtau_wt(1:i_right))
+!!$ !call dwrite_ascfile_f(trim(filename)//'.dlnA_full',f0,df,i_right,dlnA_wt(1:i_right))
+!!$ !call dwrite_ascfile_f(trim(filename)//'.transfer_full',f0,df,i_right,abs(trans(1:i_right)))
+!!$ endif
end subroutine write_trans
! --------------------------------------------------------------------
- subroutine reconstruct_syn_cc(datafile,filename,dzr,dzr2,tstart,dt,nlen,ishift,tshift_xc,dlnA, tseis_recon_cc, tseis_recon_cc_dt)
+ 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.
+
+ character(len=*), intent(in) :: filename
+ double precision, dimension(NPT), intent(in) :: dat_dtw
+ integer, intent(in) :: ishift, nlen
+ double precision, intent(in) :: tshift, dlnA, tstart, dt
+ double precision, dimension(NPT), intent(out) :: dat_dtw_cc
+ integer i
+
+ ! apply time shift (-dT) to OBSERVED seismogram
+ dat_dtw_cc(:) = dat_dtw(:)
+ do i = 1, nlen
+ if ((ishift+i) > 1 .and. (ishift+i) < nlen ) dat_dtw_cc(i) = dat_dtw(i+ishift)
+ enddo
+ ! fill the missing time window with the endpoint value
+ if (ishift < 0) dat_dtw_cc(1:-ishift+1) = dat_dtw_cc(-ishift+2)
+ if (ishift > 0) dat_dtw_cc(nlen-ishift:nlen) = dat_dtw_cc(nlen-ishift-1)
+
+ ! apply cross-correlation amplitude measurement (-DlnA) to OBSERVED seismogram
+ dat_dtw_cc(:) = dat_dtw_cc(:) * exp( -dlnA )
+
+ !if (DISPLAY_DETAILS) then
+ ! call dwrite_sacfile_f(datafile,'windowed_shifted_data.sac',tstart,nlen,dat_dtw0)
+ !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)
+
! reconstruct the synthetics using cross-correlation measurements:
- ! (1) apply dT to get tseis_recon_cc_dt
- ! (2) apply dT and dlnA to get tseis_recon_cc
+ ! (1) apply dT to get syn_dtw_cc_dt
+ ! (2) apply dT and dlnA to get syn_dtw_cc
- character(len=*), intent(in) :: datafile, filename
- double precision, dimension(NPT), intent(in) :: dzr, dzr2
+ double precision, dimension(NPT), intent(in) :: syn_dtw
integer, intent(in) :: ishift, nlen
- double precision, intent(in) :: tshift_xc, dlnA, tstart, dt
- double precision, dimension(NPT), intent(out) :: tseis_recon_cc, tseis_recon_cc_dt
+ double precision, intent(in) :: tshift, dlnA, tstart, dt
+ double precision, dimension(NPT), intent(out) :: syn_dtw_cc, syn_dtw_cc_dt
integer i
- ! shift synthetics by tshift_xc (in the main program, we shift the data instead)
- ! ishift = tshift_xc * dt
- tseis_recon_cc_dt(:) = dzr(:)
+ ! shift synthetics by tshift (in the main program, we shift the data instead)
+ ! ishift = tshift * dt
+ syn_dtw_cc_dt(:) = syn_dtw(:)
do i = 1, nlen
- if ((i-ishift) > 1 .and. (i-ishift) < nlen ) tseis_recon_cc_dt(i) = dzr(i-ishift)
+ if ((i-ishift) > 1 .and. (i-ishift) < nlen ) syn_dtw_cc_dt(i) = syn_dtw(i-ishift)
enddo
! fill the missing time window with the endpoint value
- if (ishift > 0) tseis_recon_cc_dt(1:ishift+1) = tseis_recon_cc_dt(ishift+2)
- if (ishift < 0) tseis_recon_cc_dt(nlen+ishift:nlen) = tseis_recon_cc_dt(nlen+ishift-1)
+ if (ishift > 0) syn_dtw_cc_dt(1:ishift+1) = syn_dtw_cc_dt(ishift+2)
+ if (ishift < 0) syn_dtw_cc_dt(nlen+ishift:nlen) = syn_dtw_cc_dt(nlen+ishift-1)
! apply cross-correlation amplitude measurement
- tseis_recon_cc(:) = 0.
- !tseis_recon_cc(:) = tseis_recon_cc_dt * (1. + dlnA) ! based on first-order approximation of dlnA
- tseis_recon_cc(:) = tseis_recon_cc_dt * exp( dlnA ) ! based on dlnA = ln(Aobs/Asyn)
+ syn_dtw_cc(:) = 0.
+ syn_dtw_cc(:) = syn_dtw_cc_dt * exp( dlnA ) ! based on dlnA = ln(Aobs/Asyn)
+ !syn_dtw_cc(:) = syn_dtw_cc_dt * (1. + dlnA) ! based on first-order approximation of dlnA
- if (OUTPUT_MEASUREMENT_FILES) then
- call dwsac1(trim(filename)//'.recon_syn_cc.sac',tseis_recon_cc,nlen,tstart,dt)
- call dwsac1(trim(filename)//'.recon_syn_cc_dt.sac',tseis_recon_cc_dt,nlen,tstart,dt)
-!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_cc.sac',tstart,nlen,tseis_recon_cc)
-!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_cc_dt.sac',tstart,nlen,tseis_recon_cc_dt)
-!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn',tstart,dt,nlen,tseis_recon_cc)
-!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn_dt',tstart,dt,nlen,tseis_recon_cc_dt)
- endif
-
end subroutine reconstruct_syn_cc
! -----------------------------------------------------------------------
- subroutine reconstruct_syn(datafile,filename, wseis_syn, wvec, dtau_wt, dlnA_wt, &
- i_right, tstart, dt, nlen, tseis_recon, tseis_recon_dt)
+ subroutine reconstruct_syn(filename, syn_dtwo, wvec, dtau_wt, dlnA_wt, &
+ i_right, tstart, dt, nlen, syn_dtw_mt, syn_dtw_mt_dt)
! reconstruct the synthetics using multitaper measurements:
- ! (1) apply dtau(w) to get tseis_recon_dt
- ! (2) apply dtau(w) and dlnA(w) to get tseis_recon
- ! --> d(w) = s(w) T(w) exp[-i w dT] (note that dtau_w includes tshift_xc)
+ ! (1) apply dtau(w) and dlnA(w) to get syn_dtw_mt0
- character(len=*), intent(in) :: datafile,filename
- complex*16, dimension(:), intent(in) :: wseis_syn
+ character(len=*), intent(in) :: filename
+ complex*16, dimension(:), intent(in) :: syn_dtwo
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) :: tseis_recon, tseis_recon_dt
+ double precision, dimension(:), intent(out) :: syn_dtw_mt, syn_dtw_mt_dt
complex*16, dimension(NPT) :: wseis_recon
integer i
double precision omega
! apply transfer function to synthetics (phase and amplitude)
- tseis_recon(:) = 0.
+ syn_dtw_mt(:) = 0.
wseis_recon(:) = cmplx(0.,0.)
do i = 1,i_right
omega = wvec(i)
- ! mathematically, these are identical; numerically, they produce nearly identical results
- wseis_recon(i) = wseis_syn(i) * (1.+ dlnA_wt(i)) * exp(-CCI*omega*dtau_wt(i))
- !synw_recon(i) = wseis_syn(i) * trans_mtm(i) * exp(-CCI*omega*tshift_xc)
+ 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)
enddo
- call fftinv(LNPT,wseis_recon,REVERSE_FFT,dt,tseis_recon)
+ call fftinv(LNPT,wseis_recon,REVERSE_FFT,dt,syn_dtw_mt)
! apply phase shifts only
- tseis_recon_dt(:) = 0.
+ syn_dtw_mt_dt(:) = 0.
wseis_recon(:) = cmplx(0.,0.)
do i = 1,i_right
omega = wvec(i)
- wseis_recon(i) = wseis_syn(i) * exp(-CCI*omega*dtau_wt(i))
+ wseis_recon(i) = syn_dtwo(i) * exp(-CCI*omega*dtau_wt(i))
enddo
- call fftinv(LNPT,wseis_recon,REVERSE_FFT,dt,tseis_recon_dt)
+ call fftinv(LNPT,wseis_recon,REVERSE_FFT,dt,syn_dtw_mt_dt)
if (OUTPUT_MEASUREMENT_FILES) then
- call dwsac1(trim(filename)//'.recon_syn.sac',tseis_recon,nlen,tstart,dt)
- call dwsac1(trim(filename)//'.recon_syn_dt.sac',tseis_recon_dt,nlen,tstart,dt)
-!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn.sac',tstart,nlen,tseis_recon)
-!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_dt.sac',tstart,nlen,tseis_recon_dt)
-!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn',tstart,dt,nlen,tseis_recon)
-!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn_dt',tstart,dt,nlen,tseis_recon_dt)
+ 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
! -----------------------------------------------------------------------
-!!$ subroutine check_recon_quality(filename,dzr20,dzr0,dzr2,tseis_recon,nlen,dt,tshift_xc,tshift_f1f2,cc_max_f1f2,cc_max)
+!!$ subroutine check_recon_quality(filename,dat_dtw_cc,syn_dtw,dat_dtw,syn_dtw_mt,nlen,dt,tshift,tshift_f1f2,cc_max_f1f2,cc_max)
!!$
!!$ character(len=*), intent(in) :: filename
-!!$ double precision, dimension(:), intent(in) :: dzr20, dzr0, dzr2, tseis_recon
-!!$ double precision, intent(in) :: dt, tshift_xc
+!!$ double precision, dimension(:), intent(in) :: dat_dtw_cc, syn_dtw, dat_dtw, syn_dtw_mt
+!!$ double precision, intent(in) :: dt, tshift
!!$ integer, intent(in) :: nlen
!!$ double precision, intent(out) :: tshift_f1f2, cc_max_f1f2, cc_max
!!$
@@ -1613,21 +1585,21 @@
!!$
!!$ ! Using Alessia's subroutine
!!$ ! First the shifted_obs_win vs the synthetic
-!!$ call f1f2_calc(dzr20,dzr0,nlen,1,nlen,dt, f1,f2,tshift_f1f2,cc_max_f1f2,dlnA_f1f2)
+!!$ call f1f2_calc(dat_dtw_cc,syn_dtw,nlen,1,nlen,dt, f1,f2,tshift_f1f2,cc_max_f1f2,dlnA_f1f2)
!!$
!!$ cc_max = cc_max_f1f2
!!$ if (OUTPUT_MEASUREMENT_FILES) then
!!$ open(10,file=trim(filename)//'.quality')
!!$ write(10,*) '<--------- F1 ------ F2 ---- tshift -- cc_max --- dlnA -->'
-!!$ write(10,"(a,5F10.5)") 'Before',sngl(F1),sngl(F2),sngl(tshift_xc),sngl(cc_max_f1f2),sngl(dlnA_f1f2)
+!!$ write(10,"(a,5F10.5)") 'Before',sngl(F1),sngl(F2),sngl(tshift),sngl(cc_max_f1f2),sngl(dlnA_f1f2)
!!$ endif
!!$ if (DISPLAY_DETAILS) then
!!$ write(*,*) '<--------- F1 ------ F2 ---- tshift -- cc_max --- dlnA -->'
-!!$ write(*,"(a,5F10.5)") 'Before',sngl(F1),sngl(F2),sngl(tshift_xc),sngl(cc_max_f1f2),sngl(dlnA_f1f2)
+!!$ write(*,"(a,5F10.5)") 'Before',sngl(F1),sngl(F2),sngl(tshift),sngl(cc_max_f1f2),sngl(dlnA_f1f2)
!!$ endif
!!$
!!$ ! Then the obs_win vs the reconstructed obs
-!!$ call f1f2_calc(dzr2,tseis_recon,nlen,1,nlen,dt, f1,f2,tshift_f1f2,cc_max_f1f2,dlnA_f1f2)
+!!$ call f1f2_calc(dat_dtw,syn_dtw_mt,nlen,1,nlen,dt, f1,f2,tshift_f1f2,cc_max_f1f2,dlnA_f1f2)
!!$
!!$ if (OUTPUT_MEASUREMENT_FILES) then
!!$ write(10,"(a,5F10.5)") 'After ',sngl(F1),sngl(F2),sngl(tshift_f1f2),sngl(cc_max_f1f2),sngl(dlnA_f1f2)
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2009-08-14 22:13:56 UTC (rev 15544)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2009-08-16 10:59:55 UTC (rev 15545)
@@ -83,9 +83,9 @@
}
#die("testing");
-#$imin = 1; $imax = $nevent;
+$imin = 1; $imax = $nevent;
#$imin = 100; $imax = $nevent;
-$imin = 183; $imax = $imin;
+#$imin = 183; $imax = $imin; # 183 Yorba Linda
# EXAMPLE lines from EIDs_simulation_duration
# 2 9967025 200.000 18200 0.29 4.084
@@ -159,26 +159,26 @@
print "--> run the measurement code and make adjoint sources\n";
print CSH "prepare_mt_measure_adj.pl $smodel 0 1 0 $Ts $eid\n";
-# # KEY: for T >=6, do both surface wave and body wave measurements
-# if ($Tmin <= 2 || $eid == 14095540) {
-# # special case for a particular event -- do cross-correlation only
-# print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
-# } else {
-# print CSH "run_mt_cc_plot.pl $smodel $ibp $Ts $lcut $tvec 0 $iparbools $par1 $par2 $par3\n";
-# }
+ #------------------------------------
+ # EXAMPLES for socal dataset
- # SOCAL: MT and CC, with both ajoint source plotted
+ # MT and CC, with both adjoint source plotted
#print CSH "run_mt_cc_plot.pl $smodel $ibp $Ts $lcut $tvec 0 $iparbools $par1 $par2 $par3\n";
- # SOCAL: CC only
+ # CC only
#print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
- # SOCAL: multitaper
+ # multitaper only
print CSH "run_mt_measure_adj.pl $smodel 1 1 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
- # SOCAL: coverage kernels for m16 (CC banana-doughnut adjoint sources)
+ # coverage kernels for m16 (CC banana-doughnut adjoint sources)
#print CSH "run_mt_measure_adj.pl $smodel 3 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+ # waveform difference
+ #print CSH "run_mt_measure_adj.pl $smodel 0 0 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+
+ #------------------------------------
+
# copy run files into RUN directory
$infile1 = "window_chi"; $outfile1 = "${dir_meas}/${ftag}_${infile1}";
$infile2 = "run_file"; $outfile2 = "${dir_meas}/${ftag}_${infile2}";
More information about the CIG-COMMITS
mailing list