[cig-commits] r18179 - in seismo/3D/SPECFEM3D/trunk: examples/homogeneous_halfspace utils utils/Cluster/lsf utils/adjoint_sources utils/adjoint_sources/amplitude utils/adjoint_sources/traveltime utils/lib
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Tue Apr 5 11:13:02 PDT 2011
Author: danielpeter
Date: 2011-04-05 11:13:02 -0700 (Tue, 05 Apr 2011)
New Revision: 18179
Added:
seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/
seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/README
seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/
seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/Makefile
seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/compile_cut
seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/create_adjsrc_amplitude.f90
seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/
seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/Makefile
seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/compile_cut
seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/create_adjsrc_traveltime.f90
seismo/3D/SPECFEM3D/trunk/utils/lib/rw_fortran_wrapper.c
Removed:
seismo/3D/SPECFEM3D/trunk/utils/cut_velocity/
Modified:
seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/README_kernel
seismo/3D/SPECFEM3D/trunk/utils/Cluster/lsf/go_mesher_solver_lsf_basin.kernel
seismo/3D/SPECFEM3D/trunk/utils/lib/drw_ascfile.h
seismo/3D/SPECFEM3D/trunk/utils/lib/rw_ascfile_c.c
seismo/3D/SPECFEM3D/trunk/utils/lib/rw_ascfile_f.f90
Log:
adds directory utils/adjoint_sources to create traveltime and amplitude adjoint sources
Modified: seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/README_kernel
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/README_kernel 2011-04-05 18:02:33 UTC (rev 18178)
+++ seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace/README_kernel 2011-04-05 18:13:02 UTC (rev 18179)
@@ -43,15 +43,15 @@
2. create adjoint source files:
- - compile the utility xcut_velocity:
- > cd /SPECFEM3D/UTILS/cut_velocity
+ - compile the utility xcreate_adjsrc_traveltime:
+ > cd /SPECFEM3D/utils/adjoint_sources/traveltime
> make
- specify which receiver station becomes an adjoint source,
e.g. using the seismograms from station 20, and create the
corresponding adjoint source files:
> cd /SPECFEM3D
- > ./UTILS/cut_velocity/xcut_velocity 10. 25. 0 in_out_files/OUTPUT_FILES/X20.DB.BX*.semd
+ > ./utils/adjoint_sources/traveltime/xcreate_adjsrc_traveltime 10. 25. 0 in_out_files/OUTPUT_FILES/X20.DB.BX*.semd
- make designated directory for adjoint sources:
> mkdir in_out_files/SEM/
Modified: seismo/3D/SPECFEM3D/trunk/utils/Cluster/lsf/go_mesher_solver_lsf_basin.kernel
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/Cluster/lsf/go_mesher_solver_lsf_basin.kernel 2011-04-05 18:02:33 UTC (rev 18178)
+++ seismo/3D/SPECFEM3D/trunk/utils/Cluster/lsf/go_mesher_solver_lsf_basin.kernel 2011-04-05 18:13:02 UTC (rev 18179)
@@ -36,7 +36,8 @@
mkdir -p $current_pwd/in_out_files/SEM
cd $current_pwd/in_out_files/SEM
collect_seismo_lsf_multi.pl $current_pwd/in_out_files/OUTPUT_FILES/lsf_machines $current_pwd/in_data_files/Par_file
-xcut_velocity 29.0 32.5 3 GSC*.semd
+#xcut_velocity 29.0 32.5 3 GSC*.semd
+xcreate_adjsrc_traveltime 29.0 32.5 3 GSC*.semd
rename .semd.adj .adj *.semd.ad
cd $current_pwd
##########################################################################
Added: seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/README
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/README (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/README 2011-04-05 18:13:02 UTC (rev 18179)
@@ -0,0 +1,21 @@
+------------------------------------------------------------
+README
+------------------------------------------------------------
+
+This folder contains routines to create adjoint sources for:
+
+- traveltime/ subdirectory: cross-correlation traveltime adjoint sources
+
+ this program cuts certain portion of the seismograms and converts them into
+ the adjoint sources for generating banana-dougnut kernels.
+ (see Tromp et al. (2005), eq. 45)
+
+- amplitude/ subdirectory: amplitude adjoint sources
+
+ this program cuts a certain portion of displacement seismograms and
+ converts them into adjoint sources for generating classical amplitude
+ kernels following Tromp et al. (2005) eq. 67.
+
+
+
+
Added: seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/Makefile
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/Makefile (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/Makefile 2011-04-05 18:13:02 UTC (rev 18179)
@@ -0,0 +1,29 @@
+# Makefile
+
+#############################################################
+## modify to match your compiler defaults
+
+# compilers
+F90 = gfortran
+CC = gcc
+
+#############################################################
+
+
+LIB_OBJS = create_adjsrc_amplitude.o rw_ascfile_c.o
+
+# targets
+all: xcreate_adjsrc_amplitude
+
+xcreate_adjsrc_amplitude: $(LIB_OBJS)
+ ${F90} -Wall -o xcreate_adjsrc_amplitude $(LIB_OBJS)
+
+create_adjsrc_amplitude.o: create_adjsrc_amplitude.f90
+ ${F90} -Wall -c create_adjsrc_amplitude.f90
+
+rw_ascfile_c.o: ../../lib/rw_ascfile_c.c
+ ${CC} -c -o rw_ascfile_c.o ../../lib/rw_ascfile_c.c
+
+
+clean:
+ rm -f *.o xcreate_adjsrc_amplitude
Added: seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/compile_cut
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/compile_cut (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/compile_cut 2011-04-05 18:13:02 UTC (rev 18179)
@@ -0,0 +1,12 @@
+# use make
+echo "xcreate_adjsrc_amplitude:"
+echo
+echo "please,in this directory UTILS/adjoint_sources/amplitude"
+echo "add your compiler specifics to: Makefile"
+echo
+echo "and run in this directory UTILS/adjoint_sources/amplitude:"
+echo
+echo "> make"
+echo
+echo
+
Property changes on: seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/compile_cut
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/create_adjsrc_amplitude.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/create_adjsrc_amplitude.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/amplitude/create_adjsrc_amplitude.f90 2011-04-05 18:13:02 UTC (rev 18179)
@@ -0,0 +1,173 @@
+program create_adjsrc_amplitude
+
+! this program cuts a certain portion of displacement seismograms and
+! converts them into adjoint sources for generating classical amplitude
+! kernels following Tromp et al. (2005) eq.67.
+! Modified from cut_velocity.f90 (Qinya Liu, Caltech, May 2007)
+! by Ebru, Princeton, March 2011.
+!
+! call by: ./xcreate_adjsrc_amplitude t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]
+!
+ implicit none
+
+ integer :: i, is, ie, nstep, j, itime ,ifile,ios, i1, i2, nstep_old
+ character(len=256) :: arg(100), file(100)
+ character(len=256) :: filename
+ integer,parameter :: NMAX = 30000
+ real*8, parameter :: EPS = 1.0d-17
+ real*8, parameter :: PI = 3.1415926d0
+ real*8 :: ts, te, data(5,NMAX), out(NMAX), adj(NMAX), tw(NMAX), norm
+ real*8 :: dt, t0, t0_old, dt_old, costh, sinth, th, baz
+ logical :: lrot
+
+ i = 1
+ lrot = .false.
+
+ ! reads in file arguments
+ do while (1 == 1)
+ call getarg(i,arg(i))
+ if (i < 6 .and. trim(arg(i)) == '') then
+ print*,'Usage: '
+ print*,' xcreate_adjsrc_amplitude t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ print*,'with'
+ print*,' t1: window start time'
+ print*,' t2: window end time'
+ print*,' ifile: 0 = adjoint source calculated for each seismogram component'
+ print*,' ifile: 1 = adjoint source given by East component only'
+ print*,' ifile: 2 = adjoint source given by North component'
+ print*,' ifile: 3 = adjoint source given by Z component'
+ print*,' ifile: 4 = adjoint source given by rotated transversal component (requires baz)'
+ print*,' ifile: 5 = adjoint source given by rotated radial component (requires baz)'
+ print*,' E/N/Z-ascii-files : displacement traces stored as ascii files'
+ print*,' [baz]: (optional) back-azimuth, requires ifile = 4 or ifile = 5'
+ stop 'create_adjsrc_amplitude t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ endif
+ if (trim(arg(i)) == '') exit
+ if (i == 1) then
+ read(arg(i),*,iostat=ios) ts
+ if (ios /= 0) stop 'Error reading ts'
+ else if (i == 2) then
+ read(arg(i),*,iostat=ios) te
+ if (ios /= 0) stop 'Error reading te'
+ else if (i == 3) then
+ read(arg(i),*) ifile
+ if (ios /= 0) stop 'Error reading ifile'
+ else if (i == 4 .or. i == 5 .or. i == 6) then
+ file(i-3) = trim(arg(i))
+ else if (i == 7) then
+ read(arg(i),*,iostat=ios) baz
+ if (ios /= 0) stop 'Error reading baz'
+ lrot = .true.
+ else if (i > 7) then
+ stop 'Error: create_adjsrc_amplitude t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ endif
+ i = i + 1
+ enddo
+
+ ! checks rotation baz and ifile parameter
+ i = i - 1
+ if (lrot) then
+ if (i /= 7) stop 'create_adjsrc_amplitude t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ if (ifile /= 4 .and. ifile /= 5) stop 'ifile = 4 or 5 when baz is present'
+ th = (baz - 180.0) / 180.0 * PI
+ costh = cos(th)
+ sinth = sin(th)
+ else
+ if (ifile > 3 .or. ifile < 0) stop 'Error ifile should be between 0 - 3 when baz is not present'
+ if (i /= 6) stop 'create_adjsrc_amplitude t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ endif
+
+ ! user output
+ print *, 'ifile = ', ifile, ' lrot = ', lrot
+ print *, ' '
+
+ ! reads seismograms (ascii format)
+ do i = 1, 3
+ filename = trim(file(i))
+ print *, 'reading asc file '//trim(filename)//' ...'
+ call dread_ascfile_c(trim(filename)//char(0),t0,dt,nstep,data(i,:))
+ if (nstep > NMAX) stop 'Change the data array range limit'
+ if (i == 1) then
+ t0_old = t0; dt_old = dt; nstep_old = nstep
+ else
+ if (i > 1 .and. abs(t0_old - t0) > EPS &
+ .and. abs(dt_old - dt) > EPS &
+ .and. nstep_old /= nstep) &
+ stop 'Error different t0, dt, nstep'
+ endif
+ enddo
+ print *, ' '
+ print *, 'start time:',t0
+ print *, 'time step:',dt
+ print *, 'number of steps:',nstep
+ print *, ' '
+
+ ! component rotation
+ if (lrot) then
+ data(4,:) = costh * data(1,:) - sinth * data(2,:)
+ data(5,:) = sinth * data(1,:) + costh * data(2,:)
+ call dwrite_ascfile_c(trim('t.txt')//char(0),t0,dt,nstep,data(4,:))
+ call dwrite_ascfile_c(trim('r.txt')//char(0),t0,dt,nstep,data(5,:))
+ i1 = 3; i2 = 5
+ else
+ i1 = 1; i2 = 3
+ endif
+
+ ! loops over seismogram components
+ do i = i1, i2
+ ! start and end index
+ is = (ts - t0) / dt + 1
+ ie = (te - t0) / dt + 1
+ if (is < 1 .or. ie <= is .or. ie > nstep) then
+ print *, 'Error in ts, te'; stop
+ endif
+
+ ! time window (parabola shaped)
+ tw(1:nstep) = 0.
+ if( i == i1 ) open(44,file='plot_time_window.txt',status='unknown')
+ do j = is, ie
+ tw(j) = 1 - (2 * (dble(j) - is)/(ie - is) - 1) ** 2
+ if( i == i1 ) write(44,*) j,tw(j)
+ enddo
+ if( i == i1 ) close(44)
+
+ ! displacement array
+ do itime = 1, nstep
+ out(itime) = data(i,itime)
+ enddo
+
+ ! normalization factor
+ norm = dt * sum( tw(1:nstep) * out(1:nstep) * out(1:nstep))
+ print *, 'i = ', i, 'norm = ', norm
+ if (ifile /= 0 .and. ifile /= i) norm = 0.0
+
+ ! adjoint source
+ if (abs(norm) > EPS) then
+ adj(1:nstep) = out(1:nstep) * tw(1:nstep) / norm
+ else
+ print *, 'norm < EPS for file '//trim(file(i))
+ adj(:) = 0.
+ endif
+ data(i,:) = adj(:)
+
+ enddo
+ print *, ' '
+
+ ! component rotation back to cartesian x-y-z
+ if (lrot) then
+ call dwrite_ascfile_c(trim('t-cut.txt')//char(0),t0,dt,nstep,data(4,:))
+ call dwrite_ascfile_c(trim('r-cut.txt')//char(0),t0,dt,nstep,data(5,:))
+ data(1,:) = costh * data(4,:) + sinth * data(5,:)
+ data(2,:) = -sinth * data(4,:) + costh * data(5,:)
+ endif
+
+ ! file output for component BHE/BHN/BHZ
+ do i = 1, 3
+ filename = trim(file(i))//'.adj'
+ print *, 'write to asc file '//trim(filename)
+ call dwrite_ascfile_c(trim(filename)//char(0),t0,dt,nstep,data(i,:))
+ enddo
+
+end program create_adjsrc_amplitude
+
+
Added: seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/Makefile
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/Makefile (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/Makefile 2011-04-05 18:13:02 UTC (rev 18179)
@@ -0,0 +1,29 @@
+# Makefile
+
+#############################################################
+## modify to match your compiler defaults
+
+# compilers
+F90 = gfortran
+CC = gcc
+
+#############################################################
+
+
+LIB_OBJS = create_adjsrc_traveltime.o rw_ascfile_c.o
+
+# targets
+all: xcreate_adjsrc_traveltime
+
+xcreate_adjsrc_traveltime: $(LIB_OBJS)
+ ${F90} -Wall -o xcreate_adjsrc_traveltime $(LIB_OBJS)
+
+create_adjsrc_traveltime.o: create_adjsrc_traveltime.f90
+ ${F90} -Wall -c create_adjsrc_traveltime.f90
+
+rw_ascfile_c.o: ../../lib/rw_ascfile_c.c
+ ${CC} -c -o rw_ascfile_c.o ../../lib/rw_ascfile_c.c
+
+
+clean:
+ rm -f *.o xcreate_adjsrc_traveltime
Added: seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/compile_cut
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/compile_cut (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/compile_cut 2011-04-05 18:13:02 UTC (rev 18179)
@@ -0,0 +1,12 @@
+# use make
+echo "xcreate_adjsrc_traveltime:"
+echo
+echo "please,in this directory UTILS/adjoint_sources/traveltime"
+echo "add your compiler specifics to: Makefile"
+echo
+echo "and run in this directory UTILS/adjoint_sources/traveltime:"
+echo
+echo "> make"
+echo
+echo
+
Property changes on: seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/compile_cut
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/create_adjsrc_traveltime.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/create_adjsrc_traveltime.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/adjoint_sources/traveltime/create_adjsrc_traveltime.f90 2011-04-05 18:13:02 UTC (rev 18179)
@@ -0,0 +1,174 @@
+program create_adjsrc_traveltime
+
+! this program cuts certain portion of the seismograms and converts them into
+! the adjoint sources for generating banana-dougnut kernels.
+! Qinya Liu, Caltech, May 2007
+! (renamed from cut_velocity.f90 (2011)
+!
+! call by: ./xcreate_adjsrc_traveltime t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]
+!
+ implicit none
+
+ integer :: i, is, ie, nstep, j, itime ,ifile,ios, i1, i2, nstep_old
+ character(len=256) :: arg(100), file(100)
+ character(len=256) :: filename
+ integer,parameter :: NMAX = 30000
+ real*8, parameter :: EPS = 1.0d-17
+ real*8, parameter :: PI = 3.1415926d0
+ real*8 :: ts, te, data(5,NMAX), out(NMAX), adj(NMAX), tw(NMAX), norm
+ real*8 :: dt, t0, t0_old, dt_old, costh, sinth, th, baz
+ logical :: lrot
+
+ i = 1
+ lrot = .false.
+
+ ! reads in file arguments
+ do while (1 == 1)
+ call getarg(i,arg(i))
+ if (i < 6 .and. trim(arg(i)) == '') then
+ print*,'Usage: '
+ print*,' xcreate_adjsrc_traveltime t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ print*,'with'
+ print*,' t1: window start time'
+ print*,' t2: window end time'
+ print*,' ifile: 0 = adjoint source calculated for each seismogram component'
+ print*,' ifile: 1 = adjoint source given by East component only'
+ print*,' ifile: 2 = adjoint source given by North component'
+ print*,' ifile: 3 = adjoint source given by Z component'
+ print*,' ifile: 4 = adjoint source given by rotated transversal component (requires baz)'
+ print*,' ifile: 5 = adjoint source given by rotated radial component (requires baz)'
+ print*,' E/N/Z-ascii-files : displacement traces stored as ascii files'
+ print*,' [baz]: (optional) back-azimuth, requires ifile = 4 or ifile = 5'
+ stop 'create_adjsrc_traveltime t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ endif
+ if (trim(arg(i)) == '') exit
+ if (i == 1) then
+ read(arg(i),*,iostat=ios) ts
+ if (ios /= 0) stop 'Error reading ts'
+ else if (i == 2) then
+ read(arg(i),*,iostat=ios) te
+ if (ios /= 0) stop 'Error reading te'
+ else if (i == 3) then
+ read(arg(i),*) ifile
+ if (ios /= 0) stop 'Error reading ifile'
+ else if (i == 4 .or. i == 5 .or. i == 6) then
+ file(i-3) = trim(arg(i))
+ else if (i == 7) then
+ read(arg(i),*,iostat=ios) baz
+ if (ios /= 0) stop 'Error reading baz'
+ lrot = .true.
+ else if (i > 7) then
+ stop 'Error: create_adjsrc_traveltime t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ endif
+ i = i + 1
+ enddo
+
+ ! checks rotation baz and ifile parameter
+ i = i - 1
+ if (lrot) then
+ if (i /= 7) stop 'create_adjsrc_traveltime t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ if (ifile /= 4 .and. ifile /= 5) stop 'ifile = 4 or 5 when baz is present'
+ th = (baz - 180.0) / 180.0 * PI
+ costh = cos(th)
+ sinth = sin(th)
+ else
+ if (ifile > 3 .or. ifile < 0) stop 'Error ifile should be between 0 - 3 when baz is not present'
+ if (i /= 6) stop 'create_adjsrc_traveltime t1 t2 ifile[0-5] E/N/Z-ascii-files [baz]'
+ endif
+
+ ! user output
+ print *, 'ifile = ', ifile, ' lrot = ', lrot
+ print *, ' '
+
+ ! reads seismograms (ascii format)
+ do i = 1, 3
+ filename = trim(file(i))
+ print *, 'reading asc file '//trim(filename)//' ...'
+ call dread_ascfile_c(trim(filename)//char(0),t0,dt,nstep,data(i,:))
+ if (nstep > NMAX) stop 'Change the data array range limit'
+ if (i == 1) then
+ t0_old = t0; dt_old = dt; nstep_old = nstep
+ else
+ if (i > 1 .and. abs(t0_old - t0) > EPS &
+ .and. abs(dt_old - dt) > EPS &
+ .and. nstep_old /= nstep) &
+ stop 'Error different t0, dt, nstep'
+ endif
+ enddo
+ print *, ' '
+ print *, 'start time:',t0
+ print *, 'time step:',dt
+ print *, 'number of steps:',nstep
+ print *, ' '
+
+ ! component rotation
+ if (lrot) then
+ data(4,:) = costh * data(1,:) - sinth * data(2,:)
+ data(5,:) = sinth * data(1,:) + costh * data(2,:)
+ call dwrite_ascfile_c(trim('t.txt')//char(0),t0,dt,nstep,data(4,:))
+ call dwrite_ascfile_c(trim('r.txt')//char(0),t0,dt,nstep,data(5,:))
+ i1 = 3; i2 = 5
+ else
+ i1 = 1; i2 = 3
+ endif
+
+ ! loops over seismogram components
+ do i = i1, i2
+ ! start and end index
+ is = (ts - t0) / dt + 1
+ ie = (te - t0) / dt + 1
+ if (is < 1 .or. ie <= is .or. ie > nstep) then
+ print *, 'Error in ts, te'; stop
+ endif
+
+ ! time window (parabola shaped)
+ tw(1:nstep) = 0.
+ if( i == i1 ) open(44,file='plot_time_window.txt',status='unknown')
+ do j = is, ie
+ tw(j) = 1 - (2 * (dble(j) - is)/(ie - is) - 1) ** 2
+ if( i == i1 ) write(44,*) j,tw(j)
+ enddo
+ if( i == i1 ) close(44)
+
+ ! calculates velocity (by finite-differences)
+ do itime = 2, nstep-1
+ out(itime) = (data(i,itime+1) - data(i,itime-1)) / (2 * dt)
+ enddo
+ out(1) = (data(i,2) - data(i,1)) / dt
+ out(nstep) = (data(i,nstep) - data(i,nstep-1)) /dt
+
+ ! normalization factor
+ norm = dt * sum( tw(1:nstep) * out(1:nstep) * out(1:nstep))
+ print *, 'i = ', i, 'norm = ', norm
+ if (ifile /= 0 .and. ifile /= i) norm = 0.0
+
+ ! adjoint source
+ if (abs(norm) > EPS) then
+ adj(1:nstep) = - out(1:nstep) * tw(1:nstep) / norm
+ else
+ print *, 'norm < EPS for file '//trim(file(i))
+ adj(:) = 0.
+ endif
+ data(i,:) = adj(:)
+
+ enddo
+ print *, ' '
+
+ ! component rotation back to cartesian x-y-z
+ if (lrot) then
+ call dwrite_ascfile_c(trim('t-cut.txt')//char(0),t0,dt,nstep,data(4,:))
+ call dwrite_ascfile_c(trim('r-cut.txt')//char(0),t0,dt,nstep,data(5,:))
+ data(1,:) = costh * data(4,:) + sinth * data(5,:)
+ data(2,:) = -sinth * data(4,:) + costh * data(5,:)
+ endif
+
+ ! file output for component BHE/BHN/BHZ
+ do i = 1, 3
+ filename = trim(file(i))//'.adj'
+ print *, 'write to asc file '//trim(filename)
+ call dwrite_ascfile_c(trim(filename)//char(0),t0,dt,nstep,data(i,:))
+ enddo
+
+end program create_adjsrc_traveltime
+
+
Modified: seismo/3D/SPECFEM3D/trunk/utils/lib/drw_ascfile.h
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/lib/drw_ascfile.h 2011-04-05 18:02:33 UTC (rev 18178)
+++ seismo/3D/SPECFEM3D/trunk/utils/lib/drw_ascfile.h 2011-04-05 18:13:02 UTC (rev 18179)
@@ -1,4 +1,3 @@
-// Qinya Liu, May 2007, Caltech
#ifndef _drw_ascfile_h
#define _drw_ascfile_h
@@ -6,7 +5,7 @@
double *t0, double *dt, int *n,
double *data);
void dwrite_ascfile(const char *ascfile,
- double *t0, double *dt, int *n,
+ double t0, double dt, int n,
const double *data);
void dread_ascfile_c(const char *ascfile,
Modified: seismo/3D/SPECFEM3D/trunk/utils/lib/rw_ascfile_c.c
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/lib/rw_ascfile_c.c 2011-04-05 18:02:33 UTC (rev 18178)
+++ seismo/3D/SPECFEM3D/trunk/utils/lib/rw_ascfile_c.c 2011-04-05 18:13:02 UTC (rev 18179)
@@ -5,8 +5,8 @@
void dread_ascfile_c_(const char *ascfile,
- double *t0, double *dt, int *n,
- double *data)
+ double *t0, double *dt, int *n,
+ double *data)
{
FILE *fd;
@@ -33,16 +33,16 @@
}
void dwrite_ascfile_c_(const char *ascfile,
- double *t0, double *dt, int *n,
- const double *data)
+ double *t0, double *dt, int *n,
+ const double *data)
{
FILE *fd;
int i;
double tmp,tmp0;
-
+
//printf("t0: %f ,dt: %f ,n: %i \n",*t0,*dt,*n);
-
+
if ((fd = fopen(ascfile,"w")) == NULL) {
printf(" file %s cannot be opened to write\n",ascfile);
exit(1);
@@ -50,7 +50,7 @@
i = 0;
tmp = *dt;
tmp0 = *t0;
- for (i=0; i< *n; i++) {
+ for (i=0; i< *n; i++) {
fprintf(fd,"%14.7g %18.7g\n", tmp0+i*tmp, data[i]);
}
if (fclose(fd) != 0) {
Modified: seismo/3D/SPECFEM3D/trunk/utils/lib/rw_ascfile_f.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/lib/rw_ascfile_f.f90 2011-04-05 18:02:33 UTC (rev 18178)
+++ seismo/3D/SPECFEM3D/trunk/utils/lib/rw_ascfile_f.f90 2011-04-05 18:13:02 UTC (rev 18179)
@@ -1,4 +1,3 @@
-! Qinya Liu, Caltech, May 2007
! these functions can be used directly in fortran
@@ -24,8 +23,6 @@
real*8 :: t0, dt, data(*)
integer :: n
- print *,'asc_f:',t0,dt,n
-
call dwrite_ascfile_c(trim(name)//char(0), t0, dt, n, data)
end subroutine dwrite_ascfile_f
Added: seismo/3D/SPECFEM3D/trunk/utils/lib/rw_fortran_wrapper.c
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/lib/rw_fortran_wrapper.c (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/lib/rw_fortran_wrapper.c 2011-04-05 18:13:02 UTC (rev 18179)
@@ -0,0 +1,33 @@
+#include "drw_ascfile.h"
+
+void dread_ascfile_c(const char *name, double *t0, double *dt, int * n, double *data) {
+
+ dread_ascfile(name,t0,dt,n,data);}
+
+void DREAD_ASCFILE_C(const char *name, double *t0, double *dt, int * n, double *data) {
+
+ dread_ascfile(name,t0,dt,n,data);}
+
+void dread_ascfile_c_(const char *name, double *t0, double *dt, int * n, double *data) {
+
+ dread_ascfile(name,t0,dt,n,data);}
+
+void DREAD_ASCFILE_C_(const char *name, double *t0, double *dt, int * n, double *data) {
+
+ dread_ascfile(name,t0,dt,n,data);}
+
+void dwrite_ascfile_c(const char *name, const double *t0, const double *dt, const int *n, const double *data) {
+
+ dwrite_ascfile(name,*t0,*dt,*n,data);}
+
+void DWRITE_ASCFILE_C(const char *name, const double *t0, const double *dt, const int *n, const double *data) {
+
+ dwrite_ascfile(name,*t0,*dt,*n,data);}
+
+void dwrite_ascfile_c_(const char *name, const double *t0, const double *dt, const int *n, const double *data) {
+
+ dwrite_ascfile(name,*t0,*dt,*n,data);}
+
+void DWRITE_ASCFILE_C_(const char *name, const double *t0, const double *dt, const int *n, const double *data) {
+
+ dwrite_ascfile(name,*t0,*dt,*n,data);}
More information about the CIG-COMMITS
mailing list