[cig-commits] r21221 - seismo/3D/ADJOINT_TOMO/measure_adj
hejunzhu at geodynamics.org
hejunzhu at geodynamics.org
Thu Jan 10 15:03:54 PST 2013
Author: hejunzhu
Date: 2013-01-10 15:03:54 -0800 (Thu, 10 Jan 2013)
New Revision: 21221
Modified:
seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR
seismo/3D/ADJOINT_TOMO/measure_adj/ma_sub.f90
seismo/3D/ADJOINT_TOMO/measure_adj/measure_adj.f90
Log:
Add physical dispersion
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR 2013-01-10 21:29:11 UTC (rev 21220)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR 2013-01-10 23:03:54 UTC (rev 21221)
@@ -18,3 +18,4 @@
2.500 # ERR_FAC
3.500 # DT_MAX_SCALE
1.500 # NCYCLE_IN_WINDOW
+ .false. # USE_PHYSICAL_DISPERSION
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/ma_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/ma_sub.f90 2013-01-10 21:29:11 UTC (rev 21220)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/ma_sub.f90 2013-01-10 23:03:54 UTC (rev 21221)
@@ -9,8 +9,8 @@
contains
- subroutine mt_measure(datafile,filename,dat_dt,syn_dt,t0,dt,npts,tstart,tend, &
- istart,dat_dtw,syn_dtw,nlen,tshift,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,syn_dtw_cc, &
+ subroutine mt_measure(datafile,filename,dat_dt,syn_dt,syn_dt_phydisp,t0,dt,npts,tstart,tend, &
+ istart,dat_dtw,syn_dtw,syn_dtw_phydisp,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,syn_dtw_mt, &
err_dt,err_dlnA)
@@ -48,14 +48,14 @@
implicit none
character(len=150), intent(in) :: filename,datafile
- double precision, dimension(:), intent(in) :: dat_dt,syn_dt
+ double precision, dimension(:), intent(in) :: dat_dt,syn_dt, syn_dt_phydisp
double precision, intent(in) :: t0,dt,tstart,tend
integer, intent(in) :: npts ! rarely used
double precision, intent(out) :: tshift,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,sigma_dt,sigma_dlnA
integer, intent(out) :: istart,nlen,i_pmax_dat,i_pmax_syn,i_right
- double precision, dimension(:), intent(out) :: syn_dtw,dat_dtw,syn_dtw_cc,syn_dtw_mt,dtau_w,dlnA_w
+ double precision, dimension(:), intent(out) :: syn_dtw,dat_dtw,syn_dtw_phydisp, syn_dtw_cc,syn_dtw_mt,dtau_w,dlnA_w
complex*16, dimension(:), intent(out) :: trans_w
double precision, dimension(:), intent(out), optional :: err_dt,err_dlnA
@@ -99,7 +99,7 @@
! window and interpolate data and synthetics
!--------------------------------------------------------------------------
- call interpolate_dat_and_syn(dat_dt,syn_dt,tstart,tend,t0,dt,NPT,dat_dtw,syn_dtw,nlen,istart)
+ call interpolate_dat_and_syn(dat_dt,syn_dt,syn_dt_phydisp,tstart,tend,t0,dt,NPT,dat_dtw,syn_dtw,syn_dtw_phydisp,nlen,istart)
! some constants
sfac1 = (2./dble(nlen))**2 ! for Welch window
@@ -113,6 +113,9 @@
syn_dtw(i) = syn_dtw(i) * fac ! syn, windowed
dat_dtw(i) = dat_dtw(i) * fac ! dat, windowed
+ if (USE_PHYSICAL_DISPERSION) then
+ syn_dtw_phydisp(i)=syn_dtw_phydisp(i)*fac
+ end if
enddo
if (DISPLAY_DETAILS) then
@@ -556,13 +559,13 @@
! original coding by Carl Tape and finalized by Qinya Liu
! ======================================================================================================
- subroutine mt_adj(istart,dat_dtw,syn_dtw,nlen,dt,tshift,dlnA,sigma_dt_cc,sigma_dlnA_cc, &
+ subroutine mt_adj(istart,dat_dtw,syn_dtw,syn_dtw_phydisp,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, &
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, dimension(:), intent(in) :: dat_dtw,syn_dtw,syn_dtw_phydisp
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
@@ -624,20 +627,37 @@
if( (imeas >= 3).and.(imeas <= 6) ) then
! compute synthetic velocity
- do i = 2, nlen-1
- syn_vtw(i) = (syn_dtw(i+1) - syn_dtw(i-1)) / (2.0*dt)
- enddo
- syn_vtw(1) = (syn_dtw(2) - syn_dtw(1)) / dt
- syn_vtw(nlen) = (syn_dtw(nlen) - syn_dtw(nlen-1)) / dt
+ if (USE_PHYSICAL_DISPERSION) then
+ do i = 2, nlen-1
+ syn_vtw(i) = (syn_dtw_phydisp(i+1) - syn_dtw_phydisp(i-1)) / (2.0*dt)
+ enddo
+ syn_vtw(1) = (syn_dtw_phydisp(2) - syn_dtw_phydisp(1)) / dt
+ syn_vtw(nlen) = (syn_dtw_phydisp(nlen) - syn_dtw_phydisp(nlen-1)) / dt
- ! compute CC traveltime and amplitude banana-dougnut kernels
- ft_bar_t = 0.
- Nnorm = dt * sum( syn_vtw(1:nlen) * syn_vtw(1:nlen) )
- ft_bar_t(1:nlen) = -syn_vtw(1:nlen) / Nnorm
+ ! compute CC traveltime and amplitude banana-dougnut kernels
+ ft_bar_t = 0.
+ Nnorm = dt * sum( syn_vtw(1:nlen) * syn_vtw(1:nlen) )
+ ft_bar_t(1:nlen) = -syn_vtw(1:nlen) / Nnorm
- fa_bar_t = 0.
- Mnorm = dt * sum( syn_dtw(1:nlen) * syn_dtw(1:nlen) )
- fa_bar_t(1:nlen) = syn_dtw(1:nlen) / Mnorm
+ fa_bar_t = 0.
+ Mnorm = dt * sum( syn_dtw_phydisp(1:nlen) * syn_dtw_phydisp(1:nlen) )
+ fa_bar_t(1:nlen) = syn_dtw_phydisp(1:nlen) / Mnorm
+ else
+ do i = 2, nlen-1
+ syn_vtw(i) = (syn_dtw(i+1) - syn_dtw(i-1)) / (2.0*dt)
+ enddo
+ syn_vtw(1) = (syn_dtw(2) - syn_dtw(1)) / dt
+ syn_vtw(nlen) = (syn_dtw(nlen) - syn_dtw(nlen-1)) / dt
+
+ ! compute CC traveltime and amplitude banana-dougnut kernels
+ ft_bar_t = 0.
+ Nnorm = dt * sum( syn_vtw(1:nlen) * syn_vtw(1:nlen) )
+ ft_bar_t(1:nlen) = -syn_vtw(1:nlen) / Nnorm
+
+ fa_bar_t = 0.
+ Mnorm = dt * sum( syn_dtw(1:nlen) * syn_dtw(1:nlen) )
+ fa_bar_t(1:nlen) = syn_dtw(1:nlen) / Mnorm
+ end if
endif
! ----------------------------------------------
@@ -716,8 +736,11 @@
do ictaper = 1,ntaper
! tapered synthetic displacement
- syn_dtw_h(1:nlen) = syn_dtw(1:nlen) * tas(1:nlen,ictaper)
-
+ if (USE_PHYSICAL_DISPERSION) then
+ syn_dtw_h(1:nlen) = syn_dtw_phydisp(1:nlen) * tas(1:nlen,ictaper)
+ else
+ syn_dtw_h(1:nlen) = syn_dtw(1:nlen) * tas(1:nlen,ictaper)
+ end if
! compute velocity of tapered syn
do i = 2, nlen-1
syn_vtw_h(i) = (syn_dtw_h(i+1) - syn_dtw_h(i-1)) / (2.0*dt)
@@ -1205,19 +1228,19 @@
! subroutines used in mtm_measure() and mtm_adj()
!==============================================================================
- subroutine interpolate_dat_and_syn(data, syn, tstart, tend, t0, dt, NPT, &
- dat_win, syn_win, nlen, istart)
+ subroutine interpolate_dat_and_syn(data, syn, syn_phydisp, tstart, tend, t0, dt, NPT, &
+ dat_win, syn_win, syn_win_phydisp, nlen, istart)
! extract windowed portion from original data/syn traces
implicit none
- double precision, dimension(NPT), intent(in) :: data, syn
+ double precision, dimension(NPT), intent(in) :: data, syn, syn_phydisp
double precision, intent(in) :: tstart
double precision, intent(in) :: tend, t0, dt
integer, intent(in) :: NPT
- double precision, dimension(NPT), intent(out) :: dat_win, syn_win
+ double precision, dimension(NPT), intent(out) :: dat_win, syn_win, syn_win_phydisp
integer, intent(out) :: nlen, istart
integer :: ii, i
@@ -1245,8 +1268,10 @@
t1 = floor((time-t0)/dt) * dt + t0
dat_win(i) = data(ii) + (data(ii+1)-data(ii)) * (time-t1) / dt
- syn_win(i) = syn(ii) + (syn(ii+1)-syn(ii)) * (time-t1) /dt
-
+ syn_win(i) = syn(ii) + (syn(ii+1)-syn(ii)) * (time-t1) /dt
+ if (USE_PHYSICAL_DISPERSION) then
+ syn_win_phydisp(i) = syn_phydisp(ii) + ( syn_phydisp(ii+1) - syn_phydisp(ii) ) * (time-t1)/dt
+ end if
enddo
end subroutine interpolate_dat_and_syn
@@ -1828,6 +1853,7 @@
read(10,*) ERR_FAC
read(10,*) DT_MAX_SCALE
read(10,*) NCYCLE_IN_WINDOW
+ read(10,*) USE_PHYSICAL_DISPERSION
close(10)
imeas = imeas0
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/measure_adj.f90 2013-01-10 21:29:11 UTC (rev 21220)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/measure_adj.f90 2013-01-10 23:03:54 UTC (rev 21221)
@@ -30,10 +30,11 @@
implicit none
- character(len=150) :: datafile,synfile,file_prefix,file_prefix0,file_prefix2,measure_file_prefix,adj_file_prefix
- integer :: num_meas, j, ios, npt1, npt2, npts, nn
- double precision, dimension(NDIM) :: data, syn, adj_syn_all, tr_adj_src, am_adj_src, recon_cc_all, syn_dtw_cc, syn_dtw_mt
- double precision :: t01, dt1, t02, dt2, t0, dt, tstart, tend, tt, dtt, df
+ character(len=150) :: datafile,synfile,synfile_phydisp,file_prefix,file_prefix0,file_prefix2,measure_file_prefix,adj_file_prefix
+ integer :: num_meas, j, ios, npt1, npt2,npt3, npts, nn
+ double precision, dimension(NDIM) :: data, syn, syn_phydisp, adj_syn_all, &
+ tr_adj_src, am_adj_src, recon_cc_all, syn_dtw_cc, syn_dtw_mt
+ double precision :: t01, dt1, t02, dt2, t03, dt3, t0, dt, tstart, tend, tt, dtt, df
double precision, dimension(NCHI) :: window_chi
double precision :: fend0, fstart0, fend, fstart
@@ -44,7 +45,7 @@
double precision :: tshift, sigma_dt_cc, dlnA, sigma_dlnA_cc, sigma_dt, sigma_dlnA
double precision :: all_chi, tr_chi, am_chi, cc_max, T_pmax_dat, T_pmax_syn
!double precision :: tshift_f1f2, cc_max_f1f2
- double precision, dimension(NPT) :: dtau_w, dlnA_w, err_dt, err_dlnA, syn_dtw, data_dtw
+ double precision, dimension(NPT) :: dtau_w, dlnA_w, err_dt, err_dlnA, syn_dtw, data_dtw,syn_dtw_phydisp
complex*16, dimension(NPT) :: trans_mtm
integer :: nlen, i_left, i_pmax_dat, i_pmax_syn, i_right, i_right0, istart, &
ipair, npairs, nwin, itmax
@@ -89,12 +90,19 @@
if (ios /= 0) stop 'Error reading windows file'
read(11,'(a)',iostat=ios) synfile
if (ios /= 0) stop 'Error reading windows file'
+ if (USE_PHYSICAL_DISPERSION) then
+ synfile_phydisp=trim(synfile)//'.phydisp'
+ end if
+
! read data and syn (in double precision)
! LQY: data is read last to be stored in memory for header retrieval later
call drsac1(datafile,data,npt1,t01,dt1)
call drsac1(synfile,syn,npt2,t02,dt2)
+ if (USE_PHYSICAL_DISPERSION) then
+ call drsac1(synfile_phydisp,syn_phydisp,npt3,t03,dt3)
+ end if
if (DISPLAY_DETAILS) then
print *
@@ -128,6 +136,9 @@
if(RUN_BANDPASS) then
call bandpass(data,npts,dt,fstart0,fend0)
call bandpass(syn,npts,dt,fstart0,fend0)
+ if (USE_PHYSICAL_DISPERSION) then
+ call bandpass(syn_phydisp,npts,dt,fstart0,fend0)
+ end if
endif
! find out station/network/comp names,etc from synthetics
@@ -197,8 +208,8 @@
! make measurements
! also compute reconstructed synthetics for CC (and MT, if specified) measurements
- call mt_measure(datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,&
- istart,data_dtw,syn_dtw,nlen,tshift,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,syn_dtw_cc,&
+ call mt_measure(datafile,measure_file_prefix,data,syn,syn_phydisp,t0,dt,npts,tstart,tend,&
+ istart,data_dtw,syn_dtw,syn_dtw_phydisp,nlen,tshift,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,syn_dtw_cc,&
i_pmax_dat,i_pmax_syn,i_right0,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,syn_dtw_mt,err_dt,err_dlnA)
i_right = i_right0
i_left = 1 ! LQY: is it feasible that i_left is not 1? mt_adj() inherently assumes it.
@@ -223,8 +234,8 @@
print *, ' reverting from MT measurement to CC measurement...'
imeas = imeas0 - 2
is_mtm = 3 ! LQY: WHY not is_mtm = 2?
- call mt_measure(datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,&
- istart,data_dtw,syn_dtw,nlen,&
+ call mt_measure(datafile,measure_file_prefix,data,syn,syn_phydisp,t0,dt,npts,tstart,tend,&
+ istart,data_dtw,syn_dtw,syn_dtw_phydisp,nlen,&
tshift,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,syn_dtw_cc,&
i_pmax_dat,i_pmax_syn,i_right,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,syn_dtw_mt)
else
@@ -260,7 +271,7 @@
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,&
+ call mt_adj(istart,data_dtw,syn_dtw,syn_dtw_phydisp,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,&
window_chi,tr_adj_src,tr_chi,am_adj_src,am_chi)
More information about the CIG-COMMITS
mailing list