[cig-commits] r15479 - in seismo/3D/ADJOINT_TOMO/measure_adj: . PLOTS scripts_tomo
carltape at geodynamics.org
carltape at geodynamics.org
Mon Jul 27 10:49:08 PDT 2009
Author: carltape
Date: 2009-07-27 10:49:07 -0700 (Mon, 27 Jul 2009)
New Revision: 15479
Added:
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/9818433_T006_T030_MPM_CI_m16_mtm_win_adj.pdf
seismo/3D/ADJOINT_TOMO/measure_adj/ascii_rw.f90
seismo/3D/ADJOINT_TOMO/measure_adj/run_mt_measure_adj.csh
Modified:
seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR
seismo/3D/ADJOINT_TOMO/measure_adj/Makefile
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats_all.pl
seismo/3D/ADJOINT_TOMO/measure_adj/README
seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub2.f90
seismo/3D/ADJOINT_TOMO/measure_adj/rotate_adj_src.f90
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl
Log:
Modifications so that the package no longer relies on Caltech-specific SAC libraries; also, gfortran, which is free, is now the default compiler.
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/MEASUREMENT.PAR 2009-07-27 17:49:07 UTC (rev 15479)
@@ -1,7 +1,7 @@
OUTPUT_FILES
- 3 # is_mtm -- taper type: 0 none; 1 multi-taper; 2 cosine; 3 boxcar
+ 1 # is_mtm -- taper type: 0 none; 1 multi-taper; 2 cosine; 3 boxcar
0.020 2.50 # wtr, npi (ntaper = 2*npi)
- 2 # iker -- 0 waveform; 1 multi-taper; 2 cross-correlation;
+ 1 # iker -- 0 waveform; 1 multi-taper; 2 cross-correlation;
.false. # RUN_BANDPASS: use band-pass on records
30.000 6.000 # TLONG and TSHORT: band-pass periods for records
-0.585 0.0110 18200 # tstart, DT, npts: time vector for simulations
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/Makefile
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/Makefile 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/Makefile 2009-07-27 17:49:07 UTC (rev 15479)
@@ -1,20 +1,19 @@
F90 = gfortran
-F90_FLAGS = -O2
+F90_FLAGS = -O2 -fbounds-check
-# NEW VERSION WITH DEFAULT SAC LIBRARIES
-#SACLIBDIR=${SACHOME}/lib
-#LIB = -lsacio -lsac -ltau -lm -L/opt/seismo/lib -lDRWFiles -lf90recipes
+# new version with default sac libraries
+SACLIBDIR = ${SACHOME}/lib
+LIB = -lsacio -lsac
+#LIB = -L/opt/seismo/lib -lDRWFiles -lf90recipes -lDSacio -lDSacLib -lSacTools -lm
-# ORIGINAL VERSION
-LIB = -L/opt/seismo/lib -lDRWFiles -lf90recipes -lDSacio -lDSacLib -lSacTools -lm
+# NOTE: order matters if modules depend on other modules
+MOD = mt_constants mt_variables ascii_rw mt_sub2 mt_sub
-MOD = mt_constants mt_variables mt_sub2 mt_sub
-
SRC_DIR = .
MOD_DIR = mod
OBJ_DIR = obj
BIN_DIR = .
-MAIN = mt_measure_adj
+MAIN = mt_measure_adj rotate_adj_src
MOD_FLAG = module
MOD_OBJ = $(patsubst %,$(OBJ_DIR)/%.o,$(MOD))
@@ -24,7 +23,7 @@
all : mt_measure_adj rotate_adj_src
$(MAIN) : % : $(SRC_DIR)/%.f90 $(F90_OBJ) $(MOD_OBJ)
- $(F90) -o $(BIN_DIR)/$* $(F90_FLAGS) $(SRC_DIR)/$*.f90 -I$(MOD_DIR) -J$(MOD_DIR) $(OBJ) $(LIB)
+ $(F90) -o $(BIN_DIR)/$* $(F90_FLAGS) $(SRC_DIR)/$*.f90 -I$(MOD_DIR) -J$(MOD_DIR) $(OBJ) -L${SACLIBDIR} $(LIB)
$(F90_OBJ): $(OBJ_DIR)/%.o : $(SRC_DIR)/%.f90
$(F90) -o $@ $(F90_FLAGS) -c $(SRC_DIR)/$*.f90
@@ -32,9 +31,8 @@
$(MOD_OBJ): $(OBJ_DIR)/%.o : $(SRC_DIR)/%.f90
$(F90) -o $@ $(F90_FLAGS) -c $(SRC_DIR)/$*.f90 -I$(MOD_DIR) -J$(MOD_DIR)
-# Note call to libraries
-rotate_adj_src: rotate_adj_src.f90
- $(F90) -o rotate_adj_src rotate_adj_src.f90 -L/opt/seismo/lib -lDRWFiles
+#rotate_adj_src: rotate_adj_src.f90
+# $(F90) -o rotate_adj_src rotate_adj_src.f90 -L/opt/seismo/lib -lDRWFiles
.PHONY : clean
Added: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/9818433_T006_T030_MPM_CI_m16_mtm_win_adj.pdf
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/9818433_T006_T030_MPM_CI_m16_mtm_win_adj.pdf
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2009-07-27 17:49:07 UTC (rev 15479)
@@ -97,7 +97,7 @@
}
# width of seismograms depends on whether you plot adjoint sources
-$iplot_adj = 0; # plot adjoint sources
+$iplot_adj = 1; # plot adjoint sources
$iplot_win = 1; # plot windows
$iplot_CClabel = 1; # plot CC measurement labels for each window
$iplot_recon = 1; # plot reconstructed synthetics
@@ -128,14 +128,15 @@
$fontno = 1; # font number to use for text
# set GMT defaults
-`gmtset BASEMAP_TYPE plain PAPER_MEDIA letter ANNOT_FONT_SIZE = 12p LABEL_FONT_SIZE = 12p PAGE_ORIENTATION = portrait PLOT_DEGREE_FORMAT D`;
+`gmtset BASEMAP_TYPE plain PAPER_MEDIA letter MEASURE_UNIT inch ANNOT_FONT_SIZE = 12p LABEL_FONT_SIZE = 12p PAGE_ORIENTATION = portrait PLOT_DEGREE_FORMAT D`;
#$name = "${sta}_${net}_win_adj_${klab}";
$name = "${eid}_${Ttag}_${sta}_${net}_${smodel}_${klab}_win_adj";
$ps_file = "$name.ps";
$jpg_file = "$name.jpg";
-$saclst="/opt/seismo-util/bin/saclst";
+#$saclst="/opt/seismo-util/bin/saclst";
+$saclst="/opt/sac/bin/saclst";
#------------
$xmin = -122; $xmax = -114;
@@ -205,7 +206,7 @@
#$Zsyn = "${opt_s}/${Zlab}.semd.sac${suffix}";
# RECONSTRUCTED SYNTHETICS
-$suffix = ".recon.cc"; # USER MUST CHANGE
+$suffix = ".recon.cc.sac"; # USER MUST CHANGE
$Trecon = "${opt_c}/${Tlab}${suffix}";
$Rrecon = "${opt_c}/${Rlab}${suffix}";
$Zrecon = "${opt_c}/${Zlab}${suffix}";
@@ -391,6 +392,7 @@
$maxD = max($maxT,max($maxR,$maxZ));
$maxS = max($maxT2,max($maxR2,$maxZ2));
$max = max($maxD,$maxS);
+#print "\n-- $Tsyn -- $Tdat\n -- $Rsyn -- $Rdat\n -- $Zsyn -- $Zdat\n";
#print "\n Overall max values: $maxD (data), $maxS (syn) $max (overall)\n";
$bounds = "-R$tmin/$tmax/$Smin/$Smax";
@@ -453,7 +455,7 @@
#==========================================
# TRANSVERSE component: data, synthetics, and windows
-`pssac2 -X-2 -Y4.5 $proj $Tsyn $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O $tick >> $ps_file`; # synthetics
+`pssac2 $Tsyn -X-2 -Y4.5 $proj $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O $tick >> $ps_file`; # synthetics
if ($ntline) {
if ($iplot_win==1) {
for ($i=0; $i<$nT; $i++) {
@@ -480,14 +482,14 @@
}
}
}
-if ($Tflag==1) {`pssac2 $proj $Tdat $bounds -Ent-3 -M${sizeT} $black_pen -N -K -O >> $ps_file`;} # data
-`pssac2 $proj $Tsyn $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O >> $ps_file`; # synthetics
-if ($ntline && ($iplot_recon==1)) {`pssac2 $proj $Trecon $bounds -Ent-3 -M${sizeT3} $recon_pen -N -K -O >> $ps_file`;} # reconstructed
+if ($Tflag==1) {`pssac2 $Tdat $proj $bounds -Ent-3 -M${sizeT} $black_pen -N -K -O >> $ps_file`;} # data
+`pssac2 $Tsyn $proj $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O >> $ps_file`; # synthetics
+if ($ntline && ($iplot_recon==1)) {`pssac2 $Trecon $proj $bounds -Ent-3 -M${sizeT3} $recon_pen -N -K -O >> $ps_file`;} # reconstructed
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC T\" | pstext $proj $bounds -N -O -K >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strT\" | pstext $proj $bounds -N -O -K >> $ps_file`;
# RADIAL component: data, synthetics, and windows
-`pssac2 -Y$dY $proj $Rsyn $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O $tick2 >> $ps_file`;
+`pssac2 $Rsyn -Y$dY $proj $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O $tick2 >> $ps_file`;
if ($nrline) {
if ($iplot_win==1) {
for ($i=0; $i<$nR; $i++) {
@@ -510,14 +512,14 @@
}
}
}
-if ($Rflag==1) {`pssac2 $proj $Rdat $bounds -Ent-3 -M${sizeR} $black_pen -N -K -O >> $ps_file`;}
-`pssac2 $proj $Rsyn $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O >> $ps_file`;
-if ($nrline && ($iplot_recon==1)) {`pssac2 $proj $Rrecon $bounds -Ent-3 -M${sizeR3} $recon_pen -N -K -O >> $ps_file`;}
+if ($Rflag==1) {`pssac2 $Rdat $proj $bounds -Ent-3 -M${sizeR} $black_pen -N -K -O >> $ps_file`;}
+`pssac2 $Rsyn $proj $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O >> $ps_file`;
+if ($nrline && ($iplot_recon==1)) {`pssac2 $Rrecon $proj $bounds -Ent-3 -M${sizeR3} $recon_pen -N -K -O >> $ps_file`;}
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC R\" | pstext $proj $bounds -N -O -K >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strR\" | pstext $proj $bounds -N -O -K >> $ps_file`;
# VERTICAL component: data, synthetics, and windows
-`pssac2 -Y$dY $proj $Zsyn $bounds -Ent-3 -M${sizeZ2} $red_pen -N -K -O $tick2 >> $ps_file`;
+`pssac2 $Zsyn -Y$dY $proj $bounds -Ent-3 -M${sizeZ2} $red_pen -N -K -O $tick2 >> $ps_file`;
if ($nzline) {
if ($iplot_win==1) {
for ($i=0; $i<$nZ; $i++) {
@@ -540,9 +542,9 @@
}
}
}
-if ($Zflag==1) {`pssac2 $proj $Zdat $bounds -Ent-3 -M${sizeZ} $black_pen -N -K -O $tick2 >> $ps_file`;}
-`pssac2 $proj $Zsyn $bounds -Ent-3 -M${sizeZ2} $red_pen -N -O -K >> $ps_file`;
-if ($nzline && ($iplot_recon==1)) {`pssac2 $proj $Zrecon $bounds -Ent-3 -M${sizeZ3} $recon_pen -N -K -O >> $ps_file`;}
+if ($Zflag==1) {`pssac2 $Zdat $proj $bounds -Ent-3 -M${sizeZ} $black_pen -N -K -O $tick2 >> $ps_file`;}
+`pssac2 $Zsyn $proj $bounds -Ent-3 -M${sizeZ2} $red_pen -N -O -K >> $ps_file`;
+if ($nzline && ($iplot_recon==1)) {`pssac2 $Zrecon $proj $bounds -Ent-3 -M${sizeZ3} $recon_pen -N -K -O >> $ps_file`;}
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC Z\" | pstext $proj $bounds -N -K -O >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strZ\" | pstext $proj $bounds -N -O -K >> $ps_file`;
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj_all.pl 2009-07-27 17:49:07 UTC (rev 15479)
@@ -51,5 +51,5 @@
}
# sleep command to allow the last file to convert to PDF
-`sleep 2s`;
+`sleep 1s`;
#----------------
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats_all.pl 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats_all.pl 2009-07-27 17:49:07 UTC (rev 15479)
@@ -152,7 +152,7 @@
open(CSH,">$cshfile");
# set GMT defaults
-print CSH "gmtset MEASURE_UNIT inch TICK_LENGTH 5p BASEMAP_TYPE plain PAPER_MEDIA letter ANNOT_FONT_SIZE = 10p LABEL_FONT_SIZE = 12p PLOT_DEGREE_FORMAT D\n";
+print CSH "gmtset BASEMAP_TYPE plain PAPER_MEDIA letter MEASURE_UNIT inch TICK_LENGTH 5p ANNOT_FONT_SIZE = 10p LABEL_FONT_SIZE = 12p PLOT_DEGREE_FORMAT D\n";
#-----------------------------------------
# plot histograms
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/README 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/README 2009-07-27 17:49:07 UTC (rev 15479)
@@ -1,4 +1,21 @@
==============================
+Measurement and adjoint source code stored in the SVN repository of the Computational Infrastructure for Geodynamics (CIG).
+ /cig/seismo/3D/ADJOINT_TOMO/measure_adj/
+
+This package is complementary to FLEXWIN, the automated window-picking algorithm that is also stored at CIG.
+
+Last tested by Carl Tape on 7-27-09 using Ubuntu 9.04 Linux 64-bit with gfortran 4.3.3 and 64-bit SAC beta version Version 101.3.
+
+For a test to see if things are working, try this:
+ make clean
+ make
+ mt_measure_adj
+ csh -f run_mt_measure_adj.csh
+If all goes well, you should reproduce the figure
+ PLOTS/9818433_T006_T030_MPM_CI_m16_mtm_win_adj.pdf
+
+Please email Carl Tape (carltape at fas.harvard.edu) with questions, comments, or bug reports.
+==============================
Qinya Liu, 12/2/06
Carl Tape, 2/7/08
@@ -51,7 +68,7 @@
combine_2_adj_src.pl -- combines two sets of adjoint sources and creates a new STATIONS_ADJOINT file
combine_2_adj_src_all.pl -- calls combine_2_adj_src.pl repeatedly to make the final set of adjoint sources
-These must be COPIED into the main directory to use them.
+These must be copied into the main directory to use them.
Finally, to make summary histogram plots for the measurements for each event,
as well as a single histogram for all measurements in total, see
Added: seismo/3D/ADJOINT_TOMO/measure_adj/ascii_rw.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/ascii_rw.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/ascii_rw.f90 2009-07-27 17:49:07 UTC (rev 15479)
@@ -0,0 +1,105 @@
+module ascii_rw
+
+ ! This module replaces calls to previous libraries
+ ! for reading and writing ascii files.
+ ! All operations are in double precision.
+
+ implicit none
+
+contains
+
+ !-----------------------------------------------------------------------------
+
+ subroutine dwascii(fname,data,npt,b,dt)
+
+ character(len=*),intent(in) :: fname
+ integer, intent(in) :: npt
+ double precision, intent(in) :: data(*)
+ double precision, intent(in) :: b,dt
+ integer :: ios,i
+
+ open(9,file=trim(fname),iostat=ios,status='unknown')
+ if (ios /= 0) then
+ print *,'Error opening ascii file to write: ',trim(fname)
+ stop
+ endif
+ do i = 1,npt
+ write(9,*) b + dt*(i-1), data(i)
+ enddo
+ close(9)
+
+ end subroutine dwascii
+
+ !-----------------------------------------------------------------------------
+
+ subroutine drascii(fname,data,npt,b,dt)
+
+ character(len=*), intent(in) :: fname
+ integer, intent(out) :: npt
+ double precision, intent(out) :: b,dt
+ double precision :: data(*)
+ double precision :: dtemp,ttemp
+ double precision, dimension(:), allocatable :: ti
+ integer :: ios,i
+
+ ! get number of lines in the input file
+ npt=0
+ open(10,file=trim(fname),iostat=ios,status='old')
+ if (ios /= 0) then
+ print *,'Error opening ascii file to read: ',trim(fname)
+ stop
+ endif
+ do while (1==1)
+ read(10,*,iostat=ios) ttemp, dtemp
+ if (ios /= 0) exit
+ npt=npt+1
+ enddo
+ close(10)
+
+ allocate(ti(npt))
+ open(10,file=trim(fname),iostat=ios,status='old')
+ do i=1,npt
+ read(10,*,iostat=ios) ttemp, dtemp
+ ti(i) = ttemp
+ data(i) = dtemp
+ enddo
+ close(10)
+
+ ! assign b and dt -- assume time increment is constant
+ b = ti(1)
+ dt = (ti(npt) - ti(1)) / (npt-1)
+
+ deallocate(ti)
+
+ end subroutine drascii
+
+ !-----------------------------------------------------------------------------
+
+!!$! 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
+
+end module ascii_rw
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-07-27 17:49:07 UTC (rev 15479)
@@ -5,14 +5,16 @@
use mt_variables
use mt_constants
+ use ascii_rw ! dwascii()
+ use mt_sub2 ! fft(), fftinv()
use mt_sub ! mt_measure(), mt_adj()
- use mt_sub2 ! fft(), fftinv()
implicit none
character(len=150) :: out_dir,datafile,synfile,file_prefix,file_prefix0,file_prefix2,measure_file_prefix,adj_file_prefix
- integer :: num_files, num_meas, i, j, k, iker, iker0, ios, is_mtm, is_mtm0, npt1, npt2, npts, nn
+ integer :: num_files, num_meas, i, j, k, iker, iker0, &
+ ios, is_mtm, is_mtm0, npt1, npt2, npts, nn, nerr
double precision, dimension(NDIM) :: data, syn, adj_syn_all, tr_adj_src, am_adj_src, recon_cc_all, tseis_recon_cc
double precision :: t01, dt1, t02, dt2, t0, dt, tstart, tend, chi, tt, dtt, df
double precision, dimension(NCHI) :: window_chi
@@ -138,13 +140,18 @@
read(11,'(a)',iostat=ios) synfile
if (ios /= 0) stop 'Error reading synfile'
- ! read in data (SAC format)
- !call dread_ascfile_f(datafile,t01,dt1,npt1,data) ! ASCII format
- call dread_sacdata_f(datafile,t01,dt1,npt1,data)
+!!$ ! read in data (SAC format)
+!!$ !call dread_ascfile_f(datafile,t01,dt1,npt1,data) ! ASCII format
+!!$ call dread_sacdata_f(datafile,t01,dt1,npt1,data)
+!!$
+!!$ ! read in synthetics (SAC format)
+!!$ !call dread_ascfile_f(synfile,t02,dt2,npt2,syn) ! ASCII format
+!!$ call dread_sacdata_f(synfile,t02,dt2,npt2,syn)
- ! read in synthetics (SAC format)
- !call dread_ascfile_f(synfile,t02,dt2,npt2,syn) ! ASCII format
- call dread_sacdata_f(synfile,t02,dt2,npt2,syn)
+ ! read data and syn
+ call drsac1(datafile,data,npt1,t01,dt1)
+ call drsac1(synfile,syn,npt2,t02,dt2)
+ !print *, npt1,t01,dt1,NDIM ! check: double precision
if (max(npt1,npt2) > NDIM) stop 'Too many number of points in data or syn'
@@ -154,15 +161,17 @@
npts = min(npt1,npt2)
! apply bandpass filter to data and synthetics, if desired
+ ! http://www.iris.washington.edu/pipermail/sac-help/2008-March/000376.html
+ ! Access to the kidate, xapiir, and getfil is not simple and not
+ ! supported under the current state of the SAC code base.
+
+ t0 = t01
if(.not. RUN_BANDPASS) then
t0 = t01
else
- call interpolate_syn(syn,t02,dt,npt2,t01,dt,npts)
+ !call interpolate_syn(syn,t02,dt,npt2,t01,dt,npts)
t0 = t01
-
- ! filter the data and synthetics
- call xapiir(data,npts,'BU',TRBDNDW,APARM,IORD,'BP',fstart0,fend0,dt,PASSES)
- call xapiir(syn,npts,'BU',TRBDNDW,APARM,IORD,'BP',fstart0,fend0,dt,PASSES)
+ call bandpass(data,npts,dt,fstart0,fend0)
endif
!!$ ! figure out station name, network name, component name
@@ -365,11 +374,12 @@
! write out the adjoint source
! filter the adjoint source (higher frequencies could enter from the tapering operations)
- ! (See also time_window in mt_adj.f90, which tapers the windows.)
- call xapiir(adj_syn_all,npts,'BU',TRBDNDW,APARM,IORD,'BP',fstart0,fend0,dt,PASSES)
+ ! Note: time_window in mt_adj.f90 tapers the windows.
+ call bandpass(adj_syn_all,npts,dt,fstart0,fend0)
! cut and interpolate to match time-stepping for SEM
- ! NOTE: This can leave a non-zero value to start the record, which is NOT GOOD for the SEM simulation.
+ ! NOTE: This can leave a non-zero value to start the record,
+ ! which is NOT GOOD for the SEM simulation.
call interpolate_syn(adj_syn_all,t0,dt,npts,tt,dtt,nn)
! Taper the start of the adjoint source, since the cutting of the record
@@ -378,8 +388,10 @@
call taper_start(adj_syn_all,nn,itmax)
! we can have choices of outputing the adjoint source as ASCII or SAC format
- call dwrite_ascfile_f(trim(adj_file_prefix)//'.adj',tt,dtt,nn,adj_syn_all)
- !call dwrite_sacfile_f(trim(datafile),trim(adj_file_prefix)//'.adj',tt,nn,adj_syn_all)
+ call dwascii(trim(adj_file_prefix)//'.adj',adj_syn_all,nn,tt,dtt)
+ !call dwsac1(trim(adj_file_prefix)//'.adj.sac',adj_syn_all,nn,tt,dtt)
+!!$ call dwrite_ascfile_f(trim(adj_file_prefix)//'.adj',tt,dtt,nn,adj_syn_all)
+!!$ !call dwrite_sacfile_f(trim(datafile),trim(adj_file_prefix)//'.adj',tt,nn,adj_syn_all)
!----------------------------
! write out the CC-reconstructed data from synthetics
@@ -389,7 +401,9 @@
!call dwrite_ascfile_f(trim(file_prefix2)//'.recon.cc',tt,dtt,nn,recon_cc_all)
! write SAC file
- call dwrite_sacfile_f(trim(datafile),trim(file_prefix2)//'.recon.cc',t0,npts,recon_cc_all)
+ call dwsac1(trim(file_prefix2)//'.recon.cc.sac',recon_cc_all,npts,t0,dt)
+ if (nerr > 0) stop 'Error writing reconstructed CC file'
+ !call dwrite_sacfile_f(trim(datafile),trim(file_prefix2)//'.recon.cc',t0,npts,recon_cc_all)
enddo ! npairs
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-07-27 17:49:07 UTC (rev 15479)
@@ -3,6 +3,7 @@
use mt_constants
use mt_variables
use mt_sub2
+ use ascii_rw
implicit none
@@ -67,19 +68,30 @@
character(len=150) :: filename
logical :: output_logical, display_logical
- if ( tstart < t0 .or. tend > t0+(npts-1)*dt .or. tstart >= tend) stop 'Check tstart and tend'
+ !-------------------------------------------------------------
+ if ( tstart < t0 .or. tend > t0+(npts-1)*dt .or. tstart >= tend) then
+ print *, 'tstart, t0, tend, t0+(npts-1)*dt:'
+ print *, tstart, t0, tend, t0+(npts-1)*dt
+ stop 'Check tstart and tend'
+ endif
+
! LQY -- is this too small ???
wtr_mtm = 1.e-10
filename = trim(file_prefix)
if (DISPLAY_DETAILS) then
- call dwrite_sacfile_f(datafile,trim(file_prefix)//'_data.sac',t0,npts,data)
- call dwrite_sacfile_f(datafile,trim(file_prefix)//'_syn.sac',t0,npts,syn)
- !call dwrite_sacfile_f(datafile,trim(file_prefix)//'_diff_dms.sac',t0,npts,data-syn)
- !call dwrite_ascfile_f(trim(file_prefix)//'_data.txt',t0,dt,npts,data)
- !call dwrite_ascfile_f(trim(file_prefix)//'_syn.txt',t0,dt,npts,syn)
+ call dwascii(trim(file_prefix)//'_data',data,npts,t0,dt)
+ call dwascii(trim(file_prefix)//'_syn',syn,npts,t0,dt)
+ !call dwsac1(trim(file_prefix)//'_data.sac',data,npts,t0,dt)
+ !call dwsac1(trim(file_prefix)//'_syn.sac',syn,npts,t0,dt)
+
+!!$ call dwrite_sacfile_f(datafile,trim(file_prefix)//'_data.sac',t0,npts,data)
+!!$ call dwrite_sacfile_f(datafile,trim(file_prefix)//'_syn.sac',t0,npts,syn)
+!!$ !call dwrite_sacfile_f(datafile,trim(file_prefix)//'_diff_dms.sac',t0,npts,data-syn)
+!!$ !call dwrite_ascfile_f(trim(file_prefix)//'_data.txt',t0,dt,npts,data)
+!!$ !call dwrite_ascfile_f(trim(file_prefix)//'_syn.txt',t0,dt,npts,syn)
endif
!--------------------------------------------------------------------------
@@ -107,9 +119,11 @@
enddo
if (DISPLAY_DETAILS) then
- print *, ' NPTs = ', NPT, ' new nlen = ', nlen
- call dwrite_sacfile_f(datafile,trim(file_prefix)//'.obs.sac',tstart,nlen,dzr2)
- call dwrite_sacfile_f(datafile,trim(file_prefix)//'.syn.sac',tstart,nlen,dzr)
+ print *, ' NPTs = ', NPT, ' new nlen = ', nlen
+ call dwsac1(trim(file_prefix)//'.obs.sac',dzr2,nlen,tstart,dt)
+ call dwsac1(trim(file_prefix)//'.syn.sac',dzr,nlen,tstart,dt)
+!!$ call dwrite_sacfile_f(datafile,trim(file_prefix)//'.obs.sac',tstart,nlen,dzr2)
+!!$ call dwrite_sacfile_f(datafile,trim(file_prefix)//'.syn.sac',tstart,nlen,dzr)
endif
! save a copy of the windowed data
@@ -197,15 +211,15 @@
enddo
! create complex synthetic seismogram and shifted-data seismogram
- wseis_syn(:) = cmplx(0.,0.)
- wseis_dat(:) = cmplx(0.,0.)
- wseis_syn(1:nlen) = dzr0(1:nlen)
- wseis_dat(1:nlen) = dzr20(1:nlen)
+ wseis_syn = cmplx(0.,0.)
+ wseis_dat = cmplx(0.,0.)
+ !wseis_syn(1:nlen) = dzr0(1:nlen)
+ !wseis_dat(1:nlen) = dzr20(1:nlen)
+ wseis_syn(1:nlen) = cmplx(dzr0(1:nlen))
+ wseis_dat(1:nlen) = cmplx(dzr20(1:nlen))
- print *, 'computing each time window in the frequency domain, calling fft'
call fft(LNPT,wseis_syn,FORWARD_FFT,dt)
- call fft(LNPT,wseis_dat,FORWARD_FFT,dt)
- print *, 'fft successful'
+ call fft(LNPT,wseis_dat,FORWARD_FFT,dt)
! index of the freq of the max power in the windowed data
ampmax_unw = 0.
@@ -252,8 +266,12 @@
print *, ' i_right = ', i_right, ', stopping freq = ', sngl(i_right * df)
! write out power for each signal
- call dwrite_ascfile_f(trim(file_prefix)//'.obs.power',df,df,i_right,abs(wseis_dat(1:i_right)) )
- call dwrite_ascfile_f(trim(file_prefix)//'.syn.power',df,df,i_right,abs(wseis_syn(1:i_right)) )
+ call dwascii(trim(file_prefix)//'.obs.power',abs(wseis_dat(1:i_right)),i_right,df,df)
+ call dwascii(trim(file_prefix)//'.syn.power',abs(wseis_syn(1:i_right)),i_right,df,df)
+ !call dwsac1(trim(file_prefix)//'.obs.power.sac',abs(wseis_dat(1:i_right)),i_right,df,df)
+ !call dwsac1(trim(file_prefix)//'.syn.power.sac',abs(wseis_syn(1:i_right)),i_right,df,df)
+!!$ call dwrite_ascfile_f(trim(file_prefix)//'.obs.power',df,df,i_right,abs(wseis_dat(1:i_right)) )
+!!$ call dwrite_ascfile_f(trim(file_prefix)//'.syn.power',df,df,i_right,abs(wseis_syn(1:i_right)) )
endif
@@ -888,9 +906,106 @@
!==============================================================================
!==============================================================================
+
+ !----------------------------------------------------------------------
+
+ subroutine bandpass(x,n,delta_t,f1,f2)
+ ! modified from FLEXWIN subroutines on 26-July-2009
+ integer, intent(in) :: n
+ double precision, intent(inout), dimension(*) :: x
+ double precision, intent(in) :: delta_t,f1,f2
+ real, dimension(:), allocatable :: x_sngl
+
+ allocate(x_sngl(n))
+
+ x_sngl(1:n) = sngl(x(1:n))
+ ! delta_t_sngl = sngl(delta_t)
+
+ ! old version - uses old SacLib
+ ! does band-pass filter
+ !call xapiir(x_sngl,n,'BU',sngl(TRBDNDW),sngl(APARM),IORD,'BP',sngl(FSTART),sngl(FEND),delta_t_sngl,PASSES)
+
+ ! new version, uses subroutines in libsac.a
+ ! does band-pass filter
+ call xapiir(x_sngl,n,'BU',TRBDNDW,APARM,IORD,'BP',f1,f2,delta_t,PASSES)
+
+ x(1:n) = dble(x_sngl(1:n))
+
+ deallocate(x_sngl)
+
+ end subroutine bandpass
+
!-----------------------------------------------------------------------------
+ subroutine drsac1(datafile,data,npt1,b1,dt1)
+ ! read sac file and convert to double precision
+
+ character(len=*),intent(in) :: datafile
+ real, dimension(NDIM) :: data_sngl
+ double precision, dimension(NDIM), intent(out) :: data
+ integer :: npt1, nerr
+ real :: b1_sngl,dt1_sngl
+ double precision :: b1,dt1
+
+ ! read file as single precision
+ call rsac1(datafile,data_sngl,npt1,b1_sngl,dt1_sngl,NDIM,nerr)
+ if (nerr > 0) stop ' Error reading sac file'
+
+ ! return double precision quantities
+ b1 = dble(b1_sngl)
+ dt1 = dble(dt1_sngl)
+ data = dble(data_sngl)
+
+ end subroutine drsac1
+
+ !-----------------------------------------------------------------------------
+
+ subroutine dwsac1(datafile,data,npt1,b1,dt1)
+ ! convert to single precision, then write sac file
+ ! --> includes an option to add minmax values to sac file,
+ ! which are used in the plotting scripts
+
+ character(len=*),intent(in) :: datafile
+ integer, intent(in) :: npt1
+ double precision, dimension(npt1), intent(in) :: data
+ double precision, intent(in) :: b1,dt1
+ logical, parameter :: minmax_header = .true.
+
+ real, dimension(npt1) :: data_sngl,ti_sngl
+ real :: b1_sngl,dt1_sngl,xmin_sngl,xmax_sngl
+ integer :: nerr,i
+
+ ! convert to single precision
+ b1_sngl = real(b1)
+ dt1_sngl = real(dt1)
+ data_sngl = real(data)
+
+ if (minmax_header) then
+ ! get time vector
+ ti_sngl = 0.
+ do i = 1,npt1
+ ti_sngl(i) = b1_sngl + (i-1)*dt1_sngl
+ enddo
+
+ ! set minmax values in sac file
+ xmin_sngl = minval(data_sngl)
+ xmax_sngl = maxval(data_sngl)
+ call setfhv('depmin',xmin_sngl,nerr)
+ call setfhv('depmax',xmax_sngl,nerr)
+
+ ! write file with headers
+ call wsac0(datafile,ti_sngl,data_sngl,nerr)
+
+ else
+ call wsac1(datafile,data_sngl,npt1,b1_sngl,dt1_sngl,nerr)
+ endif
+ if (nerr > 0) stop ' Error writing sac file'
+
+ end subroutine dwsac1
+
+ !-----------------------------------------------------------------------------
+
subroutine mt_measure_select(nlen,tshift_xc,i_pmax_syn,dtau_w,dlnA_w,err_dt,err_dlnA,dt,i_left,i_right,fstart,fend,use_trace)
integer, intent(in) :: nlen, i_pmax_syn
@@ -1422,10 +1537,12 @@
tseis_recon_cc(:) = tseis_recon_cc_dt * exp( dlnA ) ! based on dlnA = ln(Aobs/Asyn)
if (OUTPUT_MEASUREMENT_FILES) then
- call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_cc.sac',tstart,nlen,tseis_recon_cc)
- call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_cc_dt.sac',tstart,nlen,tseis_recon_cc_dt)
- !call dwrite_ascfile_f(trim(filename)//'.recon_syn',tstart,dt,nlen,tseis_recon_cc)
- !call dwrite_ascfile_f(trim(filename)//'.recon_syn_dt',tstart,dt,nlen,tseis_recon_cc_dt)
+ call dwsac1(trim(filename)//'.recon_syn_cc.sac',tseis_recon_cc,nlen,tstart,dt)
+ call dwsac1(trim(filename)//'.recon_syn_cc_dt.sac',tseis_recon_cc_dt,nlen,tstart,dt)
+!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_cc.sac',tstart,nlen,tseis_recon_cc)
+!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_cc_dt.sac',tstart,nlen,tseis_recon_cc_dt)
+!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn',tstart,dt,nlen,tseis_recon_cc)
+!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn_dt',tstart,dt,nlen,tseis_recon_cc_dt)
endif
end subroutine reconstruct_syn_cc
@@ -1472,10 +1589,12 @@
call fftinv(LNPT,wseis_recon,REVERSE_FFT,dt,tseis_recon_dt)
if (OUTPUT_MEASUREMENT_FILES) then
- call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn.sac',tstart,nlen,tseis_recon)
- call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_dt.sac',tstart,nlen,tseis_recon_dt)
- !call dwrite_ascfile_f(trim(filename)//'.recon_syn',tstart,dt,nlen,tseis_recon)
- !call dwrite_ascfile_f(trim(filename)//'.recon_syn_dt',tstart,dt,nlen,tseis_recon_dt)
+ call dwsac1(trim(filename)//'.recon_syn.sac',tseis_recon,nlen,tstart,dt)
+ call dwsac1(trim(filename)//'.recon_syn_dt.sac',tseis_recon_dt,nlen,tstart,dt)
+!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn.sac',tstart,nlen,tseis_recon)
+!!$ call dwrite_sacfile_f(datafile,trim(filename)//'.recon_syn_dt.sac',tstart,nlen,tseis_recon_dt)
+!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn',tstart,dt,nlen,tseis_recon)
+!!$ !call dwrite_ascfile_f(trim(filename)//'.recon_syn_dt',tstart,dt,nlen,tseis_recon_dt)
endif
end subroutine reconstruct_syn
@@ -1579,28 +1698,50 @@
character(len=*),intent(out) :: ntw,sta,comp
integer :: nsec,msec,i
- integer :: klen,nerror
-
+ integer :: klen,nerr
+
! integer header variables
- call saclst_iheader_f(data_file,'nzyear', yr)
- call saclst_iheader_f(data_file,'nzjday', jda)
- call saclst_iheader_f(data_file,'nzhour', ho)
- call saclst_iheader_f(data_file,'nzmin', mi)
- call saclst_iheader_f(data_file,'nzsec', nsec)
- call saclst_iheader_f(data_file,'nzmsec', msec)
-
+ call getnhv('nzyear',yr,nerr)
+ call getnhv('nzjday',jda,nerr)
+ call getnhv('nzhour',ho,nerr)
+ call getnhv('nzhour',mi,nerr)
+ call getnhv('nzmin',nsec,nerr)
+ call getnhv('nzmsec',msec,nerr)
+
sec=nsec+msec/1000.0
- call saclst_kheader_f(data_file,'knetwk',ntw,klen)
- call saclst_kheader_f(data_file,'kstnm', sta,klen)
- call saclst_kheader_f(data_file,'kcmpnm',comp,klen)
-
- call dsaclst_fheader_f(data_file,'dist',dist)
- call dsaclst_fheader_f(data_file,'az', az)
- call dsaclst_fheader_f(data_file,'baz', baz)
- call dsaclst_fheader_f(data_file,'stlo',slon)
- call dsaclst_fheader_f(data_file,'stla',slat)
+ ! string headers
+ call getkhv('knetwk',ntw,nerr)
+ call getkhv('kstnm',sta,nerr)
+ call getkhv('kcmpnm',comp,nerr)
+ ! decimal headers
+ call getfhv('dist',dist,nerr)
+ call getfhv('az',az,nerr)
+ call getfhv('baz',baz,nerr)
+ call getfhv('stlo',slon,nerr)
+ call getfhv('stla',slat,nerr)
+
+!!$ ! integer header variables
+!!$ call saclst_iheader_f(data_file,'nzyear', yr)
+!!$ call saclst_iheader_f(data_file,'nzjday', jda)
+!!$ call saclst_iheader_f(data_file,'nzhour', ho)
+!!$ call saclst_iheader_f(data_file,'nzmin', mi)
+!!$ call saclst_iheader_f(data_file,'nzsec', nsec)
+!!$ call saclst_iheader_f(data_file,'nzmsec', msec)
+!!$
+!!$ sec=nsec+msec/1000.0
+!!$
+!!$ call saclst_kheader_f(data_file,'knetwk',ntw,klen)
+!!$ call saclst_kheader_f(data_file,'kstnm', sta,klen)
+!!$ call saclst_kheader_f(data_file,'kcmpnm',comp,klen)
+!!$
+!!$ call dsaclst_fheader_f(data_file,'dist',dist)
+!!$ call dsaclst_fheader_f(data_file,'az', az)
+!!$ call dsaclst_fheader_f(data_file,'baz', baz)
+!!$ call dsaclst_fheader_f(data_file,'stlo',slon)
+!!$ call dsaclst_fheader_f(data_file,'stla',slat)
+
end subroutine get_sacfile_header
end module mt_sub
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub2.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub2.f90 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub2.f90 2009-07-27 17:49:07 UTC (rev 15479)
@@ -8,7 +8,7 @@
contains
!------------------------------------------------------------------
- subroutine fft(n,xi,zign,dt)
+ subroutine fft(n,xi,zzign,dt)
! Fourier transform
! This inputs AND outputs a complex function.
! The convention is FFT --> e^(-iwt)
@@ -21,12 +21,12 @@
double precision, parameter :: PI = 3.141592653589793d+00
complex*16 :: wk, hold, q
double precision :: m(25)
- double precision :: zign,flx,v
+ 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(zign >= 0.) then
+ if(zzign >= 0.) then
zign = 1.
else
zign = -1.
@@ -90,7 +90,7 @@
end subroutine fft
!------------------------------------------------------------------
- subroutine fftinv(npow,s,zign,dt,r)
+ subroutine fftinv(npow,s,zzign,dt,r)
! inverse Fourier transform -- calls fft
!------------------------------------------------------------------
@@ -101,13 +101,14 @@
complex*16, intent(in) :: s(*)
double precision, intent(out) :: r(*) ! note this is REAL
- double precision :: dt,zign
+ 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
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/rotate_adj_src.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/rotate_adj_src.f90 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/rotate_adj_src.f90 2009-07-27 17:49:07 UTC (rev 15479)
@@ -1,5 +1,7 @@
program rotate_adj_src
+ use ascii_rw ! ascii read and write module
+
implicit none
character(len=150) :: bazch, zfile,tfile,rfile,efile,nfile
@@ -11,7 +13,6 @@
integer, parameter :: NDIM = 40000 ! check mt_constants.f90
double precision, dimension(NDIM) :: rdata, tdata, edata, ndata
-
call getarg(1,bazch)
call getarg(2,zfile)
call getarg(3,tfile)
@@ -32,17 +33,21 @@
if (.not. t_exist .and. .not. r_exist) then
if (.not. z_exist) stop 'At least one file should exist: zfile, tfile, rfile'
- call dread_ascfile_f(zfile,t0,dt,npts,tdata)
+
+ call drascii(zfile,tdata,npts,t0,dt)
+ !call dread_ascfile_f(zfile,t0,dt,npts,tdata)
tdata = 0.
rdata = 0.
endif
if (t_exist) then
- call dread_ascfile_f(tfile,t0t,dtt,nptst,tdata)
+ call drascii(tfile,tdata,nptst,t0t,dtt)
+ !call dread_ascfile_f(tfile,t0t,dtt,nptst,tdata)
if (.not. r_exist) rdata = 0.
endif
if (r_exist) then
- call dread_ascfile_f(rfile,t0r,dtr,nptsr,rdata)
+ call drascii(rfile,rdata,nptsr,t0r,dtr)
+ !call dread_ascfile_f(rfile,t0r,dtr,nptsr,rdata)
if (.not. t_exist) tdata = 0.
endif
if (t_exist .and. r_exist) then
@@ -62,13 +67,15 @@
edata(1:npts) = -costh * tdata(1:npts) - sinth * rdata(1:npts)
ndata(1:npts) = sinth * tdata(1:npts) - costh * rdata(1:npts)
- call dwrite_ascfile_f(efile,t0,dt,npts,edata)
- call dwrite_ascfile_f(nfile,t0,dt,npts,ndata)
-
+ call dwascii(efile,edata,npts,t0,dt)
+ call dwascii(nfile,ndata,npts,t0,dt)
+ !call dwrite_ascfile_f(efile,t0,dt,npts,edata)
+ !call dwrite_ascfile_f(nfile,t0,dt,npts,ndata)
if (.not. z_exist) then
tdata = 0.
- call dwrite_ascfile_f(zfile,t0,dt,npts,tdata)
+ call dwascii(zfile,tdata,npts,t0,dt)
+ !call dwrite_ascfile_f(zfile,t0,dt,npts,tdata)
endif
end program rotate_adj_src
Added: seismo/3D/ADJOINT_TOMO/measure_adj/run_mt_measure_adj.csh
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/run_mt_measure_adj.csh (rev 0)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/run_mt_measure_adj.csh 2009-07-27 17:49:07 UTC (rev 15479)
@@ -0,0 +1,13 @@
+\rm -rf PLOTS/ADJOINT_SOURCES
+ mkdir PLOTS/ADJOINT_SOURCES
+cp OUTPUT_FILES/*adj PLOTS/ADJOINT_SOURCES
+cp OUTPUT_FILES/*recon.cc* PLOTS/RECON
+cp window_chi PLOTS
+\rm -rf ADJOINT_SOURCES
+ mkdir ADJOINT_SOURCES
+prepare_adj_src.pl -m CMTSOLUTION_9818433 -s PLOTS/STATIONS_TOMO -o ADJOINT_SOURCES OUTPUT_FILES/*adj
+cp STATIONS_ADJOINT ADJOINT_SOURCES
+\mv STATIONS_ADJOINT PLOTS
+cd PLOTS
+plot_win_adj_all.pl -l -10/200 -m ../CMTSOLUTION_9818433 -b 0 -k 1 -a STATIONS_ADJOINT -d DATA -s SYN -c RECON -w MEASUREMENT.WINDOWS -i m16 -j 6/30
+cd -
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl 2009-07-27 17:49:07 UTC (rev 15479)
@@ -4,15 +4,15 @@
#
# run_mt_cc_plot.pl
# Carl Tape
-# 17-June-2008
+# 27-July-2009
#
-# This runs mt_measure_adj.f90 first in cross-correlation mode, and then in cross-correlation mode.
-# The purpose is to have a directo comparison between the two types of adjoint sources.
-# In principle, we do not need both for the adjoint tomography inversion.
+# This runs mt_measure_adj.f90 first in cross-correlation mode,
+# and then in multitaper mode. The purpose is to have
+# a direct comparison between the two types of adjoint sources.
+# We do not need both for the adjoint tomography inversion.
#
# EXAMPLE:
-# run_mt_cc_plot.pl m00 0 6/40 -10/180 -0.6/0.011/18200 1 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.1 (DATA_TEST/SYN_TEST)
-# run_mt_cc_plot.pl m00 0 6/40 -10/180 -0.6/0.011/18200 0 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.1 (DATA/SYN)
+# run_mt_cc_plot.pl m16 0 6/30 -10/200 -0.585/0.011/18200 0 1/1/1 2.0/2.5/3.5/1.5 0.7/0.7/5.0/1.0 1.0/0.5 (DATA/SYN)
#
#==========================================================
@@ -62,7 +62,8 @@
open(CSH,">$cshfile");
# empty the output files directory and compile the measurement code
-print CSH "\\rm -rf OUTPUT_FILES ; mkdir OUTPUT_FILES ; make clean\n make\n";
+# (this is done in prepare_mt_measure_adj.pl)
+#print CSH "\\rm -rf OUTPUT_FILES ; mkdir OUTPUT_FILES ; make clean\n make\n";
# create MEASUREMENT.PAR file
$wtr = 0.02; $npi = 2.5; # "default" values
@@ -70,9 +71,6 @@
#---------------------------------------------
-# empty the output files directory and compile the measurement code
-print CSH "\\rm -rf OUTPUT_FILES ; mkdir OUTPUT_FILES ; make clean\n make\n";
-
# run the measurement code
print CSH "mt_measure_adj\n";
@@ -81,7 +79,7 @@
print CSH "cp OUTPUT_FILES/*adj ${plot_adj_dir}\n";
# copy reconstructed records into PLOTS
- print CSH "cp OUTPUT_FILES/*recon.cc ${plot_recon_dir}\n";
+ print CSH "cp OUTPUT_FILES/*recon.cc* ${plot_recon_dir}\n";
#////////////////////////////////////////////////////////
@@ -89,8 +87,8 @@
$iker = 1;
$itaper = 1;
-# empty the output files directory and compile the measurement code
-print CSH "\\rm -rf OUTPUT_FILES ; mkdir OUTPUT_FILES ; make clean\n make\n";
+# empty the output files directory
+print CSH "\\rm -rf OUTPUT_FILES ; mkdir OUTPUT_FILES\n";
# copy chi values into PLOTS
print CSH "cp window_chi PLOTS\n";
@@ -99,6 +97,7 @@
print CSH "write_par_file.pl OUTPUT_FILES $ibp $Ts $tvec $itaper $iker $wtr/$npi $iparbools $par1 $par2 $par3\n";
# run the measurement code and save output file
+# (re-compilation is not necessary)
print CSH "mt_measure_adj > run_file\n";
print CSH "cp OUTPUT_FILES/*adj ${plot_adj_dir}\n";
@@ -121,7 +120,7 @@
#print CSH "make_pdf_by_event.pl $isort $otag\n";
#print CSH "/home/carltape/bin/pdcat -r *.pdf mt_cc_all.pdf\n"; # replace with make_pdf_by_event.pl
-print CSH "cd -\n";
+print CSH "cd $pwd\n";
#-----------------------------------------
close(CSH);
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl 2009-07-27 17:49:07 UTC (rev 15479)
@@ -42,6 +42,8 @@
if (@ARGV < 14) {die("Usage: run_mt_measure_adj.pl smodel iker itaper ibp Tmin/Tmax tstart/dt/npt itest iplot iboth iparbools par1 par2 par3\n")}
($smodel,$iker,$itaper,$ibp,$Ts,$lcut,$tvec,$itest,$iplot,$iboth,$iparbools,$par1,$par2,$par3) = @ARGV;
+$pwd = $ENV{PWD};
+
$plot_adj_dir = "PLOTS/ADJOINT_SOURCES";
#$iboth = 0; # option to plot both cross-correlation and multitaper adjoint sources
@@ -71,7 +73,8 @@
open(CSH,">$cshfile");
# empty the output files directory and compile the measurement code
-print CSH "\\rm -rf OUTPUT_FILES ; mkdir OUTPUT_FILES ; make clean\n make\n";
+# (this is done in prepare_mt_measure_adj.pl)
+#print CSH "\\rm -rf OUTPUT_FILES ; mkdir OUTPUT_FILES ; make clean\n make\n";
#----------------------------
# choose ipar2 based on the values used in the windowing code PAR_FILE (Feb 2009)
@@ -119,7 +122,7 @@
print CSH "cp OUTPUT_FILES/*adj ${plot_adj_dir}\n";
# copy reconstructed records into PLOTS
- print CSH "cp OUTPUT_FILES/*recon.cc ${plot_recon_dir}\n";
+ print CSH "cp OUTPUT_FILES/*recon.cc* ${plot_recon_dir}\n";
# copy chi values into PLOTS
print CSH "cp window_chi PLOTS\n";
@@ -150,7 +153,7 @@
# print CSH "make_pdf_by_event.pl $isort $otag\n";
#print CSH "/home/carltape/bin/pdcat -r *.pdf mt_cc_all.pdf\n"; # replace with make_pdf_by_event.pl
- print CSH "cd -\n";
+ print CSH "cd $pwd\n";
}
#-----------------------------------------
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2009-07-25 02:48:38 UTC (rev 15478)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2009-07-27 17:49:07 UTC (rev 15479)
@@ -37,6 +37,8 @@
if (@ARGV < 8) {die("Usage: run_tomo.pl smodel ibandpass Tmin/Tmax dt iparbools par1 par2 par3 lcut_frac\n")}
($smodel,$ibp,$Ts,$dt,$iparbools,$par1,$par2,$par3,$lcut_frac) = @ARGV;
+$pwd = $ENV{PWD};
+
($Tmin,$Tmax) = split("/",$Ts);
$sTmin = sprintf("T%3.3i",$Tmin);
$sTmax = sprintf("T%3.3i",$Tmax);
@@ -49,6 +51,7 @@
# USER INPUT
$dir0 = "/home/carltape";
+#$dir0 = "/media/raid/carltape/SOCAL_ADJOINT";
$dir_out = "$dir0/RUNS";
$dir_pdf_eid = "$dir0/results/MEASUREMENTS/${smodel}/PDF_EID_${Ttag}";
$dir_pdf_all = "$dir0/results/MEASUREMENTS/${smodel}/PDF_ALL_${Ttag}";
@@ -80,9 +83,9 @@
}
#die("testing");
-$imin = 1; $imax = $nevent;
+#$imin = 1; $imax = $nevent;
#$imin = 100; $imax = $nevent;
-#$imin = 2; $imax = $imin;
+$imin = 183; $imax = $imin;
# EXAMPLE lines from EIDs_simulation_duration
# 2 9967025 200.000 18200 0.29 4.084
@@ -152,6 +155,7 @@
print CSH "rm -rf ${dir_pdf_all}/*${eid}*\n";
print CSH "sleep 3s\n";
+ # this also compiles the measurement code, which is not needed for each new event
print "--> run the measurement code and make adjoint sources\n";
print CSH "prepare_mt_measure_adj.pl $smodel 0 1 0 $Ts $eid\n";
@@ -163,15 +167,14 @@
# print CSH "run_mt_cc_plot.pl $smodel $ibp $Ts $lcut $tvec 0 $iparbools $par1 $par2 $par3\n";
# }
-# # SOCAL: as of model m04, we are no longer making MT measurements
-# print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+ # SOCAL: MT and CC, with both ajoint source plotted
+ #print CSH "run_mt_cc_plot.pl $smodel $ibp $Ts $lcut $tvec 0 $iparbools $par1 $par2 $par3\n";
-
# SOCAL: CC only
- print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+ #print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
# SOCAL: multitaper
- #print CSH "run_mt_measure_adj.pl $smodel 1 1 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+ print CSH "run_mt_measure_adj.pl $smodel 1 1 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
# SOCAL: coverage kernels for m16 (CC banana-doughnut adjoint sources)
#print CSH "run_mt_measure_adj.pl $smodel 3 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
@@ -194,14 +197,15 @@
$infile5 = "${ftag}_ALL.pdf";
$outfile5 = "${dir_meas}/${infile5}";
print CSH "cp $infile5 $outfile5\n"; # composite PDF into RUN
- print CSH "mv $infile5 ${dir_pdf_eid}\n"; # composite PDF into single directory
- print CSH "mv *.pdf ${dir_pdf_all}\n"; # individual PDFs into single directory
- print CSH "cd -\n";
+ if (0==1) {
+ print CSH "mv $infile5 ${dir_pdf_eid}\n"; # composite PDF into single directory
+ print CSH "mv *.pdf ${dir_pdf_all}\n"; # individual PDFs into single directory
+ }
+ print CSH "cd $pwd\n";
# move adjoint sources into RUN directory
- print CSH "rm -rf $dir_adj\n";
- #print CSH "cp -r ADJOINT_SOURCES $dir_adj\n";
- print CSH "mv ADJOINT_SOURCES $dir_adj\nmkdir ADJOINT_SOURCES\n";
+ print CSH "rm -rf ${dir_adj}\nsleep 5s\n";
+ print CSH "mv ADJOINT_SOURCES ${dir_adj}; mkdir ADJOINT_SOURCES\n";
print CSH "echo done with Event $eid -- $ievent out of $imax/$nevent\n";
}
More information about the CIG-COMMITS
mailing list