[cig-commits] r15703 - seismo/3D/ADJOINT_TOMO/measure_adj
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Sep 28 13:57:41 PDT 2009
Author: danielpeter
Date: 2009-09-28 13:57:41 -0700 (Mon, 28 Sep 2009)
New Revision: 15703
Modified:
seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
Log:
bug fixed for waveform misfit measurements in mt_sub.f90
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-09-28 20:08:25 UTC (rev 15702)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-09-28 20:57:41 UTC (rev 15703)
@@ -79,6 +79,9 @@
stop 'Check tstart and tend'
endif
+ ! initializes i_right
+ i_right = 0
+
! LQY -- is this too small ???
wtr_mtm = 1.e-10
@@ -629,50 +632,50 @@
! ----------------------
! TAPERS
! ----------------------
+ if( iker == 1 ) then
+ ! frequency-domain tapers
+ ! 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
+ ! type of filter in the freq domain
+ !w_taper(i) = 1. ! boxcar
+ !w_taper(i) = 1. - (2.0/nw)**2 * ((i-1) - nw/2.0)**2 ! welch
+ w_taper(i) = 1. - cos(PI*(i-i_left)/(i_right-i_left))**ipwr_w ! cosine
+ enddo
- ! frequency-domain tapers
- ! 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
- ! type of filter in the freq domain
- !w_taper(i) = 1. ! boxcar
- !w_taper(i) = 1. - (2.0/nw)**2 * ((i-1) - nw/2.0)**2 ! welch
- w_taper(i) = 1. - cos(PI*(i-i_left)/(i_right-i_left))**ipwr_w ! cosine
- enddo
+ ! compute normalization factor for w_taper
+ ! note: 2 is needed for the integration from -inf to inf
+ df = 1. /(NPT*dt)
+ ffac = 2.0 * df * sum(w_taper(i_left:i_right) ) ! CHT: 1 --> i_left
+ if (DISPLAY_DETAILS) print *, 'Taper normalization factor, ffac = ', ffac
- ! compute normalization factor for w_taper
- ! note: 2 is needed for the integration from -inf to inf
- df = 1. /(NPT*dt)
- ffac = 2.0 * df * sum(w_taper(i_left:i_right) ) ! CHT: 1 --> i_left
- if (DISPLAY_DETAILS) print *, 'Taper normalization factor, ffac = ', ffac
+ ! wp_taper and wq_taper are modified frequency-domain tapers
+ ! Notice the option to include the frequency-dependent error.
+ wp_taper(:) = 0.
+ wq_taper(:) = 0.
+ dtau_wtr = wtr * sum(abs(dtau_w(i_left:i_right)))/(i_right-i_left) ! CHT i_left
+ dlnA_wtr = wtr * sum(abs(dlnA_w(i_left:i_right)))/(i_right-i_left) ! CHT i_left
- ! wp_taper and wq_taper are modified frequency-domain tapers
- ! Notice the option to include the frequency-dependent error.
- wp_taper(:) = 0.
- wq_taper(:) = 0.
- dtau_wtr = wtr * sum(abs(dtau_w(i_left:i_right)))/(i_right-i_left) ! CHT i_left
- dlnA_wtr = wtr * sum(abs(dlnA_w(i_left:i_right)))/(i_right-i_left) ! CHT i_left
+ do i = i_left, i_right ! CHT: 1 --> i_left
- do i = i_left, i_right ! CHT: 1 --> i_left
-
- if (INCLUDE_ERROR == 0) then ! no error estimate
+ if (INCLUDE_ERROR == 0) then ! no error estimate
wp_taper(i) = w_taper(i) / ffac
wq_taper(i) = w_taper(i) / ffac
- elseif (INCLUDE_ERROR == 1) then ! MT error estimate is assigned the CC error estimate
+ elseif (INCLUDE_ERROR == 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 (INCLUDE_ERROR == 2) then ! MT jack-knife error estimate
+ elseif (INCLUDE_ERROR == 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)
if (err_dlnA(i) < dlnA_wtr) err_A = err_A + dlnA_wtr
wp_taper(i) = w_taper(i) / ffac / (err_t ** 2)
wq_taper(i) = w_taper(i) / ffac / (err_A ** 2)
- endif
- enddo
+ endif
+ enddo
!!$ open(88,file='ftaper.dat')
!!$ do i = 1,i_right
@@ -680,6 +683,9 @@
!!$ enddo
!!$ close(88)
+ endif ! iker == 1
+
+
! post-processing time-domain taper
! NOTE: If the adjoint sources will be band-pass filtered at the end,
! then perhaps time_window is not necessary (i.e., use boxcar).
@@ -698,23 +704,26 @@
! ----------------------------------
! CROSS CORRELATION ADJOINT SOURCES
! ----------------------------------
+ if( iker == 2 .or. iker == 3 ) then
- ! compute synthetic velocity
- do i = 2, nlen-1
- syn_vtw(i) = (syn_dtw(i+1) - syn_dtw(i-1)) / (2.0*dt)
- enddo
- syn_vtw(1) = (syn_dtw(2) - syn_dtw(1)) / dt
- syn_vtw(nlen) = (syn_dtw(nlen) - syn_dtw(nlen-1)) / dt
+ ! compute synthetic velocity
+ do i = 2, nlen-1
+ syn_vtw(i) = (syn_dtw(i+1) - syn_dtw(i-1)) / (2.0*dt)
+ enddo
+ syn_vtw(1) = (syn_dtw(2) - syn_dtw(1)) / dt
+ syn_vtw(nlen) = (syn_dtw(nlen) - syn_dtw(nlen-1)) / dt
- ! compute CC traveltime and amplitude banana-dougnut kernels
- ft_bar_t = 0.
- Nnorm = dt * sum( syn_vtw(1:nlen) * syn_vtw(1:nlen) )
- ft_bar_t(1:nlen) = -syn_vtw(1:nlen) / Nnorm
+ ! 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(1:nlen) * syn_dtw(1:nlen) )
+ fa_bar_t(1:nlen) = syn_dtw(1:nlen) / Mnorm
+ endif
+
! -------------------------------
! MULTITAPER ADJOINT SOURCES
! -------------------------------
More information about the CIG-COMMITS
mailing list