[cig-commits] r14107 - in seismo/3D/ADJOINT_TOMO: . mtadj mtadj/dataset mtadj/utils
liuqy at geodynamics.org
liuqy at geodynamics.org
Thu Feb 19 16:01:36 PST 2009
Author: liuqy
Date: 2009-02-19 16:01:36 -0800 (Thu, 19 Feb 2009)
New Revision: 14107
Added:
seismo/3D/ADJOINT_TOMO/mtadj/
seismo/3D/ADJOINT_TOMO/mtadj/MEASUREMENT.WINDOWS
seismo/3D/ADJOINT_TOMO/mtadj/MTADJ.PAR
seismo/3D/ADJOINT_TOMO/mtadj/Makefile
seismo/3D/ADJOINT_TOMO/mtadj/OUTPUT/
seismo/3D/ADJOINT_TOMO/mtadj/call-graph
seismo/3D/ADJOINT_TOMO/mtadj/dataset/
seismo/3D/ADJOINT_TOMO/mtadj/dataset/test_data.sac
seismo/3D/ADJOINT_TOMO/mtadj/dataset/test_syn.sac
seismo/3D/ADJOINT_TOMO/mtadj/mod/
seismo/3D/ADJOINT_TOMO/mtadj/mtadj.f90
seismo/3D/ADJOINT_TOMO/mtadj/mtadj_constants.f90
seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub.f90
seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub2.f90
seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub3.f90
seismo/3D/ADJOINT_TOMO/mtadj/mtadj_variables.f90
seismo/3D/ADJOINT_TOMO/mtadj/obj/
seismo/3D/ADJOINT_TOMO/mtadj/utils/
seismo/3D/ADJOINT_TOMO/mtadj/utils/compile
seismo/3D/ADJOINT_TOMO/mtadj/utils/readme
seismo/3D/ADJOINT_TOMO/mtadj/utils/rotate_adj_src.f90
seismo/3D/ADJOINT_TOMO/mtadj/utils/rotate_adj_src.pl
Log:
New version of the multi-taper measurements and adjoint source calculation code
Added: seismo/3D/ADJOINT_TOMO/mtadj/MEASUREMENT.WINDOWS
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/MEASUREMENT.WINDOWS (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/MEASUREMENT.WINDOWS 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,5 @@
+1
+dataset/test_data.sac
+dataset/test_syn.sac
+1
+1200 1900
Added: seismo/3D/ADJOINT_TOMO/mtadj/MTADJ.PAR
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/MTADJ.PAR (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/MTADJ.PAR 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,15 @@
+2 %kernel type, WF=0 (waveform), CC=1 (cross-corr), FD=2 (freq-dep)
+2 %taper type for FD: BC=0 (boxcar), CS=1 (cosine-taper), MT=2 (multi-taper)
+OUTPUT
+OUTPUT
+.false. % bandpass or not
+50 500 % tshort, tlong
+40.0 % cc time shift max (before_shift)
+2.5 0.020 0.020 % npi, wtr, wtr_mtm
+.true. % select window after measurement
+1 % number of cycles in a window
+0.250 0.500 3.500 % dt_fac, err_fac, dt_max_scale
+.false. % include_error
+0.2 0.05 % min_dt_sigma, min_dlnA_sigma
+-1.0 1.0 10000 % b_adj,dt_adj, npts_adj
+.false. % bandpass adjoint src between [tstart,tend]
Added: seismo/3D/ADJOINT_TOMO/mtadj/Makefile
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/Makefile (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/Makefile 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,29 @@
+F90 = gfortran
+F90_FLAGS = -fbounds-check
+F77 = gfortran
+F77_FLAGS = -ffixed-line-length-132
+
+LIB = -L/opt/sac/lib/ -lsacio -lsac
+
+MT_F90_SRC = mtadj_constants mtadj_variables mtadj_sub3 mtadj_sub2 mtadj_sub
+MT_SRC = mtadj
+
+MOD_DIR = mod
+OBJ_DIR = obj
+BIN_DIR = .
+
+
+MT_F90_OBJ = $(patsubst %,$(OBJ_DIR)/%.o,$(MT_F90_SRC))
+MT_OBJ = $(MT_F77_OBJ) $(MT_F90_OBJ)
+
+$(MT_SRC) : $(MT_SRC).f90 $(MT_OBJ)
+ $(F90) -o $(BIN_DIR)/$@ $(F90_FLAGS) $@.f90 -I $(MOD_DIR) $(MT_OBJ) $(LIB)
+
+$(MT_F90_OBJ): $(OBJ_DIR)/%.o : %.f90
+ $(F90) -o $@ $(F90_FLAGS) -c $*.f90 -I $(MOD_DIR) -J$(MOD_DIR)
+
+.PHONY : clean
+
+clean:
+ \rm -f *.o *.mod *~ $(OBJ_DIR)/*.o $(MOD_DIR)/*.mod $(MT_SRC) OUTPUT/*
+
Added: seismo/3D/ADJOINT_TOMO/mtadj/call-graph
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/call-graph (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/call-graph 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,100 @@
+This package is rewritten from the old version (named measure_adj), mainly
+to reflect the updated sequence of processing,
+i.e. we now pick time windows using FLEXWIN package before
+running the measurement code, which should have eliminated the bad picks.
+I also try to take advantage of the modules in Fortran 90 and instead of
+passing everything through the subroutine arguments, passing them through
+the global variables (since we have one module dedicated for this anyways).
+
+Concerns:
+
+1. Extra manipulations to the adjoint source, such as time-domain taper
+and band-passing (interpolation not included), which are not reflected
+on evaluation of the misfit function, need to be done with care.
+I basically set time windowing to be boxcar, and let bandpassing the signal
+to be an option.
+
+2. Presumably windowing code should have been run before making the measure-
+ments, therefore, BEFORE_SHIFT may not be useful, or you could probably get
+rid of the selection process. I took out the selection criteria before_quality
+after_shift, after_quality. They may turn out to be useful at some point, and
+they are not very hard to add back.
+
+3. I am pretty clueless when it comes to how to set MIN_SIGMA_DT and MIN_SIGMA_DTAU, just to make the adjoint source stable. I am also not sure how to compute the errors for tshift_cc and dlnA, of course misfit is one factor, also the frequency content should also be taken into account.
+
+4. comparing CC with MT(or FD in general), one can set INCLUDE_ERROR = false
+
+===================================================
+
+mtadj_sub.f90 =====
+ subroutine read_mtadj_par(par_file)
+ subroutine read_data_syn(datafile,synfile,sta,chan,net)
+ subroutine cc_fd_measure(file_prefix,tstart,tend)
+ subroutine select_cc_fd_measure(tstart, tend, use_window)
+ subroutine mt_adj_src(dt_adj_src,amp_adj_src,dt_chi,amp_chi)
+ subroutine adjust_adj_src(dt_adj_src,amp_adj_src,nlen,tstart,dt,&
+ dt_adj_src_win,amp_adj_src_win,npts_adj,b_adj,dt_adj)
+
+mtadj_sub2.f90 ====
+ call rsac1(datafile,data,npts1,b1,dt1,NDIM,nerr)
+ call getkhv('kstnm', sta, nerr)
+ call wsac1(trim(file_prefix)//'.obs.sac',dataw,nlen,tstart,dt,nerr)
+ call compute_time_shift(synw,dataw,nlen,dt,ishift,tshift_cc,BEFORE_SHIFT)
+ call reconstruct_syn_cc(synw,nlen,dt,ishift,dlnA,synw_rc_cc)
+ call compute_cc_error(dataw,synw_rc_cc,nlen,dt,i_pmax,dlnA,sigma_tshift_cc,sigma_dlnA_cc,MIN_DT_SIGMA,MIN_DlnA_SIGMA)
+ call fft(LNPT,cdataw,FORWARD_FFT,dt)
+ call staper(nlen, NPI, NTAPER, tas, NPT, ey1, ey2)
+ call costaper(nlen, NPT, tas)
+ call boxcar(nlen, NPT, tas)
+ call reconstruct_syn_fd(csynwt,dtau_fdm,dlnA_fdm,i_right,synw_rc_fd,dt,nlen)
+ call compute_dtau_dlnA(trans_fdm,tshift_cc,dtau_fdm,dlnA_fdm,i_right)
+ call compute_mt_error(ntaper,dataw,synw,tas,dt,wtr,i_right,tshift_cc, &
+ dtau_fdm,dlnA_fdm,sigma_dtau_fdm, sigma_dlnA_fdm)
+ call compute_fd_error(npi,nlen,i_right,dtau_fdm,dlnA_fdm, &
+ sigma_dtau_fdm,sigma_dlnA_fdm)
+ call compute_veloc_from_displ(synw,nlen,dt,synw_veloc)
+ call interp_adj_src(dt_adj_src,nlen,tstart,dt, &
+ dt_adj_src_win,npts_adj,b_adj,dt_adj)
+
+
+====================================================
+
+-- NDIM used only for data, syn, dt_adj_src_all, amp_adj_src_all,
+ dt_adj_src_win, amp_adj_src_win
+ most of other arrays actually use NPT as the dimension in accordance
+ with nlen
+
+-- fft and fftinv convention:
+ fft(nn,cdata,forward_fft,dt) outputs the true Fourier spectrum as
+ defined in the continuos FT.
+ same as fftinv(nn,cdata,reverse_fft,dt,data)
+ although only the first half of cdata is used, and the
+ second half mirror imaged to produce real 'data'.
+
+-- xapiir from libsac.a
+ xapiir(data,npts,'BU',trbddw,apara,iord,'bp',flo,fhi,dt,passes)
+ where trbddw,apara, flo,fhi and dt have to be double, while data is float
+
+========================================================
+
+input file:
+
+2 %kernel type, WF=0 (waveform), CC=1 (cross-corr), FD=2 (freq-dep)
+2 %taper type for FD: BC=0 (boxcar), CS=1 (cosine-taper), MT=2 (multi-taper)
+OUTPUT_MT
+OUTPUT_MT
+.false. % bandpass or not
+50 500 % tshort, tlong
+40.0 % cc time shift max (before_shift)
+2.5 0.020 0.020 % npi, wtr, wtr_mtm
+.true. % select window after measurement
+1 % number of cycles in a window
+0.250 0.500 3.500 % dt_fac, err_fac, dt_max_scale
+.false. % include_error
+0.2 0.05 % min_dt_sigma, min_dlnA_sigma
+-1.0 1.0 10000 % b_adj,dt_adj, npts_adj
+.false. % bandpass adjoint src between [tstart,tend]
+
+
+Qinya Liu
+Feb 19, 2009
Added: seismo/3D/ADJOINT_TOMO/mtadj/dataset/test_data.sac
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/mtadj/dataset/test_data.sac
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: seismo/3D/ADJOINT_TOMO/mtadj/dataset/test_syn.sac
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/mtadj/dataset/test_syn.sac
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: seismo/3D/ADJOINT_TOMO/mtadj/mtadj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/mtadj.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/mtadj.f90 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,134 @@
+! this program computes the cross-correlation/transfer function measurements between
+! data and synthetics in a WINDOWING file
+! and output the corresponding adjoint source for tomography inversions
+
+program mtadj
+
+ use mtadj_constants
+ use mtadj_variables
+ use mtadj_sub
+ ! use mtadj_sub2
+
+ implicit none
+
+ ! sac header information
+ character(len=10) :: net,sta,chan
+
+ ! file prefix
+ character(len=150) :: datafile, synfile
+ character(len=150) :: file_prefix, adj_prefix, meas_prefix
+
+ ! npairs of data and syn, number of windows
+ integer :: npairs, nwin, nwin_total, nadj_src
+ integer :: ios, ipair, j, nerr
+
+ ! data and syn measurements/adjoint outpu
+ real :: tstart, tend
+ logical :: use_adj_trace, use_window
+ real*8 :: fstart_dble,fend_dble,dt_dble,dt_adj_dble
+
+ real,dimension(NPT) :: dt_adj_src,amp_adj_src
+ real, dimension(NDIM) :: dt_adj_src_win, amp_adj_src_win, &
+ dt_adj_src_all, amp_adj_src_all
+ real :: dt_chi,amp_chi
+
+
+ ! --- PROGRAM STARTS HERE ---
+ ! read parameter file
+ call read_mtadj_par('MTADJ.PAR')
+ fstart_dble=fstart
+ fend_dble=fend_dble
+ dt_adj_dble=dt_adj
+
+ ! loop over measurement file
+ open(11,file='MEASUREMENT.WINDOWS',status='old',iostat=ios)
+ if (ios /= 0) stop 'Error opening input file: MEASUREMENT WINDOWS'
+ read(11,*,iostat=ios) npairs
+ if (ios /= 0) stop 'Error reading number of pairs of data/syn'
+
+
+ nwin_total=0
+ do ipair = 1, npairs
+
+ ! read data and syn pair
+ read(11,'(a)',iostat=ios) datafile
+ if (ios /= 0) stop 'Error reading datafile'
+ read(11,'(a)',iostat=ios) synfile
+ if (ios /= 0) stop 'Error reading synfile'
+ if (DEBUG) print *, trim(datafile), ' ', trim(synfile)
+
+ call read_data_syn(datafile,synfile,sta,net,chan)
+
+ ! output file prefix
+ file_prefix=trim(sta)//'.'//trim(net)//'.'//trim(chan)
+
+ dt_dble=dt
+ ! filter data and synthetics (check xapiir() usage in sac lib)
+ if (BANDPASS) then
+ call xapiir(data,npts,'BU',TRBDNDW,APARM,IORD,'BP',fstart_dble,fend_dble,dt_dble,PASSES)
+ call xapiir(syn,npts,'BU',TRBDNDW,APARM,IORD,'BP',fstart_dble,fend_dble,dt_dble,PASSES)
+ endif
+
+ ! loop over nwin
+ read(11,*) nwin
+ if (nwin < 0) stop 'Check nwin '
+ nwin_total = nwin_total + nwin
+
+ ! initialize
+ use_adj_trace = .false.
+ dt_adj_src_all = 0.; amp_adj_src_all = 0.; nadj_src = 0
+
+ do j = 1, nwin
+ read(11,*,iostat=ios) tstart, tend
+ if (ios /= 0) stop 'Error reading tstart and tend'
+ print *, ' Measurement No. ', j, '...', tstart, tend
+ tstart = max(tstart,b)
+ tend = min(tend, b+(npts-1)*dt)
+
+ if (iker == IKER_CC .or. iker == IKER_FD) then
+ call cc_fd_measure(file_prefix, tstart, tend)
+ if (SELECT_WINDOW) then
+ call select_cc_fd_measure(tstart, tend, use_window)
+ else
+ use_window = .true.
+ endif
+ endif
+
+ if (OUTPUT_ADJSRC .and. use_window) then
+ call mt_adj_src(file_prefix,j,tstart,dt_adj_src,amp_adj_src,dt_chi,amp_chi)
+
+ ! taper and interpolate into adj_src_all(t0,dt,npts)
+ call adjust_adj_src(dt_adj_src,amp_adj_src,nlen,tstart,dt, &
+ dt_adj_src_win,amp_adj_src_win,npts_adj,b_adj,dt_adj)
+ endif
+
+ if (use_window) use_adj_trace = .true.
+
+ enddo ! nwin
+
+ if (OUTPUT_ADJSRC .and. use_adj_trace) then
+ dt_adj_src_all = dt_adj_src_all + dt_adj_src_win
+ amp_adj_src_all = amp_adj_src_all + amp_adj_src_win
+ if (BANDPASS_ADJ) then
+ call xapiir(dt_adj_src_all,npts,'BU',TRBDNDW,APARM,IORD,'BP', &
+ fstart_dble,fend_dble,dt_adj_dble,PASSES)
+ endif
+ if (iker /= IKER_FD) then
+ adj_prefix=trim(CKER(iker+1))//'.adj'
+ else
+ adj_prefix=trim(CTAP(itap+1))//'.adj'
+ endif
+ adj_prefix=trim(adj_dir)//'/'//trim(file_prefix)//'.'//trim(adj_prefix)
+
+ call wsac1(trim(adj_prefix),dt_adj_src_all,npts_adj,b_adj,dt_adj,nerr)
+ ! write amplitude adjoint sources here
+ nadj_src = nadj_src + 1
+ if (nerr > 0) stop 'Error writing adjoint source'
+ endif
+
+ enddo ! nfiles
+
+ print *, 'Total number of windows processed ', nwin_total
+ print *, 'Total number of adjoint sources ', nadj_src
+
+end program mtadj
Added: seismo/3D/ADJOINT_TOMO/mtadj/mtadj_constants.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/mtadj_constants.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/mtadj_constants.f90 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,50 @@
+module mtadj_constants
+
+ ! adjoint source types
+
+ integer, parameter :: IKER_WF = 0
+ integer, parameter :: IKER_CC = 1
+ integer, parameter :: IKER_FD = 2
+ character(len=2), parameter :: CKER(3) = (/ 'wf', 'cc', 'fd' /)
+
+ ! taper types
+ integer, parameter :: ITAP_BC = 0
+ integer, parameter :: ITAP_CS = 1
+ integer, parameter :: ITAP_MT = 2
+ character(len=2), parameter :: CTAP(3) = (/ 'bc', 'cs', 'mt' /)
+
+ ! maximum number of tapers
+ integer, parameter :: NMAX_TAPER = 10
+
+ ! debug flags
+ logical, parameter :: DEBUG = .true.
+ logical, parameter :: OUTPUT_MEASURE = .true.
+ logical, parameter :: OUTPUT_ADJSRC = .true.
+
+ ! constants
+ real, parameter :: PI = 3.141592653
+ real, parameter :: TWOPI = 2.0 * PI
+ complex ,parameter :: CCI = cmplx(0.,1.)
+ real, parameter :: LARGE_VAL = 1.0d8
+ real, parameter :: EPS_dt=1.0e-4
+
+ ! FFT parameters
+ integer, parameter :: LNPT = 13, NPT = 2**LNPT, NDIM = 40000
+ real*8, parameter :: FORWARD_FFT = 1.0
+ real*8, parameter :: REVERSE_FFT = -1.0
+
+ ! phase correction control parameters, set this between (PI, 2PI),
+ ! use a higher value for conservative phase wrapping
+ real, parameter :: PHASE_STEP = 1.5 * PI
+
+ ! filter parameters for xapiir subroutine (filter type is BP)
+ ! TRBDNDW and APARM are used only for chebyshev filters
+ ! 2-passes guarantees a zero-phase filter (non-causal,which is ok)
+ real*8, parameter :: TRBDNDW = 0.3
+ real*8, parameter :: APARM = 30.
+ integer, parameter :: IORD = 4
+ integer, parameter :: PASSES = 2
+
+
+end module mtadj_constants
+
Added: seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub.f90 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,644 @@
+module mtadj_sub
+
+ use mtadj_constants
+ use mtadj_variables
+ use mtadj_sub2
+ use mtadj_sub3
+
+ implicit none
+
+contains
+
+! =========================================================
+
+ subroutine read_mtadj_par(par_file)
+
+ character(len=*) par_file
+ integer :: ios
+
+ open(10,file=trim(par_file),status='old',iostat=ios)
+ if (ios /= 0) stop 'Error opening parameter file'
+
+ ! kernel type, WF=0 (waveform), CC=1 (cross-corr), FD=2 (freq-dep)
+ read(10,*) iker
+ ! taper type for FD meas: BC=0 (boxcar), CS=1 (cosine-taper), MT=2 (multi-taper)
+ read(10,*) itap
+ read(10,'(a)') meas_dir
+ read(10,'(a)') adj_dir
+
+ ! bandpass data and syn
+ read(10,*) BANDPASS
+ read(10,*) tshort, tlong
+
+ ! cross-correlation time shift maximum
+ read(10,*) before_shift
+
+ ! multi-taper setting
+ read(10,*) npi,wtr,wtr_mtm
+
+ ! select windows
+ read(10,*) SELECT_WINDOW
+ ! FD meas window length, window shorter than this length
+ ! will resort to CC measurements and kernels
+ read(10,*) ncycle_in_window
+ read(10,*) dt_fac, err_fac, dt_max_scale ! ???
+ ! read(10,*) after_quality, after_tshift, dlna_sigma_min
+
+ ! interpolation of adjoint source
+ read(10,*) INCLUDE_ERROR
+ read(10,*) MIN_DT_SIGMA,MIN_DlnA_SIGMA
+ read(10,*) b_adj, dt_adj, npts_adj
+ read(10,*) BANDPASS_ADJ
+
+ close(10)
+
+ ! list input parameters:
+ write(*,*) '========= INPUTS FROM MEASUREMENT.PAR ============'
+ if (iker == IKER_WF) then
+ write(*,*) 'Adjoint source type: Waveform '
+ else if (iker == IKER_CC) then
+ write(*,*) 'Adjoint source type: Cross-correlation (T & A)'
+ else if (iker == IKER_FD) then
+ write(*,*) 'Adjoint source type: Frequency-dependent (T & A)'
+ else
+ stop 'Only allow iker = IKER_WF/CC/FD'
+ endif
+
+ if (iker == IKER_FD) then
+ if (itap == ITAP_BC) then
+ write(*,*) 'Freq dep measurements use taper: boxcar'
+ ntaper = 1
+ else if (itap == ITAP_CS) then
+ write(*,*) 'Freq dep measurements use taper: cosine taper'
+ ntaper = 1
+ else if (itap == ITAP_MT) then
+ write(*,*) 'Freq dep measurements use taper: multi-taper'
+ ntaper = 2*npi
+ write(*,*) ' only for windows with at least ', ncycle_in_window, ' number of cycles'
+ write(*,*) ' number of multi-tapers = ', ntaper
+ endif
+ write(*,*) ' only output results above water level ', wtr
+ if (ntaper > NMAX_TAPER) stop 'ntaper should be <= NMAX_TAPER'
+ endif
+
+ if (BANDPASS) then
+ write(*,*) 'Run bandpass before making measurements between',tshort, tlong
+ fstart = 1./tlong; fend = 1./ tshort
+ else
+ write(*,*) 'No filtering before making measurements'
+ endif
+ write(*,*) 'cross-correlation time shift only allowed up to ', before_shift, ' secs'
+
+ if (SELECT_WINDOW) then
+ write(*,*) 'Only output adjoint source for selected windows based on:'
+ write(*,*) dt_fac,dt_max_scale,err_fac
+ else
+ write(*,*) 'Adjoint sources for all windows will be used'
+ endif
+
+ if (npts_adj > NDIM) stop 'Npts for the adjoint source exceeds the array limit'
+ if (NPT > NDIM) stop 'Check NPT > NDIM '
+ write(*,*) 'Interpolate output adjoint source onto: ', b_adj, dt_adj, npts_adj
+ if (INCLUDE_ERROR) then
+ write(*,*) 'Errors are included in the adjoint sources'
+ else
+ write(*,*) 'Errors are not included in the adjoint sources -- test only !!'
+ endif
+ if (iker == IKER_CC) then
+ write(*,*) ' CC adjoint source minimum error: ', MIN_DT_SIGMA,MIN_DlnA_SIGMA
+ endif
+ write(*,*) ' '
+
+ ! frequency-domain taper (power of cosine) for adjoint source calculation
+ ipwr_w = 4
+ ! time domain taper for adjoint source
+ ipwr_t = 10
+
+ end subroutine read_mtadj_par
+
+ !========================================================
+
+ subroutine read_data_syn(datafile,synfile,sta,chan,net)
+
+ character(len=*),intent(in) :: datafile,synfile
+ character(len=*),intent(out) :: sta,chan,net
+
+ integer :: npts1, npts2, nerr, j
+ real :: b1, b2, dt1, dt2
+
+ ! read data and syn
+ call rsac1(datafile,data,npts1,b1,dt1,NDIM,nerr)
+ if (nerr > 0) stop ' Error reading synthetic file'
+ call rsac1(synfile,syn,npts2,b2,dt2,NDIM,nerr)
+ if (nerr > 0) stop ' Error reading data file '
+ if (npts1 /= npts2 .or. abs(b1-b2) > dt1 .or. abs(dt1-dt2) > EPS_dt) &
+ stop 'check if npts1=npts2, b1=b2 and dt1=dt2'
+
+ npts=npts1; b=b1; dt=dt1
+
+ ! sac header sta.net.chan, and file prefixes
+ call getkhv('kstnm', sta, nerr)
+ call getkhv('kcmpnm', chan, nerr)
+ call getkhv('knetwk', net, nerr)
+
+ ! add comparsion for date and time to make sure they correspond
+ ! to each other
+
+ end subroutine read_data_syn
+
+!============================================================
+
+ subroutine cc_fd_measure(file_prefix,tstart,tend)
+
+ character(len=*) file_prefix
+ real :: tstart, tend
+ real, dimension(NPT) :: dataw_save, synw_rc_cc, synw_rc_fd
+ real*8, dimension(NPT) :: ey1, ey2
+ integer :: ishift, i, nerr
+ real :: f0, df, df_new, ampmax_syn, wtr_amp_syn, ampmax_bot, wtr_amp_bot
+ integer :: nf, idf_new, iampmax_syn, i_right_stop, iampmax_bot, ictaper
+ integer :: i_ampmax_syn(1)
+ complex*16, dimension(NPT) :: cdataw, csynw, cdatawt,csynwt
+ complex,dimension(NPT) :: top_fdm, bot_fdm, trans_fdm,csynw_sngl
+ character(len=150) :: meas_prefix
+ real*8,dimension(NPT,NMAX_TAPER) :: tas_dp
+ real*8 :: dt_dble
+
+
+ ! set measurement file prefix
+ if (iker /= IKER_FD) then
+ meas_prefix=trim(CKER(iker+1))
+ else
+ meas_prefix=trim(CTAP(itap+1))
+ endif
+ meas_prefix=trim(meas_dir)//'/'//trim(file_prefix)//'.'//trim(meas_prefix)
+
+ ! at least 10 points within the window
+ if (tend-tstart < 10*dt) stop 'Check if tend > tstart+10*dt'
+ nstart = floor((tstart-b)/dt)+1
+ nend = floor((tend-b)/dt)
+ nlen = nend-nstart+1
+ if (nlen > NPT) stop 'Check if the window length is over NPT'
+
+ ! window data and syn (data -> dataw; syn -> synw)
+ dataw(1:nlen) = data(nstart:nend)
+ synw(1:nlen) = syn(nstart:nend)
+
+ ! for BC or CS tapers, rmean and rtrend
+ if (iker == IKER_FD .and. itap == ITAP_BC) then
+ call rmean(dataw,nlen); call rmean(synw,nlen)
+ call rtrend(dataw,nlen); call rtrend(synw,nlen)
+ endif
+
+ if (DEBUG) then
+ ! write windowed data and synthetics as sac files
+ call wsac1(trim(meas_prefix)//'.obs.sac',dataw,nlen,tstart,dt,nerr)
+ if (nerr > 0) stop 'Error writing obs data'
+ call wsac1(trim(meas_prefix)//'.syn.sac',synw,nlen,tstart,dt,nerr)
+ if (nerr > 0) stop 'Error writing synthetics'
+ endif
+
+ ! do not need t-domain tapers for windowed data and syn at this point
+
+ ! ============== CC measurements =======================
+ ! save a copy of observed data
+ dataw_save(1:nlen) = dataw(1:nlen)
+
+ ! tshift_cc: synw(t) * dataw(t+tshift_cc)
+ call compute_time_shift(synw,dataw,nlen,dt,ishift,tshift_cc)
+ if (abs(tshift_cc) > BEFORE_SHIFT) &
+ stop 'Check if BEFORE_SHIFT is too small for tshift'
+
+ ! align data and synthetics according to CC
+ do i = 1, nlen
+ if ( (i+nstart-1+ishift) >= 1 .and. (i+nstart-1+ishift) <= npts) then
+ dataw(i) = data(i+ishift-1+nstart)
+ else
+ dataw(i) = 0. ! may be should be replaced by 1st and last value of data(:)
+ endif
+ enddo
+ if (DEBUG) then
+ call wsac1(trim(meas_prefix)//'.obs_shifted.sac',dataw,nlen,tstart,dt,nerr)
+ if (nerr > 0) stop 'Error writing obs shifted data'
+ endif
+
+ ! CC measurements - amplitude dlnA
+ dlnA = 0.5 * log( sum(dataw(1:nlen)**2) / sum(synw(1:nlen)**2) )
+
+ ! reconstruct best-fitting syn according to tshift_cc and dlnA
+ call reconstruct_syn_cc(synw,nlen,dt,ishift,dlnA,synw_rc_cc)
+ if (DEBUG) then
+ call wsac1(trim(meas_prefix)//'.syn_rc_cc.sac',synw_rc_cc,nlen,tstart,dt,nerr)
+ if (nerr > 0) stop 'Error writing CC reconstructed syn'
+ endif
+
+ ! -------
+ ! frequency vector from fft
+ f0=0; df=1./(NPT*dt); nf=floor(NPT/2.)+1
+
+ ! true independent frequency spacing for nlen
+ df_new=1./(tend-tstart); idf_new = int(df_new/df)
+
+ ! FFT windowed data (shifted) and syn
+ cdataw = cmplx(0.,0.); csynw = cmplx(0.,0.)
+ cdataw(1:nlen)=cmplx(dataw(1:nlen)); csynw(1:nlen)=cmplx(synw(1:nlen))
+
+ call fft(LNPT,cdataw,FORWARD_FFT,dble(dt))
+ call fft(LNPT,csynw,FORWARD_FFT,dble(dt))
+
+ ! check the highest trustable frequency according to synthetics water level
+ ampmax_syn = maxval(abs(csynw(1:nf)))
+ i_ampmax_syn = maxloc(abs(csynw(1:nf)))
+ wtr_amp_syn = cmplx(ampmax_syn * wtr, 0.)
+ i_pmax=i_ampmax_syn(1)
+
+ ! estimate tshift, dlnA uncertainties
+ ! according to misfit between shifted data and reconstructed syn
+ call compute_cc_error(dataw_save,synw_rc_cc,nlen,dt,i_pmax,dlnA,&
+ sigma_tshift_cc,sigma_dlnA_cc,MIN_DT_SIGMA,MIN_DlnA_SIGMA)
+
+ ! write measurement file
+ open(30,file=trim(meas_prefix)//'.dt_dlnA',status='unknown')
+ write(30,*) tshift_cc, sigma_tshift_cc
+ write(30,*) dlnA, sigma_dlnA_cc
+ close(30)
+
+ ! DONE if just cross-correlation measurements
+ if (iker /= IKER_FD) then
+ dataw(1:nlen) = dataw_save(1:nlen)
+ return
+ endif
+
+ ! ============= FD measurements ====================
+
+ ! corresponding index to wtr_use_unw -> i_right
+ i_right = nf; i_right_stop = 0
+ do i = 1,nf
+ if (abs(csynw(i))<=abs(wtr_amp_syn) .and. i_right_stop==0 .and. i>i_pmax ) then
+ i_right_stop = 1; i_right = i
+ endif
+ if (abs(csynw(i))>=10.*abs(wtr_amp_syn) .and. i_right_stop==1 .and. i>i_pmax) then
+ i_right_stop = 0; i_right = i
+ endif
+ enddo
+ f_right = df*i_right
+
+ if (itap == ITAP_MT) then
+ df_fd = npi*df_new; idf_fd = floor(df_fd/df); ; i_left=idf_fd
+ else ! ignore first few points for cosine and boxcar tapers
+ df_fd = df_new; idf_fd = floor(df_fd/df); i_left= floor(npi*df_fd/df)
+ endif
+ f_left=df_fd
+
+ if (DEBUG) then
+ print *, 'Frequency of max power in windowed synthetic:'
+ print *, 'f_pmax = ', i_pmax * df, 'Hz, T_pmax = ', 1./(i_pmax*df), 'secs'
+ print *, 'Frequency spacing df_new = ', df_new, ' Hz'
+ print *, ' Longest T = ', 1/df_new, 's; shortest T = ', 1/(i_right * df),'secs'
+ print *, ' i_right = ', i_right
+ print *, 'Independent MTM spacing df_meas = ', df_new * npi, ' Hz'
+ endif
+
+ ! define tapers for FD measurements: tas(1:nlen,1:ntaper)
+ if (itap == ITAP_MT) then
+ call staper(nlen, dble(NPI), ntaper, tas_dp, NPT, ey1, ey2)
+ elseif (itap == ITAP_CS) then
+ call costaper(nlen, NPT, tas_dp)
+ elseif (itap == ITAP_BC) then
+ call boxcar(nlen, NPT, tas_dp)
+ endif
+ tas=tas_dp
+
+ ! compute transfer function for freq-dep measurements
+ top_fdm(:) = cmplx(0.,0.)
+ bot_fdm(:) = cmplx(0.,0.)
+ trans_fdm(:) = cmplx(0.,0.)
+
+ do ictaper = 1, ntaper
+
+ ! tapered data (win+shifted) and syn (win)
+ cdatawt(:)=cmplx(0.,0.); csynwt(:) = cmplx(0.,0.)
+ cdatawt(1:nlen)=cmplx(dataw(1:nlen)*tas(1:nlen,ictaper))
+ csynwt(1:nlen)=cmplx(synw(1:nlen)*tas(1:nlen,ictaper))
+
+ ! complex sepctra
+ call fft(LNPT,cdatawt,FORWARD_FFT,dble(dt)) ! syn
+ call fft(LNPT,csynwt,FORWARD_FFT,dble(dt)) ! data
+
+ ! top and bottom of transfer function
+ top_fdm(1:nf) = top_fdm(1:nf) + cdatawt(1:nf) * conjg(csynwt(1:nf))
+ bot_fdm(1:nf) = bot_fdm(1:nf) + csynwt(1:nf) * conjg(csynwt(1:nf))
+
+ enddo
+
+ ! water level of bottom
+ ampmax_bot = maxval(abs(bot_fdm(1:nf)))
+ if (itap == ITAP_MT) then
+ wtr_amp_bot = ampmax_bot * (wtr_mtm ** 2)
+ else
+ wtr_amp_bot = ampmax_bot * (wtr ** 2)
+ endif
+ do i = 1, nf
+ if(abs(bot_fdm(i)) > abs(wtr_amp_bot)) then
+ trans_fdm(i) = top_fdm(i) / bot_fdm(i)
+ else if(abs(bot_fdm(i)) < abs(wtr_amp_bot)) then
+ trans_fdm(i) = top_fdm(i) / (bot_fdm(i)+wtr_amp_bot)
+ endif
+ enddo
+
+ ! compute dtau_fdm and dlnA_fdm
+ call compute_dtau_dlnA(trans_fdm,dt,tshift_cc,dtau_fdm,dlnA_fdm,i_right)
+
+ ! reconstruct syn with transfer function
+ csynw_sngl=csynw
+ call reconstruct_syn_fd(csynw_sngl,dtau_fdm,dlnA_fdm,i_right,synw_rc_fd,dt,nlen)
+
+ if (DEBUG) then
+ call wsac1(trim(meas_prefix)//'.syn_rc_fd.sac',synw_rc_fd,nlen,tstart,dt,nerr)
+ if (nerr > 0) stop 'Error writing FD reconstructed syn'
+ endif
+
+ ! estimate transfer function error
+ if (itap == ITAP_MT) then
+ ! for multi-taper
+ call compute_mt_error(ntaper,dataw,synw,tas,&
+ nlen,dt,wtr_mtm,i_right,tshift_cc, &
+ dtau_fdm,dlnA_fdm,sigma_dtau_fdm, sigma_dlnA_fdm)
+ else
+ ! for single tapers (BC or CS), use generic error estimation technique
+ call compute_fd_error(npi,nlen,i_right,dt,dtau_fdm,dlnA_fdm,&
+ sigma_dtau_fdm,sigma_dlnA_fdm)
+ endif
+
+ ! write transfer function measurements
+ open(30,file=trim(meas_prefix)//'.dtau',status='unknown')
+ open(40,file=trim(meas_prefix)//'.dlnA',status='unknown')
+
+ do i = i_left, i_right, idf_new ! not all measurements are indep.
+ write(30,'(3f14.4,i6)') df*(i-1), dtau_fdm(i), sigma_dtau_fdm(i), (i-i_left)/idf_fd
+ write(40,'(3f14.4,i6)') df*(i-1), dlnA_fdm(i), sigma_dlnA_fdm(i), (i-i_left)/idf_fd
+ enddo
+ close(30)
+ close(40)
+
+ ! adjust errors according to min_dt_sigma, min_dlnA_sigma
+ do i = 1, i_right
+ sigma_dtau_fdm(i) = max(sigma_dtau_fdm(i),min_dt_sigma)
+ sigma_dlnA_fdm(i) = max(sigma_dlnA_fdm(i),min_dlnA_sigma)
+ enddo
+
+ dataw(1:nlen) = dataw_save(1:nlen)
+
+ end subroutine cc_fd_measure
+
+! ======================================================================
+
+ subroutine select_cc_fd_measure(tstart, tend, use_window)
+
+ real :: tstart, tend
+ logical :: use_window
+ real :: df, f_pmax, T_pmax, wlen
+ integer :: nf, j
+ real, dimension(NPT) :: fvec
+
+
+ use_window=.true.
+ df = 1./(dt*NPT); nf=floor(NPT/2.)+1
+ do j = 1, nf
+ fvec(j) = df*(j-1)
+ enddo
+
+ ! at least N freq points between maximum power and f=0
+ f_pmax = df * i_pmax
+ T_pmax = 1./ f_pmax
+ wlen = dt*nlen
+ if( ncycle_in_window * T_pmax > wlen ) then
+ use_window = .false.
+ print *, 'rejecting window [', tstart, tend, ']'
+ if (DEBUG) print *, ' wlen > ncycle * t_pmax= ', wlen, T_pmax
+ return
+ endif
+
+ ! at least 2 independent measurements between f_left and f_right
+ if ( (f_right-f_left) < df_fd) then
+ use_window = .false.
+ print *, 'rejecting window [',tstart, tend, ']'
+ print *, ' --> considering CC adjoint source for this window'
+ if (DEBUG) print *, 'f_left,f_right = ', f_left, f_right
+ return
+ endif
+
+ ! require all FD measurements to be sensible
+ if (iker == IKER_FD) then
+ do j = i_left, i_right, idf_fd
+ if ( abs(dtau_fdm(j)) > 1./(dt_fac*fvec(j)) .or. &
+ abs(dtau_fdm(j)) > dt_max_scale*abs(tshift_cc) .or. &
+ sigma_dtau_fdm(j) > 1./(err_fac*fvec(j)) ) then
+ if (DEBUG) then
+ write(*,*) 'Rejecting window ...'
+ write(*,*) j, dtau_fdm(j), 1./(dt_fac*fvec(j)), &
+ dt_max_scale*abs(tshift_cc)
+ write(*,*) ' ', sigma_dtau_fdm(j), 1./(err_fac*fvec(j))
+ endif
+ use_window = .false.; return
+ endif
+ enddo
+ endif
+
+ end subroutine select_cc_fd_measure
+
+
+! =============================================================================
+
+
+ subroutine mt_adj_src(file_prefix,iwin,tstart,dt_adj_src,amp_adj_src,dt_chi,amp_chi)
+
+ character(len=*) :: file_prefix
+ integer, intent(in) :: iwin
+ real,intent(in) :: tstart
+ real, dimension(NPT) :: dt_adj_src, amp_adj_src
+ real :: dt_chi,amp_chi
+
+ real, dimension(NPT) :: synw_veloc, synwt, synwt_veloc,&
+ ft_bar, fa_bar, wf_taper, wp_taper, wq_taper, &
+ d_bot_mtm, v_bot_mtm
+ real*8,dimension(NPT) :: dtau_pj_t, dlnA_qj_t
+ real :: ffac, dtau_wtr, dlnA_wtr, err_t, err_a, wtr_v_bot, wtr_d_bot,df
+ complex*16, dimension(NPT) :: csynwt, csynwt_veloc, dtau_pj, dlnA_qj
+ complex,dimension(NPT,NMAX_TAPER) :: pwc_adj, qwc_adj
+ character(len=150) :: file
+ integer :: ictaper, i, nerr
+
+ if (iker /= IKER_FD) then
+ file=trim(CKER(iker+1))
+ else
+ file=trim(CTAP(itap+1))
+ endif
+ file=trim(adj_dir)//'/'//trim(file_prefix)//'.'//trim(file)// '.'//char(iwin+48)
+
+ ! IKER_WF
+ if (iker == IKER_WF) then
+ dt_adj_src(1:nlen) = synw(1:nlen)-dataw(1:nlen)
+ dt_chi = 0.5 * sum(synw(1:nlen)-dataw(1:nlen) ** 2) * dt
+
+ ! IKER_CC
+ else if (iker == IKER_CC) then
+ call compute_veloc_from_displ(synw,nlen,dt,synw_veloc)
+ ! we could have taken away - sign from the measurements, and put them here
+ ! just to be consistent with the FD measurements
+ ft_bar = 0.; fa_bar = 0.
+ ft_bar(1:nlen) = - synw_veloc(1:nlen) / ( sum(synw_veloc(1:nlen)**2) * dt )
+ fa_bar(1:nlen) = synw(1:nlen) / ( sum(synw(1:nlen)**2) * dt )
+ if (INCLUDE_ERROR) then
+ err_t = sigma_tshift_cc; err_a = sigma_dlnA_cc
+ else
+ err_t = 1; err_a = 1
+ endif
+ ! include CC measurements (and error)
+ dt_adj_src(1:nlen) = - (tshift_cc / err_t**2) * ft_bar(1:nlen)
+ amp_adj_src(1:nlen) = - (dlnA / err_a**2 ) * fa_bar(1:nlen)
+
+ dt_chi = 0.5 * (tshift_cc/err_t) ** 2
+ amp_chi = 0.5 * (dlnA / err_a) ** 2
+
+ ! IKER_FD
+ else if (iker == IKER_FD) then
+ df = 1./(dt*NPT)
+ ! define frequency-domain taper W(f)
+ wf_taper(:) = 0.
+ i_left = 1 !! a better choice??
+ do i = i_left, i_right
+ !wf_taper(i) = 1. ! boxcar
+ !wf_taper(i) = 1. - (2.0/nw)**2 * ((i-1) - nw/2.0)**2 ! welch
+ wf_taper(i) = 1. - cos(PI*(i-i_left)/(i_right-i_left))**ipwr_w ! cosine
+ enddo
+ ! normalize freq taper
+ ffac=2*df*sum(wf_taper(i_left:i_right))
+ wf_taper(i_left:i_right) = wf_taper(i_left:i_right) / ffac
+
+ ! water level of FD measurements based on average (LQY:can be smaller for mtm!!)
+ dtau_wtr = wtr * sum(abs(dtau_fdm(i_left:i_right)))/(i_right-i_left)
+ dlnA_wtr = wtr * sum(abs(dlnA_fdm(i_left:i_right)))/(i_right-i_left)
+
+ ! include errors in the W(f) taper
+ wp_taper(:) = 0.; wq_taper(:) = 0.
+ do i = i_left, i_right
+ wp_taper(i) = wf_taper(i); wq_taper(i) = wf_taper(i)
+ if (INCLUDE_ERROR) then
+ err_t = max(sigma_dtau_fdm(i),dtau_wtr)
+ err_a = max(sigma_dlnA_fdm(i),dlnA_wtr)
+ wp_taper(i) = wp_taper(i) / (err_t ** 2)
+ wq_taper(i) = wq_taper(i) / (err_A ** 2)
+ endif
+ enddo
+
+ d_bot_mtm = 0. ! should take only 1:i_right values, all positive numbers
+ v_bot_mtm = 0.
+
+ pwc_adj = cmplx(0.,0.); qwc_adj = cmplx(0., 0.)
+
+ ! compute bottom (f) of p_j(f) and q_j(f)
+ do ictaper = 1, ntaper
+
+ synwt(1:nlen) = synw(1:nlen) * tas(1:nlen,ictaper)
+ call compute_veloc_from_displ(synwt,nlen,dt,synwt_veloc)
+ ! FFT
+ csynwt(:) = 0.; csynwt_veloc = 0.
+ csynwt(1:nlen) = cmplx(synwt(1:nlen))
+ csynwt_veloc(1:nlen) = cmplx(synwt_veloc(1:nlen))
+
+ call fft(LNPT,csynwt,FORWARD_FFT,dble(dt))
+ call fft(LNPT,csynwt_veloc,FORWARD_FFT,dble(dt))
+
+ d_bot_mtm = d_bot_mtm + real(csynwt * conjg(csynwt))
+ v_bot_mtm = v_bot_mtm + real(csynwt_veloc * conjg(csynwt_veloc))
+
+ pwc_adj(:,ictaper) = csynwt_veloc
+ qwc_adj(:,ictaper) = - csynwt
+
+ enddo
+
+ ! add water level to bottom of p_j(f) and q_j(f)
+ wtr_v_bot = wtr * sum(abs(v_bot_mtm(1:i_right)))/i_right
+ wtr_d_bot = wtr * sum(abs(d_bot_mtm(1:i_right)))/i_right
+ do i = 1, i_right
+ if (v_bot_mtm(i) < wtr_v_bot) v_bot_mtm(i) = v_bot_mtm(i)+wtr_v_bot
+ if (d_bot_mtm(i) < wtr_d_bot) d_bot_mtm(i) = d_bot_mtm(i)+wtr_d_bot
+ enddo
+
+ ! initialize adjoint source
+ dt_adj_src = 0. ; amp_adj_src = 0.
+ do ictaper = 1, ntaper
+ dtau_pj = pwc_adj(:,ictaper)
+ dlnA_qj = qwc_adj(:,ictaper)
+ dtau_pj(1:i_right) = dtau_pj(1:i_right)/ v_bot_mtm(1:i_right) &
+ * cmplx(dtau_fdm(1:i_right), 0.) * cmplx(wp_taper(1:i_right),0.)
+ dlnA_qj(1:i_right) = dlnA_qj(1:i_right) / d_bot_mtm(1:i_right) &
+ * cmplx(dlnA_fdm(1:i_right), 0.) * cmplx(wp_taper(1:i_right),0.)
+
+ ! IFFT into the time domain
+ call fftinv(LNPT,dtau_pj,REVERSE_FFT,dble(dt),dtau_pj_t)
+ call fftinv(LNPT,dlnA_qj,REVERSE_FFT,dble(dt),dlnA_qj_t)
+
+ ! add up adjoint source
+ dt_adj_src(:) = dt_adj_src(:) + tas(:,ictaper) * dtau_pj_t(:)
+ amp_adj_src(:) = amp_adj_src(:) + tas(:,ictaper) * dlnA_qj_t(:)
+
+ enddo
+
+ dt_chi = 0.5 * 2.0 * df * sum( (dtau_fdm(1:i_right))**2 * wp_taper(1:i_right) )
+ amp_chi = 0.5 * 2.0 * df * sum( (dlnA_fdm(1:i_right))**2 * wq_taper(1:i_right) )
+
+ endif
+
+ if (DEBUG) then
+ print *, 'Writing adjoint file: ', trim(file)//'.adj.sac'
+ call wsac1(trim(file)//'.adj.sac',dt_adj_src,nlen,tstart,dt,nerr)
+ if (nerr > 0) stop 'Error writing windowed adjoint file'
+ endif
+ open(40,file=trim(file)//'.chi',status='unknown')
+ write(40,*) dt_chi
+ write(40,*) amp_chi
+ close(40)
+
+ end subroutine mt_adj_src
+
+
+! ================================================
+
+ subroutine adjust_adj_src(dt_adj_src,amp_adj_src,nlen,tstart,dt,&
+ dt_adj_src_win,amp_adj_src_win,npts_adj,b_adj,dt_adj)
+
+ real,dimension(:) :: dt_adj_src_win, amp_adj_src_win
+ real, dimension(:) :: dt_adj_src, amp_adj_src
+ real :: tstart, dt, b_adj, dt_adj
+ integer :: nlen, npts_adj
+
+ real,dimension(NPT) :: time_window
+ real :: fac
+ integer :: i
+
+ ! time-domain taper
+ time_window(1:nlen) = 0.
+ do i = 1,nlen
+ fac = 1. ! boxcar window
+ !fac = 1 - sfac2*((i-1) - dble(nlen1)/2.0)**2 ! welch window
+ !fac = 1. - cos(PI*(i-1)/(nlen-1))**ipwr_t ! cosine window
+ time_window(i) = fac
+ enddo
+ dt_adj_src = dt_adj_src * time_window
+ dt_adj_src = dt_adj_src * time_window
+
+ ! interpolate adjoint source
+ call interp_adj_src(dt_adj_src,nlen,tstart,dt, &
+ dt_adj_src_win,npts_adj,b_adj,dt_adj)
+
+ call interp_adj_src(amp_adj_src,nlen,tstart,dt, &
+ amp_adj_src_win,npts_adj,b_adj,dt_adj)
+
+
+ end subroutine adjust_adj_src
+
+
+! ================================================
+end module mtadj_sub
Added: seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub2.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub2.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub2.f90 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,359 @@
+module mtadj_sub2
+
+ use mtadj_constants
+ use mtadj_sub3
+
+ implicit none
+
+contains
+
+ subroutine compute_time_shift(synw,dataw,nlen,dt,ishift,tshift_cc)
+
+ real, intent(in),dimension(:) :: synw, dataw
+ integer, intent(in) :: nlen
+ real, intent(in) :: dt
+ real, intent(out) :: tshift_cc
+ integer, intent(out) :: ishift
+
+ real :: cc,cc_max
+ integer :: i,j
+
+ ! these choices will slide the entire windowed record past the other
+ if (nlen <= 0) stop 'Error compute_time_shift(): nlen has to be > 0'
+
+ ! cross-correlation (i should be between -nlen+1 to nlen+1)
+ ishift = 0
+ do i = -nlen, nlen, 1
+ cc = 0.
+ do j = 1, nlen
+ if((j+i) >= 1 .and. (j+i) <= nlen) cc = cc + synw(j) * dataw(j+i)
+ enddo
+ if(cc > cc_max) then
+ cc_max = cc
+ ishift = i
+ endif
+
+ enddo
+ tshift_cc = ishift*dt
+
+ end subroutine compute_time_shift
+
+ ! ===========================================================================
+
+ subroutine reconstruct_syn_cc(synw,nlen,dt,ishift,dlnA,synw_rc_cc)
+
+ real,intent(in),dimension(:) :: synw
+ integer, intent(in) :: nlen, ishift
+ real, intent(in) :: dt, dlnA
+ real,intent(out),dimension(:) :: synw_rc_cc
+
+ real,dimension(NPT) :: synw_shifted
+ integer :: i
+
+ ! shift synthetics
+ synw_shifted(:) = synw(:)
+ do i = 1, nlen
+ if ((i-ishift) >= 1 .and. (i-ishift) <= nlen ) synw_shifted(i) = synw(i-ishift)
+ enddo
+ ! fill the missing time window with the endpoint value
+ if (ishift > 0) synw_shifted(1:ishift) = synw(ishift+1)
+ if (ishift < 0) synw_shifted(nlen+ishift:nlen) = synw(nlen+ishift-1)
+
+ ! multiplying amplitude factor
+ synw_rc_cc(1:nlen) = synw_shifted(1:nlen) * exp( dlnA )
+
+ end subroutine reconstruct_syn_cc
+
+ ! ======================================================================
+
+ ! this subroutine computes the cross-correlation error
+ ! sigma_tshift_cc and sigma_dlnA_cc to take into account
+ ! the data and synthetics misfit
+
+ subroutine compute_cc_error(dataw,synw,nlen,dt,i_pmax,dlnA,&
+ sigma_tshift_cc,sigma_dlnA_cc, &
+ MIN_SIGMA_TSHIFT_CC,MIN_SIGMA_DLNA_CC)
+
+ real,dimension(:),intent(in) :: dataw, synw
+ integer, intent(in) :: nlen
+ real, intent(in) :: dt, dlnA, MIN_SIGMA_TSHIFT_CC, MIN_SIGMA_DLNA_CC
+ integer, intent(in) :: i_pmax
+ real, intent(out) :: sigma_tshift_cc, sigma_dlnA_cc
+
+ ! these two parameters seem to be arbitrary, but I guess it is better
+ ! to have some estimates than not at all!
+ real, parameter :: SIGMA_TSHIFT_CC_SCALE = 0.001
+ real, parameter :: SIGMA_DLNA_CC_SCALE = 2.0
+ real :: df,Tmax,misfit
+
+ df=1./(NPT*dt)
+ if (i_pmax <= 0) stop 'Check if the maximum f has been identified'
+ Tmax=1./(df*i_pmax)
+
+ misfit=sum((dataw(1:nlen)-synw(1:nlen))**2)/sum(dataw(1:nlen)**2)
+
+ sigma_tshift_cc = max(MIN_SIGMA_TSHIFT_CC,Tmax*misfit*SIGMA_TSHIFT_CC_SCALE)
+ sigma_dlnA_cc = max(MIN_SIGMA_dlnA_CC,abs(dlnA)*misfit*SIGMA_TSHIFT_CC_SCALE)
+
+ print *,
+
+ end subroutine compute_cc_error
+
+ ! ======================================================================
+
+ subroutine reconstruct_syn_fd(csynw,dtau_fdm,dlnA_fdm,i_right,&
+ synw_rc_fd,dt,nlen)
+
+ complex, dimension(:), intent(in) :: csynw
+ real, dimension(:), intent(in) :: dtau_fdm, dlnA_fdm
+ integer, intent(in) :: i_right, nlen
+ real,dimension(:), intent(out) :: synw_rc_fd
+ real, intent(in) :: dt
+
+ real :: df,fvec(NPT)
+ integer :: nf,j
+ complex*16,dimension(NPT) :: csynw_recon
+ real*8,dimension(NPT) :: synw_rc_fd_dp
+
+ df = 1./(dt*NPT); nf=floor(NPT/2.)+1
+ do j = 1, nf
+ fvec(j) = df*(j-1)
+ enddo
+ csynw_recon = cmplx(0.,0.)
+ csynw_recon(1:i_right)=csynw(1:i_right) * (1.+ dlnA_fdm(1:i_right))&
+ * exp(-CCI*2*pi*fvec(1:i_right)*dtau_fdm(1:i_right))
+
+ call fftinv(LNPT,csynw_recon,REVERSE_FFT,dble(dt),synw_rc_fd_dp)
+ synw_rc_fd(1:nlen) = synw_rc_fd_dp(1:nlen)
+ synw_rc_fd(nlen+1:NPT) = 0.
+
+
+ end subroutine reconstruct_syn_fd
+
+ ! ======================================================================
+
+ subroutine compute_dtau_dlnA(trans_fdm,dt,tshift_cc,dtau_fdm,dlnA_fdm,i_right)
+
+ complex, dimension(:),intent(in) :: trans_fdm
+ real,intent(in):: tshift_cc,dt
+ real,dimension(:), intent(out) :: dtau_fdm, dlnA_fdm
+ integer, intent(in) :: i_right
+
+ real :: df, smth, smth1, smth2
+ real, dimension(NPT) :: fvec, phi_wt, abs_wt
+ integer :: nf, j, i
+
+ df = 1./(dt*NPT); nf=floor(NPT/2.)+1
+ do j = 1, nf
+ fvec(j) = df*(j-1)
+ enddo
+ ! obtain amplitude and phase first
+ do i = 1, i_right
+ phi_wt(i) = atan2( aimag(trans_fdm(i)) , real(trans_fdm(i)) )
+ abs_wt(i) = abs(trans_fdm(i))
+ enddo
+
+ ! convert phase to travel time delay
+ dtau_fdm(1) = 0. ! tshift_cc is probably not a good choice
+ do i = 1, i_right
+
+ dlnA_fdm(i) = abs_wt(i) - 1.
+ !
+ if (i > 1) dtau_fdm(i) = (-1./(TWOPI*fvec(i))) * phi_wt(i) + tshift_cc
+
+ ! right now phi_wt [-180,180], adjust abrupt phase discontinuities [-90,90]
+ 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)
+ 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 (DEBUG) print *, '2 pi phase correction:', fvec(i), phi_wt(i) - phi_wt(i+1)
+ do j = i+1, i_right
+ phi_wt(j) = phi_wt(j) + TWOPI
+ enddo
+ endif
+ if(abs(smth2).lt.abs(smth).and.abs(smth2).lt.abs(smth1).and. abs(phi_wt(i) - phi_wt(i+1)) > PHASE_STEP)then
+ if (DEBUG) print *, '-2 pi phase correction:', fvec(i), phi_wt(i) - phi_wt(i+1)
+ do j = i+1, i_right
+ phi_wt(j) = phi_wt(j) - TWOPI
+ enddo
+ endif
+ endif
+! write(21,*) fvec(i), aimag(trans_fdm(i)),real(trans_fdm(i)),phi_wt(i)/PI*180
+ enddo
+
+ end subroutine compute_dtau_dlnA
+
+ ! ======================================================================
+ subroutine compute_mt_error(ntaper,dataw,synw,tas,&
+ nlen,dt,wtr_mtm,i_right,tshift_cc, &
+ dtau_fdm,dlnA_fdm, sigma_dtau_fdm,sigma_dlnA_fdm)
+
+ integer, intent(in) :: ntaper, nlen, i_right
+ real,dimension(:),intent(in) :: dataw, synw
+ real,dimension(:,:),intent(in) :: tas
+ real,intent(in) :: dt,tshift_cc,wtr_mtm
+ real,dimension(:),intent(in) :: dtau_fdm, dlnA_fdm
+ real,dimension(:),intent(out) :: sigma_dtau_fdm, sigma_dlnA_fdm
+
+ integer :: iom, i, nf, ictaper
+ complex,dimension(NPT) :: top_mtm, bot_mtm, trans_mtm
+ complex*16,dimension(NPT) :: cdatawt, csynwt
+ real,dimension(NPT,NMAX_TAPER) :: dtau_mtm, dlnA_mtm
+ real :: ampmax_bot, wtr_amp_bot, edt_ave, eabs2_ave, edt_iom, eabs2_iom
+ real,dimension(NPT)::err_dt,err_dlnA ,datawt,synwt
+
+
+ nf=floor(NPT/2.)+1
+
+ do iom = 1, ntaper
+
+ top_mtm(:) = cmplx(0.,0.)
+ bot_mtm(:) = cmplx(0.,0.)
+
+ do ictaper = 1, ntaper
+ if(ictaper.eq.iom) cycle
+
+ ! apply ictaper'th taper
+ datawt(1:nlen) = dataw(1:nlen) * tas(1:nlen,ictaper)
+ synwt(1:nlen) = synw(1:nlen) * tas(1:nlen,ictaper)
+
+ ! complex tapered series
+ cdatawt(:) = cmplx(0.,0.)
+ csynwt(:) = cmplx(0.,0.)
+ cdatawt(1:nlen) = cmplx(datawt(1:nlen),0)
+ csynwt(1:nlen) = cmplx(synwt(1:nlen),0)
+
+ ! apply f.t. to get complex spectra
+ call fft(LNPT,cdatawt,FORWARD_FFT,dble(dt))
+ call fft(LNPT,csynwt,FORWARD_FFT,dble(dt))
+
+ ! calculate top and bottom of Jacknife transfer function
+ do i = 1, nf
+ top_mtm(i) = top_mtm(i) + cdatawt(i) * conjg(csynwt(i))
+ bot_mtm(i) = bot_mtm(i) + csynwt(i) * conjg(csynwt(i))
+ enddo
+ enddo ! ictaper
+
+ ! water level
+ ampmax_bot = maxval(abs(bot_mtm(1:nf)))
+ wtr_amp_bot = ampmax_bot * wtr_mtm ** 2
+
+ ! calculate transfer function using water level
+ do i = 1, i_right
+ if(abs(bot_mtm(i)) >= abs(wtr_amp_bot)) trans_mtm(i) = top_mtm(i) / bot_mtm(i)
+ if(abs(bot_mtm(i)) < abs(wtr_amp_bot)) trans_mtm(i) = top_mtm(i) /(bot_mtm(i)+wtr_amp_bot)
+ enddo
+ call compute_dtau_dlnA(trans_mtm,dt,tshift_cc,dtau_mtm(:,iom),dlnA_mtm(:,iom),i_right)
+ enddo ! iom
+
+ err_dt = 0.
+ err_dlnA = 0.
+ do i = 1, i_right
+
+ edt_ave = 0.
+ eabs2_ave = 0.
+
+ do iom = 1, ntaper
+ edt_iom = ntaper*dtau_fdm(i) - (ntaper-1)*dtau_mtm(i,iom)
+ edt_ave = edt_ave + edt_iom
+
+ eabs2_iom = ntaper*dlnA_fdm(i) - (ntaper-1)*dlnA_mtm(i,iom)
+ eabs2_ave = eabs2_ave + eabs2_iom
+ enddo
+
+ edt_ave = edt_ave / (ntaper)
+ eabs2_ave = eabs2_ave / (ntaper)
+
+ do iom = 1, ntaper
+ err_dt(i) = err_dt(i) + (dtau_mtm(i,iom) - edt_ave)**2
+ err_dlnA(i) = err_dlnA(i)+ (dlnA_mtm(i,iom) - eabs2_ave)**2
+ enddo
+
+ err_dt(i) = sqrt( err_dt(i) / (ntaper * (ntaper-1) ) )
+ ! set the error bar for the first point corresponding to
+ ! static offset to be large, which makes no contribution to
+ ! the adjoint source
+ if (i == 1) err_dt(i) = LARGE_VAL
+ err_dlnA(i) = sqrt( err_dlnA(i) / (ntaper * (ntaper-1) ) )
+
+ enddo ! i_right
+ sigma_dtau_fdm(1:i_right) = err_dt(1:i_right)
+ sigma_dlnA_fdm(1:i_right) = err_dlnA(1:i_right)
+
+ end subroutine compute_mt_error
+
+ ! ======================================================================
+ subroutine compute_fd_error(npi,nlen,i_right,dt,dtau_fdm,dlnA_fdm, &
+ sigma_dtau_fdm,sigma_dlnA_fdm)
+
+ real,intent(in) :: npi
+ integer,intent(in) :: nlen, i_right
+ real,intent(in) :: dt
+ real,dimension(:),intent(in) :: dtau_fdm,dlnA_fdm
+ real,dimension(:),intent(inout) :: sigma_dtau_fdm, sigma_dlnA_fdm
+
+ real :: df
+ integer :: idf_new,i
+
+ df = 1./(dt*NPT) * npi
+ idf_new = floor(1./(nlen)/df)
+ do i=1, i_right-idf_new
+ sigma_dtau_fdm(i) = maxval(abs(dtau_fdm(i)-dtau_fdm(i:i+idf_new)))
+ sigma_dlnA_fdm(i) = maxval(abs(dlnA_fdm(i)-dlnA_fdm(i:i+idf_new)))
+ enddo
+ sigma_dtau_fdm(i_right-idf_new+1:i_right) = sigma_dtau_fdm(i_right-idf_new)
+ sigma_dlnA_fdm(i_right-idf_new+1:i_right) = sigma_dlnA_fdm(i_right-idf_new)
+
+ end subroutine compute_fd_error
+
+ ! ======================================================================
+
+ subroutine compute_veloc_from_displ(synw,nlen,dt,synw_veloc)
+
+ real, dimension(:), intent(in) :: synw
+ integer, intent(in) :: nlen
+ real, intent(in) :: dt
+ real,dimension(:),intent(out) :: synw_veloc
+
+ integer i
+
+ do i = 2, nlen-1
+ synw_veloc(i) = (synw(i+1) - synw(i-1)) / (2.0*dt)
+ enddo
+ synw_veloc(1) = (synw(2) - synw(1)) / dt
+ synw_veloc(nlen) = (synw(nlen) - synw(nlen-1)) /dt
+
+ end subroutine compute_veloc_from_displ
+
+ ! ======================================================================
+
+ subroutine interp_adj_src(dt_adj_src,nlen,tstart,dt, &
+ dt_adj_src_win,npts_adj,b_adj,dt_adj)
+
+ real,intent(in) :: tstart,dt,b_adj,dt_adj
+ real,intent(in),dimension(:) :: dt_adj_src
+ integer,intent(in) :: nlen,npts_adj
+ real,intent(out),dimension(:) :: dt_adj_src_win
+
+ integer :: i,ii
+ real :: time, tt
+
+ ! initialize to zero
+ dt_adj_src_win(1:npts_adj) = 0.
+
+ do i = 1, npts_adj
+ time = b_adj + (i-1) * dt_adj
+ if (time > tstart .and. time < tstart + (nlen-1)*dt) then
+ ii = floor((time-tstart)/dt) + 1
+ tt = time - ((ii-1)*dt + tstart)
+ dt_adj_src_win(i) = (dt_adj_src(ii+1)-dt_adj_src(ii)) * tt/dt + dt_adj_src(ii)
+ endif
+ enddo
+
+ end subroutine interp_adj_src
+
+ ! ======================================================================
+
+end module mtadj_sub2
Added: seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub3.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub3.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/mtadj_sub3.f90 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,526 @@
+! this module contains the following subroutine/functions:
+! subroutine fft(n,xi,zign,dt)
+! subroutine fftinv(npow,s,zign,dt,r)
+
+! subroutine staper(nt, fw, nev, v, ndim, a, w)
+! subroutine tsturm(nt,n,a,b,w,nev,r,ndim,ev,ipar)
+! subroutine root(u,el,elam,a,bb,w,n,ik)
+! subroutine costaper(ipoint, ndata, tas)
+! subroutine boxcar(ipoint, ndata, tas)
+
+module mtadj_sub3
+
+ implicit none
+
+! TOLERRANCE CONTROL
+ double precision, parameter :: TOL=1e-7
+
+contains
+
+!------------------------------------------------------------------
+ subroutine fft(n,xi,zzign,dt)
+! Fourier transform
+! This inputs AND outputs a complex function.
+! The convention is FFT --> e^(-iwt)
+! numerical factor for Plancherel theorem: planch_fac = dble(NPT * dt * dt)
+!------------------------------------------------------------------
+ complex*16, dimension(*) :: xi
+ integer :: n
+ double precision :: dt
+
+ double precision, parameter :: PI = 3.141592653589793d+00
+ complex*16 :: wk, hold, q
+ double precision :: m(25)
+ double precision :: zzign,zign,flx,v
+ integer :: lblock,k,fk,jh,ii,istart
+ integer :: l,iblock,nblock,i,lbhalf,j,lx
+
+ ! sign must be +1. or -1.
+ if(zzign >= 0.) then
+ zign = 1.
+ else
+ zign = -1.
+ endif
+
+ lx = 2**n
+ do 1 i=1,n
+ 1 m(i) = 2**(n-i)
+ do 4 l=1,n
+ nblock = 2**(l-1)
+ lblock = lx/nblock
+ lbhalf = lblock/2
+ k = 0
+ do 4 iblock=1,nblock
+ fk = k
+ flx = lx
+
+ v = zign*2.*PI*fk/flx ! Fourier convention
+
+ wk = cmplx(cos(v),-sin(v)) ! sign change to -sin(v) 17-Nov-2006
+ istart = lblock*(iblock-1)
+
+ do 2 i=1,lbhalf
+ j = istart+i
+ jh = j+lbhalf
+ q = xi(jh)*wk
+ xi(jh) = xi(j)-q
+ xi(j) = xi(j)+q
+ 2 continue
+
+ do 3 i=2,n
+ ii = i
+ if(k < m(i)) go to 4
+ 3 k = k-m(i)
+ 4 k = k+m(ii)
+ k = 0
+ do 7 j=1,lx
+ if(k < j) go to 5
+ hold = xi(j)
+ xi(j) = xi(k+1)
+ xi(k+1) = hold
+ 5 do 6 i=1,n
+ ii = i
+ if(k < m(i)) go to 7
+ 6 k = k-m(i)
+ 7 k = k+m(ii)
+
+ ! final steps deal with dt factors
+ if(zign > 0.) then ! FORWARD FFT
+ do i = 1,lx
+ xi(i) = xi(i)*dt ! multiplication by dt
+ enddo
+
+ else ! REVERSE FFT
+ flx = flx*dt
+ do i = 1,lx
+ xi(i) = xi(i)/flx ! division by dt
+ enddo
+ endif
+
+ end subroutine fft
+
+!------------------------------------------------------------------
+ subroutine fftinv(npow,s,zzign,dt,r)
+! inverse Fourier transform -- calls fft
+!------------------------------------------------------------------
+
+ !implicit real*8(a-h,o-z)
+ !dimension r(4096*4)
+ !complex s(4096*4)
+
+ complex*16, intent(in) :: s(*)
+ double precision, intent(out) :: r(*) ! note this is REAL
+
+ double precision :: dt,zzign,zign
+ integer :: npow, nsmp, nhalf, i
+
+ nsmp = 2**npow
+ nhalf = nsmp/2
+ call rspec(s,nhalf) ! re-structuring
+
+ zign=zzign
+ call fft(npow,s,zign,dt) ! Fourier transform
+
+ do i = 1,nsmp
+ r(i) = real(s(i)) ! REAL part
+ enddo
+
+ end subroutine fftinv
+
+!------------------------------------------------------------------
+ subroutine rspec(s,np2)
+!------------------------------------------------------------------
+
+ !implicit real*8(a-h,o-z)
+ !complex s(4096*4)
+
+ complex*16 :: s(*)
+ integer :: np2,n,n1,i
+
+ n = 2*np2
+ n1 = np2+1
+
+ s(n1) = 0.
+! s(1) = 0.
+ s(1) = cmplx( real(s(1)),0.)
+
+ do i = 1,np2
+ s(np2+i) = conjg(s(np2+2-i))
+ enddo
+
+ end subroutine rspec
+
+!------------------------------------------------------------------
+ subroutine staper(nt, fw, nev, v, ndim, a, w)
+!------------------------------------------------------------------
+!$$$$ calls tsturm, root
+! Slepian - Thomson multi-taper procedure
+! Slepian, D. 1978 Bell Sys Tech J v57 n5 1371-1430
+! Thomson, D. J. 1982 Proc IEEE v70 n9 1055-1096
+! nt the number of points in the series
+! fw the time-bandwidth product (number of Rayleigh bins)
+! nev the desired number of tapers
+! v the eigenvectors (tapers) are returned in v(.,nev)
+! a, w work arrays dimensioned at least nt long (nt+1, nt odd)
+! a(1..nev) contains bandwidth retention factors on output.
+! The tapers are the eigenvectors of the tridiagonal matrix sigma(i,j)
+! [see Slepian(1978) eq 14 and 25.] They are also the eigenvectors of
+! the Toeplitz matrix eq. 18. We solve the tridiagonal system in
+! tsturm for the tapers and use them in Slepians eq 18 to get the
+! bandwidth retention factors (i.e. the eigenvalues) Thomson's
+! normalisation is used with no attention to sign.
+ !implicit real*8(a-h,o-z)
+ !dimension a(*),w(*),v(ndim,*)
+ !parameter (pi=3.14159265358979d0,r2=1.414213562373095d0)
+
+ integer :: nt, nev, ndim
+ double precision :: fw
+ double precision :: v(ndim,*), a(*), w(*)
+
+ double precision, parameter :: PI = 3.141592653589793d+00
+ integer :: i,j,k,m
+ integer :: nxi, lh, lp1, neven, nodd, ntot, kk, kmax, nlow, nup
+ double precision :: r2,om,com,hn,asav,rbd,dc,sm,s,sn,vmax
+
+ !-------------------------
+
+ r2 = sqrt(2.)
+
+ if(nt < 2) return
+ nxi=mod(nt,2)
+ lh=(nt/2)+nxi
+ lp1=nt+1
+ om=2.*PI*fw/nt
+ com=cos(om)
+ hn=0.5*dble(lp1)
+ do 10 i=1,lh
+ a(i)=com*(i-hn)**2
+ 10 w(i)=0.5*dble(i*(nt-i))
+ if(nxi == 0) then
+ asav=a(lh)-w(lh)
+ a(lh)=a(lh)+w(lh)
+ rbd=1./(a(lh)+w(lh-1))
+ else
+ asav=w(lh-1)
+ rbd=1./(w(lh)+w(lh-1))
+ w(lh-1)=r2*w(lh-1)
+ endif
+ do 15 i=1,lh
+ a(i+lh)=w(i)*rbd
+ w(i)=a(i+lh)**2
+ 15 a(i)=a(i)*rbd
+ neven=max0((nev+1)/2,1)
+ nodd=nev-neven
+! Do the even tapers
+ call tsturm(nt,lh,a,a(lh+1),w,neven,v,ndim,w(lh+1),0)
+ do 20 i=1,neven
+ k=2*i-1
+ if(nxi == 1) v(lh,k)=r2*v(lh,k)
+ do 20 j=1,lh
+ 20 v(lp1-j,k)=v(j,k)
+ if(nodd <= 0) goto 34
+! Do the odd tapers
+ if(nxi == 0) then
+ a(lh)=asav*rbd
+ else
+ a(nt)=asav*rbd
+ w(lh-1)=asav*asav
+ endif
+ call tsturm(nt,lh-nxi,a,a(lh+1),w,nodd,v,ndim,w(lh+1),1)
+ do 30 i=1,nodd
+ k=2*i
+ if(nxi == 1) v(lh,k)=0.
+ do 30 j=1,lh
+ 30 v(lp1-j,k)=-v(j,k)
+ 34 ntot=neven+nodd
+! Calculate bandwidth retention parameters
+ dc=2.*com
+ sm=0.
+ s=sin(om)
+ w(1)=om/PI
+ w(2)=s/PI
+ do 35 j=3,nt
+ sn=dc*s-sm
+ sm=s
+ s=sn
+ 35 w(j)=s/(PI*(j-1))
+ do 55 m=1,ntot
+ vmax=abs(v(1,m))
+ kmax=1
+ do 40 kk=2,lh
+ if(abs(v(kk,m)) <= vmax) goto 40
+ kmax=kk
+ vmax=abs(v(kk,m))
+ 40 continue
+ a(m)=0.
+ nlow=kmax-1
+ do 45 j=1,nlow
+ 45 a(m)=a(m)+w(j+1)*v(nlow+1-j,m)
+ nup=nt-nlow
+ do 50 j=1,nup
+ 50 a(m)=a(m)+w(j)*v(nlow+j,m)
+ 55 a(m)=a(m)/v(kmax,m)
+ return
+
+ end subroutine staper
+
+!------------------------------------------------------------------
+ subroutine tsturm(nt,n,a,b,w,nev,r,ndim,ev,ipar)
+!------------------------------------------------------------------
+!$$$$ calls root
+! Uses bisection and Sturm counting to isolate the eigenvalues of the
+! symmetric tridiagonal matrix with main diagonal a(.) and sub/super
+! diagonal b(.). Newton's method is used to refine the eigenvalue in
+! subroutine root then direct recursion is used to get the eigenvector
+! as this is always stable. Note ipar=0 for even tapers =1 for odd
+! tapers
+ !implicit real*8(a-h,o-z)
+ !parameter (epsi=1.d-15,epsi1=5.d-15)
+ !dimension a(*),b(*),ev(*),w(*),r(ndim,*)
+
+ double precision, parameter :: epsi = 1.d-15, epsi1 = 5.d-15
+
+ double precision, dimension(ndim) :: a, b, w, ev
+ double precision, dimension(ndim,*) :: r
+ integer :: nt,n,ndim,nev,ipar
+
+ double precision, dimension(ndim) :: bb
+ double precision :: q,el,elam,u,umeps,x,ddot,rnorm
+ integer :: i,j,ik,iag,m,jk,jm1
+
+ !-------------------------
+
+ if(n <= 0.or.nev <= 0) return
+ umeps=1.-epsi
+ do 5 i=1,nev
+ 5 ev(i)=-1.
+ u=1.
+ do 1000 ik=1,nev
+ if(ik > 1) u=ev(ik-1)*umeps
+ el=min(ev(ik),u)
+ 10 elam=0.5*(u+el)
+ if(abs(u-el) <= epsi1) goto 35
+ iag=0
+ q=a(1)-elam
+ if(q >= 0.) iag=iag+1
+ do 15 i=2,n
+ if(q == 0.) x=abs(b(i-1))/epsi
+ if(q /= 0.) x=w(i-1)/q
+ q=a(i)-elam-x
+ if(q >= 0.) iag=iag+1
+ if(iag > nev) goto 20
+ 15 continue
+ if(iag >= ik) go to 20
+ u=elam
+ go to 10
+ 20 if(iag == ik) go to 30
+ m=ik+1
+ do 25 i=m,iag
+ 25 ev(i)=elam
+ el=elam
+ go to 10
+ 30 el=elam
+ call root(u,el,elam,a,b,w,n,ik)
+ 35 ev(ik)=elam
+ jk=2*ik+ipar-1
+ r(1,jk)=1.
+ r(2,jk)=-(a(1)-ev(ik))/b(1)
+ ddot=1.+r(2,jk)*r(2,jk)
+ jm1=2
+ do 45 j=3,n
+ r(j,jk)=-((a(jm1)-ev(ik))*r(jm1,jk)+b(j-2)*r(j-2,jk))/b(jm1)
+ ddot=ddot+r(j,jk)*r(j,jk)
+ 45 jm1=j
+ rnorm=sqrt(nt/(2.*ddot))
+ do 50 j=1,n
+ 50 r(j,jk)=r(j,jk)*rnorm
+ 1000 continue
+ return
+
+ end subroutine tsturm
+
+!------------------------------------------------------------------
+ subroutine root(u,el,elam,a,bb,w,n,ik)
+!------------------------------------------------------------------
+
+ !implicit real*8(a-h,o-z)
+ !parameter (epsi = 1.d-15, epsi1 = 5.d-15)
+ !dimension a(*),bb(*),w(*)
+
+ double precision, parameter :: epsi = 1.d-15, epsi1 = 5.d-15
+ double precision :: u,el,elam
+ double precision, dimension(*) :: a,bb,w
+ integer :: n,ik
+
+ double precision :: an,b,bm,bn,del,x
+ integer :: i,iag
+
+ !----------------------
+
+ 5 elam=0.5*(u+el)
+ 10 if(abs(u-el) <= 1.5*epsi1) return
+ an=a(1)-elam
+ b=0.
+ bn=-1./an
+ iag=0
+ if(an >= 0.) iag=iag+1
+ do 20 i=2,n
+ if(an == 0.) x=abs(bb(i-1))/epsi
+ if(an /= 0.) x=w(i-1)/an
+ an=a(i)-elam-x
+ if(an == 0.) an=epsi
+ bm=b
+ b=bn
+ bn=((a(i)-elam)*b-bm*x-1.)/an
+ if(an >= 0.) iag=iag+1
+ 20 continue
+ if(iag == ik) goto 25
+ u=elam
+ goto 30
+ 25 el=elam
+ 30 del=1./bn
+ if(abs(del) <= epsi1) del=sign(epsi1,del)
+ elam=elam-del
+ if(elam >= u.or.elam <= el) goto 5
+ goto 10
+
+ end subroutine root
+
+! ------------------------------------------------------------------
+! subroutine costaper(ipoint, ndata, tas)
+! ------------------------------------------------------------------
+ subroutine costaper(ipoint, ndata, tas)
+ implicit none
+
+ integer ipoint, ndata
+ double precision tas(ndata,*)
+ double precision sum, pi
+ integer i
+
+ pi = asin(1.0d0)*2
+ sum = 0.
+ do i =1,ipoint
+ tas(i,1) = 1 - cos( 2*pi*i/ipoint)
+ tas(i,1) = tas(i,1) / sqrt(1.5)
+ enddo
+ return
+ end subroutine costaper
+
+! ------------------------------------------------------------------
+! subroutine boxcar(ipoint, ndata, tas)
+! ------------------------------------------------------------------
+ subroutine boxcar(ipoint, ndata, tas)
+
+ integer ipoint, ndata
+ double precision tas(ndata,*)
+ integer i
+
+ do i =1,ipoint
+ tas(i,1) = 1.0
+ enddo
+ return
+ end subroutine boxcar
+
+! ----------------------------
+ subroutine rmean(dat,npts)
+! ---------------------------
+
+! remove mean from data
+
+ implicit none
+
+ integer npts
+ real dat(npts),sum
+ integer i
+
+
+ sum=0
+ do i = 1, npts
+ sum = sum + dat(i)
+ enddo
+ do i = 1, npts
+ dat(i) = dat(i) - sum/npts
+ enddo
+
+ return
+
+ end subroutine rmean
+
+! ----------------------------
+ subroutine rtrend(y, n)
+! ---------------------------
+
+! remove trend a*i + b
+! revised from lupei's c code
+
+ implicit none
+
+ integer n
+ real y(n)
+ integer i
+ real*8 a, b, a11, a12, a22, y1, y2
+
+ y1 = 0.
+ y2 = 0.
+ do i = 1, n
+ y1 = y1 + (i-1)*y(i)
+ y2 = y2 + y(i);
+ enddo
+ a12 = 0.5*n*(n-1)
+ a11 = a12*(2*n-1)/3.
+ a22 = n;
+ b = a11*a22-a12*a12
+ a = (a22*y1-a12*y2)/b
+ b = (a11*y2-a12*y1)/b
+ do i = 1, n
+ y(i) = y(i) - a * (i-1) - b
+ enddo
+
+ return
+ end subroutine rtrend
+
+! ---------------------------------------
+ subroutine taper (a, n, start, end, b)
+! ---------------------------------------
+ integer n
+ real start, end
+ real a(*), b(*)
+
+! taper a(1:n) to b(1:n)
+! cosine taper the beginning portion 'start' and end portion 'end' between [0,1]
+
+ real :: an,ang,xi,cs
+ integer :: i, m1,m2,m3,m4,m5
+
+ an = n
+ m1 = an*start+0.5
+ m2 = m1 + 1
+ if(m1 > 0) then
+ ang = 3.1415926 / float(m1)
+ do i=1, m1
+ xi = i
+ cs = (1.0-cos(xi*ang))/2.0
+ b(i) = a(i)*cs
+ enddo
+ endif
+
+ m3 = an*end+0.5
+ m5 = n-m3
+ m4 = m5 + 1
+
+ if(m3 > 0) then
+ ang = 3.1415926 / float (m3)
+ do i=m4,n
+ xi = i-n-1
+ cs = (1.0-cos(xi*ang))/2.0
+ b(i) = a(i)*cs
+ enddo
+ endif
+ do i=m2, m5
+ b(i) = a(i)
+ enddo
+
+ return
+ end subroutine taper
+
+ end module mtadj_sub3
Added: seismo/3D/ADJOINT_TOMO/mtadj/mtadj_variables.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/mtadj_variables.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/mtadj_variables.f90 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,74 @@
+module mtadj_variables
+
+ use mtadj_constants
+
+ ! ------------------------------------------------------------
+ ! parameter file
+ integer :: iker ! kernel type
+ integer :: itap ! taper type
+ character(len=150) :: meas_dir ! output dir for measurements
+ character(len=150) :: adj_dir ! output dir for adj src
+
+ ! preprocessing
+ logical :: BANDPASS ! band-pass option
+ real :: tshort, tlong ! band-pass range
+ real :: fstart, fend
+
+ ! measurements
+ real :: BEFORE_SHIFT ! tshift_cc max control
+ real :: wtr ! water level for estimating f_right through amp. of s(f)
+ ! multi-taper measurements
+ real :: npi
+ integer :: ntaper ! ntaper = 2*npi
+ real :: wtr_mtm ! water level for bottom of MTM
+ real :: MIN_DT_SIGMA,MIN_DlnA_SIGMA ! minimum dt,dlnA standard error
+
+ ! compute error ????
+
+ ! window selection after measurements
+ logical :: SELECT_WINDOW
+ ! selection criteria
+ integer :: ncycle_in_window
+ ! dtau_i < min(T_i/dt_fac,dt_max_scale*tshift), err_dt_i < T_i/err_fac
+ real :: dt_fac, dt_max_scale, err_fac
+ ! real :: after_quality,after_tshift ! ??????
+
+ ! write adjoint source
+ logical :: INCLUDE_ERROR, BANDPASS_ADJ
+ real :: b_adj, dt_adj
+ integer :: npts_adj !interpolation of adjoint source
+
+ ! -----------------------------------------------------
+ ! global parameters
+
+ ! data and syn array
+ real, dimension(NDIM) :: data, syn
+ real :: b, dt
+ integer :: npts
+
+ ! windowed data and syn
+ real, dimension(NPT) :: dataw, synw
+ integer :: nstart, nend, nlen
+
+ ! cross-correlation measurements
+ real :: tshift_cc, dlnA
+ real :: sigma_tshift_cc, sigma_dlnA_cc ! error estimates
+
+ ! freq-dependent measurements
+ real,dimension(NPT,NMAX_TAPER) :: tas
+
+ ! left and right of frequency range to output measurements
+ integer :: i_left, i_right, idf_fd
+ real :: f_left, f_right, df_fd ! freq-dep indep df spacing
+ ! maximum power index
+ integer :: i_pmax
+
+ ! complex transfer function from freq-dep measurements
+ complex, dimension(NPT) :: trans_fdm
+ real, dimension(NPT) :: dtau_fdm, dlnA_fdm ! dtau and dlnA
+ real, dimension(NPT) :: sigma_dtau_fdm, sigma_dlnA_fdm ! error bars
+
+ ! adjoint source
+ integer :: ipwr_t, ipwr_w
+
+end module mtadj_variables
Added: seismo/3D/ADJOINT_TOMO/mtadj/utils/compile
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/utils/compile (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/utils/compile 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1 @@
+gfortran -o rotate_adj_src rotate_adj_src.f90 -L/opt/sac/lib -lsacio
Property changes on: seismo/3D/ADJOINT_TOMO/mtadj/utils/compile
___________________________________________________________________
Name: svn:executable
+
Added: seismo/3D/ADJOINT_TOMO/mtadj/utils/readme
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/utils/readme (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/utils/readme 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,3 @@
+rotate_adj_src utility helps to rotate the adjoint sources
+calculated from the measurements to E, N, Z components which
+is ready to go directly into SEM adjoint code.
Added: seismo/3D/ADJOINT_TOMO/mtadj/utils/rotate_adj_src.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/utils/rotate_adj_src.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/utils/rotate_adj_src.f90 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,88 @@
+program rotate_adj_src
+
+ ! this program rotates the three component seismograms in the sac format
+ ! from (R, T) to (N, E)
+
+ implicit none
+
+ integer, parameter :: NDIM = 40000 ! from mt_constants.f90
+ real, parameter :: EPS = 1.0e-2 ! check constistency
+
+ character(len=150) :: bazch, zfile,tfile,rfile,efile,nfile
+ integer :: nptst, nptsr, npts, nerr
+ logical :: r_exist, t_exist, z_exist
+ real :: baz,t0t,dtt,t0r,dtr, t0,dt, costh, sinth
+ real, dimension(NDIM) :: rdata, tdata, edata, ndata
+
+
+ ! read in arguments
+ call getarg(1,bazch)
+ call getarg(2,zfile)
+ call getarg(3,tfile)
+ call getarg(4,rfile)
+ call getarg(5,efile)
+ call getarg(6,nfile)
+
+ if (trim(bazch)=='' .or. trim(zfile) =='' .or. trim(tfile)=='' .or. &
+ trim(rfile) == '' .or. trim(efile) =='' .or. trim(nfile) =='') then
+ stop 'rotate_adj_src baz(radians) zfile tfile rfile efile nfile'
+ endif
+
+ read(bazch,*) baz
+
+ ! check existence of Z, T, R files
+ inquire(file=trim(zfile),exist=z_exist)
+ inquire(file=trim(tfile),exist=t_exist)
+ inquire(file=trim(rfile),exist=r_exist)
+
+ ! only Z component -> t0, dt, npts
+ if (.not. t_exist .and. .not. r_exist) then
+ if (.not. z_exist) stop 'At least one file should exist: zfile, tfile, rfile'
+ call rsac1(zfile,tdata,npts,t0,dt,NDIM,nerr)
+ if (nerr > 0) stop 'Error reading Z file'
+ tdata = 0.
+ rdata = 0.
+ endif
+
+ ! read R/T components
+ if (t_exist) then
+ call rsac1(tfile,tdata,nptst,t0t,dtt,NDIM,nerr)
+ if (nerr > 0) stop 'Error reading T file'
+ if (.not. r_exist) rdata = 0.
+ endif
+ if (r_exist) then
+ call rsac1(rfile,rdata,nptsr,t0r,dtr,NDIM,nerr)
+ if (nerr > 0) stop 'Error reading R file'
+ if (.not. t_exist) tdata = 0.
+ endif
+ ! check consistency of T and R components -> t0,dt,npts
+ if (t_exist .and. r_exist) then
+ if (abs(t0t-t0r)>EPS .or. abs(dtt-dtr)>EPS .or. nptst /= nptsr) &
+ stop 'check t0 and npts'
+ endif
+ if (t_exist) then
+ t0 =t0t; dt = dtt; npts = nptst
+ else if (r_exist) then
+ t0 =t0r; dt = dtr; npts = nptsr
+ endif
+ if (NDIM < npts) stop 'Increase NDIM'
+
+ ! [N E]' = [ -costh sinth; -sinth costh ] [R T]'
+ costh = cos(baz)
+ sinth = sin(baz)
+ ndata(1:npts) = - costh * rdata(1:npts) + sinth * tdata(1:npts)
+ edata(1:npts) = - sinth * rdata(1:npts) - costh * tdata(1:npts)
+
+ ! write N/E adjoint source
+ call wsac1(nfile,ndata,npts,t0,dt,nerr)
+ if (nerr > 0) stop 'Error writing N file'
+ call wsac1(efile,edata,npts,t0,dt,nerr)
+ if (nerr > 0) stop 'Error writing E file'
+
+ if (.not. z_exist) then
+ tdata = 0.
+ call wsac1(zfile,tdata,npts,t0,dt,nerr)
+ if (nerr > 0) stop 'Error writing Z file'
+ endif
+
+end program rotate_adj_src
Added: seismo/3D/ADJOINT_TOMO/mtadj/utils/rotate_adj_src.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/mtadj/utils/rotate_adj_src.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/mtadj/utils/rotate_adj_src.pl 2009-02-20 00:01:36 UTC (rev 14107)
@@ -0,0 +1,72 @@
+#!/usr/bin/perl -w
+
+use File::Basename;
+use Getopt::Std;
+use lib '/opt/seismo/lib/perl';
+use CMT_TOOLS;
+use DELAZ5;
+
+# works only on BHT/R/E components
+
+#Reading input arguments:
+ at ARGV>0 || die("rotate_adj_src.pl -m CMT -s STATION -o OUTDIR -n [BH/LH] all_adj_files\n");
+
+if (!getopts('m:s:o:n:')) {die("Check arguments\n");}
+
+if ($opt_m) {$cmt = $opt_m;} else {$cmt = "CMTSOLUTION";}
+if ($opt_s) {$stafile = $opt_s;} else {$stafile = "STATIONS";}
+if ($opt_o) {$outdir = $opt_o;} else {$outdir = "ADJ_SRC";}
+if ($opt_n) {$name = $opt_n;} else {$name="BH";}
+$tname="${name}T"; $rname="${name}R"; $zname="${name}Z";
+$ename="${name}E"; $nname="${name}N";
+
+system("rm -f $outdir/*");
+if (not -f $cmt or not -f $stafile or not -d $outdir) {die("Check files/dirs\n");}
+
+($elat,$elon) = get_cmt_location($cmt);
+$bin=".";
+
+# fill up the hash table that stores the adj src existence information
+foreach $file (@ARGV) {
+ if (not -f $file) {die("Check file $file\n");}
+ ($basefile) = basename($file);
+ ($sta,undef,$cmp) = split(/\./,$basefile);
+ $adj{$sta}{$cmp} = $file;
+ if (defined $adj{$sta}{all}) {$adj{$sta}{all} ++ ;}
+ else {$adj{$sta}{all} = 1;}
+# print "$basefile -- $sta\n";
+}
+
+system("rm -f station.tmp"); $nstation = 0;
+print "Copying stations adjoint source file to: $outdir ...\n";
+foreach $sta (keys(%adj)) {
+ if (defined $adj{$sta}{$tname} or defined $adj{$sta}{$rname} or defined $adj{$sta}{$zname}) {
+ if (defined $adj{$sta}{$tname} ) {
+ $tcomp = $adj{$sta}{$tname}; $rcomp = $tcomp; $zcomp = $tcomp;
+ $ncomp = $rcomp; $ecomp = $tcomp; $zcomp=~s/$tname/$zname/;
+ $rcomp =~s/$tname/$rname/; $ncomp=~s/$tname/$nname/; $ecomp=~s/$tname/$ename/;}
+ elsif (defined $adj{$sta}{$rname}) {
+ $rcomp = $adj{$sta}{$rname}; $tcomp = $rcomp; $zcomp = $rcomp;
+ $ncomp = $rcomp; $ecomp = $rcomp; $zcomp=~s/$rname/$zname/;
+ $tcomp =~s/$rname/$tname/; $ncomp=~s/$rname/$nname/; $ecomp=~s/$rname/$ename/;}
+ else {
+ $zcomp = $adj{$sta}{$zname}; $tcomp = $zcomp; $rcomp = $zcomp;
+ $ncomp = $rcomp; $ecomp = $rcomp; $rcomp=~s/$zname/$rname/;
+ $tcomp =~s/$zname/$tname/; $ncomp=~s/$zname/$nname/; $ecomp=~s/$zname/$ename/;}
+ (undef,undef,$slat,$slon) = split(" ",`grep "$sta " $stafile | head -n 1`);
+ print "*** $sta ***\n";
+ `grep "$sta " $stafile | head -n 1 >> station.tmp`; $nstation ++ ;
+ # station/event locs in degrees (input i = 0) and output baz in radian
+ (undef,undef,undef,$baz) = delaz5($slat,$slon,$elat,$elon,0);
+ print "rotate_adj_src $baz $zcomp $tcomp $rcomp $ecomp $ncomp\n\n";
+ system("$bin/rotate_adj_src $baz $zcomp $tcomp $rcomp $ecomp $ncomp");
+ if ($? != 0) {die("Error: rotate_adj_src $baz $zcomp $tcomp $rcomp $ecomp $ncomp\n");}
+ system("cp -f $ecomp $ncomp $zcomp $outdir");
+ }
+}
+
+system("echo $nstation > STATIONS_ADJOINT; cat station.tmp >> STATIONS_ADJOINT; mv STATIONS_ADJOINT $outdir/");
+#system("mv STATIONS_ADJOINT $outdir/STATIONS_ADJOINT");
+
+system("sac2ascii.pl $outdir/*");
+
Property changes on: seismo/3D/ADJOINT_TOMO/mtadj/utils/rotate_adj_src.pl
___________________________________________________________________
Name: svn:executable
+
More information about the CIG-COMMITS
mailing list