[cig-commits] r18177 - in seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS: . adjoint_sources adjoint_sources/amplitude adjoint_sources/traveltime
bozdag at geodynamics.org
bozdag at geodynamics.org
Tue Apr 5 09:25:05 PDT 2011
Author: bozdag
Date: 2011-04-05 09:25:05 -0700 (Tue, 05 Apr 2011)
New Revision: 18177
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/README
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/Makefile
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/compile_cut
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/create_adjsrc_amplitude.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/drw_ascfile.h
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_ascfile_c.c
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_ascfile_f.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_fortran_wrapper.c
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/create_adjsrc_traveltime.f90
Removed:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/cut_velocity.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/Makefile
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/compile_cut
Log:
adds UTILS/adjoint_sources directory with subdirectories to create traveltime and amplitude adjoint sources
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/README
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/README (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/README 2011-04-05 16:25:05 UTC (rev 18177)
@@ -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_GLOBE/trunk/UTILS/adjoint_sources/amplitude/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/Makefile (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/Makefile 2011-04-05 16:25:05 UTC (rev 18177)
@@ -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_GLOBE/trunk/UTILS/adjoint_sources/amplitude/compile_cut
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/compile_cut (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/compile_cut 2011-04-05 16:25:05 UTC (rev 18177)
@@ -0,0 +1,16 @@
+# 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
+
+# obsolete...
+#gcc -c -o rw_ascfile_c.o rw_ascfile_c.c
+#gcc -c -o rw_fortran_wrapper.o rw_fortran_wrapper.c
+#ifort -O3 -o xcreate_adjsrc_amplitude create_adjsrc_amplitude.f90 rw_ascfile_f.f90 rw_ascfile_c.o rw_fortran_wrapper.o
Property changes on: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/compile_cut
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/create_adjsrc_amplitude.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/create_adjsrc_amplitude.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/create_adjsrc_amplitude.f90 2011-04-05 16:25:05 UTC (rev 18177)
@@ -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_GLOBE/trunk/UTILS/adjoint_sources/amplitude/drw_ascfile.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/drw_ascfile.h (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/drw_ascfile.h 2011-04-05 16:25:05 UTC (rev 18177)
@@ -0,0 +1,29 @@
+#ifndef _drw_ascfile_h
+#define _drw_ascfile_h
+
+void dread_ascfile(const char *ascfile,
+ double *t0, double *dt, int *n,
+ double *data);
+void dwrite_ascfile(const char *ascfile,
+ double t0, double dt, int n,
+ const double *data);
+
+void dread_ascfile_c(const char *ascfile,
+ double *t0, double *dt, int *n,
+ double *data);
+void dwrite_ascfile_c(const char *ascfile,
+ const double *t0, const double *dt, const int *n,
+ const double *data);
+
+void DREAD_ASCFILE_C(const char *ascfile,
+ double *t0, double *dt, int *n,
+ double *data);
+void DWRITE_ASCFILE_C(const char *ascfile,
+ const double *t0, const double *dt, const int *n,
+ const double *data);
+
+
+
+
+
+#endif
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_ascfile_c.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_ascfile_c.c (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_ascfile_c.c 2011-04-05 16:25:05 UTC (rev 18177)
@@ -0,0 +1,50 @@
+#include <stdio.h>
+#include "drw_ascfile.h"
+
+void dread_ascfile(const char *ascfile,
+ double *t0, double *dt, int *n,
+ double *data)
+
+{
+ FILE *fd;
+ double junk,junk1;
+ int i;
+
+ if ((fd = fopen(ascfile,"r")) == NULL) {
+ printf(" file %s cannot be opened\n",ascfile);
+ exit(1);
+ }
+ i = 0;
+ while ( fscanf(fd,"%lf %lf\n",&junk, data+i) != EOF ) {
+ if (i == 0) junk1 = junk;
+ if (i == 1) *dt = junk - junk1;
+ i++;}
+ *t0 = junk1;
+ *n = i;
+ if (fclose(fd) != 0) {
+ printf(" file %s cannot be closed\n",ascfile);
+ exit(1);}
+
+}
+
+void dwrite_ascfile(const char *ascfile,
+ double t0, double dt, int n,
+ const double *data)
+
+{
+ FILE *fd;
+ int i;
+
+ if ((fd = fopen(ascfile,"w")) == NULL) {
+ printf(" file %s cannot be opened to write\n",ascfile);
+ exit(1);
+ }
+ i = 0;
+ for (i=0; i<n; i++) {
+ fprintf(fd,"%14.7g %18.7g\n", t0+i*dt, data[i]);
+ }
+ if (fclose(fd) != 0) {
+ printf("file %s cannot be closed\n",ascfile);
+ exit(1);}
+
+}
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_ascfile_f.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_ascfile_f.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_ascfile_f.f90 2011-04-05 16:25:05 UTC (rev 18177)
@@ -0,0 +1,29 @@
+
+
+! these functions can be used directly in fortran
+
+subroutine dread_ascfile_f(name,t0,dt,n,data)
+
+ implicit none
+
+ character(len=*) :: name
+ real*8 :: t0, dt, data(*)
+ integer :: n
+
+ call dread_ascfile_c(trim(name)//char(0), t0, dt, n, data)
+
+end subroutine dread_ascfile_f
+
+
+subroutine dwrite_ascfile_f(name,t0,dt,n,data)
+
+ implicit none
+
+ character(len=*) :: name
+ real*8 :: t0, dt, data(*)
+ integer :: n
+
+ call dwrite_ascfile_c(trim(name)//char(0), t0, dt, n, data)
+
+end subroutine dwrite_ascfile_f
+
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_fortran_wrapper.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_fortran_wrapper.c (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/amplitude/rw_fortran_wrapper.c 2011-04-05 16:25:05 UTC (rev 18177)
@@ -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);}
Copied: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime (from rev 18176, seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/Makefile 2011-04-05 14:28:05 UTC (rev 18176)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/Makefile 2011-04-05 16:25:05 UTC (rev 18177)
@@ -10,20 +10,20 @@
#############################################################
-LIB_OBJS = cut_velocity.o rw_ascfile_c.o
+LIB_OBJS = create_adjsrc_traveltime.o rw_ascfile_c.o
# targets
-all: xcut_velocity
+all: xcreate_adjsrc_traveltime
-xcut_velocity: $(LIB_OBJS)
- ${F90} -Wall -o xcut_velocity $(LIB_OBJS)
+xcreate_adjsrc_traveltime: $(LIB_OBJS)
+ ${F90} -Wall -o xcreate_adjsrc_traveltime $(LIB_OBJS)
-cut_velocity.o: cut_velocity.f90
- ${F90} -Wall -c cut_velocity.f90
+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 xcut_velocity
+ rm -f *.o xcreate_adjsrc_traveltime
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/compile_cut
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/compile_cut 2011-04-05 14:28:05 UTC (rev 18176)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/compile_cut 2011-04-05 16:25:05 UTC (rev 18177)
@@ -1,10 +1,10 @@
# use make
-echo "xcut_velocity:"
+echo "xcreate_adjsrc_traveltime:"
echo
-echo "please,in this directory UTILS/cut_velocity"
+echo "please,in this directory UTILS/adjoint_sources/traveltime"
echo "add your compiler specifics to: Makefile"
echo
-echo "and run in this directory UTILS/cut_velocity:"
+echo "and run in this directory UTILS/adjoint_sources/traveltime:"
echo
echo "> make"
echo
@@ -13,4 +13,4 @@
# obsolete...
#gcc -c -o rw_ascfile_c.o rw_ascfile_c.c
#gcc -c -o rw_fortran_wrapper.o rw_fortran_wrapper.c
-#ifort -O3 -o xcut_velocity cut_velocity.f90 rw_ascfile_f.f90 rw_ascfile_c.o rw_fortran_wrapper.o
+#ifort -O3 -o xcreate_adjsrc_traveltime create_adjsrc_traveltime.f90 rw_ascfile_f.f90 rw_ascfile_c.o rw_fortran_wrapper.o
Copied: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/create_adjsrc_traveltime.f90 (from rev 18176, seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/cut_velocity.f90)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/create_adjsrc_traveltime.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/create_adjsrc_traveltime.f90 2011-04-05 16:25:05 UTC (rev 18177)
@@ -0,0 +1,173 @@
+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
+!
+! 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
+
+
Deleted: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/cut_velocity.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/cut_velocity/cut_velocity.f90 2011-04-05 14:28:05 UTC (rev 18176)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/adjoint_sources/traveltime/cut_velocity.f90 2011-04-05 16:25:05 UTC (rev 18177)
@@ -1,173 +0,0 @@
-program cut_velocity
-
-! 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
-!
-! call by: ./xcut_velocity 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*,' xcut_velocity 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 'cut_velocity 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: cut_velocity 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 'cut_velocity 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 'cut_velocity 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('t.txt',t0,dt,nstep,data(4,:))
- call dwrite_ascfile_c('r.txt',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('t-cut.txt',t0,dt,nstep,data(4,:))
- call dwrite_ascfile_c('r-cut.txt',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 cut_velocity
-
-
More information about the CIG-COMMITS
mailing list