[cig-commits] r15004 - seismo/3D/ADJOINT_TOMO/measure_adj
carltape at geodynamics.org
carltape at geodynamics.org
Sun May 17 17:26:41 PDT 2009
Author: carltape
Date: 2009-05-17 17:26:40 -0700 (Sun, 17 May 2009)
New Revision: 15004
Modified:
seismo/3D/ADJOINT_TOMO/measure_adj/Makefile
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
Log:
In preparation for switching from ifort to gfortran, I removed the f1f2 subroutines, which are not available using gfortran. f1f2 was not essential for the measurement package. The present version of mt_measure_adj.f90 compiles using Makefile_ifort and the ifort compiler at Caltech, but not using 64-bit Linux Ubuntu with gfortran at Harvard. I modified cross-correlation measurement subroutine to match the one used in FLEXWIN.
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/Makefile
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/Makefile 2009-05-17 18:31:16 UTC (rev 15003)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/Makefile 2009-05-18 00:26:40 UTC (rev 15004)
@@ -1,8 +1,9 @@
F90 = gfortran
-F90_FLAGS = -O2 -ffixed-line-length-132
+#F90_FLAGS = -O2 -ffixed-line-length-132
+F90_FLAGS = -O2
-LIB = -L/opt/seismo/lib -lDRWFiles -lf90recipes -lDSacio -lDSacLib -lSacTools -lm
+LIB = -L/opt/seismo/lib -lDRWFiles -lf90recipes_gfortran -lDSacio -lDSacLib -lSacTools -lm
MOD = mt_constants mt_variables mt_sub2 mt_sub
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-05-17 18:31:16 UTC (rev 15003)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-05-18 00:26:40 UTC (rev 15004)
@@ -27,7 +27,8 @@
double precision :: BEFORE_DLNA ! temporary: this should be an input parameter
double precision :: tshift_xc, sigma_dt_cc, dlnA, sigma_dlnA_cc, sigma_dt, sigma_dlnA
- double precision :: tr_chi, am_chi, tshift_f1f2, cc_max_f1f2, cc_max, T_pmax_dat, T_pmax_syn
+ 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
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
@@ -250,7 +251,7 @@
! 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,tshift_f1f2,cc_max_f1f2,tseis_recon_cc,&
+ 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
@@ -267,17 +268,17 @@
print *, 'fstart0/fend0 :', fstart0, fend0
print *, 'fstart/fend :', fstart, fend
print *, 'Tpmax :', T_pmax_dat, T_pmax_syn
- if (use_trace .and. (cc_max < BEFORE_QUALITY .or. cc_max_f1f2 < AFTER_QUALITY .or. abs(tshift_xc) > BEFORE_TSHIFT .or. abs(tshift_f1f2) > AFTER_TSHIFT)) then
- use_trace = .false.
- if (DISPLAY_DETAILS) then
- print *, 'Fail if ANY of these is true :'
- print *, ' cc_max : ', cc_max, BEFORE_QUALITY, cc_max < BEFORE_QUALITY
- print *, ' cc_max_f1f2 : ', cc_max_f1f2, AFTER_QUALITY, cc_max_f1f2 < AFTER_QUALITY
- print *, ' tshift_xc : ', tshift_xc, BEFORE_TSHIFT, abs(tshift_xc) > BEFORE_TSHIFT
- print *, ' tshift_f1f2 : ', tshift_f1f2, AFTER_TSHIFT, abs(tshift_f1f2) > AFTER_TSHIFT
- endif
- endif
- print *, ' Status: ', trim(measure_file_prefix), use_trace
+!!$ if (use_trace .and. (cc_max < BEFORE_QUALITY .or. cc_max_f1f2 < AFTER_QUALITY .or. abs(tshift_xc) > BEFORE_TSHIFT .or. abs(tshift_f1f2) > AFTER_TSHIFT)) then
+!!$ use_trace = .false.
+!!$ if (DISPLAY_DETAILS) then
+!!$ print *, 'Fail if ANY of these is true :'
+!!$ print *, ' cc_max : ', cc_max, BEFORE_QUALITY, cc_max < BEFORE_QUALITY
+!!$ print *, ' cc_max_f1f2 : ', cc_max_f1f2, AFTER_QUALITY, cc_max_f1f2 < AFTER_QUALITY
+!!$ print *, ' tshift_xc : ', tshift_xc, BEFORE_TSHIFT, abs(tshift_xc) > BEFORE_TSHIFT
+!!$ print *, ' tshift_f1f2 : ', tshift_f1f2, AFTER_TSHIFT, abs(tshift_f1f2) > AFTER_TSHIFT
+!!$ endif
+!!$ endif
+!!$ print *, ' Status: ', trim(measure_file_prefix), use_trace
! CHT: If MT measurement window is rejected by mt_measure_select, then use a CC measurement.
! --> May also want to consider quality checks for CC as well.
@@ -287,7 +288,7 @@
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,&
- tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tshift_f1f2,cc_max_f1f2,tseis_recon_cc,&
+ 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.
endif
@@ -299,7 +300,7 @@
! 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 BEFORE_TSHIFT in subroutine compute_time_shift in mt_sub.f90.
+ ! See BEFORE_TSHIFT in subroutine compute_cc in mt_sub.f90.
! NOTE: BEFORE_DLNA should be included in the PARFILE
if ( abs(tshift_xc) > BEFORE_TSHIFT ) then
if (DISPLAY_DETAILS) then
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-05-17 18:31:16 UTC (rev 15003)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-05-18 00:26:40 UTC (rev 15004)
@@ -29,8 +29,8 @@
!
! original coding in Fortran 77 by Ying Zhou
! upgraded to Fortran 90 by Alessia Maggi
- ! minor modifications by Vala Hjorleifsdottir and Carl Tape
- ! finalized by Qinya Liu
+ ! 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
@@ -40,7 +40,7 @@
! =================================================================================================
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,tshift_f1f2,cc_max_f1f2,tseis_recon_cc, &
+ istart,dzr20,dzr0,nlen,tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tseis_recon_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
@@ -48,7 +48,7 @@
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, tshift_f1f2,cc_max_f1f2, sigma_dt,sigma_dlnA
+ double precision, intent(out) :: tshift_xc, 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
@@ -119,8 +119,9 @@
! cross-correlation traveltime and amplitude measurements
!------------------------------------------------------------------
+ ! compute cross-correltaion time shift and also amplitude measurmement
! LQY: Ying suggested to align them at relatively long periods
- call compute_time_shift(dzr, dzr2, nlen, dt, ishift, tshift_xc)
+ call compute_cc(dzr, dzr2, nlen, dt, ishift, tshift_xc, dlnA, cc_max)
! apply time shift to OBSERVED seismogram
dzr0(:) = dzr(:)
@@ -155,17 +156,6 @@
syn_accel(1) = (syn_veloc(2) - syn_veloc(1)) / dt
syn_accel(nlen) = (syn_veloc(nlen) - syn_veloc(nlen-1)) / 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(1:nlen) * dat_displ(1:nlen) )) &
- ! / (dt * sum( syn_displ(1:nlen) * syn_displ(1:nlen) )) ) - 1.
-
- ! CHT revision 15-Oct-2008
- ! The previously used expression for dlnA is a first-order perturbation of ln(A1/A2) = (A1-A2)/A2
- ! The new expression is better suited to getting values between -1 and 1,
- ! with dlnA = 0 indicating perfect fit, as before.
- dlnA = 0.5 * log( sum(dat_displ(1:nlen) * dat_displ(1:nlen)) / sum(syn_displ(1:nlen) * syn_displ(1:nlen)) )
-
! 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)
@@ -340,7 +330,7 @@
! 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 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_average_meas(filename, IKER_MT, dtau_wa, dlnA_wa, sigma_dt, sigma_dlnA)
endif
@@ -378,7 +368,7 @@
call reconstruct_syn(datafile,filename,wseis_syn,wvec,dtau_mtm,dlnA_mtm,i_right,tstart,dt,nlen,tseis_recon,tseis_recon_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,dzr20,dzr0,dzr3,tseis_recon,nlen,dt,tshift_xc, 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.
@@ -1051,48 +1041,104 @@
!-----------------------------------------------------------------------------
- subroutine compute_time_shift(dzr, dzr2, nlen, dt, ishift, tshift_xc)
+ subroutine compute_cc(syn, data, nlen, dt, ishift, tshift_xc, dlnA, cc_max)
- ! time shift MEASUREMENT
- ! If the input are windowed, then the information outside the window is not used.
- ! WARNING: If the data are NOT windowed, then this will use information that is
- ! a full window's width outside the window of interest.
+ ! time shift MEASUREMENT between data (data) and synthetics (syn)
+ ! CHT: modified the subroutine to resemble the one used in FLEXWIN
- double precision, dimension(*), intent(in) :: dzr, dzr2
+ double precision, dimension(*), intent(in) :: syn, data
integer, intent(in) :: nlen
double precision, intent(in) :: dt
- double precision, intent(out) :: tshift_xc
+ double precision, intent(out) :: tshift_xc, dlnA, cc_max
integer, intent(out) :: ishift
- double precision cr_shift,cc,cc_max
- integer n_left, n_right, i, j
+ double precision :: cr_shift, cc, norm_s, norm
+ integer i1, i2, i, j, i_left, i_right, id_left, id_right
- ! these choices will slide the entire windowed record past the other
- cr_shift = nlen*dt
- n_left = ceiling( -1.0 * cr_shift / dt )
- n_right = floor( cr_shift / dt )
+!!$ ! these choices will slide the entire windowed record past the other
+!!$ cr_shift = nlen*dt
+!!$ i_left = ceiling( -1.0 * cr_shift / dt )
+!!$ i_right = floor( cr_shift / dt )
+!!$
+!!$ ! cross-correlation
+!!$ 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)
+!!$ enddo
+!!$
+!!$ !if(cc > cc_max) then
+!!$ ! CHT, 07-Sept-2008: Do not allow time shifts larger than the specified input
+!!$ if(cc > cc_max .and. abs(i*dt) <= BEFORE_TSHIFT ) then
+!!$ cc_max = cc
+!!$ ishift = i
+!!$ endif
+!!$
+!!$ enddo
+!!$ tshift_xc = ishift*dt
- ! cross-correlation
- ishift = 0
- do i = n_left, n_right, 1
+ ! initialise shift and cross correlation to zero
+ ishift = 0
+ cc_max = 0.0
- cc = 0.
- do j = 1, nlen
- if((j+i) > 1 .and. (j+i) < nlen) cc = cc + dzr(j) * dzr2(j+i)
- enddo
+ ! index of window limits
+ i1 = 1
+ i2 = nlen
- !if(cc > cc_max) then
- ! CHT, 07-Sept-2008: Do not allow time shifts larger than the specified input
- if(cc > cc_max .and. abs(i*dt) <= BEFORE_TSHIFT ) then
- cc_max = cc
- ishift = i
- endif
+ ! length of window (number of points, including ends)
+ !nlen = i2 - i1 + 1
+ ! power of synthetic signal in window
+ norm_s = sqrt(sum(syn(i1:i2)*syn(i1:i2)))
+
+ ! left and right limits of index (time) shift search
+ ! NOTE: This looks OUTSIDE the time window of interest to compute TSHIFT and CC.
+ ! How far to look outside, in theory, should be another parameter.
+ ! However, it does not matter as much if the data and synthetics are
+ ! zeroed outside the windows.
+ i_left = -1*int(nlen/2.0)
+ i_right = int(nlen/2.0)
+
+ ! i is the index to shift to be applied to DATA (data)
+ do i = i_left, i_right
+
+ ! normalization factor varies as you take different windows of data
+ id_left = max(1,i1+i) ! left index for data window
+ id_right = min(nlen,i2+i) ! right index for data window
+ norm = norm_s * sqrt(sum(data(id_left:id_right)*(data(id_left:id_right))))
+
+ ! cc as a function of i
+ cc = 0.
+ do j = i1, i2 ! loop over full window length
+ if((j+i).ge.1 .and. (j+i).le.nlen) cc = cc + syn(j)*data(j+i) ! d is shifted by i
+ enddo
+ cc = cc/norm
+
+ !if(cc > cc_max) then
+ ! CHT, 07-Sept-2008: Do not allow time shifts larger than the specified input
+ if(cc .gt. cc_max .and. abs(i*dt) <= BEFORE_TSHIFT ) then
+ cc_max = cc
+ ishift = i
+ endif
+
enddo
- tshift_xc = ishift*dt
+ tshift_xc = ishift*dt
- end subroutine compute_time_shift
+ ! 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.
+ ! CHT revision 15-Oct-2008
+ ! The previously used expression for dlnA is a first-order perturbation of ln(A1/A2) = (A1-A2)/A2
+ ! The new expression is better suited to getting 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
+
! ---------------------------------------------------------------------------
!!$ subroutine compute_cc_error(dzr,dzr2,dzr20,nlen,dt,tshift_xc,dlnA,sigma_dt_cc,sigma_dlnA_cc)
@@ -1428,43 +1474,43 @@
! -----------------------------------------------------------------------
- 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,dzr20,dzr0,dzr2,tseis_recon,nlen,dt,tshift_xc,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
+!!$ integer, intent(in) :: nlen
+!!$ double precision, intent(out) :: tshift_f1f2, cc_max_f1f2, cc_max
+!!$
+!!$ double precision :: f1,f2, dlnA_f1f2
+!!$
+!!$ ! 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)
+!!$
+!!$ 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)
+!!$ 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)
+!!$ 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)
+!!$
+!!$ 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)
+!!$ close(10)
+!!$ endif
+!!$
+!!$ if (DISPLAY_DETAILS) write(*,"(a,5F10.5)") 'After ',sngl(F1),sngl(F2),sngl(tshift_f1f2),sngl(cc_max_f1f2),sngl(dlnA_f1f2)
+!!$
+!!$ end subroutine check_recon_quality
- character(len=*), intent(in) :: filename
- double precision, dimension(:), intent(in) :: dzr20, dzr0, dzr2, tseis_recon
- double precision, intent(in) :: dt, tshift_xc
- integer, intent(in) :: nlen
- double precision, intent(out) :: tshift_f1f2, cc_max_f1f2, cc_max
-
- double precision :: f1,f2, dlnA_f1f2
-
- ! 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)
-
- 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)
- 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)
- 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)
-
- 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)
- close(10)
- endif
-
- if (DISPLAY_DETAILS) write(*,"(a,5F10.5)") 'After ',sngl(F1),sngl(F2),sngl(tshift_f1f2),sngl(cc_max_f1f2),sngl(dlnA_f1f2)
-
- end subroutine check_recon_quality
-
!-------------------------------------------------------------------
subroutine interpolate_syn(syn,t1,dt1,npt1,t2,dt2,npt2)
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub2.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub2.f90 2009-05-17 18:31:16 UTC (rev 15003)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub2.f90 2009-05-18 00:26:40 UTC (rev 15004)
@@ -410,199 +410,199 @@
!
! -----------------------------------------------------------------
- subroutine F1F2_calc(d,s,npts,i1,i2,dt,F1,F2,tshift,cc_max,dlnA)
+!!$ subroutine F1F2_calc(d,s,npts,i1,i2,dt,F1,F2,tshift,cc_max,dlnA)
+!!$
+!!$
+!!$ double precision, dimension(*), intent(in) :: d, s
+!!$ integer, intent(in) :: npts,i1,i2
+!!$ double precision, intent (in) :: dt
+!!$ double precision, intent(out) :: F1,F2,tshift,cc_max,dlnA
+!!$
+!!$ double precision, dimension(:), allocatable :: s_cor,d_loc
+!!$
+!!$ double precision :: cr_shift, cc
+!!$ integer :: n_left,n_right,ishift,npts_win, i, j
+!!$
+!!$ real ax,bx,cx,fa,fb,fc,f1_min,f2_min,f1_top,f1_bot,A1,A2
+!!$ real golden !f1,f2
+!!$
+!!$ npts_win=i2-i1+1
+!!$
+!!$! allocate memory for s_cor (the corrected synthetic)
+!!$ allocate(s_cor(npts_win))
+!!$ allocate(d_loc(npts_win))
+!!$
+!!$ d_loc(1:npts_win)=d(i1:i2)
+!!$
+!!$! do cross-correlation:
+!!$ call xcorr_calc(d,s,npts,i1,i2,ishift,cc_max)
+!!$! n_left = int((-1.0) * cr_shift / dt)
+!!$! n_right = int(cr_shift / dt)
+!!$! ishift=0
+!!$! cc_max=0.
+!!$! do i = n_left, n_right
+!!$! cc = 0
+!!$! do j = 1, npts
+!!$! if((j+i).gt.1.and.(j+i).lt.npts) cc = cc + s(j) * d(j+i)
+!!$! enddo
+!!$! if( cc .gt. cc_max) then
+!!$! cc_max = cc
+!!$! ishift = i
+!!$! endif
+!!$! enddo
+!!$ tshift=ishift*dt
+!!$
+!!$! apply time shift to synthetic seismogram
+!!$! write(*,*)'shift synth seismogram by ', tshift, 'seconds'
+!!$ do i = 1, npts_win
+!!$ s_cor(i) = 0
+!!$ if( (i1-1+i-ishift) .gt. 1 .and. (i1-1+i-ishift) .lt.npts ) s_cor(i) = s(i1-1+i-ishift)
+!!$ enddo
+!!$
+!!$! DEBUG: output
+!!$! open(unit=11, file='DEBUG_calcF1F2.dat')
+!!$! do i = 1, npts_win
+!!$! write(11,'(4(e12.4,1x))') b+(i-1)*dt, s_cor(i), s(i1-1+i), d(i1-1+i)
+!!$! enddo
+!!$! close(11)
+!!$
+!!$! calculate dlnA
+!!$ dlnA = sqrt( ( sum( d(i1:i2) * d(i1:i2) )) / (sum( s_cor(1:npts_win) * s_cor(1:npts_win) )) ) - 1
+!!$
+!!$
+!!$! calculate F1, the least squares misfit
+!!$ f1_top=0.0
+!!$ f1_bot=0.0
+!!$ do i = 1,npts_win
+!!$ f1_top=f1_top+(sngl(d_loc(i))-sngl(s_cor(i)))**2
+!!$! f1_bot=f1_bot+sqrt(sngl(d_loc(i))**2*sngl(s_cor(i))**2)
+!!$ f1_bot=f1_bot+sngl(d_loc(i))**2
+!!$ enddo
+!!$ if ( f1_bot .gt. 0.0 ) then
+!!$ F1 = dble(f1_top / f1_bot)
+!!$ else
+!!$ write(*,*) 'Sum d(t)**2 = 0 : empty observed seismogram.'
+!!$ F1=0
+!!$ F2=0
+!!$ return
+!!$ endif
+!!$
+!!$! do fa1 minimization to find A1
+!!$ ax=1e-3
+!!$ bx=1e3
+!!$ call mnbrak(ax,bx,cx,fa,fb,fc,fa1)
+!!$ f1_min=golden(ax,bx,cx,fa1,sngl(tol),A1)
+!!$
+!!$! do fa2 minimization to find A2
+!!$ ax=1e-3
+!!$ bx=1e3
+!!$ call mnbrak(ax,bx,cx,fa,fb,fc,fa2)
+!!$ f2_min=golden(ax,bx,cx,fa2,sngl(TOL),A2)
+!!$
+!!$! calculate F2
+!!$ F2=dble(min(A1,A2)/max(A1,A2))
+!!$
+!!$! Turn F1 around
+!!$ F1=1-F1
+!!$
+!!$ deallocate(s_cor)
+!!$ deallocate(d_loc)
+!!$
+!!$ contains
+!!$
+!!$! -----------------------------------------------------------------
+!!$
+!!$ real function fa1(a1)
+!!$ real a1
+!!$
+!!$ if (abs(a1).lt.TOL) then
+!!$ write(*,*) 'value of a1 close to zero : ', a1
+!!$ stop
+!!$ endif
+!!$
+!!$ fa1=0.0
+!!$ do i = 1,npts_win
+!!$ fa1=fa1+(sngl(d_loc(i))-a1*sngl(s_cor(i)))**2
+!!$ enddo
+!!$
+!!$ end function
+!!$
+!!$! -----------------------------------------------------------------
+!!$
+!!$ real function fa2(a2)
+!!$ real a2
+!!$
+!!$ if (abs(a2).lt.TOL) then
+!!$ write(*,*) 'value of a2 close to zero : ', a2
+!!$ stop
+!!$ endif
+!!$
+!!$ fa2=0.0
+!!$ do i = 1,npts_win
+!!$ fa2=fa2+((1/a2)*sngl(d_loc(i))-sngl(s_cor(i)))**2
+!!$ enddo
+!!$
+!!$ end function
+!!$
+!!$ end subroutine F1F2_calc
+!!$
+!!$! --------------------------------------------------------------------
+!!$
+!!$ subroutine xcorr_calc(d,s,npts,i1,i2,ishift,cc_max)
+!!$
+!!$ ! inputs:
+!!$ ! s(npts) = synthetic
+!!$ ! d(npts) = data (or observed)
+!!$ ! i1, i2 = start and stop indexes of window within s and d
+!!$
+!!$ double precision, dimension(*), intent(in) :: s,d
+!!$ integer, intent(in) :: npts, i1, i2
+!!$
+!!$ ! outputs:
+!!$ ! ishift = index lag (d-s) for max cross correlation
+!!$ ! cc_max = maximum of cross correlation (normalised by sqrt(synthetic*data))
+!!$ integer, intent(out) :: ishift
+!!$ double precision, intent(out) :: cc_max
+!!$
+!!$ ! local variables
+!!$ integer :: nlen
+!!$ integer :: i_left, i_right, i, j
+!!$ double precision :: cc
+!!$
+!!$ ! initialise shift and cross correlation to zero
+!!$ ishift=0
+!!$ cc_max=0
+!!$
+!!$ if (i1.gt.i2 .or. i2.gt.npts) then
+!!$ write(*,*) 'Error with window limits: i1, i2, npts ', i1, i2, npts
+!!$ return
+!!$ endif
+!!$
+!!$ ! length of window (number of points including ends)
+!!$ nlen = i2 - i1 + 1
+!!$
+!!$ ! left and right limits of index (time) shift search
+!!$ i_left=-1*int(nlen/2)
+!!$ i_right=int(nlen/2)
+!!$
+!!$
+!!$ ! i -> shift (to be applied to d in cc search)
+!!$ do i = i_left, i_right
+!!$ cc=0
+!!$ do j = i1, i2
+!!$ if((j+i).ge.1 .and. (j+i).le.npts) cc = cc + s(j)*d(j+i)
+!!$ enddo
+!!$ if (cc .gt. cc_max) then
+!!$ cc_max=cc
+!!$ ishift=i
+!!$ endif
+!!$ enddo
+!!$
+!!$ cc_max=cc_max / sqrt(sum(s(i1:i2)*s(i1:i2)) * sum(d(i1:i2)*(d(i1:i2))))
+!!$
+!!$end subroutine xcorr_calc
- double precision, dimension(*), intent(in) :: d, s
- integer, intent(in) :: npts,i1,i2
- double precision, intent (in) :: dt
- double precision, intent(out) :: F1,F2,tshift,cc_max,dlnA
-
- double precision, dimension(:), allocatable :: s_cor,d_loc
-
- double precision :: cr_shift, cc
- integer :: n_left,n_right,ishift,npts_win, i, j
-
- real ax,bx,cx,fa,fb,fc,f1_min,f2_min,f1_top,f1_bot,A1,A2
- real golden !f1,f2
-
- npts_win=i2-i1+1
-
-! allocate memory for s_cor (the corrected synthetic)
- allocate(s_cor(npts_win))
- allocate(d_loc(npts_win))
-
- d_loc(1:npts_win)=d(i1:i2)
-
-! do cross-correlation:
- call xcorr_calc(d,s,npts,i1,i2,ishift,cc_max)
-! n_left = int((-1.0) * cr_shift / dt)
-! n_right = int(cr_shift / dt)
-! ishift=0
-! cc_max=0.
-! do i = n_left, n_right
-! cc = 0
-! do j = 1, npts
-! if((j+i).gt.1.and.(j+i).lt.npts) cc = cc + s(j) * d(j+i)
-! enddo
-! if( cc .gt. cc_max) then
-! cc_max = cc
-! ishift = i
-! endif
-! enddo
- tshift=ishift*dt
-
-! apply time shift to synthetic seismogram
-! write(*,*)'shift synth seismogram by ', tshift, 'seconds'
- do i = 1, npts_win
- s_cor(i) = 0
- if( (i1-1+i-ishift) .gt. 1 .and. (i1-1+i-ishift) .lt.npts ) s_cor(i) = s(i1-1+i-ishift)
- enddo
-
-! DEBUG: output
-! open(unit=11, file='DEBUG_calcF1F2.dat')
-! do i = 1, npts_win
-! write(11,'(4(e12.4,1x))') b+(i-1)*dt, s_cor(i), s(i1-1+i), d(i1-1+i)
-! enddo
-! close(11)
-
-! calculate dlnA
- dlnA = sqrt( ( sum( d(i1:i2) * d(i1:i2) )) / (sum( s_cor(1:npts_win) * s_cor(1:npts_win) )) ) - 1
-
-
-! calculate F1, the least squares misfit
- f1_top=0.0
- f1_bot=0.0
- do i = 1,npts_win
- f1_top=f1_top+(sngl(d_loc(i))-sngl(s_cor(i)))**2
-! f1_bot=f1_bot+sqrt(sngl(d_loc(i))**2*sngl(s_cor(i))**2)
- f1_bot=f1_bot+sngl(d_loc(i))**2
- enddo
- if ( f1_bot .gt. 0.0 ) then
- F1 = dble(f1_top / f1_bot)
- else
- write(*,*) 'Sum d(t)**2 = 0 : empty observed seismogram.'
- F1=0
- F2=0
- return
- endif
-
-! do fa1 minimization to find A1
- ax=1e-3
- bx=1e3
- call mnbrak(ax,bx,cx,fa,fb,fc,fa1)
- f1_min=golden(ax,bx,cx,fa1,sngl(tol),A1)
-
-! do fa2 minimization to find A2
- ax=1e-3
- bx=1e3
- call mnbrak(ax,bx,cx,fa,fb,fc,fa2)
- f2_min=golden(ax,bx,cx,fa2,sngl(TOL),A2)
-
-! calculate F2
- F2=dble(min(A1,A2)/max(A1,A2))
-
-! Turn F1 around
- F1=1-F1
-
- deallocate(s_cor)
- deallocate(d_loc)
-
- contains
-
-! -----------------------------------------------------------------
-
- real function fa1(a1)
- real a1
-
- if (abs(a1).lt.TOL) then
- write(*,*) 'value of a1 close to zero : ', a1
- stop
- endif
-
- fa1=0.0
- do i = 1,npts_win
- fa1=fa1+(sngl(d_loc(i))-a1*sngl(s_cor(i)))**2
- enddo
-
- end function
-
-! -----------------------------------------------------------------
-
- real function fa2(a2)
- real a2
-
- if (abs(a2).lt.TOL) then
- write(*,*) 'value of a2 close to zero : ', a2
- stop
- endif
-
- fa2=0.0
- do i = 1,npts_win
- fa2=fa2+((1/a2)*sngl(d_loc(i))-sngl(s_cor(i)))**2
- enddo
-
- end function
-
- end subroutine F1F2_calc
-
-! --------------------------------------------------------------------
-
- subroutine xcorr_calc(d,s,npts,i1,i2,ishift,cc_max)
-
- ! inputs:
- ! s(npts) = synthetic
- ! d(npts) = data (or observed)
- ! i1, i2 = start and stop indexes of window within s and d
-
- double precision, dimension(*), intent(in) :: s,d
- integer, intent(in) :: npts, i1, i2
-
- ! outputs:
- ! ishift = index lag (d-s) for max cross correlation
- ! cc_max = maximum of cross correlation (normalised by sqrt(synthetic*data))
- integer, intent(out) :: ishift
- double precision, intent(out) :: cc_max
-
- ! local variables
- integer :: nlen
- integer :: i_left, i_right, i, j
- double precision :: cc
-
- ! initialise shift and cross correlation to zero
- ishift=0
- cc_max=0
-
- if (i1.gt.i2 .or. i2.gt.npts) then
- write(*,*) 'Error with window limits: i1, i2, npts ', i1, i2, npts
- return
- endif
-
- ! length of window (number of points including ends)
- nlen = i2 - i1 + 1
-
- ! left and right limits of index (time) shift search
- i_left=-1*int(nlen/2)
- i_right=int(nlen/2)
-
-
- ! i -> shift (to be applied to d in cc search)
- do i = i_left, i_right
- cc=0
- do j = i1, i2
- if((j+i).ge.1 .and. (j+i).le.npts) cc = cc + s(j)*d(j+i)
- enddo
- if (cc .gt. cc_max) then
- cc_max=cc
- ishift=i
- endif
- enddo
-
- cc_max=cc_max / sqrt(sum(s(i1:i2)*s(i1:i2)) * sum(d(i1:i2)*(d(i1:i2))))
-
-end subroutine xcorr_calc
-
-
! ------------------------------------------------------------------
! subroutine costaper(ipoint, ndata, tas)
! ------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list