[cig-commits] r16168 - in seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate: . INPUT PLOTTING matlab src
carltape at geodynamics.org
carltape at geodynamics.org
Sun Jan 24 15:32:09 PST 2010
Author: carltape
Date: 2010-01-24 15:32:07 -0800 (Sun, 24 Jan 2010)
New Revision: 16168
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xy_pert_sigmas.dat
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xyt_pert.dat
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xyt_pert_sigmas.dat
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_test.m
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/Makefile
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex02.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex03.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex04.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex05.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex06.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex07.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub.f90
Log:
Eliminated the second structural variable (bulk sound) for the membrane wave inversions.
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xy_pert_sigmas.dat
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xy_pert_sigmas.dat (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xy_pert_sigmas.dat 2010-01-24 23:32:07 UTC (rev 16168)
@@ -0,0 +1 @@
+ 2.000000000000e+03 2.000000000000e+03 5.000000000000e-01
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xyt_pert.dat
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xyt_pert.dat (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xyt_pert.dat 2010-01-24 23:32:07 UTC (rev 16168)
@@ -0,0 +1,25 @@
+ 1.267383299789e+02 -9.570977571789e+01 -6.286303233874e-01
+ 1.969883956224e+03 1.446792358321e+03 3.182555521271e-01
+ -8.144705530881e+02 -2.390402904540e+02 -1.996440975760e-01
+ -2.194790660469e+03 -8.608647334048e+02 -1.466697292799e-01
+ -2.668343248704e+03 3.093533072563e+03 -2.523697752226e-01
+ -4.824024741007e+00 -2.437783873262e+02 -3.765361466706e-01
+ 4.641535642308e+03 1.129355023043e+03 4.294870695967e-01
+ 1.143152245849e+03 -1.105850697102e+03 -2.812512701569e-01
+ -1.836714742299e+03 -4.758128723658e+02 -1.845912730572e-01
+ -3.384600275225e+02 2.300381319903e+03 -2.725172520900e-01
+ 2.262996471575e+03 -1.603336211662e+03 3.827927475534e-02
+ -1.324845384413e+03 1.155804447578e+03 -1.373179179902e-01
+ -7.985609180655e+02 3.949525656308e+02 -2.273747165639e-01
+ -1.877787227075e+03 1.132380269566e+03 -2.667791187071e-01
+ -6.944090254420e+01 2.677108493911e+03 -2.776452522937e-01
+ -4.416364228716e+02 -8.376638282015e+02 2.963486554201e-01
+ 1.385871535299e+03 -8.417102927914e+02 1.778779430494e-01
+ 2.643637190779e+03 5.021008993197e+02 -2.882550971269e-01
+ 1.671586741220e+03 -9.327837754046e+02 -5.723542443655e-02
+ -1.328798763517e+02 -3.268021867436e+03 -1.799233935729e-01
+ 5.237739337400e+02 3.376722966085e+03 -5.322138237609e-01
+ 2.994986799913e+03 1.976546294905e+03 1.252838013974e-01
+ 3.111593417954e+03 1.490199531609e+03 4.590259909292e-01
+ -3.967165351755e+03 2.703428143245e+03 2.042804539675e-01
+ -4.780331544454e+02 -4.979069835893e+03 -3.947758168940e-01
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xyt_pert_sigmas.dat
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xyt_pert_sigmas.dat (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/INPUT/events_xyt_pert_sigmas.dat 2010-01-24 23:32:07 UTC (rev 16168)
@@ -0,0 +1 @@
+ 2.000000000000e+03 2.000000000000e+03 3.000000000000e-01
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/Makefile
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/Makefile 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/Makefile 2010-01-24 23:32:07 UTC (rev 16168)
@@ -17,7 +17,7 @@
OBJ_DIR = obj
BIN_DIR = .
#MAIN = wave2d wave2d_2 wave2d_3 wave2d_4 test_smooth
-MAIN = wave2d
+MAIN = wave2d wave2d_2
MOD_FLAG = module
MOD_OBJ = $(patsubst %,$(OBJ_DIR)/%.o,$(MOD))
@@ -25,7 +25,7 @@
OBJ = $(F90_OBJ) $(MOD_OBJ)
#all : wave2d wave2d_2 wave2d_3 wave2d_4 test_smooth
-all : wave2d
+all : wave2d wave2d_2
$(MAIN) : % : $(SRC_DIR)/%.f90 $(F90_OBJ) $(MOD_OBJ)
$(F90) -o $(BIN_DIR)/$* $(F90_FLAGS) $(SRC_DIR)/$*.f90 -I$(MOD_FLAG) -J$(MOD_DIR) $(OBJ)
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl 2010-01-24 23:32:07 UTC (rev 16168)
@@ -20,6 +20,10 @@
# ../joint_inversion.pl -6/3.0/0/80/0.08 1 600 1 16 0 # 25 sources, structure only
# ../joint_inversion.pl -6/3.0/0/80/0.08 1 700 1 2 1 # 25 sources, structure only, SUBSPACE
#
+# ../joint_inversion.pl -6/3.0/0/80/0.08 1 800 1 16 0 # 1 sources, structure
+# ../joint_inversion.pl -6/3.0/0/80/0.08 1 XXX 1 16 0 # 5 sources, joint
+# ../joint_inversion.pl -6/3.0/0/80/0.08 1 XXX 1 16 0 # 25 sources, joint
+#
#==========================================================
if (@ARGV < 6) {die("Usage: plot_kernels.pl colors iker irun0 ipoly qmax HESSIAN\n");}
@@ -452,9 +456,12 @@
}
if ($icolor==1) {
#print CSH "awk '{print \$1,\$2,\$6}' $str_files[$k] | pscontour $R $J -A- -C$cpt_vel -I -O -K -V >> $psfile\n";
- print CSH "awk '{print \$1,\$2,\$7}' $str_files[$k] | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "awk '{print \$1,\$2,\$6}' $str_files[$k] | nearneighbor -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C${cpt_vel} $J -K -O -V -Q >> $psfile\n";
}
+
+ print "\n $str_files[$k]\n";
+
print CSH "pscoast $J $R $coast_info -K -O -V >> $psfile\n";
print CSH "awk '{print \$1,\$2}' $evefile | psxy -N $J $R -K -O -V $src >> $psfile\n";
print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $recfile | psxy -N $J $R -K -O -V $rec >> $psfile\n";
@@ -502,7 +509,7 @@
print CSH "psbasemap $B $R $J -K -O -V $shift >> $psfile\n";
if ($icolor==1) {
- print CSH "awk '{print \$1,\$2,\$7}' $file1dat_str | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "awk '{print \$1,\$2,\$6}' $file1dat_str | nearneighbor -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C${cpt_vel} $J -K -O -V -Q >> $psfile\n";
}
print CSH "pscoast $J $R $coast_info -K -O -V >> $psfile\n";
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl 2010-01-24 23:32:07 UTC (rev 16168)
@@ -185,8 +185,8 @@
# phase velocity map
print CSH "psbasemap $B1syn $R $J1 -P -K -V $origin > $psfile\n"; # START
- #print CSH "awk '{print \$1,\$2,\$7}' $file1syn | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile\n";
- print CSH "awk '{print \$1,\$2,\$7}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
+ #print CSH "awk '{print \$1,\$2,\$6}' $file1syn | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$1,\$2,\$6}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
print CSH "pscoast $J1 $R $B1syn -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
#print CSH "awk '{print \$2,\$1}' $idir/BOUNDARIES/oms_shelf |psxy $J1 $R $Wshelf -K -O -V >> $psfile\n";
@@ -206,8 +206,8 @@
# phase velocity map -- data
print CSH "psbasemap $B1dat $R $J1 -P -K -V $origin > $psfile\n"; # START
- #print CSH "awk '{print \$1,\$2,\$7}' $file1dat | pscontour $R $J1 -A- -C$cptfile1 -I -V -K -O >> $psfile\n";
- print CSH "awk '{print \$1,\$2,\$7}' $file1dat | nearneighbor -G$grdfile $R $interp\n";
+ #print CSH "awk '{print \$1,\$2,\$6}' $file1dat | pscontour $R $J1 -A- -C$cptfile1 -I -V -K -O >> $psfile\n";
+ print CSH "awk '{print \$1,\$2,\$6}' $file1dat | nearneighbor -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
print CSH "pscoast $J1 $R $B1dat -W1p -Na/1p -Dh -V -K -O >> $psfile\n";
#print CSH "awk '{print \$2,\$1}' $idir/BOUNDARIES/oms_shelf |psxy $J1 $R $Wshelf -V -K -O >> $psfile\n";
@@ -219,8 +219,8 @@
# phase velocity map -- synthetics
print CSH "psbasemap $B1syn $R $J1 -K -O -V $shift >> $psfile\n";
- #print CSH "awk '{print \$1,\$2,\$7}' $file1syn | pscontour $R $J1 -A- -C$cptfile1 -I -V -K -O >> $psfile\n";
- print CSH "awk '{print \$1,\$2,\$7}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
+ #print CSH "awk '{print \$1,\$2,\$6}' $file1syn | pscontour $R $J1 -A- -C$cptfile1 -I -V -K -O >> $psfile\n";
+ print CSH "awk '{print \$1,\$2,\$6}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
print CSH "pscoast $J1 $R $B1syn -W1p -Na/1p -Dh -V -K -O >> $psfile\n";
print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -V -K -O >> $psfile \n";
@@ -264,8 +264,8 @@
# phase velocity map
print CSH "psbasemap $B1dat $R $J1 -K -O -V $shift >> $psfile\n";
- #print CSH "awk '{print \$1,\$2,\$7}' $file1dat | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile\n";
- print CSH "awk '{print \$1,\$2,\$7}' $file1dat | nearneighbor -G$grdfile $R $interp\n";
+ #print CSH "awk '{print \$1,\$2,\$6}' $file1dat | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$1,\$2,\$6}' $file1dat | nearneighbor -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
print CSH "pscoast $J1 $R $B1dat -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
@@ -292,8 +292,8 @@
# phase velocity map
print CSH "psbasemap $B1syn $R $J1 -K -O -V $shift >> $psfile\n";
- #print CSH "awk '{print \$1,\$2,\$7}' $file1syn | pscontour $R $J1 -A- -C$cptfile1 $shift -I -O -K -V >> $psfile\n";
- print CSH "awk '{print \$1,\$2,\$7}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
+ #print CSH "awk '{print \$1,\$2,\$6}' $file1syn | pscontour $R $J1 -A- -C$cptfile1 $shift -I -O -K -V >> $psfile\n";
+ print CSH "awk '{print \$1,\$2,\$6}' $file1syn | nearneighbor -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
print CSH "pscoast $J1 $R $B1syn -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model_all.pl 2010-01-24 23:32:07 UTC (rev 16168)
@@ -6,7 +6,9 @@
# Carl Tape
# 15-Jan-2010
#
-# Plot a 6-panel figure for each model iteration.
+# Full display of models and model perturbations.
+# FIGURE 1: 6-panel figure
+# FIGURE 2: 9-panel figure (emphasizing CG test model)
# Input allows for a range of models as well.
#
# Example:
@@ -18,6 +20,11 @@
($irun0,$imodels,$ievents,$pbar,$ifinite) = @ARGV;
$comp = 1;
+# USER INPUT
+$ifig1 = 1;
+$ifig2 = 1;
+$itest = 1; # plot CG test models or not
+
$icolor = 1;
$ijpg = 0;
$ipdf = 1;
@@ -100,23 +107,12 @@
#-------------------------
-# write plotting scripts
-$wid = 2.5;
-$J1 = "-JM${wid}i"; # in lat-lon
-$origin = "-X1.5 -Y7.25";
+$wid1 = 2.5;
+$wid2 = 2.0;
$B1dat = "-B1:.\" \":WeSN";
$B1syn = "-B1:.\" \":wESN";
-$Dlen = 2.0;
-$Dx = $wid/2;
-$Dy = -0.35;
-
-$Dscale1 = "-D$Dx/$Dy/$Dlen/0.15h";
-
-$Bscale1 = sprintf("-B%2.2e:\" Perturbation from %3.3f km s\@+-1\@+ \": -E10p",$ptick,$beta0/1000);
-$Bscale2 = sprintf("-B%2.2e:\" Perturbation from m00 \": -E10p",$ptick);
-
#===============================================
$cshfile = "plot_surf_model_all.csh";
print "\nWriting to CSH file ${cshfile}...\n";
@@ -137,17 +133,22 @@
# sides to show the lat-lon tick marks
$B0 = "-B1:.\" \":";
- at Bopts = ("WesN","wEsN","Wesn","wEsn","Wesn","wEsn");
+ at Bopts1 = ("WesN","wEsN","Wesn","wEsn","Wesn","wEsn");
+ at Bopts2 = ("WesN","wesN","wEsN","Wesn","wesn","wEsn","Wesn","wesn","wEsn");
-$dX = -1.1*$wid; $dY = -1.1*$wid; $shift1 = "-X$dX -Y$dY";
-$dX = 1.1*$wid; $dY = 0; $shift2 = "-X$dX -Y$dY";
- at shifts = ($origin,$shift2,$shift1,$shift2,$shift1,$shift2);
+# figure 1
+$dX = 1.1*$wid1; $dY = 0; $shift1 = "-X$dX -Y$dY";
+$dX = -1.1*$wid1; $dY = -1.1*$wid1; $shift2 = "-X$dX -Y$dY";
+ at shifts1 = ($origin,$shift1,$shift2,$shift1,$shift2,$shift1);
- at finds = (7,7,7,7,7,7,7);
+# figure 2
+$dX = 1.1*$wid2; $dY = 0; $shift1 = "-X$dX -Y$dY";
+$dX = -2.2*$wid2; $dY = -1.1*$wid2; $shift2 = "-X$dX -Y$dY";
+ at shifts2 = ($origin,$shift1,$shift1,$shift2,$shift1,$shift1,$shift2,$shift1,$shift1);
-$xlab = $wid/2;
-$ylab = -0.02*$wid;
-$olab = "-Xa$xlab -Ya$ylab";
+$fcol = 6;
+#@finds = (7,7,7,7,7,7,7,7,7);
+
$J_title = "-JM1";
$R_title = "-R0/1/0/1";
@@ -156,111 +157,325 @@
# loop over models
for ($m = $imodelmin; $m <= $imodelmax; $m = $m+1) {
- $imodel = $m;
+ $imodel = $m; # current model
+ $mtest = $m % 2; # =1 for a test model
+
$stirun = sprintf("%4.4i",$irun0+$imodel);
- $stmod = sprintf("m%2.2i",$imodel/2);
+ if ($mtest==1) {
+ $stmod = sprintf("m%2.2it",$imodel/2);
+ } else {
+ $stmod = sprintf("m%2.2i",$imodel/2);
+ }
- $dir1 = "$odir/run_${stirun}";
+ # directories (current model and previos model)
+ $stirun = sprintf("%4.4i",$irun0+$imodel);
+ $dir1 = "$odir/run_${stirun}";
$name = "model_all_${stirun}";
$psfile = "$name.ps";
$jpgfile = "$name.jpg";
- @files = ("$dir0/${mfile_syn}","$dir0/${mfile_syn}","$dir1/${mfile_syn}","$dir1/${mfile_syn}","$dir0/${mfile_dat}","$dir0/${mfile_dat}");
- @stlabs = ("m00","ln(m00/m00)",$stmod,"ln($stmod/m00)","mtar","ln(mtar/m00)");
+ if ($ifig1 == 1) {
- # only plot even models (m04 = iteration 8)
- $kmax = 6*(1 - $m % 2);
+# write plotting scripts
+$wid = $wid1;
+$J1 = "-JM${wid}i"; # in lat-lon
+$origin = "-X1.5 -Y7.25";
+$Dlen = 2.0; $Dx = $wid/2; $Dy = -0.35;
+$Dscale1 = "-D$Dx/$Dy/$Dlen/0.15h";
+$Bscale1 = sprintf("-B%2.2e:\" Perturbation from %3.3f km s\@+-1\@+ \": -E10p",$ptick,$beta0/1000);
+$Bscale2 = sprintf("-B%2.2e:\" Perturbation from m00 \": -E10p",$ptick);
+$xlab = $wid/2;
+$ylab = -0.02*$wid;
+$olab = "-Xa$xlab -Ya$ylab";
- # loop over subplots -- only for even m
- for ($k = 0; $k < $kmax; $k = $k+1) {
+ print "\n----FIGURE 1--model $stmod----\n";
- $iB = $Bopts[$k];
- $shift = $shifts[$k];
- $B = "${B0}${iB}";
- $find = $finds[$k];
- $find2 = 7 + $find;
+ @files = ("$dir0/${mfile_syn}","$dir0/${mfile_syn}",
+ "$dir1/${mfile_syn}","$dir1/${mfile_syn}",
+ "$dir0/${mfile_dat}","$dir0/${mfile_dat}");
+ @stlabs = ("m00","ln(m00/m00)",$stmod,"ln($stmod/m00)","mtar","ln(mtar/m00)");
- $file0 = $files[0];
- $file1 = $files[$k];
- if (not -f $file1) {die("Check if $file1 exist or not\n");}
-
- if ($k % 2 == 1) {
- $file2 = "ftemp";
- print CSH "paste $file1 $file0 > $file2\n";
+ # only plot even models (m04 = iteration 8)
+ $kmaxt = 6;
+ $kmaxm = 6*(1 - $mtest); # =0 for test model
+ if ($itest == 1) {
+ $kmax = $kmaxt;
+ } else {
+ $kmax = $kmaxm;
}
- # phase velocity map
- if ($k==0) {
- # set bounds for the plotting
- #$xmin = -121; $xmax = -115.0; $zmin = 31.75; $zmax = 36.5;
- ($xmin,$xmax,$zmin,$zmax,$smin,$smax,$tmin,$tmax) = split(" ",`minmax -C $file1`);
- $dinc = 0.25; # buffer around min/max
- $xmin = $xmin-$dinc; $xmax = $xmax+$dinc;
- $zmin = $zmin-$dinc; $zmax = $zmax+$dinc;
- $R = "-R$xmin/$xmax/$zmin/$zmax";
- print "\n$R\n";
+ # loop over subplots
+ for ($k = 0; $k < $kmax; $k = $k+1) {
- print CSH "psbasemap $B $R $J1 -K -V -P $origin > $psfile\n"; # START
- } else {
- print CSH "psbasemap $B $R $J1 -K -O -V $shift >> $psfile\n";
- }
+ $iB = $Bopts1[$k];
+ $shift = $shifts1[$k];
+ $B = "${B0}${iB}";
+ $find = $fcol;
+ $find2 = 2*$fcol;
- if ($icolor == 1) {
- #print CSH "awk '{print \$1,\$2,\$7}' $file1 | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile\n";
- if ($k % 2 == 0) {
- print CSH "awk '{print \$1,\$2,\$${find}}' $file1 | nearneighbor -G$grdfile $R $interp\n";
+ $file0 = $files[0]; # initial model
+ $file1 = $files[$k];
+ if (not -f $file1) {
+ die("Check if $file1 exist or not\n");
+ }
+
+ # column-concatenate current model and initial model
+ if ($k % 2 == 1) {
+ $file2 = "ftemp";
+ print CSH "paste $file1 $file0 > $file2\n";
+ }
+
+ # phase velocity map
+ if ($k==0 && $m==$imodelmin) {
+ # set bounds for the plotting
+ #$xmin = -121; $xmax = -115.0; $zmin = 31.75; $zmax = 36.5;
+ ($xmin,$xmax,$zmin,$zmax,$smin,$smax,$tmin,$tmax) = split(" ",`minmax -C $file1`);
+ $dinc = 0.25; # buffer around min/max
+ $xmin = $xmin-$dinc; $xmax = $xmax+$dinc;
+ $zmin = $zmin-$dinc; $zmax = $zmax+$dinc;
+ $R = "-R$xmin/$xmax/$zmin/$zmax";
+ print "\n$R\n";
+ }
+
+ if ($k==0) {
+ print CSH "psbasemap $B $R $J1 -K -V -P $origin > $psfile\n"; # START
} else {
- print CSH "awk '{print \$1,\$2,\$${find}-\$${find2}}' $file2 | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "psbasemap $B $R $J1 -K -O -V $shift >> $psfile\n";
}
- print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
- }
- print CSH "pscoast $J1 $R $B -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
- # plot receivers
- if (0==1) { # numbered large circles
- print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N $J1 $R -K -O -V $rec0 >> $psfile\n";
- $rec_file2 = text_rec; $angle = 0; $just = "CM";
- print CSH "awk '\$1 == \"R\" {print \$2,\$3,$fsize3,$angle,$fontno,\"$just\",\$4}' $rec_file > $rec_file2\n";
- print CSH "pstext $rec_file2 -N $J1 $R -K -O -V >> $psfile\n";
+ if ($icolor == 1) {
+ #print CSH "awk '{print \$1,\$2,\$7}' $file1 | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile\n";
+ if ($k % 2 == 0) {
+ print CSH "awk '{print \$1,\$2,\$${find}}' $file1 | nearneighbor -G$grdfile $R $interp\n";
+ } else {
+ print CSH "awk '{print \$1,\$2,\$${find}-\$${find2}}' $file2 | nearneighbor -G$grdfile $R $interp\n";
+ }
+ print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile\n";
+ }
+ print CSH "pscoast $J1 $R $B -W1p -Na/1p -Dh -K -O -V >> $psfile\n";
- } else { # small circles
- print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N -Sc5p -W0.5p $J1 $R -K -O -V >> $psfile\n";
+ # plot receivers
+ if (0==1) { # numbered large circles
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N $J1 $R -K -O -V $rec0 >> $psfile\n";
+ $rec_file2 = text_rec; $angle = 0; $just = "CM";
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3,$fsize3,$angle,$fontno,\"$just\",\$4}' $rec_file > $rec_file2\n";
+ print CSH "pstext $rec_file2 -N $J1 $R -K -O -V >> $psfile\n";
- }
+ } else { # small circles
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N -Sc5p -W0.5p $J1 $R -K -O -V >> $psfile\n";
- # plot sources
- if ($ievent_one) {print CSH "awk '\$1 == \"S\" {print \$2,\$3}' $rec_file | psxy $src0 -N $J1 $R -K -O -V >> $psfile\n";}
- if ($ievent_all) {print CSH "awk '{print \$1,\$2}' $evefile | psxy $src -N $J1 $R -K -O -V >> $psfile\n";}
+ }
- # plot label
- $stlab = $stlabs[$k];
- $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
- print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 12 0 $fontno CM $stlab\nEOF\n";
+ # plot sources
+ if ($ievent_one) {
+ print CSH "awk '\$1 == \"S\" {print \$2,\$3}' $rec_file | psxy $src0 -N $J1 $R -K -O -V >> $psfile\n";
+ }
+ if ($ievent_all) {
+ print CSH "awk '{print \$1,\$2}' $evefile | psxy $src -N $J1 $R -K -O -V >> $psfile\n";
+ }
- # colorscale bars
- if ($k==4) {print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile \n";}
- if ($k==5) {print CSH "psscale -C$cptfile1 $Dscale1 $Bscale2 -K -O -V >> $psfile \n";}
+ # plot label
+ $stlab = $stlabs[$k];
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 12 0 $fontno CM $stlab\nEOF\n";
- # plot title and GMT header
- if ($k==5) {
- $plot_title = "Model m00, model $stmod, and target model (run_${stirun})";
- $Utag = "-U/0/0.25/$plabel";
- $shift = "-Xa-$dX -Ya8.5";
- print CSH "pstext -N $J1 $R $Utag -O -V $shift >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno LM $plot_title\nEOF\n"; # FINISH
+ # colorscale bars
+ if ($k==4) {
+ print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile\n";
+ }
+ if ($k==5) {
+ print CSH "psscale -C$cptfile1 $Dscale1 $Bscale2 -K -O -V >> $psfile\n";
+ }
- if($ijpg == 1) {print CSH "convert $psfile $jpgfile\n";}
- if($ipdf == 1) {print CSH "ps2pdf $psfile\n";}
+ # plot title and GMT header
+ if ($k==5) {
+ $plot_title = "Model m00, model $stmod, and target model (run_${stirun})";
+ $Utag = "-U/0/0.25/$plabel";
+ $shift = "-Xa-3.5 -Ya8.5";
+ print CSH "pstext -N $J1 $R $Utag -O -V $shift >>$psfile<<EOF\n $xmin $zmin $fsize0 0 $fontno LM $plot_title\nEOF\n"; # FINISH
- }
- } # for loop
+ if ($ijpg == 1) {
+ print CSH "convert $psfile $jpgfile\n";
+ }
+ if ($ipdf == 1) {
+ print CSH "ps2pdf $psfile\n";
+ }
+ }
+ } # loop over subplots
+ }
+
#-----------------------------
+ # FIGURE 2
+ # if you want fig 2 and the model is m01,m02,m03,etc (iterations 2,3,4,etc)
+ if ( ($ifig2 == 1) && (($imodel >= 1) && ($mtest==0)) ) {
+
+ # previous test model
+ $imodelt = $imodel - 1;
+ $stmodt = sprintf("m%2.2it",$imodelt/2);
+ $stirunt = sprintf("%4.4i",$irun0+$imodelt);
+ $dir1t = "$odir/run_${stirunt}";
+
+ # previous model
+ $imodelp = $imodel - 2;
+ $stmodp = sprintf("m%2.2i",$imodelp/2);
+ $stirunp = sprintf("%4.4i",$irun0+$imodelp);
+ $dir1p = "$odir/run_${stirunp}";
+
+ print "\n--------FIGURE 2---models $stmodp, $stmodt, $stmod---\n";
+
+ $name2 = "model_pert_${stirun}";
+ $psfile2 = "$name2.ps";
+ $jpgfile2 = "$name2.jpg";
+
+# write plotting scripts
+$wid = $wid2;
+$J1 = "-JM${wid}i"; # in lat-lon
+$origin = "-X1 -Y7.25";
+$Dlen = 1.5; $Dx = $wid/2; $Dy = -0.35;
+$Dscale1 = "-D$Dx/$Dy/$Dlen/0.15h";
+$Bscale1 = sprintf("-B%2.2e:\" Perturbation from %3.3f km s\@+-1\@+ \": -E10p",$ptick,$beta0/1000);
+$Bscale2 = sprintf("-B%2.2e:\" Perturbation from $stmodp \": -E10p",$ptick);
+$Bscale3 = sprintf("-B%2.2e:\" Perturbation from mtar\": -E10p",$ptick);
+$xlab = $wid/2;
+$ylab = -0.02*$wid;
+$olab = "-Xa$xlab -Ya$ylab";
+
+ # nine subplots
+ @files0 = ("$dir1p/${mfile_syn}","$dir1p/${mfile_syn}","$dir1p/${mfile_syn}",
+ "$dir1t/${mfile_syn}","$dir1p/${mfile_syn}","$dir1t/${mfile_syn}",
+ "$dir1/${mfile_syn}","$dir1p/${mfile_syn}","$dir1/${mfile_syn}");
+ @files1 = ("$dir1p/${mfile_syn}","$dir1p/${mfile_syn}","$dir0/${mfile_dat}",
+ "$dir1t/${mfile_syn}","$dir1t/${mfile_syn}","$dir0/${mfile_dat}",
+ "$dir1/${mfile_syn}","$dir1/${mfile_syn}","$dir0/${mfile_dat}");
+ @stlabs = ($stmodp,"ln($stmodp/$stmodp)","ln(mtar/$stmodp)",
+ $stmodt,"ln($stmodt/$stmodp)","ln(mtar/$stmodt)",
+ $stmod,"ln($stmod/$stmodp)","ln(mtar/$stmod)");
+
+ $kmax = 9;
+
+ # loop over subplots
+ for ($k = 0; $k < $kmax; $k = $k+1) {
+
+ $iB = $Bopts2[$k];
+ $shift = $shifts2[$k];
+ $B = "${B0}${iB}";
+ $find = $fcol;
+ $find2 = 2*$fcol;
+
+ $kmod = $k % 3;
+ #print "\n -- k = $k -- kmod = $kmod --\n";
+
+ # two files to concatenate
+ $file0 = $files0[$k];
+ $file1 = $files1[$k];
+ if (not -f $file0) {
+ die("Check if file0 $file0 exist or not\n");
+ }
+
+ # column-concatenate two files (columns 2 and 3)
+ if ($kmod > 0) {
+ if (not -f $file1) {
+ die("Check if file1 $file1 exist or not\n");
+ }
+ $file2 = "ftemp";
+ print "concatenating two files:\n $file0\n $file1\n";
+ print CSH "paste $file1 $file0 > $file2\n";
+ }
+
+ # phase velocity map
+ if ($k==0) {
+ # set bounds for the plotting
+ #$xmin = -121; $xmax = -115.0; $zmin = 31.75; $zmax = 36.5;
+ ($xmin,$xmax,$zmin,$zmax,$smin,$smax,$tmin,$tmax) = split(" ",`minmax -C $file1`);
+ $dinc = 0.25; # buffer around min/max
+ $xmin = $xmin-$dinc; $xmax = $xmax+$dinc;
+ $zmin = $zmin-$dinc; $zmax = $zmax+$dinc;
+ $R = "-R$xmin/$xmax/$zmin/$zmax";
+ print "\n$R\n";
+ }
+
+ if ($k==0) {
+ print CSH "psbasemap $B $R $J1 -K -V -P $origin > $psfile2\n"; # START
+ } else {
+ print CSH "psbasemap $B $R $J1 -K -O -V $shift >> $psfile2\n";
+ }
+
+ if ($icolor == 1) {
+ #print CSH "awk '{print \$1,\$2,\$7}' $file1 | pscontour $R $J1 -A- -C$cptfile1 -I -K -O -V >> $psfile2\n";
+ if ($kmod == 0) {
+ print CSH "awk '{print \$1,\$2,\$${find}}' $file1 | nearneighbor -G$grdfile $R $interp\n";
+ } else {
+ print CSH "awk '{print \$1,\$2,\$${find}-\$${find2}}' $file2 | nearneighbor -G$grdfile $R $interp\n";
+ }
+ print CSH "grdimage $grdfile -C$cptfile1 $J1 -K -O -V -Q >> $psfile2\n";
+ }
+ print CSH "pscoast $J1 $R $B -W1p -Na/1p -Dh -K -O -V >> $psfile2\n";
+
+ # plot receivers
+ if (0==1) { # numbered large circles
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N $J1 $R -K -O -V $rec0 >> $psfile2\n";
+ $rec_file2 = text_rec; $angle = 0; $just = "CM";
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3,$fsize3,$angle,$fontno,\"$just\",\$4}' $rec_file > $rec_file2\n";
+ print CSH "pstext $rec_file2 -N $J1 $R -K -O -V >> $psfile2\n";
+
+ } else { # small circles
+ print CSH "awk '\$1 == \"R\" {print \$2,\$3}' $rec_file | psxy -N -Sc5p -W0.5p $J1 $R -K -O -V >> $psfile2\n";
+
+ }
+
+ # plot sources
+ if ($ievent_one) {
+ print CSH "awk '\$1 == \"S\" {print \$2,\$3}' $rec_file | psxy $src0 -N $J1 $R -K -O -V >> $psfile2\n";
+ }
+ if ($ievent_all) {
+ print CSH "awk '{print \$1,\$2}' $evefile | psxy $src -N $J1 $R -K -O -V >> $psfile2\n";
+ }
+
+ # plot label
+ $stlab = $stlabs[$k];
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile2<<EOF\n0 0 10 0 $fontno CM $stlab\nEOF\n";
+
+ # colorscale bars
+ if ($k==6) {
+ print CSH "psscale -C$cptfile1 $Dscale1 $Bscale1 -K -O -V >> $psfile2\n";
+ }
+ if ($k==7) {
+ print CSH "psscale -C$cptfile1 $Dscale1 $Bscale2 -K -O -V >> $psfile2\n";
+ }
+ if ($k==8) {
+ print CSH "psscale -C$cptfile1 $Dscale1 $Bscale3 -K -O -V >> $psfile2\n";
+ }
+
+ # plot title and GMT header
+ if ($k==8) {
+ $plot_title = "Models $stmodp, $stmodt, $stmod (run_${stirun})";
+ $Utag = "-U/0/0.25/$plabel";
+ $shift = "-Xa-4.5 -Ya7";
+ print CSH "pstext -N $J1 $R $Utag -O -V $shift >>$psfile2<<EOF\n $xmin $zmin $fsize0 0 $fontno LM $plot_title\nEOF\n"; # FINISH
+
+ if ($ijpg == 1) {
+ print CSH "convert $psfile2 $jpgfile\n";
+ }
+ if ($ipdf == 1) {
+ print CSH "ps2pdf $psfile2\n";
+ }
+
+ }
+ } # loop over subplots
+
+ }
+
+ #-----------------------------
+
} # loop over models
close (CSH);
system("csh -f $cshfile");
-system("gv $psfile &")
+if($ifig1==1) {system("gv $psfile &")}
+if($ifig2==1) {system("gv $psfile2 &")}
#=================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README 2010-01-24 23:32:07 UTC (rev 16168)
@@ -1,4 +1,4 @@
-Carl Tape, 01-Dec-2009
+Carl Tape, 24-Jan-2010
carltape at fas.harvard.edu
/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README
@@ -76,6 +76,13 @@
Example 5: 25 sources, joint inversion using CG algorithm
Example 6: 25 sources, structure inversion using CG algorithm
Example 7: 25 sources, structure inversion using subspace method
+------------
+Using 2 Gaussian fields:
+Example 8: 1 source, structure inversion
+Example 9: 1 source, source inversion
+Example 10: 1 source, joint inversion
+Example 11: 5 sources, joint inversion
+Example 12: 25 sources, joint inversion
To get a sense for the output from these examples, open the composite PDF files
/PLOTTING/FIGURES/Example1_run_0100.pdf
@@ -117,13 +124,13 @@
------------------
EXAMPLE 2 (run_0200)
Make the following modifications to src/wave2d_constants.f90
- IRUNZ = 100 --> IRUNZ = 200
- PERT_STRUCT = 1 --> PERT_STRUCT = 0
- PERT_SOURCE_T = 0 --> PERT_SOURCE_T = 1
- PERT_SOURCE_X = 0 --> PERT_SOURCE_X = 1
- INV_STRUCT = 1 --> INV_STRUCT = 0
- INV_SOURCE_T = 0 --> INV_SOURCE_T = 1
- INV_SOURCE_X = 0 --> INV_SOURCE_X = 1
+ IRUNZ = 100 --> IRUNZ = 200
+ PERT_STRUCT_BETA = 1 --> PERT_STRUCT_BETA = 0
+ PERT_SOURCE_T = 0 --> PERT_SOURCE_T = 1
+ PERT_SOURCE_X = 0 --> PERT_SOURCE_X = 1
+ INV_STRUCT_BETA = 1 --> INV_STRUCT_BETA = 0
+ INV_SOURCE_T = 0 --> INV_SOURCE_T = 1
+ INV_SOURCE_X = 0 --> INV_SOURCE_X = 1
Then compile and execute:
make clean ; make wave2d ; wave2d
Follow steps from Example 1 to make figures.
@@ -131,9 +138,9 @@
------------------
EXAMPLE 3 (run_0300)
Make the following modifications to src/wave2d_constants.f90
- IRUNZ = 200 --> IRUNZ = 300
- PERT_STRUCT = 1 --> PERT_STRUCT = 1
- INV_STRUCT = 1 --> INV_STRUCT = 1
+ IRUNZ = 200 --> IRUNZ = 300
+ PERT_STRUCT_BETA = 1 --> PERT_STRUCT_BETA = 1
+ INV_STRUCT_BETA = 1 --> INV_STRUCT_BETA = 1
Then compile and execute:
make clean ; make wave2d ; wave2d
Follow steps from Example 1 to make figures.
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_test.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_test.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_test.m 2010-01-24 23:32:07 UTC (rev 16168)
@@ -0,0 +1,265 @@
+%
+% wave2d_cg_test.m
+% Carl Tape, 21-Jan-2010
+%
+% This program computes an updated model using a conjugate gradient
+% algorithm. This is essentially a test program in anticipation of
+% performing the inversion in Matlab OUTSIDE of wave2d.f90.
+%
+% calls xxx
+% called by xxx
+%
+
+format long
+format compact
+close all
+clear
+
+colors;
+npts = 100;
+stfm = '%4.4i';
+
+icheck = 1;
+
+stpars = {'C = ln(c/c0)','B = ln(beta/beta0)','Ts = ts = ts0','Xs = xs - xs0','Ys = ys - ys0'};
+
+%---------------------------------------------------------
+% USER INPUT
+
+% base directory
+dir0 = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/';
+if ~exist(dir0), error([dir0 ' does not exist']); end
+
+icheck = 1;
+
+parms = [850 5 5];
+irun0 = parms(1);
+ieventmin = parms(2);
+ieventmax = parms(3);
+
+nevent_run = ieventmax - ieventmin + 1;
+einds = [ieventmin:ieventmax];
+stirun0 = sprintf(stfm,irun0);
+dir1 = [dir0 'SEM2D_iterate_OUTPUT/run_' stirun0 '/'];
+if ~exist(dir1), error([dir1 ' does not exist']); end
+
+%---------------------------------------------------------
+% load files
+
+% misc constants
+fac_str = 1.0;
+fac_ts = 1.0;
+fac_xs = 1.0;
+fac_ys = 1.0;
+fac_total = 1.0;
+ugrad_str = 1.0;
+ugrad_ts = 1.0;
+ugrad_xs = 1.0;
+ugrad_ys = 1.0;
+
+coverage_str = 0.666 / 0.962;
+coverage_src = 0.946 / 1.018;
+sigma_DT = 0.1;
+nrec = 132;
+ncomp = 1;
+nmeas_run = nevent_run*nrec*ncomp;
+nparm_src_run = nevent_run;
+
+% % double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+% % double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+% % double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+% rho0 = 2.60e3;
+% kap0 = 5.20e10;
+% mu0 = 2.66e10;
+% beta0 = sqrt( mu0/rho0 )
+% bulk0 = sqrt( kap0 /rho0 )
+
+% indexing
+m_inds = load([dir1 'm_inds.dat']);
+nparm = m_inds(5,2);
+ngrid = m_inds(1,2);
+indC = m_inds(1,1):m_inds(1,2);
+indB = m_inds(2,1):m_inds(2,2);
+indT = m_inds(3,1):m_inds(3,2);
+indX = m_inds(4,1):m_inds(4,2);
+indY = m_inds(5,1):m_inds(5,2);
+indTXY = m_inds(3,1):m_inds(5,2);
+
+% local mesh
+% write(15,'(6i8,4e18.8)') k, ispec, i, j, iglob, valence(iglob), &
+% x(iglob), z(iglob), da_local(i,j,ispec), da_global(iglob)
+[i1,i2,i3,i4,i5,i6,xg,yg,da_local_vec,d10] = textread([dir1 'local_mesh.dat'],'%f%f%f%f%f%f%f%f%f%f');
+xmin = min(xg); xmax = max(xg);
+ymin = min(yg); ymax = max(yg);
+Atot = (xmax-xmin)*(ymax-ymin);
+% check
+sum(da_local_vec), (xmax-xmin)*(ymax-ymin)
+
+% plotting mesh
+xp = linspace(xmin,xmax,npts);
+yp = linspace(ymin,ymax,npts);
+[Xp,Yp] = meshgrid(xp,yp);
+
+% load current model
+mtemp = load([dir1 'structure_syn.dat']);
+m_str_lon = mtemp(:,1);
+m_str_lat = mtemp(:,2);
+m_str_kappa = mtemp(:,3);
+m_str_mu = mtemp(:,4);
+m_str_rho = mtemp(:,5);
+m_str_C = mtemp(:,6);
+m_str_B = mtemp(:,7);
+
+% model vector
+%write(88,'(2e24.12)') m0_vec(i), m0(i)
+[m0_vec, m0] = textread([dir1 'initial_model_vector.dat'],'%f%f');
+m0_C = m0(indC);
+m0_B = m0(indB);
+m0_T = m0(indT);
+m0_X = m0(indX);
+m0_Y = m0(indY);
+
+% load kernels
+kernels_all = zeros(ngrid,nevent_run);
+iplotker = 0;
+for ii=1:nevent_run
+ ievent = einds(ii)
+ dirK = [dir1 'event_' sprintf('%3.3i',ievent) '/']
+ kernel = load([dirK 'kernel_basis']); Kbeta = kernel(:,7);
+ %kernel = load([dirK 'kernel_basis_smooth']); Kbeta = kernel(:,3);
+ kernels_all(:,ii) = Kbeta;
+
+ if iplotker==1
+ Zp = griddata(xg,yg,Kbeta,Xp,Yp,'nearest');
+ figure; imagesc(Zp); set(gca,'ydir','normal');
+ caxis([-1 1]*max(abs(Kbeta))*0.8);
+ end
+end
+summed_ker = sum(kernels_all,2);
+
+% data covariance
+% nmeas_run = nevent_run * nrec * NCOMP
+% double precision, parameter :: SIGMA_DT = 0.10
+%cov_data(:) = SIGMA_DT * SIGMA_DT * nmeas_run
+covd_diag = load([dir1 'cov_data_diagonal.dat']);
+covd_const = sigma_DT^2 * nmeas_run;
+% check
+min(covd_diag), max(covd_diag), covd_const
+
+% data vector
+% chi_data(ievent,irec,icomp,1) = (tshift_xc_pert )**2 / cov_data(imeasure)
+% write(19,'(3i8,1e20.10)') ievent, irec, icomp, chi_data(ievent,irec,icomp,1)
+chi_all = load([dir1 'chi_all.dat']);
+chi_vec = chi_all(:,4);
+
+% measurement vector
+%measure_vec(imeasure,1) = tshift_xc_pert
+%measure_vec(imeasure,2) = tshift_xc
+%measure_vec(imeasure,3) = dlnA_pert
+%measure_vec(imeasure,4) = dlnA
+%measure_vec(imeasure,5) = 0.5 * DT*sum( adj_syn(:,icomp,irec)**2 )
+meas_all = load([dir1 'measure_vec.dat']);
+DTvec = meas_all(:,1);
+% check
+sum(DTvec.^2 / covd_const), sum(chi_vec)
+
+% sigma values
+% write(19,'(2e20.10)') sigma_bulk, m_scale_str(1)
+% write(19,'(2e20.10)') sigma_beta, m_scale_str(2)
+% write(19,'(2e20.10)') sigma_ts, m_scale_src(1)
+% write(19,'(2e20.10)') sigma_xs, m_scale_src(2)
+% write(19,'(2e20.10)') sigma_zs, m_scale_src(3)
+sigma_all = load([dir1 'sigma_values.dat']);
+sigma_bulk = sigma_all(1,1);
+sigma_beta = sigma_all(2,1);
+sigma_ts = sigma_all(3,1);
+sigma_xs = sigma_all(4,1);
+sigma_zs = sigma_all(5,1);
+
+%---------------------------------------------------------
+% emulate CG variables in wave2d.f90
+
+% ! structure part
+% cov_model(m_inds(1,1):m_inds(1,2)) = ( sigma_bulk )**2 / da_local_vec(:) * AREA * coverage_str
+% cov_model(m_inds(2,1):m_inds(2,2)) = ( sigma_beta )**2 / da_local_vec(:) * AREA * coverage_str
+%
+% ! source part
+% cov_model(m_inds(3,1):m_inds(3,2)) = sigma_ts**2 * dnparm_src_run * coverage_src
+% cov_model(m_inds(4,1):m_inds(4,2)) = sigma_xs**2 * dnparm_src_run * coverage_src
+% cov_model(m_inds(5,1):m_inds(5,2)) = sigma_zs**2 * dnparm_src_run * coverage_src
+%
+% ! incorporate relative weighting to make the final metric
+% ! STRUCTURE: (fac_str / ugrad_str) * ugrad_str --> fac_str
+% cov_imetric(m_inds(1,1):m_inds(1,2)) = cov_model(m_inds(1,1):m_inds(1,2)) * (fac_str / ugrad_str) * fac_total
+% cov_imetric(m_inds(2,1):m_inds(2,2)) = cov_model(m_inds(2,1):m_inds(2,2)) * (fac_str / ugrad_str) * fac_total
+% cov_imetric(m_inds(3,1):m_inds(3,2)) = cov_model(m_inds(3,1):m_inds(3,2)) * (fac_ts / ugrad_ts) * fac_total
+% cov_imetric(m_inds(4,1):m_inds(4,2)) = cov_model(m_inds(4,1):m_inds(4,2)) * (fac_xs / ugrad_xs) * fac_total
+% cov_imetric(m_inds(5,1):m_inds(5,2)) = cov_model(m_inds(5,1):m_inds(5,2)) * (fac_ys / ugrad_ys) * fac_total
+
+cov_model = zeros(nparm,1);
+cov_model(indC) = ( sigma_bulk )^2 ./ da_local_vec * Atot * coverage_str;
+cov_model(indB) = ( sigma_beta )^2 ./ da_local_vec * Atot * coverage_str;
+cov_model(indT) = sigma_ts^2 * nparm_src_run * coverage_src;
+cov_model(indX) = sigma_xs^2 * nparm_src_run * coverage_src;
+cov_model(indY) = sigma_zs^2 * nparm_src_run * coverage_src;
+
+cov_imetric = zeros(nparm,1);
+cov_imetric(indC) = cov_model(indC) * (fac_str / ugrad_str) * fac_total;
+cov_imetric(indB) = cov_model(indB) * (fac_str / ugrad_str) * fac_total;
+cov_imetric(indT) = cov_model(indT) * (fac_ts / ugrad_ts) * fac_total;
+cov_imetric(indX) = cov_model(indX) * (fac_xs / ugrad_xs) * fac_total;
+cov_imetric(indY) = cov_model(indY) * (fac_ys / ugrad_ys) * fac_total;
+
+% cov_model -- TO CHECK
+%write(19,'(4e20.10)') cov_imetric(i), icov_metric(i), cov_model(i), icov_model(i)
+cov_model_all = load([dir1 'cov_model_diagonal.dat']);
+cov_imetric0 = cov_model_all(:,1);
+% check
+for ii=1:5
+ inds = m_inds(ii,1):m_inds(ii,2);
+ disp('------');
+ disp(sprintf('%.3e %.3e %.3e %.3e',...
+ norm(cov_imetric0(inds)), norm(cov_imetric(inds)),...
+ norm(cov_imetric0(inds)-cov_imetric(inds)),...
+ norm(cov_imetric0(inds)-cov_imetric(inds)) / norm(cov_imetric0(inds))))
+end
+
+% gradient
+%gradient_model(:) = m0(:) / cov_imetric(:)
+gradient = zeros(nparm,1);
+source_gradient = load([dir1 'source_gradient.dat']);
+gradient_model = m0 ./ cov_imetric;
+gradient_data = zeros(nparm,1);
+gradient_data(indC) = 0;
+gradient_data(indB) = summed_ker .* da_local_vec; % NOTE: summed_ker is NOT smoothed
+gradient_data(m_inds(3,1):m_inds(5,2)) = source_gradient;
+gradient = gradient_data + gradient_model;
+
+% gradient -- TO CHECK
+%write(19,'(3e20.10)') gradient(i), gradient_data(i), gradient_model(i)
+[gradient0,gradient_data0,gradient_model0] = textread([dir1 'gradient_vec.dat'],'%f%f%f');
+gradient0(indsTXY) = 0; % no source inversion
+% check
+for ii=1:5
+ inds = m_inds(ii,1):m_inds(ii,2);
+ disp('------');
+ disp(sprintf('%.3e %.3e %.3e %.3e',...
+ norm(gradient0(inds)), norm(gradient(inds)),...
+ norm(gradient0(inds)-gradient(inds)),...
+ norm(gradient0(inds)-gradient(inds)) / norm(gradient0(inds)) ))
+ figure; nr=3; nc=1;
+ subplot(nr,nc,1); plot( gradient0(inds), '.'); title(stpars{ii});
+ subplot(nr,nc,2); plot( gradient(inds), '.')
+ subplot(nr,nc,3); plot( gradient0(inds), gradient(inds), '.')
+end
+
+break
+
+%---------------------------------------------------------
+% emulate CG algorithm in wave2d.f90
+
+
+
+%/home/carltape/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/random_fields/model_001/
+
+%=============================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -114,7 +114,7 @@
! source perturbations (01-Feb-2006)
! presently only set up for 2D surface waves
!logical :: ISOURCE_LOG
- !integer :: PERT_SOURCE, PERT_STRUCT, INV_SOURCE, INV_STRUCT
+ !integer :: PERT_SOURCE, PERT_STRUCT, INV_SOURCE, INV_STRUCT_BETA
logical :: stop_cg
double precision :: src_pert_time, origin_time, origin_time_dat, origin_time_syn, ptmin, ptmax, ptime
double precision :: src_pert_dist, rand1, amin, amax, pangle, rmin, rmax, prad
@@ -170,7 +170,7 @@
character(len=200) :: last_frame_name
double precision, dimension(:,:,:,:,:), allocatable :: absorb_field
double precision, dimension(:,:,:), allocatable :: atype_kernel, btype_kernel, rtype_kernel
- double precision, dimension(NGLLX,NGLLZ,NSPEC) :: atype_kernel_sum, btype_kernel_sum
+ double precision, dimension(NGLLX,NGLLZ,NSPEC) :: btype_kernel_sum
double precision, dimension(NLOCAL) :: ktemp
double precision :: xtemp,ztemp,dx,dz,xmesh,zmesh
@@ -374,6 +374,7 @@
!src_pert_time = 1.0 ! source temporal perturbation (s)
! STANDARD DEVIATIONS of source perturbations (see wave2d_sigmas.m)
+ ! NOTE: These should be READ IN below
src_pert_dist = 2000.0 ! source spatial perturbation (m)
src_pert_time = 0.5 ! source temporal perturbation (s)
@@ -514,23 +515,23 @@
print *
write(*,*) ' PROPERTIES OF THE GRID :'
- write(*,*) ' GLOBAL :'
- write(*,'(a,1f20.10)') ' da_min (m^2) : ', minval(da_global)
- write(*,'(a,1f20.10)') ' da_mean (m^2) : ', da_global_mean
- write(*,'(a,1f20.10)') ' da_max (m^2) : ', maxval(da_global)
+ write(*,*) ' GLOBAL :'
+ write(*,'(a,1f20.10)') ' da_min (m^2) : ', minval(da_global)
+ write(*,'(a,1f20.10)') ' da_mean (m^2) : ', da_global_mean
+ write(*,'(a,1f20.10)') ' da_max (m^2) : ', maxval(da_global)
write(*,*) ' sum [ da ] : ', sum(da_global)
- write(*,*) ' LOCAL :'
- write(*,'(a,1f20.10)') ' da_min (m^2) : ', minval(da_local)
- write(*,'(a,1f20.10)') ' da_mean (m^2) : ', da_local_mean
- write(*,'(a,1f20.10)') ' da_max (m^2) : ', maxval(da_local)
- write(*,*) ' sum [ da ] : ', sum(da_local)
+ write(*,*) ' LOCAL :'
+ write(*,'(a,1f20.10)') ' da_min (m^2) : ', minval(da_local)
+ write(*,'(a,1f20.10)') ' da_mean (m^2) : ', da_local_mean
+ write(*,'(a,1f20.10)') ' da_max (m^2) : ', maxval(da_local)
+ write(*,*) ' sum [ da ] : ', sum(da_local)
print *
- write(*,*) ' LENGTH * HEIGHT : ', LENGTH * HEIGHT
+ write(*,*) ' LENGTH * HEIGHT : ', LENGTH * HEIGHT
print *
- write(*,*) ' GAMMA_SMOOTH : ', GAMMA_SMOOTH
- write(*,*) ' Nfac : ', Nfac
- write(*,*) ' irun0 : ', irun0
- write(*,*) ' IKER : ', IKER
+ write(*,*) ' GAMMA_SMOOTH_KERNEL : ', GAMMA_SMOOTH_KERNEL
+ write(*,*) ' Nfac : ', Nfac
+ write(*,*) ' irun0 : ', irun0
+ write(*,*) ' IKER : ', IKER
print *
print *, ' GLL weights:'
do i = 1,NGLLX
@@ -719,7 +720,7 @@
do ispec = 1,NSPEC
do j = 1,NGLLZ
do i = 1,NGLLX
- read(19,'(7e20.10)') temp1, temp2, temp3, temp4, temp5, temp6, temp7
+ read(19,'(6e20.10)') temp1, temp2, temp3, temp4, temp5, temp6
kappa_syn(i,j,ispec) = temp3
mu_syn(i,j,ispec) = temp4
rho_syn(i,j,ispec) = temp5
@@ -734,7 +735,7 @@
do ispec = 1,NSPEC
do j = 1,NGLLZ
do i = 1,NGLLX
- read(20,'(7e20.10)') temp1, temp2, temp3, temp4, temp5, temp6, temp7
+ read(20,'(6e20.10)') temp1, temp2, temp3, temp4, temp5, temp6
kappa_dat(i,j,ispec) = temp3
mu_dat(i,j,ispec) = temp4
rho_dat(i,j,ispec) = temp5
@@ -746,7 +747,7 @@
beta_syn = sqrt( mu_syn / rho_syn )
beta_dat = sqrt( mu_dat / rho_dat )
- ! read reference values to file
+ ! read reference values from file -- these were used to make the perturbed field
idir = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_INPUT/meshes/mesh_001/'
open(unit=16,file=trim(idir)//'reference_values.dat', status='unknown')
read(16,'(6e20.10)') alpha0, beta0, rho0, bulk0, kappa0, mu0
@@ -782,6 +783,7 @@
endif
! compute alpha0 and beta0 by finding the average of the phase velocity map
+ ! (only if you are not reading in these values)
if ( (IMODEL_SYN /= 3) .or. (IMODEL_DAT /= 3) ) then
print *, 'computing the average of the phase velocity map'
kappa0 = sum(kappa_syn * da_local) / (LENGTH * HEIGHT)
@@ -814,7 +816,7 @@
! if HESSIAN option, then READ in the structure file (local level)
! NOTE: Assume that the mesh and scaling values are IDENTICAL
! to what was used in the base directory for the CG algorithm.
- if (HESSIAN .and. INV_STRUCT == 1) then
+ if (HESSIAN .and. INV_STRUCT_BETA == 1) then
kappa_syn = 0.0 ; mu_syn = 0.0 ; rho_syn = 0.0
alpha_syn = 0.0 ; beta_syn = 0.0 ; bulk_syn = 0.0
@@ -830,7 +832,7 @@
! kappa_syn(i,j,ispec), mu_syn(i,j,ispec), rho_syn(i,j,ispec), &
! temp3, temp4
- read(19,'(7e20.10)') temp1, temp2, temp3, temp4, temp5, temp6, temp7
+ read(19,'(6e20.10)') temp1, temp2, temp3, temp4, temp5, temp6
kappa_syn(i,j,ispec) = temp3
mu_syn(i,j,ispec) = temp4
rho_syn(i,j,ispec) = temp5
@@ -858,13 +860,13 @@
endif
! unperturbed structure (set target map to be the map used to compute synthetics)
- if(PERT_STRUCT == 0) then
+ if(PERT_STRUCT_BETA == 0) then
+ mu_dat = mu_syn
+ beta_dat = beta_syn
kappa_dat = kappa_syn
- mu_dat = mu_syn
+ bulk_dat = bulk_syn
rho_dat = rho_syn
- bulk_dat = bulk_syn
alpha_dat = alpha_syn
- beta_dat = beta_syn
endif
! write velocity structure for data and synthetics to file -- LOCAL LEVEL
@@ -878,14 +880,14 @@
iglob = ibool(i,j,ispec)
! TARGET MODEL (data)
- write(18,'(7e20.10)') x_plot(iglob), z_plot(iglob), &
+ write(18,'(6e20.10)') x_plot(iglob), z_plot(iglob), &
kappa_dat(i,j,ispec), mu_dat(i,j,ispec), rho_dat(i,j,ispec), &
- log( bulk_dat(i,j,ispec) / bulk0 ), log( beta_dat(i,j,ispec) / beta0 )
+ log( beta_dat(i,j,ispec) / beta0 )
! CURRENT MODEL (synthetics)
- write(19,'(7e20.10)') x_plot(iglob), z_plot(iglob), &
+ write(19,'(6e20.10)') x_plot(iglob), z_plot(iglob), &
kappa_syn(i,j,ispec), mu_syn(i,j,ispec), rho_syn(i,j,ispec), &
- log( bulk_syn(i,j,ispec) / bulk0 ), log( beta_syn(i,j,ispec) / beta0 )
+ log( beta_syn(i,j,ispec) / beta0 )
enddo
enddo
enddo
@@ -915,6 +917,11 @@
!!$ close(19)
! write reference structure values to file
+ open(unit=16,file=trim(out_dir2)//'reference_values0.dat', status='unknown')
+ write(16,'(6e20.10)') PWAVESPEED, SWAVESPEED, DENSITY, BWAVESPEED, INCOMPRESSIBILITY, RIGIDITY
+ close(16)
+
+ ! write reference structure values to file
open(unit=16,file=trim(out_dir2)//'reference_values.dat', status='unknown')
write(16,'(6e20.10)') alpha0, beta0, rho0, bulk0, kappa0, mu0
close(16)
@@ -1090,7 +1097,7 @@
! perturbations generated from gji_src_inversion.m
! NOTE: OLD ordering scheme for source indexing
- open(19,file='INPUT/events_xy_pert.dat',status='unknown')
+ open(19,file='INPUT/events_xyt_pert.dat',status='unknown')
do ievent = 1,25
read(19,*) temp1,temp2,temp3
if( PERT_SOURCE_T == 0 ) temp3 = 0.0
@@ -1103,6 +1110,12 @@
enddo
close(19)
+ ! read sigma values used to generate the source parameter perturbations
+ ! NOTE: these over-write what is defined above
+ open(20,file='INPUT/events_xyt_pert_sigmas.dat',status='unknown')
+ read(20,*) temp1,src_pert_dist,src_pert_time
+ close(20)
+
!!$ open(19,file='OUTPUT/run_2550/events_xy.dat',status='unknown')
!!$ do ievent = 1,25
!!$ read(19,'(3f20.10,1i12)') x_eve_dat(ievent), z_eve_dat(ievent), otime_dat(ievent), itemp1
@@ -1549,10 +1562,14 @@
m_inds(:,:) = 0
m_inds(1,1) = 1 ; m_inds(1,2) = NLOCAL
- m_inds(2,1) = NLOCAL+1 ; m_inds(2,2) = nmod_str
- m_inds(3,1) = nmod_str+1 ; m_inds(3,2) = nmod_str+nevent
- m_inds(4,1) = nmod_str+nevent+1 ; m_inds(4,2) = nmod_str+2*nevent
- m_inds(5,1) = nmod_str+2*nevent+1 ; m_inds(5,2) = nmod
+ m_inds(2,1) = nmod_str+1 ; m_inds(2,2) = nmod_str+nevent
+ m_inds(3,1) = nmod_str+nevent+1 ; m_inds(3,2) = nmod_str+2*nevent
+ m_inds(4,1) = nmod_str+2*nevent+1 ; m_inds(4,2) = nmod
+!!$ m_inds(1,1) = 1 ; m_inds(1,2) = NLOCAL
+!!$ m_inds(2,1) = NLOCAL+1 ; m_inds(2,2) = nmod_str
+!!$ m_inds(3,1) = nmod_str+1 ; m_inds(3,2) = nmod_str+nevent
+!!$ m_inds(4,1) = nmod_str+nevent+1 ; m_inds(4,2) = nmod_str+2*nevent
+!!$ m_inds(5,1) = nmod_str+2*nevent+1 ; m_inds(5,2) = nmod
! write source indexing to file
open(unit=18,file=trim(out_dir2)//'m_inds.dat',status='unknown')
@@ -1644,12 +1661,11 @@
! initial model vector (physical coordinates)
! create m0 -- should be identical to what is "un-done" in the CG algorithm
- temp_local1(:,:,:) = log( bulk_syn(:,:,:) / bulk0 )
- temp_local2(:,:,:) = log( beta_syn(:,:,:) / beta0 )
- call local2mvec(temp_local1, temp_local2, nmod_src, m_src_syn, nmod, m0)
+ temp_local1(:,:,:) = log( beta_syn(:,:,:) / beta0 )
+ call local2mvec(temp_local1, nmod_src, m_src_syn, nmod, m0)
! create m0_vec
- call local2mvec(bulk_syn, beta_syn, nmod_src, m_src_syn_vec, nmod, m0_vec)
+ call local2mvec(beta_syn, nmod_src, m_src_syn_vec, nmod, m0_vec)
!m0_prior(:) = m0(:)
m0_vec_initial(:) = m0_vec(:)
@@ -1705,11 +1721,11 @@
! scaling vector for STRUCTURE parameters
sigma_checker_scale = 2.0 ! to mimic a Gaussian distibution
if ((IMODEL_SYN == 3).and.(IMODEL_DAT == 3)) then
- m_scale_str(1) = sigma_bulk0
- m_scale_str(2) = sigma_beta0
+ m_scale_str(1) = sigma_beta0
+ !m_scale_str(2) = sigma_bulk0
else
- m_scale_str(1) = (afac/sigma_checker_scale)/100.0 ! dimensionless : bulk* = ln(bulk/bulk0)
- m_scale_str(2) = (afac/sigma_checker_scale)/100.0 ! dimensionless : beta* = ln(beta/beta0)
+ m_scale_str(1) = (afac/sigma_checker_scale)/100.0 ! dimensionless : beta* = ln(beta/beta0)
+ !m_scale_str(2) = (afac/sigma_checker_scale)/100.0 ! dimensionless : bulk* = ln(bulk/bulk0)
endif
! scaling vector for SOURCE parameters
@@ -1727,8 +1743,8 @@
m_scale_src(3) = src_pert_dist
endif
- sigma_bulk = m_scale_str(1)
- sigma_beta = m_scale_str(2)
+ sigma_beta = m_scale_str(1)
+ !sigma_bulk = m_scale_str(2)
sigma_ts = m_scale_src(1)
sigma_xs = m_scale_src(2)
sigma_zs = m_scale_src(3)
@@ -1737,7 +1753,7 @@
!!$ ! A simple way is to incorporate a scaling factor in the sigma-structure terms.
!!$ ! This will control the relative weight of each set of parameters.
!!$ ! NOTE : We approximate our checkerboard-generated structure values by a Gaussian distribution.
-!!$ if( INV_SOURCE == 1 .and. INV_STRUCT == 1) then
+!!$ if( INV_SOURCE == 1 .and. INV_STRUCT_BETA == 1) then
!!$ if( HESSIAN ) then
!!$ joint_str_src_grad_ratio = 1.0
!!$ else
@@ -1767,35 +1783,40 @@
! (1) unbalanced initial gradient values (gradient_norm_all_unbalanced.dat) -- 3 x 3 checker model
! (2) how much toward the norm should each variable contribute
- if ( INV_STRUCT == 1 .and. INV_SOURCE_T == 0 .and. INV_SOURCE_X == 0 ) then ! structure
+ if ( INV_STRUCT_BETA == 1 .and. INV_SOURCE_T == 0 .and. INV_SOURCE_X == 0 ) then ! structure
fac_str = 1.00
- elseif ( INV_STRUCT == 0 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 0 ) then ! origin time
+ elseif ( INV_STRUCT_BETA == 0 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 0 ) then ! origin time
fac_ts = 1.00
- elseif ( INV_STRUCT == 0 .and. INV_SOURCE_T == 0 .and. INV_SOURCE_X == 1 ) then ! location
+ elseif ( INV_STRUCT_BETA == 0 .and. INV_SOURCE_T == 0 .and. INV_SOURCE_X == 1 ) then ! location
fac_xs = 0.50 ; fac_ys = 0.50
- elseif ( INV_STRUCT == 1 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 0 ) then ! structure + origin time
+ elseif ( INV_STRUCT_BETA == 1 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 0 ) then ! structure + origin time
fac_str = 0.50 ; fac_ts = 0.50
- elseif ( INV_STRUCT == 1 .and. INV_SOURCE_T == 0 .and. INV_SOURCE_X == 1 ) then ! structure + location
+ elseif ( INV_STRUCT_BETA == 1 .and. INV_SOURCE_T == 0 .and. INV_SOURCE_X == 1 ) then ! structure + location
fac_str = 0.50 ; fac_xs = 0.25 ; fac_ys = 0.25
- elseif ( INV_STRUCT == 0 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 1 ) then ! source
+ elseif ( INV_STRUCT_BETA == 0 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 1 ) then ! source
!ugrad_ts = 6.538943677 ; ugrad_xs = 3.288688153 ; ugrad_ys = 3.900641302 ! fac_total = 43.0
!fac_ts = 0.50 ; fac_xs = 0.25 ; fac_ys = 0.25
fac_total = (ugrad_ts/fac_ts) + (ugrad_xs/fac_xs) + (ugrad_ys/fac_ys)
- elseif ( INV_STRUCT == 1 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 1 ) then ! joint
- ! we obtain these values by looking at the initial unbalanced data gradient (gradient_norm_data_all_unbalanced.dat)
+ elseif ( INV_STRUCT_BETA == 1 .and. INV_SOURCE_T == 1 .and. INV_SOURCE_X == 1 ) then ! joint
+ ! NOTE: we obtain these values by looking at the initial unbalanced data gradient (gradient_norm_data_all_unbalanced.dat)
! TEST CASES: 1-25 events, 1-5 events, 5-5 event
!ugrad_str = 0.1296496361d6 ; ugrad_ts = 0.2706999550d4; ugrad_xs = 0.1334372694d4 ; ugrad_ys = 0.1743549898d4 ! Nfac = 2, 25 events
!ugrad_str = 0.2070999127d6 ; ugrad_ts = 0.3507740673d4 ; ugrad_xs = 0.1940368586d4 ; ugrad_ys = 0.2142118702d4 ! Nfac = 3, 25 events
- !ugrad_str = 0.4247330856d6; ugrad_ts = 0.5896527525d4; ugrad_xs = 0.2643538184d4; ugrad_ys = 0.3342815391d4 ! Nfac = 3, 5 events
- ugrad_str = 0.1205146534d6 ; ugrad_ts = 0.2926409454d3; ugrad_xs = 0.9695606936d3; ugrad_ys = 0.7224889563d3 ! Nfac = 3, 1 event
+ !ugrad_str = 0.4247330856d6 ; ugrad_ts = 0.5896527525d4 ; ugrad_xs = 0.2643538184d4 ; ugrad_ys = 0.3342815391d4 ! Nfac = 3, 5 events
+ !ugrad_str = 0.1205146534d6 ; ugrad_ts = 0.2926409454d3 ; ugrad_xs = 0.9695606936d3 ; ugrad_ys = 0.7224889563d3 ! Nfac = 3, 1 event
+ !ugrad_str = ; ugrad_ts = ; ugrad_xs = ; ugrad_ys = ! Gaussians
+ ugrad_str = 0.8251756295d5 ; ugrad_ts = 0.4212723180d3 ; ugrad_xs = 0.5753336089d2 ; ugrad_ys = 0.4077345971d2 ! Gaussians, 1 event
+ !ugrad_str = 0.3812017066d6 ; ugrad_ts = 0.2140202263d4 ; ugrad_xs = 0.2970951189d4 ; ugrad_ys = 0.4102594535d4 ! Gaussians, 5 events
+ !ugrad_str = 0.1907942047d6 ; ugrad_ts = 0.1687351289d4 ; ugrad_xs = 0.2308837150d4 ; ugrad_ys = 0.2986802305d4 ! Gaussians, 25 events
+
fac_str = 1.0/2.0 ; fac_ts = 1.0/6.0 ; fac_xs = 1.0/6.0 ; fac_ys = 1.0/6.0
fac_total = (ugrad_str/fac_str) + (ugrad_ts/fac_ts) + (ugrad_xs/fac_xs) + (ugrad_ys/fac_ys)
endif
@@ -1824,21 +1845,19 @@
!coverage_src = 1.0
! structure part
- cov_model(m_inds(1,1):m_inds(1,2)) = ( sigma_bulk )**2 / da_local_vec(:) * AREA * coverage_str
- cov_model(m_inds(2,1):m_inds(2,2)) = ( sigma_beta )**2 / da_local_vec(:) * AREA * coverage_str
+ cov_model(m_inds(1,1):m_inds(1,2)) = ( sigma_beta )**2 / da_local_vec(:) * AREA * coverage_str
! source part
- cov_model(m_inds(3,1):m_inds(3,2)) = sigma_ts**2 * dnparm_src_run * coverage_src
- cov_model(m_inds(4,1):m_inds(4,2)) = sigma_xs**2 * dnparm_src_run * coverage_src
- cov_model(m_inds(5,1):m_inds(5,2)) = sigma_zs**2 * dnparm_src_run * coverage_src
+ cov_model(m_inds(2,1):m_inds(2,2)) = sigma_ts**2 * dnparm_src_run * coverage_src
+ cov_model(m_inds(3,1):m_inds(3,2)) = sigma_xs**2 * dnparm_src_run * coverage_src
+ cov_model(m_inds(4,1):m_inds(4,2)) = sigma_zs**2 * dnparm_src_run * coverage_src
! incorporate relative weighting to make the final metric
! STRUCTURE: (fac_str / ugrad_str) * ugrad_str --> fac_str
cov_imetric(m_inds(1,1):m_inds(1,2)) = cov_model(m_inds(1,1):m_inds(1,2)) * (fac_str / ugrad_str) * fac_total
- cov_imetric(m_inds(2,1):m_inds(2,2)) = cov_model(m_inds(2,1):m_inds(2,2)) * (fac_str / ugrad_str) * fac_total
- cov_imetric(m_inds(3,1):m_inds(3,2)) = cov_model(m_inds(3,1):m_inds(3,2)) * (fac_ts / ugrad_ts) * fac_total
- cov_imetric(m_inds(4,1):m_inds(4,2)) = cov_model(m_inds(4,1):m_inds(4,2)) * (fac_xs / ugrad_xs) * fac_total
- cov_imetric(m_inds(5,1):m_inds(5,2)) = cov_model(m_inds(5,1):m_inds(5,2)) * (fac_ys / ugrad_ys) * fac_total
+ cov_imetric(m_inds(2,1):m_inds(2,2)) = cov_model(m_inds(2,1):m_inds(2,2)) * (fac_ts / ugrad_ts) * fac_total
+ cov_imetric(m_inds(3,1):m_inds(3,2)) = cov_model(m_inds(3,1):m_inds(3,2)) * (fac_xs / ugrad_xs) * fac_total
+ cov_imetric(m_inds(4,1):m_inds(4,2)) = cov_model(m_inds(4,1):m_inds(4,2)) * (fac_ys / ugrad_ys) * fac_total
! METRIC TENSOR: inverse diagponal covariance matrix
icov_metric = 1.0 / cov_imetric
@@ -1848,9 +1867,8 @@
! create mtarget -- should be identical to what is shown in the CG algorithm
allocate(mtarget(nmod))
mtarget(:) = 0.0
- temp_local1(:,:,:) = log( bulk_dat(:,:,:) / bulk0 )
- temp_local2(:,:,:) = log( beta_dat(:,:,:) / beta0 )
- call local2mvec(temp_local1, temp_local2, nmod_src, m_src_dat, nmod, mtarget)
+ temp_local1(:,:,:) = log( beta_dat(:,:,:) / beta0 )
+ call local2mvec(temp_local1, nmod_src, m_src_dat, nmod, mtarget)
! possible stopping criteria based on the target model
chi_model_stop = 0.5 * model_target_norm
@@ -1861,8 +1879,8 @@
! write model covariance values to file
open(unit=19,file=trim(out_dir2)//'sigma_values.dat',status='unknown')
- write(19,'(2e20.10)') sigma_bulk, m_scale_str(1)
- write(19,'(2e20.10)') sigma_beta, m_scale_str(2)
+ write(19,'(2e20.10)') sigma_beta, m_scale_str(1)
+ !write(19,'(2e20.10)') sigma_bulk, m_scale_str(2)
write(19,'(2e20.10)') sigma_ts, m_scale_src(1)
write(19,'(2e20.10)') sigma_xs, m_scale_src(2)
write(19,'(2e20.10)') sigma_zs, m_scale_src(3)
@@ -1952,8 +1970,6 @@
cov_data(:) = SIGMA_DLNA * SIGMA_DLNA * nmeas_run
endif
- !stop 'Should this not be SIGMA_DT SQUARED?'
-
! write data covariance matrix diagonal to file
open(unit=19,file=trim(out_dir2)//'cov_data_diagonal.dat',status='unknown')
write(19,'(1e20.10)') (cov_data(i), i = 1,nmeas)
@@ -2004,7 +2020,7 @@
irun = irun0 + istep
print *,'=============================================================='
print *,' istep, imod, irun : ', istep, imod, irun
- if(INV_STRUCT==1) print *, ' inverting for structure parameters'
+ if(INV_STRUCT_BETA==1) print *, ' inverting for structure parameters'
if(INV_SOURCE_T==1) print *, ' inverting for source origin times'
if(INV_SOURCE_X==1) print *, ' inverting for source locations'
print *,'=============================================================='
@@ -2025,15 +2041,15 @@
! structure parameters fill the top portion (bulk_syn, beta_syn);
! source parameters fill the bottom portion (m_src_syn_vec)
if(itest==0) then ! reference model (current model)
- call mvec2local(nmod, nmod_src, m0_vec, bulk_syn, beta_syn, m_src_syn_vec)
+ call mvec2local(nmod, nmod_src, m0_vec, beta_syn, m_src_syn_vec)
elseif(itest==1) then ! test model
- call mvec2local(nmod, nmod_src, mt_vec, bulk_syn, beta_syn, m_src_syn_vec)
+ call mvec2local(nmod, nmod_src, mt_vec, beta_syn, m_src_syn_vec)
endif
! for CG algorithm, update kappa_syn and mu_syn
if( .not. HESSIAN) then
- kappa_syn = rho_syn * bulk_syn * bulk_syn
+ !kappa_syn = rho_syn * bulk_syn * bulk_syn
mu_syn = rho_syn * beta_syn * beta_syn
endif
@@ -2126,7 +2142,6 @@
! initialize gradient vectors
gradient(:) = 0.0 ; gradient_data(:) = 0.0 ; gradient_model(:) = 0.0
source_gradient(:) = 0.0
- atype_kernel_sum = 0.0
btype_kernel_sum = 0.0
! initialize misfit
@@ -2627,14 +2642,14 @@
iglob = ibool(i,j,ispec)
! TARGET MODEL (data)
- write(18,'(7e20.10)') x_plot(iglob), z_plot(iglob), &
+ write(18,'(6e20.10)') x_plot(iglob), z_plot(iglob), &
kappa_dat(i,j,ispec), mu_dat(i,j,ispec), rho_dat(i,j,ispec), &
- log( bulk_dat(i,j,ispec) / bulk0 ), log( beta_dat(i,j,ispec) / beta0 )
+ log( beta_dat(i,j,ispec) / beta0 )
! CURRENT MODEL (synthetics)
- write(19,'(7e20.10)') x_plot(iglob), z_plot(iglob), &
+ write(19,'(6e20.10)') x_plot(iglob), z_plot(iglob), &
kappa_syn(i,j,ispec), mu_syn(i,j,ispec), rho_syn(i,j,ispec), &
- log( bulk_syn(i,j,ispec) / bulk0 ), log( beta_syn(i,j,ispec) / beta0 )
+ log( beta_syn(i,j,ispec) / beta0 )
enddo
enddo
enddo
@@ -2960,8 +2975,8 @@
! smooth EACH event kernel for the data subspace method
if(ISMOOTH_EVENT_KERNEL == 1) then
- call smooth_function(btype_kernel, GAMMA_SMOOTH, kbeta_smooth)
- !call smooth_function(atype_kernel, GAMMA_SMOOTH, kbulk_smooth)
+ call smooth_function(btype_kernel, GAMMA_SMOOTH_KERNEL, kbeta_smooth)
+ !call smooth_function(atype_kernel, GAMMA_SMOOTH_KERNEL, kbulk_smooth)
else
kbeta_smooth = btype_kernel
@@ -2982,7 +2997,7 @@
iglob = ibool(i,j,ispec)
write(13,'(3e16.6)') x_plot(iglob), z_plot(iglob), kbeta_smooth(i,j,ispec)
- ktemp(k) = kbeta_smooth(i,j,ispec)
+ ktemp(k) = kbeta_smooth(i,j,ispec) ! used to compute norm of gradient
k = k+1
enddo
enddo
@@ -2993,23 +3008,21 @@
itemp1 = index_source(1,ievent)
itemp2 = index_source(2,ievent)
itemp3 = index_source(3,ievent)
- gradient_data_all(ievent,1) = 0.0
- gradient_data_all(ievent,2) = sum( ktemp**2 * cov_model(m_inds(2,1):m_inds(2,2)) )
- gradient_data_all(ievent,3) = source_gradient(itemp1)**2 * cov_model(itemp1+nmod_str)
- gradient_data_all(ievent,4) = source_gradient(itemp2)**2 * cov_model(itemp2+nmod_str)
- gradient_data_all(ievent,5) = source_gradient(itemp3)**2 * cov_model(itemp3+nmod_str)
+ gradient_data_all(ievent,1) = sum( ktemp**2 * cov_model(m_inds(1,1):m_inds(1,2)) )
+ gradient_data_all(ievent,2) = source_gradient(itemp1)**2 * cov_model(itemp1+nmod_str)
+ gradient_data_all(ievent,3) = source_gradient(itemp2)**2 * cov_model(itemp2+nmod_str)
+ gradient_data_all(ievent,4) = source_gradient(itemp3)**2 * cov_model(itemp3+nmod_str)
! write the gradient parts for each event to file
open(unit=18,file=trim(out_dir)//'gradient_data_unbalanced.dat',status='unknown')
- write(18,'(5e20.10)') gradient_data_all(ievent,1), gradient_data_all(ievent,2), &
- gradient_data_all(ievent,3), gradient_data_all(ievent,4), gradient_data_all(ievent,5)
+ write(18,'(4e20.10)') gradient_data_all(ievent,1), gradient_data_all(ievent,2), &
+ gradient_data_all(ievent,3), gradient_data_all(ievent,4)
close(18)
!!$ endif
- ! summed kernel (= misfit kernel) -- UNSMOOTHED
- atype_kernel_sum = atype_kernel_sum + atype_kernel
- btype_kernel_sum = btype_kernel_sum + btype_kernel
+ ! summed kernel (= misfit kernel)
+ btype_kernel_sum = btype_kernel_sum + kbeta_smooth
deallocate(atype_kernel,btype_kernel,rtype_kernel)
@@ -3125,11 +3138,11 @@
if (itest==0) then
+ ! m_src_syn: w.r.t initial source values (m_src_syn_vec_initial)
m0(:) = 0.0
- temp_local1(:,:,:) = log( bulk_syn(:,:,:) / bulk0 )
- temp_local2(:,:,:) = log( beta_syn(:,:,:) / beta0 )
- m_src_syn(:) = m_src_syn_vec(:) - m0_vec_initial(nmod_str+1 : nmod) ! w.r.t initial source values (m_src_syn_vec_initial)
- call local2mvec(temp_local1, temp_local2, nmod_src, m_src_syn, nmod, m0)
+ temp_local1(:,:,:) = log( beta_syn(:,:,:) / beta0 )
+ m_src_syn(:) = m_src_syn_vec(:) - m0_vec_initial(nmod_str+1 : nmod)
+ call local2mvec(temp_local1, nmod_src, m_src_syn, nmod, m0)
endif
@@ -3142,7 +3155,7 @@
print *, ' Computing the gradient of the misfit function for a given model'
gradient(:) = 0.0 ; gradient_data(:) = 0.0 ; gradient_model(:) = 0.0
- if(INV_STRUCT == 1) then
+ if(INV_STRUCT_BETA == 1) then
! source gradient
open(unit=19,file=trim(out_dir1)//'source_gradient.dat',status='unknown')
@@ -3152,40 +3165,37 @@
! gradient for each event
open(unit=18,file=trim(out_dir)//'gradient_data_all_unbalanced.dat',status='unknown')
do i = 1,nevent
- write(18,'(5e20.10)') gradient_data_all(i,1), gradient_data_all(i,2), &
- gradient_data_all(i,3), gradient_data_all(i,4), gradient_data_all(i,5)
+ write(18,'(4e20.10)') gradient_data_all(i,1), gradient_data_all(i,2), &
+ gradient_data_all(i,3), gradient_data_all(i,4)
enddo
close(18)
! WRITE OUT summed event kernel (= misfit kernel)
- call local2global(atype_kernel_sum, temp_global1)
+ !call local2global(atype_kernel_sum, temp_global1)
call local2global(btype_kernel_sum, temp_global2)
open(19,file=trim(out_dir1)//'summed_ker.dat',status='unknown')
do iglob = 1,NGLOB
- write(19,'(4e16.6)') x_plot(iglob), z_plot(iglob), temp_global1(iglob), temp_global2(iglob)
+ write(19,'(3e16.6)') x_plot(iglob), z_plot(iglob), temp_global2(iglob)
enddo
close(19)
! Smooth the summed event kernel (= misfit kernel) to remove spurious src-rec effects;
! this is done via Gaussian convolution.
if( ISMOOTH_MISFIT_KERNEL == 1) then
- call smooth_function(btype_kernel_sum, GAMMA_SMOOTH, kbeta_smooth)
- !call smooth_function(atype_kernel_sum, GAMMA_SMOOTH, kbulk_smooth)
+ call smooth_function(btype_kernel_sum, GAMMA_SMOOTH_KERNEL, kbeta_smooth)
else
kbeta_smooth = btype_kernel_sum
- !kbulk_smooth = atype_kernel_sum
endif
kbulk_smooth = 0.0 ! testing for S-wave only
- endif ! INV_STRUCT
+ endif ! INV_STRUCT_BETA
!------------------------
! assemble the gradient vector
! NOTE: structure gradient g_i = K_i * dA_i
- temp_local1 = kbulk_smooth * da_local
- temp_local2 = kbeta_smooth * da_local
- call local2mvec(temp_local1, temp_local2, nmod_src, source_gradient, nmod, gradient_data)
+ temp_local1 = kbeta_smooth * da_local
+ call local2mvec(temp_local1, nmod_src, source_gradient, nmod, gradient_data)
! add the gradient term associated with the MODEL NORM in the misfit function
if( INCLUDE_MODEL_NORM ) then
@@ -3197,17 +3207,17 @@
endif
! if NOT inverting for structure or source, then set those parts of the gradient to zero
- if(INV_STRUCT == 0) then
- gradient_data(m_inds(1,1) : m_inds(2,2)) = 0.0
- gradient_model(m_inds(1,1) : m_inds(2,2)) = 0.0
+ if(INV_STRUCT_BETA == 0) then
+ gradient_data(m_inds(1,1) : m_inds(1,2)) = 0.0
+ gradient_model(m_inds(1,1) : m_inds(1,2)) = 0.0
endif
if(INV_SOURCE_T == 0) then
- gradient_data(m_inds(3,1) : m_inds(3,2)) = 0.0
- gradient_model(m_inds(3,1) : m_inds(3,2)) = 0.0
+ gradient_data(m_inds(2,1) : m_inds(2,2)) = 0.0
+ gradient_model(m_inds(2,1) : m_inds(2,2)) = 0.0
endif
if(INV_SOURCE_X == 0) then
- gradient_data( m_inds(4,1) : m_inds(5,2)) = 0.0
- gradient_model(m_inds(4,1) : m_inds(5,2)) = 0.0
+ gradient_data( m_inds(3,1) : m_inds(4,2)) = 0.0
+ gradient_model(m_inds(3,1) : m_inds(4,2)) = 0.0
endif
! overall gradient
@@ -3246,7 +3256,7 @@
call write_norm_sq(filename2,gradient_data_norm_parts)
call write_norm_sq(filename3,gradient_model_norm_parts)
- ! SAME AS ABOVE, only use the unbalanced icov_model instead of icov_metric
+ ! SAME AS ABOVE, but use the unbalanced icov_model instead of icov_metric
call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, gradient, icov_model, &
gradient_norm, gradient_norm_struct, gradient_norm_source, gradient_norm_parts)
call compute_norm_sq(ievent_min, ievent_max, nevent, index_source, nmod, gradient_data, icov_model, &
@@ -3262,56 +3272,56 @@
call write_norm_sq(filename2,gradient_data_norm_parts)
call write_norm_sq(filename3,gradient_model_norm_parts)
- ! scaling for joint inversions
- if( istep==0 .and. 0==1 ) then
+!!$ ! scaling for joint inversions
+!!$ if( istep==0 .and. 0==1 ) then
+!!$
+!!$ ! scaling for joint source-structure inversions
+!!$
+!!$ ! Formulas used for GJI paper
+!!$ !norm_bulk_2 = sum( kbulk_smooth**2 * da_local )
+!!$ !norm_beta_2 = sum( kbeta_smooth**2 * da_local )
+!!$ !norm_source_2 = sum( ( source_gradient(:) * m_scale_src_all(:) )**2 )
+!!$
+!!$ ! NORM-SQUARED of the steepest ascent vector : sqrt( ghat C ghat )
+!!$ norm_bulk_2 = sum( kbulk_smooth(:,:,:)**2 * da_local(:,:,:) * (m_scale_str(1))**2 * AREA )
+!!$ norm_beta_2 = sum( kbeta_smooth(:,:,:)**2 * da_local(:,:,:) * (m_scale_str(2))**2 * AREA )
+!!$ norm_source_2 = sum( ( source_gradient(:) * m_scale_src_all(:) )**2 * dnparm_src_run ) ! m_scale_src_all NOT ALLOCATED
+!!$
+!!$ ! This shows the balance used in the GJI-2007 paper
+!!$ scale_struct_gji = sqrt( norm_source_2 / (norm_bulk_2 + norm_beta_2) )
+!!$
+!!$ !scale_struct = scale_struct_gji
+!!$ scale_struct = 1.0
+!!$
+!!$ write(*,*) 'Scaling for joint inversions:'
+!!$ write(*,'(a24,1e14.4)') ' norm_bulk : ', sqrt( norm_bulk_2 )
+!!$ write(*,'(a24,1e14.4)') ' norm_beta : ', sqrt( norm_beta_2 )
+!!$ write(*,'(a24,1e14.4)') ' norm_source : ', sqrt( norm_source_2 )
+!!$ write(*,'(a24,1e14.4)') ' scale_struct_gji (F) : ', scale_struct_gji
+!!$ write(*,'(a24,1e14.4)') ' scale_struct (F) : ', scale_struct
+!!$
+!!$ open(unit=19,file=trim(out_dir1)//'scaling_values.dat',status='unknown')
+!!$ write(19,*) 'GJI_PAPER = ', GJI_PAPER
+!!$ write(19,*) 'IRUNZ = ', IRUNZ
+!!$ write(19,'(6i14)') 1, NLOCAL, NLOCAL+1, nmod_str, nmod_str+1, nmod
+!!$ write(19,'(2i14)') nparm_src_run, NLOCAL
+!!$ write(19,'(5e14.6)') sum(da_local(:,:,:)), sum(da_local_vec(:)), sum(da_global(:)), LENGTH*HEIGHT, AREA
+!!$ write(19,'(3e14.6)') minval(da_local_vec(:)), maxval(da_local_vec(:)), sum(da_local_vec(:))/NLOCAL
+!!$ write(19,'(3e14.6)') minval(da_global(:)), maxval(da_global(:)), sum(da_global(:))/NGLOB
+!!$ write(19,'(a20,2e14.6)') 'cov_str_fac : ', cov_str_fac, cov_str_fac**2
+!!$ write(19,'(a20,2e14.6)') 'cov_src_fac : ', cov_src_fac, cov_src_fac**2
+!!$ write(19,*) 'SIGMA VALUES :'
+!!$ write(19,'(5e14.6)') m_scale_str(1), m_scale_str(2), m_scale_src(1), m_scale_src(2), m_scale_src(3)
+!!$ write(19,'(5e14.6)') (m_scale_str(1))**2, (m_scale_str(2))**2, &
+!!$ (m_scale_src(1))**2, (m_scale_src(2))**2, (m_scale_src(3))**2
+!!$ write(19,'(5e14.6)') sigma_bulk, sigma_beta, sigma_xs, sigma_zs, sigma_ts
+!!$ write(19,'(5e14.6)') (sigma_bulk)**2, (sigma_beta)**2, (sigma_xs)**2, (sigma_zs)**2, (sigma_ts)**2
+!!$ write(19,'(a20,4e14.6)') 'scale_struct : ', scale_struct, scale_struct**2, 1.0/scale_struct, (1.0/scale_struct)**2
+!!$ write(19,'(a20,4e14.6)') 'scale_struct_gji : ', scale_struct_gji , scale_struct_gji **2, &
+!!$ 1.0/scale_struct_gji , (1.0/scale_struct_gji )**2
+!!$ close(19)
+!!$ endif
- ! scaling for joint source-structure inversions
-
- ! Formulas used for GJI paper
- !norm_bulk_2 = sum( kbulk_smooth**2 * da_local )
- !norm_beta_2 = sum( kbeta_smooth**2 * da_local )
- !norm_source_2 = sum( ( source_gradient(:) * m_scale_src_all(:) )**2 )
-
- ! NORM-SQUARED of the steepest ascent vector : sqrt( ghat C ghat )
- norm_bulk_2 = sum( kbulk_smooth(:,:,:)**2 * da_local(:,:,:) * (m_scale_str(1))**2 * AREA )
- norm_beta_2 = sum( kbeta_smooth(:,:,:)**2 * da_local(:,:,:) * (m_scale_str(2))**2 * AREA )
- norm_source_2 = sum( ( source_gradient(:) * m_scale_src_all(:) )**2 * dnparm_src_run ) ! m_scale_src_all NOT ALLOCATED
-
- ! This shows the balance used in the GJI-2007 paper
- scale_struct_gji = sqrt( norm_source_2 / (norm_bulk_2 + norm_beta_2) )
-
- !scale_struct = scale_struct_gji
- scale_struct = 1.0
-
- write(*,*) 'Scaling for joint inversions:'
- write(*,'(a24,1e14.4)') ' norm_bulk : ', sqrt( norm_bulk_2 )
- write(*,'(a24,1e14.4)') ' norm_beta : ', sqrt( norm_beta_2 )
- write(*,'(a24,1e14.4)') ' norm_source : ', sqrt( norm_source_2 )
- write(*,'(a24,1e14.4)') ' scale_struct_gji (F) : ', scale_struct_gji
- write(*,'(a24,1e14.4)') ' scale_struct (F) : ', scale_struct
-
- open(unit=19,file=trim(out_dir1)//'scaling_values.dat',status='unknown')
- write(19,*) 'GJI_PAPER = ', GJI_PAPER
- write(19,*) 'IRUNZ = ', IRUNZ
- write(19,'(6i14)') 1, NLOCAL, NLOCAL+1, nmod_str, nmod_str+1, nmod
- write(19,'(2i14)') nparm_src_run, NLOCAL
- write(19,'(5e14.6)') sum(da_local(:,:,:)), sum(da_local_vec(:)), sum(da_global(:)), LENGTH*HEIGHT, AREA
- write(19,'(3e14.6)') minval(da_local_vec(:)), maxval(da_local_vec(:)), sum(da_local_vec(:))/NLOCAL
- write(19,'(3e14.6)') minval(da_global(:)), maxval(da_global(:)), sum(da_global(:))/NGLOB
- write(19,'(a20,2e14.6)') 'cov_str_fac : ', cov_str_fac, cov_str_fac**2
- write(19,'(a20,2e14.6)') 'cov_src_fac : ', cov_src_fac, cov_src_fac**2
- write(19,*) 'SIGMA VALUES :'
- write(19,'(5e14.6)') m_scale_str(1), m_scale_str(2), m_scale_src(1), m_scale_src(2), m_scale_src(3)
- write(19,'(5e14.6)') (m_scale_str(1))**2, (m_scale_str(2))**2, &
- (m_scale_src(1))**2, (m_scale_src(2))**2, (m_scale_src(3))**2
- write(19,'(5e14.6)') sigma_bulk, sigma_beta, sigma_xs, sigma_zs, sigma_ts
- write(19,'(5e14.6)') (sigma_bulk)**2, (sigma_beta)**2, (sigma_xs)**2, (sigma_zs)**2, (sigma_ts)**2
- write(19,'(a20,4e14.6)') 'scale_struct : ', scale_struct, scale_struct**2, 1.0/scale_struct, (1.0/scale_struct)**2
- write(19,'(a20,4e14.6)') 'scale_struct_gji : ', scale_struct_gji , scale_struct_gji **2, &
- 1.0/scale_struct_gji , (1.0/scale_struct_gji )**2
- close(19)
- endif
-
endif ! itest==0 .or. POLY_ORDER==3
!====================
@@ -3378,11 +3388,12 @@
!!$ enddo
! STRUCTURE PARAMETERS (bulk, then beta)
- if(INV_STRUCT == 0) then
+ if(INV_STRUCT_BETA == 0) then
mt_vec(1 : nmod_str) = m0_vec_initial(1 : nmod_str) ! initial structure
else
- mt_vec(1 : NLOCAL) = bulk0 * exp( mt(1 : NLOCAL) )
- mt_vec(NLOCAL+1 : nmod_str) = beta0 * exp( mt(NLOCAL+1 : nmod_str) )
+ mt_vec(1 : NLOCAL) = beta0 * exp( mt(1 : NLOCAL) )
+ !mt_vec(1 : NLOCAL) = bulk0 * exp( mt(1 : NLOCAL) )
+ !mt_vec(NLOCAL+1 : nmod_str) = beta0 * exp( mt(NLOCAL+1 : nmod_str) )
endif
! SOURCE PARAMETERS
@@ -3488,11 +3499,12 @@
!!$ enddo
! STRUCTURE PARAMETERS
- if(INV_STRUCT == 0) then
+ if(INV_STRUCT_BETA == 0) then
m0_vec(1 : nmod_str) = m0_vec_initial(1 : nmod_str) ! initial structure
else
- m0_vec(1 : NLOCAL) = bulk0 * exp( m0(1 : NLOCAL) )
- m0_vec(NLOCAL+1 : nmod_str) = beta0 * exp( m0(NLOCAL+1 : nmod_str) )
+ m0_vec(1 : NLOCAL) = beta0 * exp( m0(1 : NLOCAL) )
+ !m0_vec(1 : NLOCAL) = bulk0 * exp( m0(1 : NLOCAL) )
+ !m0_vec(NLOCAL+1 : nmod_str) = beta0 * exp( m0(NLOCAL+1 : nmod_str) )
endif
! SOURCE PARAMETERS
@@ -3525,9 +3537,9 @@
close(19)
endif
- ! NOTE: if the structure parameters are unrealistic, then we will exit
- ! the program before the NEXT wave simulation. This way, we at least write
- ! the unrealistic model to file.
+ ! FUTURE WORK: if the structure or source parameters are unrealistic,
+ ! then we should exit the program before the NEXT wave simulation.
+ ! This way, we at least write the unrealistic model to file.
!====================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -1,28 +1,28 @@
module wave2d_constants
-! This file is copied into the output directory when wave2d.f90 is run
-! There are two basic options: membrane (surface) waves and body waves.
-! Several lines need to be commented/uncommented to go from one to the other.
+ ! This file is copied into the output directory when wave2d.f90 is run
+ ! There are two basic options: membrane (surface) waves and body waves.
+ ! Several lines need to be commented/uncommented to go from one to the other.
-! index of the reference directory for the simulation output
- integer, parameter :: IRUNZ = 0
+ ! index of the reference directory for the simulation output
+ integer, parameter :: IRUNZ = 700
-!========================
-! GRID, TIME-STEP, AND SOURCE PARAMETERS
+ !========================
+ ! GRID, TIME-STEP, AND SOURCE PARAMETERS
-! NFRAME : number of frames to save: membrane, 10; body, 20
-! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
-! NSTEP : number of timesteps
+ ! NFRAME : number of frames to save: membrane, 10; body, 20
+ ! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+ ! NSTEP : number of timesteps
integer, parameter :: NFRAME = 10 ! membrane
!integer, parameter :: NFRAME = 11 ! body
integer, parameter :: NSAVE = 400 ! 200,400
integer, parameter :: NSTEP = NFRAME*NSAVE
-! time step in seconds
+ ! time step in seconds
!double precision, parameter :: DT = 2.0e-2 ! body waves
double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
-! temporal properties of source (source time function)
+ ! temporal properties of source (source time function)
integer, parameter :: ISRC_TIME = 1 ! type (1)
double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
!double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
@@ -30,49 +30,49 @@
!double precision, parameter :: tshift = 8.0*hdur
logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
-! normalization factor of point source force
+ ! normalization factor of point source force
double precision, parameter :: FNORM = 1.0d10
-! forward wavefield source-time function
-! (For membrane waves, FOR_Y = 1 is all that matters.)
+ ! forward wavefield source-time function
+ ! (For membrane waves, FOR_Y = 1 is all that matters.)
integer, parameter :: FOR_X = 1 ! boolean
integer, parameter :: FOR_Y = 1 ! boolean
integer, parameter :: FOR_Z = 1 ! boolean
-! adjoint wavefield source-time function
-! (For membrane waves, REV_Y = 1 is all that matters.)
-! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ ! adjoint wavefield source-time function
+ ! (For membrane waves, REV_Y = 1 is all that matters.)
+ ! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
integer, parameter :: REV_X = 1 ! boolean
integer, parameter :: REV_Y = 1 ! boolean
integer, parameter :: REV_Z = 1 ! boolean
-! spatial properties of sources
-! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
-! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
-! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ ! spatial properties of sources
+ ! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+ ! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+ ! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
-! spatial properties of receivers
-! IREC_SPACE
-! (1) individual station(s)
-! (2) SoCal (used for GJI paper)
-! (3) regular mesh on land
-! (4) regular mesh
-! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
-! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
-! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ ! spatial properties of receivers
+ ! IREC_SPACE
+ ! (1) individual station(s)
+ ! (2) SoCal (used for GJI paper)
+ ! (3) regular mesh on land
+ ! (4) regular mesh
+ ! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+ ! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+ ! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
integer, parameter :: NMESH_REC = 17
- double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
- double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m
double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
-! lower right corner for membrane surface waves plotting grid
+ ! lower right corner for membrane surface waves plotting grid
double precision, parameter :: LAT_MIN = 32.0
double precision, parameter :: LON_MIN = -120.0
integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
-! mesh specifications: membrane surface waves
+ ! mesh specifications: membrane surface waves
double precision, parameter :: LENGTH = 480.0e3 ! m
double precision, parameter :: HEIGHT = 480.0e3 ! m
double precision, parameter :: AREA = LENGTH*HEIGHT
@@ -85,12 +85,12 @@
!!$ integer, parameter :: NEX = 80 ! 160
!!$ integer, parameter :: NEZ = 32 ! 32
-!========================
-! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+ !========================
+ ! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
-! model perturbations for HOMOGENEOUS model (or perturbation)
-! scaling from beta to alpha
-! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ ! model perturbations for HOMOGENEOUS model (or perturbation)
+ ! scaling from beta to alpha
+ ! value is from Masters et al. (2000), "The relative behavior of shear velocity..."
double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
double precision, parameter :: PBETA = 10.0
double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
@@ -103,61 +103,71 @@
! 0 1 2 3
!-----------------------------------------------------------------------------------------
! ISURFACE=1, IMODEL_SYN : homo checker het
- ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het
!-----------------------------------------------------------------------------------------
- ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
- ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker NA
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert NA
!-----------------------------------------------------------------------------------------
- ! whether you want to smooth the structure model for computing the synthetics
+ ! smooth the structure model or kernels
+ ! For the CG algorithm, ISMOOTH_EVENT_KERNEL --> use smoothed event kernels
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 0
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
- integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
- integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
- ! smoothing scalelength
- double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+ ! scalelength of smoothing Gaussian
+ ! GJI paper : 30,60,90
+ ! body waves : 10
+ double precision, parameter :: SIGMA_SMOOTH_KERNEL = 2.121320d4
+ double precision, parameter :: SIGMA_SMOOTH_MODEL = 1.060660d4
+ double precision, parameter :: GAMMA_SMOOTH_KERNEL = sqrt(8.0)*SIGMA_SMOOTH_KERNEL
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = sqrt(8.0)*SIGMA_SMOOTH_MODEL
-!========================
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
-! boolean parameters
-! IKER: (0) waveform
-! (1) traveltime, cross-correlation, misfit
-! (2) amplitude, cross-correlation, misfit
-! (3) traveltime, multitaper
-! (4) amplitude, multitaper
-! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
-! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ !========================
+
+ ! boolean parameters
+ ! IKER: (0) waveform
+ ! (1) traveltime, cross-correlation, misfit
+ ! (2) amplitude, cross-correlation, misfit
+ ! (3) traveltime, multitaper
+ ! (4) amplitude, multitaper
+ ! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+ ! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
integer, parameter :: IKER = 1
integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
-! KEY: USE ONE LINE OR THE OTHER
+ ! KEY: USE ONE LINE OR THE OTHER
integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
-! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+ ! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
-! parameters controlling what to write to file
-! NOTE: for the tomography simulations, ALL of these can be .false.
+ ! parameters controlling what to write to file
+ ! NOTE: for the tomography simulations, ALL of these can be .false.
logical, parameter :: WRITE_STF_F = .false.
logical, parameter :: WRITE_SEISMO_F = .false.
-! logical, parameter :: WRITE_SPECTRA_F = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+ ! logical, parameter :: WRITE_SPECTRA_F = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
-! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+ ! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
logical, parameter :: WRITE_STF_A = .false.
logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
-! logical, parameter :: WRITE_SPECTRA_A = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+ ! logical, parameter :: WRITE_SPECTRA_A = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
-
-!--------------------------------------
-! INVERSION PARAMETERS
+ !--------------------------------------
+ ! INVERSION PARAMETERS
+
! whether you want to compute kernels, or simply the misfit function
logical, parameter :: COMPUTE_KERNELS = .true.
@@ -169,8 +179,10 @@
! VAR_RED_MIN : minimum variance reduction (in percent)
! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
- integer, parameter :: NITERATION = 1
+ integer, parameter :: NITERATION = 16
double precision, parameter :: VAR_RED_MIN = 8.0
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
! Gaussian errors containted in input file
! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
@@ -178,9 +190,6 @@
double precision, parameter :: SIGMA_DLNA = 1.0
double precision, parameter :: SIGMA_WAVEFORM = 1.0
logical, parameter :: ADD_DATA_ERRORS = .true.
-
- !double precision, parameter :: SIGMA_FAC = 2.0
- !double precision, parameter :: CONV_STOP = 1.0e-4
! order of interpolating polynomial in conjugate gradient algorithm
! using POLY_ORDER = 3 required computing the gradient of the test model
@@ -190,114 +199,115 @@
!integer, parameter :: GJI_PAPER = 1
! what to perturb, what to invert
- integer, parameter :: PERT_STRUCT = 1
- integer, parameter :: PERT_SOURCE_T = 0
- integer, parameter :: PERT_SOURCE_X = 0
+ ! (For the inverse tests, we only allow perturbations in beta.)
+ integer, parameter :: PERT_STRUCT_BETA = 1
+ integer, parameter :: PERT_SOURCE_T = 1
+ integer, parameter :: PERT_SOURCE_X = 1
- integer, parameter :: INV_STRUCT = 1
- integer, parameter :: INV_SOURCE_T = 0
- integer, parameter :: INV_SOURCE_X = 0
+ integer, parameter :: INV_STRUCT_BETA = 1
+ integer, parameter :: INV_SOURCE_T = 1
+ integer, parameter :: INV_SOURCE_X = 1
! whether to include the model norm term in the misfit function, which acts like damping
logical, parameter :: INCLUDE_MODEL_NORM = .true.
-
+
! log file showing source loactions
logical, parameter :: ISOURCE_LOG = .true.
-
+
! DO NOT CHANGE THESE
- integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ ! integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_STRUCT = 1 ! beta only (24-Jan-2010)
integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
-! parameterization for structure in the inversion
-! 1 : kappa-mu-rho
-! 2 : alpha-beta-rho
-! 3 : c-beta-rho
+ ! parameterization for structure in the inversion
+ ! 1 : kappa-mu-rho
+ ! 2 : alpha-beta-rho
+ ! 3 : c-beta-rho
integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
-! scalelength of smoothing Gaussian
-! GJI paper : 30,60,90
-! body waves : 10
- double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+ ! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 4.50e10 ! Pa
+ double precision, parameter :: RIGIDITY = 3.185e10 ! Pa
+ ! double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ ! double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
- ! parameters for smoothing
- logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
- logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+ ! compute additional parameters
+ double precision, parameter :: PWAVESPEED = &
+ sqrt( (INCOMPRESSIBILITY + (4.0/3.0)*RIGIDITY)/DENSITY )
+ double precision, parameter :: SWAVESPEED = sqrt( RIGIDITY/DENSITY )
+ double precision, parameter :: BWAVESPEED = sqrt( INCOMPRESSIBILITY/DENSITY )
-! homogeneous background model (S.I. units)
- double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
- double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
- double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+ !---------------------------------------------------------------
+ ! measurement windows
-!---------------------------------------------------------------
-! measurement windows
-
double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
-!---------------------------------------------------------------
-! CHT: do not change these
+ !---------------------------------------------------------------
+ ! CHT: do not change these
-! UTM zone for Southern California region
-! integer, parameter :: UTM_PROJECTION_ZONE = 11
+ ! UTM zone for Southern California region
+ ! integer, parameter :: UTM_PROJECTION_ZONE = 11
-! to suppress UTM projection for SCEC benchmarks
+ ! to suppress UTM projection for SCEC benchmarks
logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
-! flag for projection from latitude/longitude to UTM, and back
+ ! flag for projection from latitude/longitude to UTM, and back
integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
-! flag for projection from latitude/longitude to mesh-UTM, and back
+ ! flag for projection from latitude/longitude to mesh-UTM, and back
integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
-! max number of fake receivers
+ ! max number of fake receivers
integer, parameter :: MAX_SR_FAKE = 1000
-! max number of events, receivers, and phass
+ ! max number of events, receivers, and phass
integer, parameter :: MAX_EVENT = 50
integer, parameter :: MAX_SR = 1400
integer, parameter :: MAX_PHASE = 1
integer, parameter :: MAX_COMP = NCOMP
-!========================
-! GRID AND GLL POINTS
+ !========================
+ ! GRID AND GLL POINTS
integer, parameter :: NELE = MAX(NEX,NEZ)
integer, parameter :: NSPEC = NEX*NEZ
-! number of GLL points (polynomial degree plus one)
+ ! number of GLL points (polynomial degree plus one)
integer, parameter :: NGLLX = 5
integer, parameter :: NGLLZ = 5
integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
-! number of global points
+ ! number of global points
integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
-! number of nodes for 2D and 3D shape functions for hexahedra
-! we use 8-node mesh bricks, which are more stable than 27-node elements
+ ! number of nodes for 2D and 3D shape functions for hexahedra
+ ! we use 8-node mesh bricks, which are more stable than 27-node elements
integer, parameter :: NGNOD = 8, NGNOD2D = 4
-
-! number of iterations to solve the system for xi and eta
+
+ ! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 1
-! very large and very small values
+ ! very large and very small values
double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
-! for the Gauss-Lobatto-Legendre points and weights
+ ! for the Gauss-Lobatto-Legendre points and weights
double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
-!
-! CONSTANTS
-!
+ !
+ ! CONSTANTS
+ !
double precision, parameter :: PI = 3.141592653589793
double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
double precision, parameter :: ONE_THIRD = 1.0 / 3.0
double precision, parameter :: ONEOVERTWO = 0.5
double precision, parameter :: DEG = 180.0/PI
-! double precision, parameter :: EPS = 1.0e-35
+ ! double precision, parameter :: EPS = 1.0e-35
!!$! parameter for FFTW
!!$ integer, parameter :: NOUT = NSTEP/2 + 1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -1,28 +1,28 @@
module wave2d_constants
-! This file is copied into the output directory when wave2d.f90 is run
-! There are two basic options: membrane (surface) waves and body waves.
-! Several lines need to be commented/uncommented to go from one to the other.
+ ! This file is copied into the output directory when wave2d.f90 is run
+ ! There are two basic options: membrane (surface) waves and body waves.
+ ! Several lines need to be commented/uncommented to go from one to the other.
-! index of the reference directory for the simulation output
+ ! index of the reference directory for the simulation output
integer, parameter :: IRUNZ = 0
-!========================
-! GRID, TIME-STEP, AND SOURCE PARAMETERS
+ !========================
+ ! GRID, TIME-STEP, AND SOURCE PARAMETERS
-! NFRAME : number of frames to save: membrane, 10; body, 20
-! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
-! NSTEP : number of timesteps
+ ! NFRAME : number of frames to save: membrane, 10; body, 20
+ ! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+ ! NSTEP : number of timesteps
integer, parameter :: NFRAME = 10 ! membrane
!integer, parameter :: NFRAME = 11 ! body
integer, parameter :: NSAVE = 400 ! 200,400
integer, parameter :: NSTEP = NFRAME*NSAVE
-! time step in seconds
+ ! time step in seconds
!double precision, parameter :: DT = 2.0e-2 ! body waves
double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
-! temporal properties of source (source time function)
+ ! temporal properties of source (source time function)
integer, parameter :: ISRC_TIME = 1 ! type (1)
double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
!double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
@@ -30,49 +30,49 @@
!double precision, parameter :: tshift = 8.0*hdur
logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
-! normalization factor of point source force
+ ! normalization factor of point source force
double precision, parameter :: FNORM = 1.0d10
-! forward wavefield source-time function
-! (For membrane waves, FOR_Y = 1 is all that matters.)
+ ! forward wavefield source-time function
+ ! (For membrane waves, FOR_Y = 1 is all that matters.)
integer, parameter :: FOR_X = 1 ! boolean
integer, parameter :: FOR_Y = 1 ! boolean
integer, parameter :: FOR_Z = 1 ! boolean
-! adjoint wavefield source-time function
-! (For membrane waves, REV_Y = 1 is all that matters.)
-! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ ! adjoint wavefield source-time function
+ ! (For membrane waves, REV_Y = 1 is all that matters.)
+ ! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
integer, parameter :: REV_X = 1 ! boolean
integer, parameter :: REV_Y = 1 ! boolean
integer, parameter :: REV_Z = 1 ! boolean
-! spatial properties of sources
-! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
-! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
-! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ ! spatial properties of sources
+ ! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+ ! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+ ! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
-! spatial properties of receivers
-! IREC_SPACE
-! (1) individual station(s)
-! (2) SoCal (used for GJI paper)
-! (3) regular mesh on land
-! (4) regular mesh
-! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
-! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
-! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ ! spatial properties of receivers
+ ! IREC_SPACE
+ ! (1) individual station(s)
+ ! (2) SoCal (used for GJI paper)
+ ! (3) regular mesh on land
+ ! (4) regular mesh
+ ! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+ ! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+ ! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
integer, parameter :: NMESH_REC = 17
- double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
- double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m
double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
-! lower right corner for membrane surface waves plotting grid
+ ! lower right corner for membrane surface waves plotting grid
double precision, parameter :: LAT_MIN = 32.0
double precision, parameter :: LON_MIN = -120.0
integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
-! mesh specifications: membrane surface waves
+ ! mesh specifications: membrane surface waves
double precision, parameter :: LENGTH = 480.0e3 ! m
double precision, parameter :: HEIGHT = 480.0e3 ! m
double precision, parameter :: AREA = LENGTH*HEIGHT
@@ -85,12 +85,12 @@
!!$ integer, parameter :: NEX = 80 ! 160
!!$ integer, parameter :: NEZ = 32 ! 32
-!========================
-! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+ !========================
+ ! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
-! model perturbations for HOMOGENEOUS model (or perturbation)
-! scaling from beta to alpha
-! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ ! model perturbations for HOMOGENEOUS model (or perturbation)
+ ! scaling from beta to alpha
+ ! value is from Masters et al. (2000), "The relative behavior of shear velocity..."
double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
double precision, parameter :: PBETA = 10.0
double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
@@ -103,61 +103,71 @@
! 0 1 2 3
!-----------------------------------------------------------------------------------------
! ISURFACE=1, IMODEL_SYN : homo checker het
- ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het
!-----------------------------------------------------------------------------------------
- ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
- ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker NA
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert NA
!-----------------------------------------------------------------------------------------
- ! whether you want to smooth the structure model for computing the synthetics
+ ! smooth the structure model or kernels
+ ! For the CG algorithm, ISMOOTH_EVENT_KERNEL --> use smoothed event kernels
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 0
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
- integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
- integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
- ! smoothing scalelength
- double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+ ! scalelength of smoothing Gaussian
+ ! GJI paper : 30,60,90
+ ! body waves : 10
+ double precision, parameter :: SIGMA_SMOOTH_KERNEL = 2.121320d4
+ double precision, parameter :: SIGMA_SMOOTH_MODEL = 1.060660d4
+ double precision, parameter :: GAMMA_SMOOTH_KERNEL = sqrt(8.0)*SIGMA_SMOOTH_KERNEL
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = sqrt(8.0)*SIGMA_SMOOTH_MODEL
-!========================
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
-! boolean parameters
-! IKER: (0) waveform
-! (1) traveltime, cross-correlation, misfit
-! (2) amplitude, cross-correlation, misfit
-! (3) traveltime, multitaper
-! (4) amplitude, multitaper
-! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
-! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ !========================
+
+ ! boolean parameters
+ ! IKER: (0) waveform
+ ! (1) traveltime, cross-correlation, misfit
+ ! (2) amplitude, cross-correlation, misfit
+ ! (3) traveltime, multitaper
+ ! (4) amplitude, multitaper
+ ! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+ ! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
integer, parameter :: IKER = 1
integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
-! KEY: USE ONE LINE OR THE OTHER
+ ! KEY: USE ONE LINE OR THE OTHER
integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
-! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+ ! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
-! parameters controlling what to write to file
-! NOTE: for the tomography simulations, ALL of these can be .false.
+ ! parameters controlling what to write to file
+ ! NOTE: for the tomography simulations, ALL of these can be .false.
logical, parameter :: WRITE_STF_F = .false.
logical, parameter :: WRITE_SEISMO_F = .false.
-! logical, parameter :: WRITE_SPECTRA_F = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+ ! logical, parameter :: WRITE_SPECTRA_F = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
-! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+ ! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
logical, parameter :: WRITE_STF_A = .false.
logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
-! logical, parameter :: WRITE_SPECTRA_A = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+ ! logical, parameter :: WRITE_SPECTRA_A = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
-
-!--------------------------------------
-! INVERSION PARAMETERS
+ !--------------------------------------
+ ! INVERSION PARAMETERS
+
! whether you want to compute kernels, or simply the misfit function
logical, parameter :: COMPUTE_KERNELS = .true.
@@ -171,6 +181,8 @@
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
integer, parameter :: NITERATION = 1
double precision, parameter :: VAR_RED_MIN = 8.0
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
! Gaussian errors containted in input file
! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
@@ -178,9 +190,6 @@
double precision, parameter :: SIGMA_DLNA = 1.0
double precision, parameter :: SIGMA_WAVEFORM = 1.0
logical, parameter :: ADD_DATA_ERRORS = .true.
-
- !double precision, parameter :: SIGMA_FAC = 2.0
- !double precision, parameter :: CONV_STOP = 1.0e-4
! order of interpolating polynomial in conjugate gradient algorithm
! using POLY_ORDER = 3 required computing the gradient of the test model
@@ -190,114 +199,115 @@
!integer, parameter :: GJI_PAPER = 1
! what to perturb, what to invert
- integer, parameter :: PERT_STRUCT = 1
+ ! (For the inverse tests, we only allow perturbations in beta.)
+ integer, parameter :: PERT_STRUCT_BETA = 1
integer, parameter :: PERT_SOURCE_T = 0
integer, parameter :: PERT_SOURCE_X = 0
- integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_STRUCT_BETA = 1
integer, parameter :: INV_SOURCE_T = 0
integer, parameter :: INV_SOURCE_X = 0
! whether to include the model norm term in the misfit function, which acts like damping
logical, parameter :: INCLUDE_MODEL_NORM = .true.
-
+
! log file showing source loactions
logical, parameter :: ISOURCE_LOG = .true.
-
+
! DO NOT CHANGE THESE
- integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ ! integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_STRUCT = 1 ! beta only (24-Jan-2010)
integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
-! parameterization for structure in the inversion
-! 1 : kappa-mu-rho
-! 2 : alpha-beta-rho
-! 3 : c-beta-rho
+ ! parameterization for structure in the inversion
+ ! 1 : kappa-mu-rho
+ ! 2 : alpha-beta-rho
+ ! 3 : c-beta-rho
integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
-! scalelength of smoothing Gaussian
-! GJI paper : 30,60,90
-! body waves : 10
- double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+ ! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 4.50e10 ! Pa
+ double precision, parameter :: RIGIDITY = 3.185e10 ! Pa
+ ! double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ ! double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
- ! parameters for smoothing
- logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
- logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+ ! compute additional parameters
+ double precision, parameter :: PWAVESPEED = &
+ sqrt( (INCOMPRESSIBILITY + (4.0/3.0)*RIGIDITY)/DENSITY )
+ double precision, parameter :: SWAVESPEED = sqrt( RIGIDITY/DENSITY )
+ double precision, parameter :: BWAVESPEED = sqrt( INCOMPRESSIBILITY/DENSITY )
-! homogeneous background model (S.I. units)
- double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
- double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
- double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+ !---------------------------------------------------------------
+ ! measurement windows
-!---------------------------------------------------------------
-! measurement windows
-
double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
-!---------------------------------------------------------------
-! CHT: do not change these
+ !---------------------------------------------------------------
+ ! CHT: do not change these
-! UTM zone for Southern California region
-! integer, parameter :: UTM_PROJECTION_ZONE = 11
+ ! UTM zone for Southern California region
+ ! integer, parameter :: UTM_PROJECTION_ZONE = 11
-! to suppress UTM projection for SCEC benchmarks
+ ! to suppress UTM projection for SCEC benchmarks
logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
-! flag for projection from latitude/longitude to UTM, and back
+ ! flag for projection from latitude/longitude to UTM, and back
integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
-! flag for projection from latitude/longitude to mesh-UTM, and back
+ ! flag for projection from latitude/longitude to mesh-UTM, and back
integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
-! max number of fake receivers
+ ! max number of fake receivers
integer, parameter :: MAX_SR_FAKE = 1000
-! max number of events, receivers, and phass
+ ! max number of events, receivers, and phass
integer, parameter :: MAX_EVENT = 50
integer, parameter :: MAX_SR = 1400
integer, parameter :: MAX_PHASE = 1
integer, parameter :: MAX_COMP = NCOMP
-!========================
-! GRID AND GLL POINTS
+ !========================
+ ! GRID AND GLL POINTS
integer, parameter :: NELE = MAX(NEX,NEZ)
integer, parameter :: NSPEC = NEX*NEZ
-! number of GLL points (polynomial degree plus one)
+ ! number of GLL points (polynomial degree plus one)
integer, parameter :: NGLLX = 5
integer, parameter :: NGLLZ = 5
integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
-! number of global points
+ ! number of global points
integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
-! number of nodes for 2D and 3D shape functions for hexahedra
-! we use 8-node mesh bricks, which are more stable than 27-node elements
+ ! number of nodes for 2D and 3D shape functions for hexahedra
+ ! we use 8-node mesh bricks, which are more stable than 27-node elements
integer, parameter :: NGNOD = 8, NGNOD2D = 4
-
-! number of iterations to solve the system for xi and eta
+
+ ! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 1
-! very large and very small values
+ ! very large and very small values
double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
-! for the Gauss-Lobatto-Legendre points and weights
+ ! for the Gauss-Lobatto-Legendre points and weights
double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
-!
-! CONSTANTS
-!
+ !
+ ! CONSTANTS
+ !
double precision, parameter :: PI = 3.141592653589793
double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
double precision, parameter :: ONE_THIRD = 1.0 / 3.0
double precision, parameter :: ONEOVERTWO = 0.5
double precision, parameter :: DEG = 180.0/PI
-! double precision, parameter :: EPS = 1.0e-35
+ ! double precision, parameter :: EPS = 1.0e-35
!!$! parameter for FFTW
!!$ integer, parameter :: NOUT = NSTEP/2 + 1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -1,28 +1,28 @@
module wave2d_constants
-! This file is copied into the output directory when wave2d.f90 is run
-! There are two basic options: membrane (surface) waves and body waves.
-! Several lines need to be commented/uncommented to go from one to the other.
+ ! This file is copied into the output directory when wave2d.f90 is run
+ ! There are two basic options: membrane (surface) waves and body waves.
+ ! Several lines need to be commented/uncommented to go from one to the other.
-! index of the reference directory for the simulation output
+ ! index of the reference directory for the simulation output
integer, parameter :: IRUNZ = 100
-!========================
-! GRID, TIME-STEP, AND SOURCE PARAMETERS
+ !========================
+ ! GRID, TIME-STEP, AND SOURCE PARAMETERS
-! NFRAME : number of frames to save: membrane, 10; body, 20
-! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
-! NSTEP : number of timesteps
+ ! NFRAME : number of frames to save: membrane, 10; body, 20
+ ! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+ ! NSTEP : number of timesteps
integer, parameter :: NFRAME = 10 ! membrane
!integer, parameter :: NFRAME = 11 ! body
integer, parameter :: NSAVE = 400 ! 200,400
integer, parameter :: NSTEP = NFRAME*NSAVE
-! time step in seconds
+ ! time step in seconds
!double precision, parameter :: DT = 2.0e-2 ! body waves
double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
-! temporal properties of source (source time function)
+ ! temporal properties of source (source time function)
integer, parameter :: ISRC_TIME = 1 ! type (1)
double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
!double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
@@ -30,49 +30,49 @@
!double precision, parameter :: tshift = 8.0*hdur
logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
-! normalization factor of point source force
+ ! normalization factor of point source force
double precision, parameter :: FNORM = 1.0d10
-! forward wavefield source-time function
-! (For membrane waves, FOR_Y = 1 is all that matters.)
+ ! forward wavefield source-time function
+ ! (For membrane waves, FOR_Y = 1 is all that matters.)
integer, parameter :: FOR_X = 1 ! boolean
integer, parameter :: FOR_Y = 1 ! boolean
integer, parameter :: FOR_Z = 1 ! boolean
-! adjoint wavefield source-time function
-! (For membrane waves, REV_Y = 1 is all that matters.)
-! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ ! adjoint wavefield source-time function
+ ! (For membrane waves, REV_Y = 1 is all that matters.)
+ ! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
integer, parameter :: REV_X = 1 ! boolean
integer, parameter :: REV_Y = 1 ! boolean
integer, parameter :: REV_Z = 1 ! boolean
-! spatial properties of sources
-! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
-! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
-! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ ! spatial properties of sources
+ ! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+ ! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+ ! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
-! spatial properties of receivers
-! IREC_SPACE
-! (1) individual station(s)
-! (2) SoCal (used for GJI paper)
-! (3) regular mesh on land
-! (4) regular mesh
-! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
-! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
-! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ ! spatial properties of receivers
+ ! IREC_SPACE
+ ! (1) individual station(s)
+ ! (2) SoCal (used for GJI paper)
+ ! (3) regular mesh on land
+ ! (4) regular mesh
+ ! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+ ! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+ ! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
integer, parameter :: NMESH_REC = 17
- double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
- double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m
double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
-! lower right corner for membrane surface waves plotting grid
+ ! lower right corner for membrane surface waves plotting grid
double precision, parameter :: LAT_MIN = 32.0
double precision, parameter :: LON_MIN = -120.0
integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
-! mesh specifications: membrane surface waves
+ ! mesh specifications: membrane surface waves
double precision, parameter :: LENGTH = 480.0e3 ! m
double precision, parameter :: HEIGHT = 480.0e3 ! m
double precision, parameter :: AREA = LENGTH*HEIGHT
@@ -85,12 +85,12 @@
!!$ integer, parameter :: NEX = 80 ! 160
!!$ integer, parameter :: NEZ = 32 ! 32
-!========================
-! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+ !========================
+ ! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
-! model perturbations for HOMOGENEOUS model (or perturbation)
-! scaling from beta to alpha
-! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ ! model perturbations for HOMOGENEOUS model (or perturbation)
+ ! scaling from beta to alpha
+ ! value is from Masters et al. (2000), "The relative behavior of shear velocity..."
double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
double precision, parameter :: PBETA = 10.0
double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
@@ -103,61 +103,71 @@
! 0 1 2 3
!-----------------------------------------------------------------------------------------
! ISURFACE=1, IMODEL_SYN : homo checker het
- ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het
!-----------------------------------------------------------------------------------------
- ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
- ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker NA
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert NA
!-----------------------------------------------------------------------------------------
- ! whether you want to smooth the structure model for computing the synthetics
+ ! smooth the structure model or kernels
+ ! For the CG algorithm, ISMOOTH_EVENT_KERNEL --> use smoothed event kernels
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 0
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
- integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
- integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
- ! smoothing scalelength
- double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+ ! scalelength of smoothing Gaussian
+ ! GJI paper : 30,60,90
+ ! body waves : 10
+ double precision, parameter :: SIGMA_SMOOTH_KERNEL = 2.121320d4
+ double precision, parameter :: SIGMA_SMOOTH_MODEL = 1.060660d4
+ double precision, parameter :: GAMMA_SMOOTH_KERNEL = sqrt(8.0)*SIGMA_SMOOTH_KERNEL
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = sqrt(8.0)*SIGMA_SMOOTH_MODEL
-!========================
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
-! boolean parameters
-! IKER: (0) waveform
-! (1) traveltime, cross-correlation, misfit
-! (2) amplitude, cross-correlation, misfit
-! (3) traveltime, multitaper
-! (4) amplitude, multitaper
-! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
-! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ !========================
+
+ ! boolean parameters
+ ! IKER: (0) waveform
+ ! (1) traveltime, cross-correlation, misfit
+ ! (2) amplitude, cross-correlation, misfit
+ ! (3) traveltime, multitaper
+ ! (4) amplitude, multitaper
+ ! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+ ! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
integer, parameter :: IKER = 1
integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
-! KEY: USE ONE LINE OR THE OTHER
+ ! KEY: USE ONE LINE OR THE OTHER
integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
-! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+ ! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
-! parameters controlling what to write to file
-! NOTE: for the tomography simulations, ALL of these can be .false.
+ ! parameters controlling what to write to file
+ ! NOTE: for the tomography simulations, ALL of these can be .false.
logical, parameter :: WRITE_STF_F = .false.
logical, parameter :: WRITE_SEISMO_F = .false.
-! logical, parameter :: WRITE_SPECTRA_F = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+ ! logical, parameter :: WRITE_SPECTRA_F = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
-! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+ ! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
logical, parameter :: WRITE_STF_A = .false.
logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
-! logical, parameter :: WRITE_SPECTRA_A = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+ ! logical, parameter :: WRITE_SPECTRA_A = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
-
-!--------------------------------------
-! INVERSION PARAMETERS
+ !--------------------------------------
+ ! INVERSION PARAMETERS
+
! whether you want to compute kernels, or simply the misfit function
logical, parameter :: COMPUTE_KERNELS = .true.
@@ -171,6 +181,8 @@
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
integer, parameter :: NITERATION = 16
double precision, parameter :: VAR_RED_MIN = 8.0
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
! Gaussian errors containted in input file
! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
@@ -178,9 +190,6 @@
double precision, parameter :: SIGMA_DLNA = 1.0
double precision, parameter :: SIGMA_WAVEFORM = 1.0
logical, parameter :: ADD_DATA_ERRORS = .true.
-
- !double precision, parameter :: SIGMA_FAC = 2.0
- !double precision, parameter :: CONV_STOP = 1.0e-4
! order of interpolating polynomial in conjugate gradient algorithm
! using POLY_ORDER = 3 required computing the gradient of the test model
@@ -190,114 +199,115 @@
!integer, parameter :: GJI_PAPER = 1
! what to perturb, what to invert
- integer, parameter :: PERT_STRUCT = 1
+ ! (For the inverse tests, we only allow perturbations in beta.)
+ integer, parameter :: PERT_STRUCT_BETA = 1
integer, parameter :: PERT_SOURCE_T = 0
integer, parameter :: PERT_SOURCE_X = 0
- integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_STRUCT_BETA = 1
integer, parameter :: INV_SOURCE_T = 0
integer, parameter :: INV_SOURCE_X = 0
! whether to include the model norm term in the misfit function, which acts like damping
logical, parameter :: INCLUDE_MODEL_NORM = .true.
-
+
! log file showing source loactions
logical, parameter :: ISOURCE_LOG = .true.
-
+
! DO NOT CHANGE THESE
- integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ ! integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_STRUCT = 1 ! beta only (24-Jan-2010)
integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
-! parameterization for structure in the inversion
-! 1 : kappa-mu-rho
-! 2 : alpha-beta-rho
-! 3 : c-beta-rho
+ ! parameterization for structure in the inversion
+ ! 1 : kappa-mu-rho
+ ! 2 : alpha-beta-rho
+ ! 3 : c-beta-rho
integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
-! scalelength of smoothing Gaussian
-! GJI paper : 30,60,90
-! body waves : 10
- double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+ ! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 4.50e10 ! Pa
+ double precision, parameter :: RIGIDITY = 3.185e10 ! Pa
+ ! double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ ! double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
- ! parameters for smoothing
- logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
- logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+ ! compute additional parameters
+ double precision, parameter :: PWAVESPEED = &
+ sqrt( (INCOMPRESSIBILITY + (4.0/3.0)*RIGIDITY)/DENSITY )
+ double precision, parameter :: SWAVESPEED = sqrt( RIGIDITY/DENSITY )
+ double precision, parameter :: BWAVESPEED = sqrt( INCOMPRESSIBILITY/DENSITY )
-! homogeneous background model (S.I. units)
- double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
- double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
- double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+ !---------------------------------------------------------------
+ ! measurement windows
-!---------------------------------------------------------------
-! measurement windows
-
double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
-!---------------------------------------------------------------
-! CHT: do not change these
+ !---------------------------------------------------------------
+ ! CHT: do not change these
-! UTM zone for Southern California region
-! integer, parameter :: UTM_PROJECTION_ZONE = 11
+ ! UTM zone for Southern California region
+ ! integer, parameter :: UTM_PROJECTION_ZONE = 11
-! to suppress UTM projection for SCEC benchmarks
+ ! to suppress UTM projection for SCEC benchmarks
logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
-! flag for projection from latitude/longitude to UTM, and back
+ ! flag for projection from latitude/longitude to UTM, and back
integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
-! flag for projection from latitude/longitude to mesh-UTM, and back
+ ! flag for projection from latitude/longitude to mesh-UTM, and back
integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
-! max number of fake receivers
+ ! max number of fake receivers
integer, parameter :: MAX_SR_FAKE = 1000
-! max number of events, receivers, and phass
+ ! max number of events, receivers, and phass
integer, parameter :: MAX_EVENT = 50
integer, parameter :: MAX_SR = 1400
integer, parameter :: MAX_PHASE = 1
integer, parameter :: MAX_COMP = NCOMP
-!========================
-! GRID AND GLL POINTS
+ !========================
+ ! GRID AND GLL POINTS
integer, parameter :: NELE = MAX(NEX,NEZ)
integer, parameter :: NSPEC = NEX*NEZ
-! number of GLL points (polynomial degree plus one)
+ ! number of GLL points (polynomial degree plus one)
integer, parameter :: NGLLX = 5
integer, parameter :: NGLLZ = 5
integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
-! number of global points
+ ! number of global points
integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
-! number of nodes for 2D and 3D shape functions for hexahedra
-! we use 8-node mesh bricks, which are more stable than 27-node elements
+ ! number of nodes for 2D and 3D shape functions for hexahedra
+ ! we use 8-node mesh bricks, which are more stable than 27-node elements
integer, parameter :: NGNOD = 8, NGNOD2D = 4
-
-! number of iterations to solve the system for xi and eta
+
+ ! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 1
-! very large and very small values
+ ! very large and very small values
double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
-! for the Gauss-Lobatto-Legendre points and weights
+ ! for the Gauss-Lobatto-Legendre points and weights
double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
-!
-! CONSTANTS
-!
+ !
+ ! CONSTANTS
+ !
double precision, parameter :: PI = 3.141592653589793
double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
double precision, parameter :: ONE_THIRD = 1.0 / 3.0
double precision, parameter :: ONEOVERTWO = 0.5
double precision, parameter :: DEG = 180.0/PI
-! double precision, parameter :: EPS = 1.0e-35
+ ! double precision, parameter :: EPS = 1.0e-35
!!$! parameter for FFTW
!!$ integer, parameter :: NOUT = NSTEP/2 + 1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex02.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex02.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex02.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -1,28 +1,28 @@
module wave2d_constants
-! This file is copied into the output directory when wave2d.f90 is run
-! There are two basic options: membrane (surface) waves and body waves.
-! Several lines need to be commented/uncommented to go from one to the other.
+ ! This file is copied into the output directory when wave2d.f90 is run
+ ! There are two basic options: membrane (surface) waves and body waves.
+ ! Several lines need to be commented/uncommented to go from one to the other.
-! index of the reference directory for the simulation output
+ ! index of the reference directory for the simulation output
integer, parameter :: IRUNZ = 200
-!========================
-! GRID, TIME-STEP, AND SOURCE PARAMETERS
+ !========================
+ ! GRID, TIME-STEP, AND SOURCE PARAMETERS
-! NFRAME : number of frames to save: membrane, 10; body, 20
-! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
-! NSTEP : number of timesteps
+ ! NFRAME : number of frames to save: membrane, 10; body, 20
+ ! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+ ! NSTEP : number of timesteps
integer, parameter :: NFRAME = 10 ! membrane
!integer, parameter :: NFRAME = 11 ! body
integer, parameter :: NSAVE = 400 ! 200,400
integer, parameter :: NSTEP = NFRAME*NSAVE
-! time step in seconds
+ ! time step in seconds
!double precision, parameter :: DT = 2.0e-2 ! body waves
double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
-! temporal properties of source (source time function)
+ ! temporal properties of source (source time function)
integer, parameter :: ISRC_TIME = 1 ! type (1)
double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
!double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
@@ -30,49 +30,49 @@
!double precision, parameter :: tshift = 8.0*hdur
logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
-! normalization factor of point source force
+ ! normalization factor of point source force
double precision, parameter :: FNORM = 1.0d10
-! forward wavefield source-time function
-! (For membrane waves, FOR_Y = 1 is all that matters.)
+ ! forward wavefield source-time function
+ ! (For membrane waves, FOR_Y = 1 is all that matters.)
integer, parameter :: FOR_X = 1 ! boolean
integer, parameter :: FOR_Y = 1 ! boolean
integer, parameter :: FOR_Z = 1 ! boolean
-! adjoint wavefield source-time function
-! (For membrane waves, REV_Y = 1 is all that matters.)
-! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ ! adjoint wavefield source-time function
+ ! (For membrane waves, REV_Y = 1 is all that matters.)
+ ! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
integer, parameter :: REV_X = 1 ! boolean
integer, parameter :: REV_Y = 1 ! boolean
integer, parameter :: REV_Z = 1 ! boolean
-! spatial properties of sources
-! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
-! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
-! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ ! spatial properties of sources
+ ! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+ ! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+ ! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
-! spatial properties of receivers
-! IREC_SPACE
-! (1) individual station(s)
-! (2) SoCal (used for GJI paper)
-! (3) regular mesh on land
-! (4) regular mesh
-! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
-! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
-! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ ! spatial properties of receivers
+ ! IREC_SPACE
+ ! (1) individual station(s)
+ ! (2) SoCal (used for GJI paper)
+ ! (3) regular mesh on land
+ ! (4) regular mesh
+ ! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+ ! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+ ! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
integer, parameter :: NMESH_REC = 17
- double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
- double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m
double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
-! lower right corner for membrane surface waves plotting grid
+ ! lower right corner for membrane surface waves plotting grid
double precision, parameter :: LAT_MIN = 32.0
double precision, parameter :: LON_MIN = -120.0
integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
-! mesh specifications: membrane surface waves
+ ! mesh specifications: membrane surface waves
double precision, parameter :: LENGTH = 480.0e3 ! m
double precision, parameter :: HEIGHT = 480.0e3 ! m
double precision, parameter :: AREA = LENGTH*HEIGHT
@@ -85,12 +85,12 @@
!!$ integer, parameter :: NEX = 80 ! 160
!!$ integer, parameter :: NEZ = 32 ! 32
-!========================
-! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+ !========================
+ ! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
-! model perturbations for HOMOGENEOUS model (or perturbation)
-! scaling from beta to alpha
-! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ ! model perturbations for HOMOGENEOUS model (or perturbation)
+ ! scaling from beta to alpha
+ ! value is from Masters et al. (2000), "The relative behavior of shear velocity..."
double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
double precision, parameter :: PBETA = 10.0
double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
@@ -103,61 +103,71 @@
! 0 1 2 3
!-----------------------------------------------------------------------------------------
! ISURFACE=1, IMODEL_SYN : homo checker het
- ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het
!-----------------------------------------------------------------------------------------
- ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
- ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker NA
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert NA
!-----------------------------------------------------------------------------------------
- ! whether you want to smooth the structure model for computing the synthetics
+ ! smooth the structure model or kernels
+ ! For the CG algorithm, ISMOOTH_EVENT_KERNEL --> use smoothed event kernels
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 0
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
- integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
- integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
- ! smoothing scalelength
- double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+ ! scalelength of smoothing Gaussian
+ ! GJI paper : 30,60,90
+ ! body waves : 10
+ double precision, parameter :: SIGMA_SMOOTH_KERNEL = 2.121320d4
+ double precision, parameter :: SIGMA_SMOOTH_MODEL = 1.060660d4
+ double precision, parameter :: GAMMA_SMOOTH_KERNEL = sqrt(8.0)*SIGMA_SMOOTH_KERNEL
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = sqrt(8.0)*SIGMA_SMOOTH_MODEL
-!========================
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
-! boolean parameters
-! IKER: (0) waveform
-! (1) traveltime, cross-correlation, misfit
-! (2) amplitude, cross-correlation, misfit
-! (3) traveltime, multitaper
-! (4) amplitude, multitaper
-! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
-! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ !========================
+
+ ! boolean parameters
+ ! IKER: (0) waveform
+ ! (1) traveltime, cross-correlation, misfit
+ ! (2) amplitude, cross-correlation, misfit
+ ! (3) traveltime, multitaper
+ ! (4) amplitude, multitaper
+ ! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+ ! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
integer, parameter :: IKER = 1
integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
-! KEY: USE ONE LINE OR THE OTHER
+ ! KEY: USE ONE LINE OR THE OTHER
integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
-! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+ ! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
-! parameters controlling what to write to file
-! NOTE: for the tomography simulations, ALL of these can be .false.
+ ! parameters controlling what to write to file
+ ! NOTE: for the tomography simulations, ALL of these can be .false.
logical, parameter :: WRITE_STF_F = .false.
logical, parameter :: WRITE_SEISMO_F = .false.
-! logical, parameter :: WRITE_SPECTRA_F = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+ ! logical, parameter :: WRITE_SPECTRA_F = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
-! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+ ! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
logical, parameter :: WRITE_STF_A = .false.
logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
-! logical, parameter :: WRITE_SPECTRA_A = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+ ! logical, parameter :: WRITE_SPECTRA_A = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
-
-!--------------------------------------
-! INVERSION PARAMETERS
+ !--------------------------------------
+ ! INVERSION PARAMETERS
+
! whether you want to compute kernels, or simply the misfit function
logical, parameter :: COMPUTE_KERNELS = .true.
@@ -171,6 +181,8 @@
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
integer, parameter :: NITERATION = 16
double precision, parameter :: VAR_RED_MIN = 8.0
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
! Gaussian errors containted in input file
! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
@@ -178,9 +190,6 @@
double precision, parameter :: SIGMA_DLNA = 1.0
double precision, parameter :: SIGMA_WAVEFORM = 1.0
logical, parameter :: ADD_DATA_ERRORS = .true.
-
- !double precision, parameter :: SIGMA_FAC = 2.0
- !double precision, parameter :: CONV_STOP = 1.0e-4
! order of interpolating polynomial in conjugate gradient algorithm
! using POLY_ORDER = 3 required computing the gradient of the test model
@@ -190,114 +199,115 @@
!integer, parameter :: GJI_PAPER = 1
! what to perturb, what to invert
- integer, parameter :: PERT_STRUCT = 0
+ ! (For the inverse tests, we only allow perturbations in beta.)
+ integer, parameter :: PERT_STRUCT_BETA = 0
integer, parameter :: PERT_SOURCE_T = 1
integer, parameter :: PERT_SOURCE_X = 1
- integer, parameter :: INV_STRUCT = 0
+ integer, parameter :: INV_STRUCT_BETA = 0
integer, parameter :: INV_SOURCE_T = 1
integer, parameter :: INV_SOURCE_X = 1
! whether to include the model norm term in the misfit function, which acts like damping
logical, parameter :: INCLUDE_MODEL_NORM = .true.
-
+
! log file showing source loactions
logical, parameter :: ISOURCE_LOG = .true.
-
+
! DO NOT CHANGE THESE
- integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ ! integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_STRUCT = 1 ! beta only (24-Jan-2010)
integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
-! parameterization for structure in the inversion
-! 1 : kappa-mu-rho
-! 2 : alpha-beta-rho
-! 3 : c-beta-rho
+ ! parameterization for structure in the inversion
+ ! 1 : kappa-mu-rho
+ ! 2 : alpha-beta-rho
+ ! 3 : c-beta-rho
integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
-! scalelength of smoothing Gaussian
-! GJI paper : 30,60,90
-! body waves : 10
- double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+ ! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 4.50e10 ! Pa
+ double precision, parameter :: RIGIDITY = 3.185e10 ! Pa
+ ! double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ ! double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
- ! parameters for smoothing
- logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
- logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+ ! compute additional parameters
+ double precision, parameter :: PWAVESPEED = &
+ sqrt( (INCOMPRESSIBILITY + (4.0/3.0)*RIGIDITY)/DENSITY )
+ double precision, parameter :: SWAVESPEED = sqrt( RIGIDITY/DENSITY )
+ double precision, parameter :: BWAVESPEED = sqrt( INCOMPRESSIBILITY/DENSITY )
-! homogeneous background model (S.I. units)
- double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
- double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
- double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+ !---------------------------------------------------------------
+ ! measurement windows
-!---------------------------------------------------------------
-! measurement windows
-
double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
-!---------------------------------------------------------------
-! CHT: do not change these
+ !---------------------------------------------------------------
+ ! CHT: do not change these
-! UTM zone for Southern California region
-! integer, parameter :: UTM_PROJECTION_ZONE = 11
+ ! UTM zone for Southern California region
+ ! integer, parameter :: UTM_PROJECTION_ZONE = 11
-! to suppress UTM projection for SCEC benchmarks
+ ! to suppress UTM projection for SCEC benchmarks
logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
-! flag for projection from latitude/longitude to UTM, and back
+ ! flag for projection from latitude/longitude to UTM, and back
integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
-! flag for projection from latitude/longitude to mesh-UTM, and back
+ ! flag for projection from latitude/longitude to mesh-UTM, and back
integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
-! max number of fake receivers
+ ! max number of fake receivers
integer, parameter :: MAX_SR_FAKE = 1000
-! max number of events, receivers, and phass
+ ! max number of events, receivers, and phass
integer, parameter :: MAX_EVENT = 50
integer, parameter :: MAX_SR = 1400
integer, parameter :: MAX_PHASE = 1
integer, parameter :: MAX_COMP = NCOMP
-!========================
-! GRID AND GLL POINTS
+ !========================
+ ! GRID AND GLL POINTS
integer, parameter :: NELE = MAX(NEX,NEZ)
integer, parameter :: NSPEC = NEX*NEZ
-! number of GLL points (polynomial degree plus one)
+ ! number of GLL points (polynomial degree plus one)
integer, parameter :: NGLLX = 5
integer, parameter :: NGLLZ = 5
integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
-! number of global points
+ ! number of global points
integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
-! number of nodes for 2D and 3D shape functions for hexahedra
-! we use 8-node mesh bricks, which are more stable than 27-node elements
+ ! number of nodes for 2D and 3D shape functions for hexahedra
+ ! we use 8-node mesh bricks, which are more stable than 27-node elements
integer, parameter :: NGNOD = 8, NGNOD2D = 4
-
-! number of iterations to solve the system for xi and eta
+
+ ! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 1
-! very large and very small values
+ ! very large and very small values
double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
-! for the Gauss-Lobatto-Legendre points and weights
+ ! for the Gauss-Lobatto-Legendre points and weights
double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
-!
-! CONSTANTS
-!
+ !
+ ! CONSTANTS
+ !
double precision, parameter :: PI = 3.141592653589793
double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
double precision, parameter :: ONE_THIRD = 1.0 / 3.0
double precision, parameter :: ONEOVERTWO = 0.5
double precision, parameter :: DEG = 180.0/PI
-! double precision, parameter :: EPS = 1.0e-35
+ ! double precision, parameter :: EPS = 1.0e-35
!!$! parameter for FFTW
!!$ integer, parameter :: NOUT = NSTEP/2 + 1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex03.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex03.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex03.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -1,28 +1,28 @@
module wave2d_constants
-! This file is copied into the output directory when wave2d.f90 is run
-! There are two basic options: membrane (surface) waves and body waves.
-! Several lines need to be commented/uncommented to go from one to the other.
+ ! This file is copied into the output directory when wave2d.f90 is run
+ ! There are two basic options: membrane (surface) waves and body waves.
+ ! Several lines need to be commented/uncommented to go from one to the other.
-! index of the reference directory for the simulation output
+ ! index of the reference directory for the simulation output
integer, parameter :: IRUNZ = 300
-!========================
-! GRID, TIME-STEP, AND SOURCE PARAMETERS
+ !========================
+ ! GRID, TIME-STEP, AND SOURCE PARAMETERS
-! NFRAME : number of frames to save: membrane, 10; body, 20
-! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
-! NSTEP : number of timesteps
+ ! NFRAME : number of frames to save: membrane, 10; body, 20
+ ! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+ ! NSTEP : number of timesteps
integer, parameter :: NFRAME = 10 ! membrane
!integer, parameter :: NFRAME = 11 ! body
integer, parameter :: NSAVE = 400 ! 200,400
integer, parameter :: NSTEP = NFRAME*NSAVE
-! time step in seconds
+ ! time step in seconds
!double precision, parameter :: DT = 2.0e-2 ! body waves
double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
-! temporal properties of source (source time function)
+ ! temporal properties of source (source time function)
integer, parameter :: ISRC_TIME = 1 ! type (1)
double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
!double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
@@ -30,49 +30,49 @@
!double precision, parameter :: tshift = 8.0*hdur
logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
-! normalization factor of point source force
+ ! normalization factor of point source force
double precision, parameter :: FNORM = 1.0d10
-! forward wavefield source-time function
-! (For membrane waves, FOR_Y = 1 is all that matters.)
+ ! forward wavefield source-time function
+ ! (For membrane waves, FOR_Y = 1 is all that matters.)
integer, parameter :: FOR_X = 1 ! boolean
integer, parameter :: FOR_Y = 1 ! boolean
integer, parameter :: FOR_Z = 1 ! boolean
-! adjoint wavefield source-time function
-! (For membrane waves, REV_Y = 1 is all that matters.)
-! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ ! adjoint wavefield source-time function
+ ! (For membrane waves, REV_Y = 1 is all that matters.)
+ ! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
integer, parameter :: REV_X = 1 ! boolean
integer, parameter :: REV_Y = 1 ! boolean
integer, parameter :: REV_Z = 1 ! boolean
-! spatial properties of sources
-! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
-! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
-! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ ! spatial properties of sources
+ ! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+ ! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+ ! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
-! spatial properties of receivers
-! IREC_SPACE
-! (1) individual station(s)
-! (2) SoCal (used for GJI paper)
-! (3) regular mesh on land
-! (4) regular mesh
-! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
-! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
-! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ ! spatial properties of receivers
+ ! IREC_SPACE
+ ! (1) individual station(s)
+ ! (2) SoCal (used for GJI paper)
+ ! (3) regular mesh on land
+ ! (4) regular mesh
+ ! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+ ! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+ ! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
integer, parameter :: NMESH_REC = 17
- double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
- double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m
double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
-! lower right corner for membrane surface waves plotting grid
+ ! lower right corner for membrane surface waves plotting grid
double precision, parameter :: LAT_MIN = 32.0
double precision, parameter :: LON_MIN = -120.0
integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
-! mesh specifications: membrane surface waves
+ ! mesh specifications: membrane surface waves
double precision, parameter :: LENGTH = 480.0e3 ! m
double precision, parameter :: HEIGHT = 480.0e3 ! m
double precision, parameter :: AREA = LENGTH*HEIGHT
@@ -85,12 +85,12 @@
!!$ integer, parameter :: NEX = 80 ! 160
!!$ integer, parameter :: NEZ = 32 ! 32
-!========================
-! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+ !========================
+ ! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
-! model perturbations for HOMOGENEOUS model (or perturbation)
-! scaling from beta to alpha
-! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ ! model perturbations for HOMOGENEOUS model (or perturbation)
+ ! scaling from beta to alpha
+ ! value is from Masters et al. (2000), "The relative behavior of shear velocity..."
double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
double precision, parameter :: PBETA = 10.0
double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
@@ -103,61 +103,71 @@
! 0 1 2 3
!-----------------------------------------------------------------------------------------
! ISURFACE=1, IMODEL_SYN : homo checker het
- ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het
!-----------------------------------------------------------------------------------------
- ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
- ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker NA
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert NA
!-----------------------------------------------------------------------------------------
- ! whether you want to smooth the structure model for computing the synthetics
+ ! smooth the structure model or kernels
+ ! For the CG algorithm, ISMOOTH_EVENT_KERNEL --> use smoothed event kernels
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 0
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
- integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
- integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
- ! smoothing scalelength
- double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+ ! scalelength of smoothing Gaussian
+ ! GJI paper : 30,60,90
+ ! body waves : 10
+ double precision, parameter :: SIGMA_SMOOTH_KERNEL = 2.121320d4
+ double precision, parameter :: SIGMA_SMOOTH_MODEL = 1.060660d4
+ double precision, parameter :: GAMMA_SMOOTH_KERNEL = sqrt(8.0)*SIGMA_SMOOTH_KERNEL
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = sqrt(8.0)*SIGMA_SMOOTH_MODEL
-!========================
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
-! boolean parameters
-! IKER: (0) waveform
-! (1) traveltime, cross-correlation, misfit
-! (2) amplitude, cross-correlation, misfit
-! (3) traveltime, multitaper
-! (4) amplitude, multitaper
-! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
-! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ !========================
+
+ ! boolean parameters
+ ! IKER: (0) waveform
+ ! (1) traveltime, cross-correlation, misfit
+ ! (2) amplitude, cross-correlation, misfit
+ ! (3) traveltime, multitaper
+ ! (4) amplitude, multitaper
+ ! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+ ! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
integer, parameter :: IKER = 1
integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
-! KEY: USE ONE LINE OR THE OTHER
+ ! KEY: USE ONE LINE OR THE OTHER
integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
-! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+ ! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
-! parameters controlling what to write to file
-! NOTE: for the tomography simulations, ALL of these can be .false.
+ ! parameters controlling what to write to file
+ ! NOTE: for the tomography simulations, ALL of these can be .false.
logical, parameter :: WRITE_STF_F = .false.
logical, parameter :: WRITE_SEISMO_F = .false.
-! logical, parameter :: WRITE_SPECTRA_F = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+ ! logical, parameter :: WRITE_SPECTRA_F = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
-! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+ ! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
logical, parameter :: WRITE_STF_A = .false.
logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
-! logical, parameter :: WRITE_SPECTRA_A = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+ ! logical, parameter :: WRITE_SPECTRA_A = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
-
-!--------------------------------------
-! INVERSION PARAMETERS
+ !--------------------------------------
+ ! INVERSION PARAMETERS
+
! whether you want to compute kernels, or simply the misfit function
logical, parameter :: COMPUTE_KERNELS = .true.
@@ -171,6 +181,8 @@
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
integer, parameter :: NITERATION = 16
double precision, parameter :: VAR_RED_MIN = 8.0
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
! Gaussian errors containted in input file
! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
@@ -178,9 +190,6 @@
double precision, parameter :: SIGMA_DLNA = 1.0
double precision, parameter :: SIGMA_WAVEFORM = 1.0
logical, parameter :: ADD_DATA_ERRORS = .true.
-
- !double precision, parameter :: SIGMA_FAC = 2.0
- !double precision, parameter :: CONV_STOP = 1.0e-4
! order of interpolating polynomial in conjugate gradient algorithm
! using POLY_ORDER = 3 required computing the gradient of the test model
@@ -190,114 +199,115 @@
!integer, parameter :: GJI_PAPER = 1
! what to perturb, what to invert
- integer, parameter :: PERT_STRUCT = 1
+ ! (For the inverse tests, we only allow perturbations in beta.)
+ integer, parameter :: PERT_STRUCT_BETA = 1
integer, parameter :: PERT_SOURCE_T = 1
integer, parameter :: PERT_SOURCE_X = 1
- integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_STRUCT_BETA = 1
integer, parameter :: INV_SOURCE_T = 1
integer, parameter :: INV_SOURCE_X = 1
! whether to include the model norm term in the misfit function, which acts like damping
logical, parameter :: INCLUDE_MODEL_NORM = .true.
-
+
! log file showing source loactions
logical, parameter :: ISOURCE_LOG = .true.
-
+
! DO NOT CHANGE THESE
- integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ ! integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_STRUCT = 1 ! beta only (24-Jan-2010)
integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
-! parameterization for structure in the inversion
-! 1 : kappa-mu-rho
-! 2 : alpha-beta-rho
-! 3 : c-beta-rho
+ ! parameterization for structure in the inversion
+ ! 1 : kappa-mu-rho
+ ! 2 : alpha-beta-rho
+ ! 3 : c-beta-rho
integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
-! scalelength of smoothing Gaussian
-! GJI paper : 30,60,90
-! body waves : 10
- double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+ ! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 4.50e10 ! Pa
+ double precision, parameter :: RIGIDITY = 3.185e10 ! Pa
+ ! double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ ! double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
- ! parameters for smoothing
- logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
- logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+ ! compute additional parameters
+ double precision, parameter :: PWAVESPEED = &
+ sqrt( (INCOMPRESSIBILITY + (4.0/3.0)*RIGIDITY)/DENSITY )
+ double precision, parameter :: SWAVESPEED = sqrt( RIGIDITY/DENSITY )
+ double precision, parameter :: BWAVESPEED = sqrt( INCOMPRESSIBILITY/DENSITY )
-! homogeneous background model (S.I. units)
- double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
- double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
- double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+ !---------------------------------------------------------------
+ ! measurement windows
-!---------------------------------------------------------------
-! measurement windows
-
double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
-!---------------------------------------------------------------
-! CHT: do not change these
+ !---------------------------------------------------------------
+ ! CHT: do not change these
-! UTM zone for Southern California region
-! integer, parameter :: UTM_PROJECTION_ZONE = 11
+ ! UTM zone for Southern California region
+ ! integer, parameter :: UTM_PROJECTION_ZONE = 11
-! to suppress UTM projection for SCEC benchmarks
+ ! to suppress UTM projection for SCEC benchmarks
logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
-! flag for projection from latitude/longitude to UTM, and back
+ ! flag for projection from latitude/longitude to UTM, and back
integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
-! flag for projection from latitude/longitude to mesh-UTM, and back
+ ! flag for projection from latitude/longitude to mesh-UTM, and back
integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
-! max number of fake receivers
+ ! max number of fake receivers
integer, parameter :: MAX_SR_FAKE = 1000
-! max number of events, receivers, and phass
+ ! max number of events, receivers, and phass
integer, parameter :: MAX_EVENT = 50
integer, parameter :: MAX_SR = 1400
integer, parameter :: MAX_PHASE = 1
integer, parameter :: MAX_COMP = NCOMP
-!========================
-! GRID AND GLL POINTS
+ !========================
+ ! GRID AND GLL POINTS
integer, parameter :: NELE = MAX(NEX,NEZ)
integer, parameter :: NSPEC = NEX*NEZ
-! number of GLL points (polynomial degree plus one)
+ ! number of GLL points (polynomial degree plus one)
integer, parameter :: NGLLX = 5
integer, parameter :: NGLLZ = 5
integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
-! number of global points
+ ! number of global points
integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
-! number of nodes for 2D and 3D shape functions for hexahedra
-! we use 8-node mesh bricks, which are more stable than 27-node elements
+ ! number of nodes for 2D and 3D shape functions for hexahedra
+ ! we use 8-node mesh bricks, which are more stable than 27-node elements
integer, parameter :: NGNOD = 8, NGNOD2D = 4
-
-! number of iterations to solve the system for xi and eta
+
+ ! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 1
-! very large and very small values
+ ! very large and very small values
double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
-! for the Gauss-Lobatto-Legendre points and weights
+ ! for the Gauss-Lobatto-Legendre points and weights
double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
-!
-! CONSTANTS
-!
+ !
+ ! CONSTANTS
+ !
double precision, parameter :: PI = 3.141592653589793
double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
double precision, parameter :: ONE_THIRD = 1.0 / 3.0
double precision, parameter :: ONEOVERTWO = 0.5
double precision, parameter :: DEG = 180.0/PI
-! double precision, parameter :: EPS = 1.0e-35
+ ! double precision, parameter :: EPS = 1.0e-35
!!$! parameter for FFTW
!!$ integer, parameter :: NOUT = NSTEP/2 + 1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex04.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex04.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex04.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -1,28 +1,28 @@
module wave2d_constants
-! This file is copied into the output directory when wave2d.f90 is run
-! There are two basic options: membrane (surface) waves and body waves.
-! Several lines need to be commented/uncommented to go from one to the other.
+ ! This file is copied into the output directory when wave2d.f90 is run
+ ! There are two basic options: membrane (surface) waves and body waves.
+ ! Several lines need to be commented/uncommented to go from one to the other.
-! index of the reference directory for the simulation output
+ ! index of the reference directory for the simulation output
integer, parameter :: IRUNZ = 400
-!========================
-! GRID, TIME-STEP, AND SOURCE PARAMETERS
+ !========================
+ ! GRID, TIME-STEP, AND SOURCE PARAMETERS
-! NFRAME : number of frames to save: membrane, 10; body, 20
-! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
-! NSTEP : number of timesteps
+ ! NFRAME : number of frames to save: membrane, 10; body, 20
+ ! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+ ! NSTEP : number of timesteps
integer, parameter :: NFRAME = 10 ! membrane
!integer, parameter :: NFRAME = 11 ! body
integer, parameter :: NSAVE = 400 ! 200,400
integer, parameter :: NSTEP = NFRAME*NSAVE
-! time step in seconds
+ ! time step in seconds
!double precision, parameter :: DT = 2.0e-2 ! body waves
double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
-! temporal properties of source (source time function)
+ ! temporal properties of source (source time function)
integer, parameter :: ISRC_TIME = 1 ! type (1)
double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
!double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
@@ -30,49 +30,49 @@
!double precision, parameter :: tshift = 8.0*hdur
logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
-! normalization factor of point source force
+ ! normalization factor of point source force
double precision, parameter :: FNORM = 1.0d10
-! forward wavefield source-time function
-! (For membrane waves, FOR_Y = 1 is all that matters.)
+ ! forward wavefield source-time function
+ ! (For membrane waves, FOR_Y = 1 is all that matters.)
integer, parameter :: FOR_X = 1 ! boolean
integer, parameter :: FOR_Y = 1 ! boolean
integer, parameter :: FOR_Z = 1 ! boolean
-! adjoint wavefield source-time function
-! (For membrane waves, REV_Y = 1 is all that matters.)
-! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ ! adjoint wavefield source-time function
+ ! (For membrane waves, REV_Y = 1 is all that matters.)
+ ! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
integer, parameter :: REV_X = 1 ! boolean
integer, parameter :: REV_Y = 1 ! boolean
integer, parameter :: REV_Z = 1 ! boolean
-! spatial properties of sources
-! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
-! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
-! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ ! spatial properties of sources
+ ! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+ ! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+ ! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
-! spatial properties of receivers
-! IREC_SPACE
-! (1) individual station(s)
-! (2) SoCal (used for GJI paper)
-! (3) regular mesh on land
-! (4) regular mesh
-! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
-! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
-! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ ! spatial properties of receivers
+ ! IREC_SPACE
+ ! (1) individual station(s)
+ ! (2) SoCal (used for GJI paper)
+ ! (3) regular mesh on land
+ ! (4) regular mesh
+ ! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+ ! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+ ! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
integer, parameter :: NMESH_REC = 17
- double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
- double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m
double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
-! lower right corner for membrane surface waves plotting grid
+ ! lower right corner for membrane surface waves plotting grid
double precision, parameter :: LAT_MIN = 32.0
double precision, parameter :: LON_MIN = -120.0
integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
-! mesh specifications: membrane surface waves
+ ! mesh specifications: membrane surface waves
double precision, parameter :: LENGTH = 480.0e3 ! m
double precision, parameter :: HEIGHT = 480.0e3 ! m
double precision, parameter :: AREA = LENGTH*HEIGHT
@@ -85,12 +85,12 @@
!!$ integer, parameter :: NEX = 80 ! 160
!!$ integer, parameter :: NEZ = 32 ! 32
-!========================
-! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+ !========================
+ ! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
-! model perturbations for HOMOGENEOUS model (or perturbation)
-! scaling from beta to alpha
-! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ ! model perturbations for HOMOGENEOUS model (or perturbation)
+ ! scaling from beta to alpha
+ ! value is from Masters et al. (2000), "The relative behavior of shear velocity..."
double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
double precision, parameter :: PBETA = 10.0
double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
@@ -103,61 +103,71 @@
! 0 1 2 3
!-----------------------------------------------------------------------------------------
! ISURFACE=1, IMODEL_SYN : homo checker het
- ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het
!-----------------------------------------------------------------------------------------
- ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
- ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker NA
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert NA
!-----------------------------------------------------------------------------------------
- ! whether you want to smooth the structure model for computing the synthetics
+ ! smooth the structure model or kernels
+ ! For the CG algorithm, ISMOOTH_EVENT_KERNEL --> use smoothed event kernels
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 0
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
- integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
- integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
- ! smoothing scalelength
- double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+ ! scalelength of smoothing Gaussian
+ ! GJI paper : 30,60,90
+ ! body waves : 10
+ double precision, parameter :: SIGMA_SMOOTH_KERNEL = 2.121320d4
+ double precision, parameter :: SIGMA_SMOOTH_MODEL = 1.060660d4
+ double precision, parameter :: GAMMA_SMOOTH_KERNEL = sqrt(8.0)*SIGMA_SMOOTH_KERNEL
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = sqrt(8.0)*SIGMA_SMOOTH_MODEL
-!========================
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
-! boolean parameters
-! IKER: (0) waveform
-! (1) traveltime, cross-correlation, misfit
-! (2) amplitude, cross-correlation, misfit
-! (3) traveltime, multitaper
-! (4) amplitude, multitaper
-! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
-! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ !========================
+
+ ! boolean parameters
+ ! IKER: (0) waveform
+ ! (1) traveltime, cross-correlation, misfit
+ ! (2) amplitude, cross-correlation, misfit
+ ! (3) traveltime, multitaper
+ ! (4) amplitude, multitaper
+ ! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+ ! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
integer, parameter :: IKER = 1
integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
-! KEY: USE ONE LINE OR THE OTHER
+ ! KEY: USE ONE LINE OR THE OTHER
integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
-! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+ ! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
-! parameters controlling what to write to file
-! NOTE: for the tomography simulations, ALL of these can be .false.
+ ! parameters controlling what to write to file
+ ! NOTE: for the tomography simulations, ALL of these can be .false.
logical, parameter :: WRITE_STF_F = .false.
logical, parameter :: WRITE_SEISMO_F = .false.
-! logical, parameter :: WRITE_SPECTRA_F = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+ ! logical, parameter :: WRITE_SPECTRA_F = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
-! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+ ! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
logical, parameter :: WRITE_STF_A = .false.
logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
-! logical, parameter :: WRITE_SPECTRA_A = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+ ! logical, parameter :: WRITE_SPECTRA_A = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
-
-!--------------------------------------
-! INVERSION PARAMETERS
+ !--------------------------------------
+ ! INVERSION PARAMETERS
+
! whether you want to compute kernels, or simply the misfit function
logical, parameter :: COMPUTE_KERNELS = .true.
@@ -171,6 +181,8 @@
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
integer, parameter :: NITERATION = 16
double precision, parameter :: VAR_RED_MIN = 8.0
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
! Gaussian errors containted in input file
! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
@@ -178,9 +190,6 @@
double precision, parameter :: SIGMA_DLNA = 1.0
double precision, parameter :: SIGMA_WAVEFORM = 1.0
logical, parameter :: ADD_DATA_ERRORS = .true.
-
- !double precision, parameter :: SIGMA_FAC = 2.0
- !double precision, parameter :: CONV_STOP = 1.0e-4
! order of interpolating polynomial in conjugate gradient algorithm
! using POLY_ORDER = 3 required computing the gradient of the test model
@@ -190,114 +199,115 @@
!integer, parameter :: GJI_PAPER = 1
! what to perturb, what to invert
- integer, parameter :: PERT_STRUCT = 1
+ ! (For the inverse tests, we only allow perturbations in beta.)
+ integer, parameter :: PERT_STRUCT_BETA = 1
integer, parameter :: PERT_SOURCE_T = 1
integer, parameter :: PERT_SOURCE_X = 1
- integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_STRUCT_BETA = 1
integer, parameter :: INV_SOURCE_T = 1
integer, parameter :: INV_SOURCE_X = 1
! whether to include the model norm term in the misfit function, which acts like damping
logical, parameter :: INCLUDE_MODEL_NORM = .true.
-
+
! log file showing source loactions
logical, parameter :: ISOURCE_LOG = .true.
-
+
! DO NOT CHANGE THESE
- integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ ! integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_STRUCT = 1 ! beta only (24-Jan-2010)
integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
-! parameterization for structure in the inversion
-! 1 : kappa-mu-rho
-! 2 : alpha-beta-rho
-! 3 : c-beta-rho
+ ! parameterization for structure in the inversion
+ ! 1 : kappa-mu-rho
+ ! 2 : alpha-beta-rho
+ ! 3 : c-beta-rho
integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
-! scalelength of smoothing Gaussian
-! GJI paper : 30,60,90
-! body waves : 10
- double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+ ! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 4.50e10 ! Pa
+ double precision, parameter :: RIGIDITY = 3.185e10 ! Pa
+ ! double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ ! double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
- ! parameters for smoothing
- logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
- logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+ ! compute additional parameters
+ double precision, parameter :: PWAVESPEED = &
+ sqrt( (INCOMPRESSIBILITY + (4.0/3.0)*RIGIDITY)/DENSITY )
+ double precision, parameter :: SWAVESPEED = sqrt( RIGIDITY/DENSITY )
+ double precision, parameter :: BWAVESPEED = sqrt( INCOMPRESSIBILITY/DENSITY )
-! homogeneous background model (S.I. units)
- double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
- double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
- double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+ !---------------------------------------------------------------
+ ! measurement windows
-!---------------------------------------------------------------
-! measurement windows
-
double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
-!---------------------------------------------------------------
-! CHT: do not change these
+ !---------------------------------------------------------------
+ ! CHT: do not change these
-! UTM zone for Southern California region
-! integer, parameter :: UTM_PROJECTION_ZONE = 11
+ ! UTM zone for Southern California region
+ ! integer, parameter :: UTM_PROJECTION_ZONE = 11
-! to suppress UTM projection for SCEC benchmarks
+ ! to suppress UTM projection for SCEC benchmarks
logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
-! flag for projection from latitude/longitude to UTM, and back
+ ! flag for projection from latitude/longitude to UTM, and back
integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
-! flag for projection from latitude/longitude to mesh-UTM, and back
+ ! flag for projection from latitude/longitude to mesh-UTM, and back
integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
-! max number of fake receivers
+ ! max number of fake receivers
integer, parameter :: MAX_SR_FAKE = 1000
-! max number of events, receivers, and phass
+ ! max number of events, receivers, and phass
integer, parameter :: MAX_EVENT = 50
integer, parameter :: MAX_SR = 1400
integer, parameter :: MAX_PHASE = 1
integer, parameter :: MAX_COMP = NCOMP
-!========================
-! GRID AND GLL POINTS
+ !========================
+ ! GRID AND GLL POINTS
integer, parameter :: NELE = MAX(NEX,NEZ)
integer, parameter :: NSPEC = NEX*NEZ
-! number of GLL points (polynomial degree plus one)
+ ! number of GLL points (polynomial degree plus one)
integer, parameter :: NGLLX = 5
integer, parameter :: NGLLZ = 5
integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
-! number of global points
+ ! number of global points
integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
-! number of nodes for 2D and 3D shape functions for hexahedra
-! we use 8-node mesh bricks, which are more stable than 27-node elements
+ ! number of nodes for 2D and 3D shape functions for hexahedra
+ ! we use 8-node mesh bricks, which are more stable than 27-node elements
integer, parameter :: NGNOD = 8, NGNOD2D = 4
-
-! number of iterations to solve the system for xi and eta
+
+ ! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 1
-! very large and very small values
+ ! very large and very small values
double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
-! for the Gauss-Lobatto-Legendre points and weights
+ ! for the Gauss-Lobatto-Legendre points and weights
double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
-!
-! CONSTANTS
-!
+ !
+ ! CONSTANTS
+ !
double precision, parameter :: PI = 3.141592653589793
double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
double precision, parameter :: ONE_THIRD = 1.0 / 3.0
double precision, parameter :: ONEOVERTWO = 0.5
double precision, parameter :: DEG = 180.0/PI
-! double precision, parameter :: EPS = 1.0e-35
+ ! double precision, parameter :: EPS = 1.0e-35
!!$! parameter for FFTW
!!$ integer, parameter :: NOUT = NSTEP/2 + 1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex05.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex05.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex05.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -1,28 +1,28 @@
module wave2d_constants
-! This file is copied into the output directory when wave2d.f90 is run
-! There are two basic options: membrane (surface) waves and body waves.
-! Several lines need to be commented/uncommented to go from one to the other.
+ ! This file is copied into the output directory when wave2d.f90 is run
+ ! There are two basic options: membrane (surface) waves and body waves.
+ ! Several lines need to be commented/uncommented to go from one to the other.
-! index of the reference directory for the simulation output
+ ! index of the reference directory for the simulation output
integer, parameter :: IRUNZ = 500
-!========================
-! GRID, TIME-STEP, AND SOURCE PARAMETERS
+ !========================
+ ! GRID, TIME-STEP, AND SOURCE PARAMETERS
-! NFRAME : number of frames to save: membrane, 10; body, 20
-! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
-! NSTEP : number of timesteps
+ ! NFRAME : number of frames to save: membrane, 10; body, 20
+ ! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+ ! NSTEP : number of timesteps
integer, parameter :: NFRAME = 10 ! membrane
!integer, parameter :: NFRAME = 11 ! body
integer, parameter :: NSAVE = 400 ! 200,400
integer, parameter :: NSTEP = NFRAME*NSAVE
-! time step in seconds
+ ! time step in seconds
!double precision, parameter :: DT = 2.0e-2 ! body waves
double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
-! temporal properties of source (source time function)
+ ! temporal properties of source (source time function)
integer, parameter :: ISRC_TIME = 1 ! type (1)
double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
!double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
@@ -30,49 +30,49 @@
!double precision, parameter :: tshift = 8.0*hdur
logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
-! normalization factor of point source force
+ ! normalization factor of point source force
double precision, parameter :: FNORM = 1.0d10
-! forward wavefield source-time function
-! (For membrane waves, FOR_Y = 1 is all that matters.)
+ ! forward wavefield source-time function
+ ! (For membrane waves, FOR_Y = 1 is all that matters.)
integer, parameter :: FOR_X = 1 ! boolean
integer, parameter :: FOR_Y = 1 ! boolean
integer, parameter :: FOR_Z = 1 ! boolean
-! adjoint wavefield source-time function
-! (For membrane waves, REV_Y = 1 is all that matters.)
-! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ ! adjoint wavefield source-time function
+ ! (For membrane waves, REV_Y = 1 is all that matters.)
+ ! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
integer, parameter :: REV_X = 1 ! boolean
integer, parameter :: REV_Y = 1 ! boolean
integer, parameter :: REV_Z = 1 ! boolean
-! spatial properties of sources
-! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
-! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
-! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ ! spatial properties of sources
+ ! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+ ! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+ ! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
-! spatial properties of receivers
-! IREC_SPACE
-! (1) individual station(s)
-! (2) SoCal (used for GJI paper)
-! (3) regular mesh on land
-! (4) regular mesh
-! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
-! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
-! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ ! spatial properties of receivers
+ ! IREC_SPACE
+ ! (1) individual station(s)
+ ! (2) SoCal (used for GJI paper)
+ ! (3) regular mesh on land
+ ! (4) regular mesh
+ ! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+ ! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+ ! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
integer, parameter :: NMESH_REC = 17
- double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
- double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m
double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
-! lower right corner for membrane surface waves plotting grid
+ ! lower right corner for membrane surface waves plotting grid
double precision, parameter :: LAT_MIN = 32.0
double precision, parameter :: LON_MIN = -120.0
integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
-! mesh specifications: membrane surface waves
+ ! mesh specifications: membrane surface waves
double precision, parameter :: LENGTH = 480.0e3 ! m
double precision, parameter :: HEIGHT = 480.0e3 ! m
double precision, parameter :: AREA = LENGTH*HEIGHT
@@ -85,12 +85,12 @@
!!$ integer, parameter :: NEX = 80 ! 160
!!$ integer, parameter :: NEZ = 32 ! 32
-!========================
-! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+ !========================
+ ! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
-! model perturbations for HOMOGENEOUS model (or perturbation)
-! scaling from beta to alpha
-! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ ! model perturbations for HOMOGENEOUS model (or perturbation)
+ ! scaling from beta to alpha
+ ! value is from Masters et al. (2000), "The relative behavior of shear velocity..."
double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
double precision, parameter :: PBETA = 10.0
double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
@@ -103,61 +103,71 @@
! 0 1 2 3
!-----------------------------------------------------------------------------------------
! ISURFACE=1, IMODEL_SYN : homo checker het
- ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het
!-----------------------------------------------------------------------------------------
- ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
- ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker NA
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert NA
!-----------------------------------------------------------------------------------------
- ! whether you want to smooth the structure model for computing the synthetics
+ ! smooth the structure model or kernels
+ ! For the CG algorithm, ISMOOTH_EVENT_KERNEL --> use smoothed event kernels
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 0
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
- integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
- integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
- ! smoothing scalelength
- double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+ ! scalelength of smoothing Gaussian
+ ! GJI paper : 30,60,90
+ ! body waves : 10
+ double precision, parameter :: SIGMA_SMOOTH_KERNEL = 2.121320d4
+ double precision, parameter :: SIGMA_SMOOTH_MODEL = 1.060660d4
+ double precision, parameter :: GAMMA_SMOOTH_KERNEL = sqrt(8.0)*SIGMA_SMOOTH_KERNEL
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = sqrt(8.0)*SIGMA_SMOOTH_MODEL
-!========================
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
-! boolean parameters
-! IKER: (0) waveform
-! (1) traveltime, cross-correlation, misfit
-! (2) amplitude, cross-correlation, misfit
-! (3) traveltime, multitaper
-! (4) amplitude, multitaper
-! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
-! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ !========================
+
+ ! boolean parameters
+ ! IKER: (0) waveform
+ ! (1) traveltime, cross-correlation, misfit
+ ! (2) amplitude, cross-correlation, misfit
+ ! (3) traveltime, multitaper
+ ! (4) amplitude, multitaper
+ ! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+ ! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
integer, parameter :: IKER = 1
integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
-! KEY: USE ONE LINE OR THE OTHER
+ ! KEY: USE ONE LINE OR THE OTHER
integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
-! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+ ! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
-! parameters controlling what to write to file
-! NOTE: for the tomography simulations, ALL of these can be .false.
+ ! parameters controlling what to write to file
+ ! NOTE: for the tomography simulations, ALL of these can be .false.
logical, parameter :: WRITE_STF_F = .false.
logical, parameter :: WRITE_SEISMO_F = .false.
-! logical, parameter :: WRITE_SPECTRA_F = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+ ! logical, parameter :: WRITE_SPECTRA_F = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
-! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+ ! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
logical, parameter :: WRITE_STF_A = .false.
logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
-! logical, parameter :: WRITE_SPECTRA_A = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+ ! logical, parameter :: WRITE_SPECTRA_A = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
-
-!--------------------------------------
-! INVERSION PARAMETERS
+ !--------------------------------------
+ ! INVERSION PARAMETERS
+
! whether you want to compute kernels, or simply the misfit function
logical, parameter :: COMPUTE_KERNELS = .true.
@@ -171,6 +181,8 @@
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
integer, parameter :: NITERATION = 16
double precision, parameter :: VAR_RED_MIN = 8.0
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
! Gaussian errors containted in input file
! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
@@ -178,9 +190,6 @@
double precision, parameter :: SIGMA_DLNA = 1.0
double precision, parameter :: SIGMA_WAVEFORM = 1.0
logical, parameter :: ADD_DATA_ERRORS = .true.
-
- !double precision, parameter :: SIGMA_FAC = 2.0
- !double precision, parameter :: CONV_STOP = 1.0e-4
! order of interpolating polynomial in conjugate gradient algorithm
! using POLY_ORDER = 3 required computing the gradient of the test model
@@ -190,114 +199,115 @@
!integer, parameter :: GJI_PAPER = 1
! what to perturb, what to invert
- integer, parameter :: PERT_STRUCT = 1
+ ! (For the inverse tests, we only allow perturbations in beta.)
+ integer, parameter :: PERT_STRUCT_BETA = 1
integer, parameter :: PERT_SOURCE_T = 1
integer, parameter :: PERT_SOURCE_X = 1
- integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_STRUCT_BETA = 1
integer, parameter :: INV_SOURCE_T = 1
integer, parameter :: INV_SOURCE_X = 1
! whether to include the model norm term in the misfit function, which acts like damping
logical, parameter :: INCLUDE_MODEL_NORM = .true.
-
+
! log file showing source loactions
logical, parameter :: ISOURCE_LOG = .true.
-
+
! DO NOT CHANGE THESE
- integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ ! integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_STRUCT = 1 ! beta only (24-Jan-2010)
integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
-! parameterization for structure in the inversion
-! 1 : kappa-mu-rho
-! 2 : alpha-beta-rho
-! 3 : c-beta-rho
+ ! parameterization for structure in the inversion
+ ! 1 : kappa-mu-rho
+ ! 2 : alpha-beta-rho
+ ! 3 : c-beta-rho
integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
-! scalelength of smoothing Gaussian
-! GJI paper : 30,60,90
-! body waves : 10
- double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+ ! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 4.50e10 ! Pa
+ double precision, parameter :: RIGIDITY = 3.185e10 ! Pa
+ ! double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ ! double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
- ! parameters for smoothing
- logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
- logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+ ! compute additional parameters
+ double precision, parameter :: PWAVESPEED = &
+ sqrt( (INCOMPRESSIBILITY + (4.0/3.0)*RIGIDITY)/DENSITY )
+ double precision, parameter :: SWAVESPEED = sqrt( RIGIDITY/DENSITY )
+ double precision, parameter :: BWAVESPEED = sqrt( INCOMPRESSIBILITY/DENSITY )
-! homogeneous background model (S.I. units)
- double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
- double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
- double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+ !---------------------------------------------------------------
+ ! measurement windows
-!---------------------------------------------------------------
-! measurement windows
-
double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
-!---------------------------------------------------------------
-! CHT: do not change these
+ !---------------------------------------------------------------
+ ! CHT: do not change these
-! UTM zone for Southern California region
-! integer, parameter :: UTM_PROJECTION_ZONE = 11
+ ! UTM zone for Southern California region
+ ! integer, parameter :: UTM_PROJECTION_ZONE = 11
-! to suppress UTM projection for SCEC benchmarks
+ ! to suppress UTM projection for SCEC benchmarks
logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
-! flag for projection from latitude/longitude to UTM, and back
+ ! flag for projection from latitude/longitude to UTM, and back
integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
-! flag for projection from latitude/longitude to mesh-UTM, and back
+ ! flag for projection from latitude/longitude to mesh-UTM, and back
integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
-! max number of fake receivers
+ ! max number of fake receivers
integer, parameter :: MAX_SR_FAKE = 1000
-! max number of events, receivers, and phass
+ ! max number of events, receivers, and phass
integer, parameter :: MAX_EVENT = 50
integer, parameter :: MAX_SR = 1400
integer, parameter :: MAX_PHASE = 1
integer, parameter :: MAX_COMP = NCOMP
-!========================
-! GRID AND GLL POINTS
+ !========================
+ ! GRID AND GLL POINTS
integer, parameter :: NELE = MAX(NEX,NEZ)
integer, parameter :: NSPEC = NEX*NEZ
-! number of GLL points (polynomial degree plus one)
+ ! number of GLL points (polynomial degree plus one)
integer, parameter :: NGLLX = 5
integer, parameter :: NGLLZ = 5
integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
-! number of global points
+ ! number of global points
integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
-! number of nodes for 2D and 3D shape functions for hexahedra
-! we use 8-node mesh bricks, which are more stable than 27-node elements
+ ! number of nodes for 2D and 3D shape functions for hexahedra
+ ! we use 8-node mesh bricks, which are more stable than 27-node elements
integer, parameter :: NGNOD = 8, NGNOD2D = 4
-
-! number of iterations to solve the system for xi and eta
+
+ ! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 1
-! very large and very small values
+ ! very large and very small values
double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
-! for the Gauss-Lobatto-Legendre points and weights
+ ! for the Gauss-Lobatto-Legendre points and weights
double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
-!
-! CONSTANTS
-!
+ !
+ ! CONSTANTS
+ !
double precision, parameter :: PI = 3.141592653589793
double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
double precision, parameter :: ONE_THIRD = 1.0 / 3.0
double precision, parameter :: ONEOVERTWO = 0.5
double precision, parameter :: DEG = 180.0/PI
-! double precision, parameter :: EPS = 1.0e-35
+ ! double precision, parameter :: EPS = 1.0e-35
!!$! parameter for FFTW
!!$ integer, parameter :: NOUT = NSTEP/2 + 1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex06.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex06.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex06.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -1,28 +1,28 @@
module wave2d_constants
-! This file is copied into the output directory when wave2d.f90 is run
-! There are two basic options: membrane (surface) waves and body waves.
-! Several lines need to be commented/uncommented to go from one to the other.
+ ! This file is copied into the output directory when wave2d.f90 is run
+ ! There are two basic options: membrane (surface) waves and body waves.
+ ! Several lines need to be commented/uncommented to go from one to the other.
-! index of the reference directory for the simulation output
+ ! index of the reference directory for the simulation output
integer, parameter :: IRUNZ = 600
-!========================
-! GRID, TIME-STEP, AND SOURCE PARAMETERS
+ !========================
+ ! GRID, TIME-STEP, AND SOURCE PARAMETERS
-! NFRAME : number of frames to save: membrane, 10; body, 20
-! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
-! NSTEP : number of timesteps
+ ! NFRAME : number of frames to save: membrane, 10; body, 20
+ ! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+ ! NSTEP : number of timesteps
integer, parameter :: NFRAME = 10 ! membrane
!integer, parameter :: NFRAME = 11 ! body
integer, parameter :: NSAVE = 400 ! 200,400
integer, parameter :: NSTEP = NFRAME*NSAVE
-! time step in seconds
+ ! time step in seconds
!double precision, parameter :: DT = 2.0e-2 ! body waves
double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
-! temporal properties of source (source time function)
+ ! temporal properties of source (source time function)
integer, parameter :: ISRC_TIME = 1 ! type (1)
double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
!double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
@@ -30,49 +30,49 @@
!double precision, parameter :: tshift = 8.0*hdur
logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
-! normalization factor of point source force
+ ! normalization factor of point source force
double precision, parameter :: FNORM = 1.0d10
-! forward wavefield source-time function
-! (For membrane waves, FOR_Y = 1 is all that matters.)
+ ! forward wavefield source-time function
+ ! (For membrane waves, FOR_Y = 1 is all that matters.)
integer, parameter :: FOR_X = 1 ! boolean
integer, parameter :: FOR_Y = 1 ! boolean
integer, parameter :: FOR_Z = 1 ! boolean
-! adjoint wavefield source-time function
-! (For membrane waves, REV_Y = 1 is all that matters.)
-! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ ! adjoint wavefield source-time function
+ ! (For membrane waves, REV_Y = 1 is all that matters.)
+ ! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
integer, parameter :: REV_X = 1 ! boolean
integer, parameter :: REV_Y = 1 ! boolean
integer, parameter :: REV_Z = 1 ! boolean
-! spatial properties of sources
-! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
-! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
-! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ ! spatial properties of sources
+ ! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+ ! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+ ! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
-! spatial properties of receivers
-! IREC_SPACE
-! (1) individual station(s)
-! (2) SoCal (used for GJI paper)
-! (3) regular mesh on land
-! (4) regular mesh
-! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
-! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
-! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ ! spatial properties of receivers
+ ! IREC_SPACE
+ ! (1) individual station(s)
+ ! (2) SoCal (used for GJI paper)
+ ! (3) regular mesh on land
+ ! (4) regular mesh
+ ! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+ ! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+ ! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
integer, parameter :: NMESH_REC = 17
- double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
- double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m
double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
-! lower right corner for membrane surface waves plotting grid
+ ! lower right corner for membrane surface waves plotting grid
double precision, parameter :: LAT_MIN = 32.0
double precision, parameter :: LON_MIN = -120.0
integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
-! mesh specifications: membrane surface waves
+ ! mesh specifications: membrane surface waves
double precision, parameter :: LENGTH = 480.0e3 ! m
double precision, parameter :: HEIGHT = 480.0e3 ! m
double precision, parameter :: AREA = LENGTH*HEIGHT
@@ -85,12 +85,12 @@
!!$ integer, parameter :: NEX = 80 ! 160
!!$ integer, parameter :: NEZ = 32 ! 32
-!========================
-! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+ !========================
+ ! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
-! model perturbations for HOMOGENEOUS model (or perturbation)
-! scaling from beta to alpha
-! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ ! model perturbations for HOMOGENEOUS model (or perturbation)
+ ! scaling from beta to alpha
+ ! value is from Masters et al. (2000), "The relative behavior of shear velocity..."
double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
double precision, parameter :: PBETA = 10.0
double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
@@ -103,61 +103,71 @@
! 0 1 2 3
!-----------------------------------------------------------------------------------------
! ISURFACE=1, IMODEL_SYN : homo checker het
- ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het
!-----------------------------------------------------------------------------------------
- ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
- ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker NA
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert NA
!-----------------------------------------------------------------------------------------
- ! whether you want to smooth the structure model for computing the synthetics
+ ! smooth the structure model or kernels
+ ! For the CG algorithm, ISMOOTH_EVENT_KERNEL --> use smoothed event kernels
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 0
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
- integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
- integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
- ! smoothing scalelength
- double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+ ! scalelength of smoothing Gaussian
+ ! GJI paper : 30,60,90
+ ! body waves : 10
+ double precision, parameter :: SIGMA_SMOOTH_KERNEL = 2.121320d4
+ double precision, parameter :: SIGMA_SMOOTH_MODEL = 1.060660d4
+ double precision, parameter :: GAMMA_SMOOTH_KERNEL = sqrt(8.0)*SIGMA_SMOOTH_KERNEL
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = sqrt(8.0)*SIGMA_SMOOTH_MODEL
-!========================
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
-! boolean parameters
-! IKER: (0) waveform
-! (1) traveltime, cross-correlation, misfit
-! (2) amplitude, cross-correlation, misfit
-! (3) traveltime, multitaper
-! (4) amplitude, multitaper
-! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
-! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ !========================
+
+ ! boolean parameters
+ ! IKER: (0) waveform
+ ! (1) traveltime, cross-correlation, misfit
+ ! (2) amplitude, cross-correlation, misfit
+ ! (3) traveltime, multitaper
+ ! (4) amplitude, multitaper
+ ! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+ ! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
integer, parameter :: IKER = 1
integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
-! KEY: USE ONE LINE OR THE OTHER
+ ! KEY: USE ONE LINE OR THE OTHER
integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
-! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+ ! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
-! parameters controlling what to write to file
-! NOTE: for the tomography simulations, ALL of these can be .false.
+ ! parameters controlling what to write to file
+ ! NOTE: for the tomography simulations, ALL of these can be .false.
logical, parameter :: WRITE_STF_F = .false.
logical, parameter :: WRITE_SEISMO_F = .false.
-! logical, parameter :: WRITE_SPECTRA_F = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+ ! logical, parameter :: WRITE_SPECTRA_F = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
-! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+ ! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
logical, parameter :: WRITE_STF_A = .false.
logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
-! logical, parameter :: WRITE_SPECTRA_A = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+ ! logical, parameter :: WRITE_SPECTRA_A = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
-
-!--------------------------------------
-! INVERSION PARAMETERS
+ !--------------------------------------
+ ! INVERSION PARAMETERS
+
! whether you want to compute kernels, or simply the misfit function
logical, parameter :: COMPUTE_KERNELS = .true.
@@ -169,8 +179,10 @@
! VAR_RED_MIN : minimum variance reduction (in percent)
! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
- integer, parameter :: NITERATION = 16
+ integer, parameter :: NITERATION = 0
double precision, parameter :: VAR_RED_MIN = 8.0
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
! Gaussian errors containted in input file
! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
@@ -178,9 +190,6 @@
double precision, parameter :: SIGMA_DLNA = 1.0
double precision, parameter :: SIGMA_WAVEFORM = 1.0
logical, parameter :: ADD_DATA_ERRORS = .true.
-
- !double precision, parameter :: SIGMA_FAC = 2.0
- !double precision, parameter :: CONV_STOP = 1.0e-4
! order of interpolating polynomial in conjugate gradient algorithm
! using POLY_ORDER = 3 required computing the gradient of the test model
@@ -190,114 +199,115 @@
!integer, parameter :: GJI_PAPER = 1
! what to perturb, what to invert
- integer, parameter :: PERT_STRUCT = 1
- integer, parameter :: PERT_SOURCE_T = 0
- integer, parameter :: PERT_SOURCE_X = 0
+ ! (For the inverse tests, we only allow perturbations in beta.)
+ integer, parameter :: PERT_STRUCT_BETA = 1
+ integer, parameter :: PERT_SOURCE_T = 1
+ integer, parameter :: PERT_SOURCE_X = 1
- integer, parameter :: INV_STRUCT = 1
- integer, parameter :: INV_SOURCE_T = 0
- integer, parameter :: INV_SOURCE_X = 0
+ integer, parameter :: INV_STRUCT_BETA = 1
+ integer, parameter :: INV_SOURCE_T = 1
+ integer, parameter :: INV_SOURCE_X = 1
! whether to include the model norm term in the misfit function, which acts like damping
logical, parameter :: INCLUDE_MODEL_NORM = .true.
-
+
! log file showing source loactions
logical, parameter :: ISOURCE_LOG = .true.
-
+
! DO NOT CHANGE THESE
- integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ ! integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_STRUCT = 1 ! beta only (24-Jan-2010)
integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
-! parameterization for structure in the inversion
-! 1 : kappa-mu-rho
-! 2 : alpha-beta-rho
-! 3 : c-beta-rho
+ ! parameterization for structure in the inversion
+ ! 1 : kappa-mu-rho
+ ! 2 : alpha-beta-rho
+ ! 3 : c-beta-rho
integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
-! scalelength of smoothing Gaussian
-! GJI paper : 30,60,90
-! body waves : 10
- double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+ ! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 4.50e10 ! Pa
+ double precision, parameter :: RIGIDITY = 3.185e10 ! Pa
+ ! double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ ! double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
- ! parameters for smoothing
- logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
- logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+ ! compute additional parameters
+ double precision, parameter :: PWAVESPEED = &
+ sqrt( (INCOMPRESSIBILITY + (4.0/3.0)*RIGIDITY)/DENSITY )
+ double precision, parameter :: SWAVESPEED = sqrt( RIGIDITY/DENSITY )
+ double precision, parameter :: BWAVESPEED = sqrt( INCOMPRESSIBILITY/DENSITY )
-! homogeneous background model (S.I. units)
- double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
- double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
- double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+ !---------------------------------------------------------------
+ ! measurement windows
-!---------------------------------------------------------------
-! measurement windows
-
double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
-!---------------------------------------------------------------
-! CHT: do not change these
+ !---------------------------------------------------------------
+ ! CHT: do not change these
-! UTM zone for Southern California region
-! integer, parameter :: UTM_PROJECTION_ZONE = 11
+ ! UTM zone for Southern California region
+ ! integer, parameter :: UTM_PROJECTION_ZONE = 11
-! to suppress UTM projection for SCEC benchmarks
+ ! to suppress UTM projection for SCEC benchmarks
logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
-! flag for projection from latitude/longitude to UTM, and back
+ ! flag for projection from latitude/longitude to UTM, and back
integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
-! flag for projection from latitude/longitude to mesh-UTM, and back
+ ! flag for projection from latitude/longitude to mesh-UTM, and back
integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
-! max number of fake receivers
+ ! max number of fake receivers
integer, parameter :: MAX_SR_FAKE = 1000
-! max number of events, receivers, and phass
+ ! max number of events, receivers, and phass
integer, parameter :: MAX_EVENT = 50
integer, parameter :: MAX_SR = 1400
integer, parameter :: MAX_PHASE = 1
integer, parameter :: MAX_COMP = NCOMP
-!========================
-! GRID AND GLL POINTS
+ !========================
+ ! GRID AND GLL POINTS
integer, parameter :: NELE = MAX(NEX,NEZ)
integer, parameter :: NSPEC = NEX*NEZ
-! number of GLL points (polynomial degree plus one)
+ ! number of GLL points (polynomial degree plus one)
integer, parameter :: NGLLX = 5
integer, parameter :: NGLLZ = 5
integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
-! number of global points
+ ! number of global points
integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
-! number of nodes for 2D and 3D shape functions for hexahedra
-! we use 8-node mesh bricks, which are more stable than 27-node elements
+ ! number of nodes for 2D and 3D shape functions for hexahedra
+ ! we use 8-node mesh bricks, which are more stable than 27-node elements
integer, parameter :: NGNOD = 8, NGNOD2D = 4
-
-! number of iterations to solve the system for xi and eta
+
+ ! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 1
-! very large and very small values
+ ! very large and very small values
double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
-! for the Gauss-Lobatto-Legendre points and weights
+ ! for the Gauss-Lobatto-Legendre points and weights
double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
-!
-! CONSTANTS
-!
+ !
+ ! CONSTANTS
+ !
double precision, parameter :: PI = 3.141592653589793
double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
double precision, parameter :: ONE_THIRD = 1.0 / 3.0
double precision, parameter :: ONEOVERTWO = 0.5
double precision, parameter :: DEG = 180.0/PI
-! double precision, parameter :: EPS = 1.0e-35
+ ! double precision, parameter :: EPS = 1.0e-35
!!$! parameter for FFTW
!!$ integer, parameter :: NOUT = NSTEP/2 + 1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex07.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex07.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex07.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -1,28 +1,28 @@
module wave2d_constants
-! This file is copied into the output directory when wave2d.f90 is run
-! There are two basic options: membrane (surface) waves and body waves.
-! Several lines need to be commented/uncommented to go from one to the other.
+ ! This file is copied into the output directory when wave2d.f90 is run
+ ! There are two basic options: membrane (surface) waves and body waves.
+ ! Several lines need to be commented/uncommented to go from one to the other.
-! index of the reference directory for the simulation output
+ ! index of the reference directory for the simulation output
integer, parameter :: IRUNZ = 700
-!========================
-! GRID, TIME-STEP, AND SOURCE PARAMETERS
+ !========================
+ ! GRID, TIME-STEP, AND SOURCE PARAMETERS
-! NFRAME : number of frames to save: membrane, 10; body, 20
-! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
-! NSTEP : number of timesteps
+ ! NFRAME : number of frames to save: membrane, 10; body, 20
+ ! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+ ! NSTEP : number of timesteps
integer, parameter :: NFRAME = 10 ! membrane
!integer, parameter :: NFRAME = 11 ! body
integer, parameter :: NSAVE = 400 ! 200,400
integer, parameter :: NSTEP = NFRAME*NSAVE
-! time step in seconds
+ ! time step in seconds
!double precision, parameter :: DT = 2.0e-2 ! body waves
double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
-! temporal properties of source (source time function)
+ ! temporal properties of source (source time function)
integer, parameter :: ISRC_TIME = 1 ! type (1)
double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
!double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
@@ -30,49 +30,49 @@
!double precision, parameter :: tshift = 8.0*hdur
logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
-! normalization factor of point source force
+ ! normalization factor of point source force
double precision, parameter :: FNORM = 1.0d10
-! forward wavefield source-time function
-! (For membrane waves, FOR_Y = 1 is all that matters.)
+ ! forward wavefield source-time function
+ ! (For membrane waves, FOR_Y = 1 is all that matters.)
integer, parameter :: FOR_X = 1 ! boolean
integer, parameter :: FOR_Y = 1 ! boolean
integer, parameter :: FOR_Z = 1 ! boolean
-! adjoint wavefield source-time function
-! (For membrane waves, REV_Y = 1 is all that matters.)
-! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ ! adjoint wavefield source-time function
+ ! (For membrane waves, REV_Y = 1 is all that matters.)
+ ! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
integer, parameter :: REV_X = 1 ! boolean
integer, parameter :: REV_Y = 1 ! boolean
integer, parameter :: REV_Z = 1 ! boolean
-! spatial properties of sources
-! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
-! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
-! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ ! spatial properties of sources
+ ! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+ ! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+ ! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
-! spatial properties of receivers
-! IREC_SPACE
-! (1) individual station(s)
-! (2) SoCal (used for GJI paper)
-! (3) regular mesh on land
-! (4) regular mesh
-! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
-! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
-! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ ! spatial properties of receivers
+ ! IREC_SPACE
+ ! (1) individual station(s)
+ ! (2) SoCal (used for GJI paper)
+ ! (3) regular mesh on land
+ ! (4) regular mesh
+ ! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+ ! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+ ! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
integer, parameter :: NMESH_REC = 17
- double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
- double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m
double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
-! lower right corner for membrane surface waves plotting grid
+ ! lower right corner for membrane surface waves plotting grid
double precision, parameter :: LAT_MIN = 32.0
double precision, parameter :: LON_MIN = -120.0
integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
-! mesh specifications: membrane surface waves
+ ! mesh specifications: membrane surface waves
double precision, parameter :: LENGTH = 480.0e3 ! m
double precision, parameter :: HEIGHT = 480.0e3 ! m
double precision, parameter :: AREA = LENGTH*HEIGHT
@@ -85,12 +85,12 @@
!!$ integer, parameter :: NEX = 80 ! 160
!!$ integer, parameter :: NEZ = 32 ! 32
-!========================
-! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+ !========================
+ ! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
-! model perturbations for HOMOGENEOUS model (or perturbation)
-! scaling from beta to alpha
-! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ ! model perturbations for HOMOGENEOUS model (or perturbation)
+ ! scaling from beta to alpha
+ ! value is from Masters et al. (2000), "The relative behavior of shear velocity..."
double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
double precision, parameter :: PBETA = 10.0
double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
@@ -103,61 +103,71 @@
! 0 1 2 3
!-----------------------------------------------------------------------------------------
! ISURFACE=1, IMODEL_SYN : homo checker het
- ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het
!-----------------------------------------------------------------------------------------
- ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
- ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker NA
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert NA
!-----------------------------------------------------------------------------------------
- ! whether you want to smooth the structure model for computing the synthetics
+ ! smooth the structure model or kernels
+ ! For the CG algorithm, ISMOOTH_EVENT_KERNEL --> use smoothed event kernels
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 0
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
- integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
- integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
- ! smoothing scalelength
- double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+ ! scalelength of smoothing Gaussian
+ ! GJI paper : 30,60,90
+ ! body waves : 10
+ double precision, parameter :: SIGMA_SMOOTH_KERNEL = 2.121320d4
+ double precision, parameter :: SIGMA_SMOOTH_MODEL = 1.060660d4
+ double precision, parameter :: GAMMA_SMOOTH_KERNEL = sqrt(8.0)*SIGMA_SMOOTH_KERNEL
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = sqrt(8.0)*SIGMA_SMOOTH_MODEL
-!========================
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
-! boolean parameters
-! IKER: (0) waveform
-! (1) traveltime, cross-correlation, misfit
-! (2) amplitude, cross-correlation, misfit
-! (3) traveltime, multitaper
-! (4) amplitude, multitaper
-! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
-! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ !========================
+
+ ! boolean parameters
+ ! IKER: (0) waveform
+ ! (1) traveltime, cross-correlation, misfit
+ ! (2) amplitude, cross-correlation, misfit
+ ! (3) traveltime, multitaper
+ ! (4) amplitude, multitaper
+ ! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+ ! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
integer, parameter :: IKER = 1
integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
-! KEY: USE ONE LINE OR THE OTHER
+ ! KEY: USE ONE LINE OR THE OTHER
integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
-! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+ ! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
-! parameters controlling what to write to file
-! NOTE: for the tomography simulations, ALL of these can be .false.
+ ! parameters controlling what to write to file
+ ! NOTE: for the tomography simulations, ALL of these can be .false.
logical, parameter :: WRITE_STF_F = .false.
logical, parameter :: WRITE_SEISMO_F = .false.
-! logical, parameter :: WRITE_SPECTRA_F = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+ ! logical, parameter :: WRITE_SPECTRA_F = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
-! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+ ! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
logical, parameter :: WRITE_STF_A = .false.
logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
-! logical, parameter :: WRITE_SPECTRA_A = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+ ! logical, parameter :: WRITE_SPECTRA_A = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
-
-!--------------------------------------
-! INVERSION PARAMETERS
+ !--------------------------------------
+ ! INVERSION PARAMETERS
+
! whether you want to compute kernels, or simply the misfit function
logical, parameter :: COMPUTE_KERNELS = .true.
@@ -169,8 +179,10 @@
! VAR_RED_MIN : minimum variance reduction (in percent)
! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
- integer, parameter :: NITERATION = 0
+ integer, parameter :: NITERATION = 16
double precision, parameter :: VAR_RED_MIN = 8.0
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
! Gaussian errors containted in input file
! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
@@ -178,9 +190,6 @@
double precision, parameter :: SIGMA_DLNA = 1.0
double precision, parameter :: SIGMA_WAVEFORM = 1.0
logical, parameter :: ADD_DATA_ERRORS = .true.
-
- !double precision, parameter :: SIGMA_FAC = 2.0
- !double precision, parameter :: CONV_STOP = 1.0e-4
! order of interpolating polynomial in conjugate gradient algorithm
! using POLY_ORDER = 3 required computing the gradient of the test model
@@ -190,114 +199,115 @@
!integer, parameter :: GJI_PAPER = 1
! what to perturb, what to invert
- integer, parameter :: PERT_STRUCT = 1
- integer, parameter :: PERT_SOURCE_T = 0
- integer, parameter :: PERT_SOURCE_X = 0
+ ! (For the inverse tests, we only allow perturbations in beta.)
+ integer, parameter :: PERT_STRUCT_BETA = 1
+ integer, parameter :: PERT_SOURCE_T = 1
+ integer, parameter :: PERT_SOURCE_X = 1
- integer, parameter :: INV_STRUCT = 1
- integer, parameter :: INV_SOURCE_T = 0
- integer, parameter :: INV_SOURCE_X = 0
+ integer, parameter :: INV_STRUCT_BETA = 1
+ integer, parameter :: INV_SOURCE_T = 1
+ integer, parameter :: INV_SOURCE_X = 1
! whether to include the model norm term in the misfit function, which acts like damping
logical, parameter :: INCLUDE_MODEL_NORM = .true.
-
+
! log file showing source loactions
logical, parameter :: ISOURCE_LOG = .true.
-
+
! DO NOT CHANGE THESE
- integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ ! integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_STRUCT = 1 ! beta only (24-Jan-2010)
integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
-! parameterization for structure in the inversion
-! 1 : kappa-mu-rho
-! 2 : alpha-beta-rho
-! 3 : c-beta-rho
+ ! parameterization for structure in the inversion
+ ! 1 : kappa-mu-rho
+ ! 2 : alpha-beta-rho
+ ! 3 : c-beta-rho
integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
-! scalelength of smoothing Gaussian
-! GJI paper : 30,60,90
-! body waves : 10
- double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+ ! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 4.50e10 ! Pa
+ double precision, parameter :: RIGIDITY = 3.185e10 ! Pa
+ ! double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ ! double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
- ! parameters for smoothing
- logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
- logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+ ! compute additional parameters
+ double precision, parameter :: PWAVESPEED = &
+ sqrt( (INCOMPRESSIBILITY + (4.0/3.0)*RIGIDITY)/DENSITY )
+ double precision, parameter :: SWAVESPEED = sqrt( RIGIDITY/DENSITY )
+ double precision, parameter :: BWAVESPEED = sqrt( INCOMPRESSIBILITY/DENSITY )
-! homogeneous background model (S.I. units)
- double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
- double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
- double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+ !---------------------------------------------------------------
+ ! measurement windows
-!---------------------------------------------------------------
-! measurement windows
-
double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
-!---------------------------------------------------------------
-! CHT: do not change these
+ !---------------------------------------------------------------
+ ! CHT: do not change these
-! UTM zone for Southern California region
-! integer, parameter :: UTM_PROJECTION_ZONE = 11
+ ! UTM zone for Southern California region
+ ! integer, parameter :: UTM_PROJECTION_ZONE = 11
-! to suppress UTM projection for SCEC benchmarks
+ ! to suppress UTM projection for SCEC benchmarks
logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
-! flag for projection from latitude/longitude to UTM, and back
+ ! flag for projection from latitude/longitude to UTM, and back
integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
-! flag for projection from latitude/longitude to mesh-UTM, and back
+ ! flag for projection from latitude/longitude to mesh-UTM, and back
integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
-! max number of fake receivers
+ ! max number of fake receivers
integer, parameter :: MAX_SR_FAKE = 1000
-! max number of events, receivers, and phass
+ ! max number of events, receivers, and phass
integer, parameter :: MAX_EVENT = 50
integer, parameter :: MAX_SR = 1400
integer, parameter :: MAX_PHASE = 1
integer, parameter :: MAX_COMP = NCOMP
-!========================
-! GRID AND GLL POINTS
+ !========================
+ ! GRID AND GLL POINTS
integer, parameter :: NELE = MAX(NEX,NEZ)
integer, parameter :: NSPEC = NEX*NEZ
-! number of GLL points (polynomial degree plus one)
+ ! number of GLL points (polynomial degree plus one)
integer, parameter :: NGLLX = 5
integer, parameter :: NGLLZ = 5
integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
-! number of global points
+ ! number of global points
integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
-! number of nodes for 2D and 3D shape functions for hexahedra
-! we use 8-node mesh bricks, which are more stable than 27-node elements
+ ! number of nodes for 2D and 3D shape functions for hexahedra
+ ! we use 8-node mesh bricks, which are more stable than 27-node elements
integer, parameter :: NGNOD = 8, NGNOD2D = 4
-
-! number of iterations to solve the system for xi and eta
+
+ ! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 1
-! very large and very small values
+ ! very large and very small values
double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
-! for the Gauss-Lobatto-Legendre points and weights
+ ! for the Gauss-Lobatto-Legendre points and weights
double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
-!
-! CONSTANTS
-!
+ !
+ ! CONSTANTS
+ !
double precision, parameter :: PI = 3.141592653589793
double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
double precision, parameter :: ONE_THIRD = 1.0 / 3.0
double precision, parameter :: ONEOVERTWO = 0.5
double precision, parameter :: DEG = 180.0/PI
-! double precision, parameter :: EPS = 1.0e-35
+ ! double precision, parameter :: EPS = 1.0e-35
!!$! parameter for FFTW
!!$ integer, parameter :: NOUT = NSTEP/2 + 1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex08.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -1,28 +1,28 @@
module wave2d_constants
-! This file is copied into the output directory when wave2d.f90 is run
-! There are two basic options: membrane (surface) waves and body waves.
-! Several lines need to be commented/uncommented to go from one to the other.
+ ! This file is copied into the output directory when wave2d.f90 is run
+ ! There are two basic options: membrane (surface) waves and body waves.
+ ! Several lines need to be commented/uncommented to go from one to the other.
-! index of the reference directory for the simulation output
+ ! index of the reference directory for the simulation output
integer, parameter :: IRUNZ = 800
-!========================
-! GRID, TIME-STEP, AND SOURCE PARAMETERS
+ !========================
+ ! GRID, TIME-STEP, AND SOURCE PARAMETERS
-! NFRAME : number of frames to save: membrane, 10; body, 20
-! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
-! NSTEP : number of timesteps
+ ! NFRAME : number of frames to save: membrane, 10; body, 20
+ ! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+ ! NSTEP : number of timesteps
integer, parameter :: NFRAME = 10 ! membrane
!integer, parameter :: NFRAME = 11 ! body
integer, parameter :: NSAVE = 400 ! 200,400
integer, parameter :: NSTEP = NFRAME*NSAVE
-! time step in seconds
+ ! time step in seconds
!double precision, parameter :: DT = 2.0e-2 ! body waves
double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
-! temporal properties of source (source time function)
+ ! temporal properties of source (source time function)
integer, parameter :: ISRC_TIME = 1 ! type (1)
double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
!double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
@@ -30,49 +30,49 @@
!double precision, parameter :: tshift = 8.0*hdur
logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
-! normalization factor of point source force
+ ! normalization factor of point source force
double precision, parameter :: FNORM = 1.0d10
-! forward wavefield source-time function
-! (For membrane waves, FOR_Y = 1 is all that matters.)
+ ! forward wavefield source-time function
+ ! (For membrane waves, FOR_Y = 1 is all that matters.)
integer, parameter :: FOR_X = 1 ! boolean
integer, parameter :: FOR_Y = 1 ! boolean
integer, parameter :: FOR_Z = 1 ! boolean
-! adjoint wavefield source-time function
-! (For membrane waves, REV_Y = 1 is all that matters.)
-! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ ! adjoint wavefield source-time function
+ ! (For membrane waves, REV_Y = 1 is all that matters.)
+ ! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
integer, parameter :: REV_X = 1 ! boolean
integer, parameter :: REV_Y = 1 ! boolean
integer, parameter :: REV_Z = 1 ! boolean
-! spatial properties of sources
-! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
-! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
-! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ ! spatial properties of sources
+ ! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+ ! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+ ! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
-! spatial properties of receivers
-! IREC_SPACE
-! (1) individual station(s)
-! (2) SoCal (used for GJI paper)
-! (3) regular mesh on land
-! (4) regular mesh
-! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
-! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
-! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ ! spatial properties of receivers
+ ! IREC_SPACE
+ ! (1) individual station(s)
+ ! (2) SoCal (used for GJI paper)
+ ! (3) regular mesh on land
+ ! (4) regular mesh
+ ! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+ ! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+ ! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
integer, parameter :: NMESH_REC = 17
- double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
- double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m
double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
-! lower right corner for membrane surface waves plotting grid
+ ! lower right corner for membrane surface waves plotting grid
double precision, parameter :: LAT_MIN = 32.0
double precision, parameter :: LON_MIN = -120.0
integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
-! mesh specifications: membrane surface waves
+ ! mesh specifications: membrane surface waves
double precision, parameter :: LENGTH = 480.0e3 ! m
double precision, parameter :: HEIGHT = 480.0e3 ! m
double precision, parameter :: AREA = LENGTH*HEIGHT
@@ -85,12 +85,12 @@
!!$ integer, parameter :: NEX = 80 ! 160
!!$ integer, parameter :: NEZ = 32 ! 32
-!========================
-! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+ !========================
+ ! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
-! model perturbations for HOMOGENEOUS model (or perturbation)
-! scaling from beta to alpha
-! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ ! model perturbations for HOMOGENEOUS model (or perturbation)
+ ! scaling from beta to alpha
+ ! value is from Masters et al. (2000), "The relative behavior of shear velocity..."
double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
double precision, parameter :: PBETA = 10.0
double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
@@ -103,61 +103,71 @@
! 0 1 2 3
!-----------------------------------------------------------------------------------------
! ISURFACE=1, IMODEL_SYN : homo checker het
- ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het
!-----------------------------------------------------------------------------------------
- ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
- ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker NA
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert NA
!-----------------------------------------------------------------------------------------
- ! whether you want to smooth the structure model for computing the synthetics
+ ! smooth the structure model or kernels
+ ! For the CG algorithm, ISMOOTH_EVENT_KERNEL --> use smoothed event kernels
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 0
integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
- integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
- integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
- ! smoothing scalelength
- double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+ ! scalelength of smoothing Gaussian
+ ! GJI paper : 30,60,90
+ ! body waves : 10
+ double precision, parameter :: SIGMA_SMOOTH_KERNEL = 2.121320d4
+ double precision, parameter :: SIGMA_SMOOTH_MODEL = 1.060660d4
+ double precision, parameter :: GAMMA_SMOOTH_KERNEL = sqrt(8.0)*SIGMA_SMOOTH_KERNEL
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = sqrt(8.0)*SIGMA_SMOOTH_MODEL
-!========================
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
-! boolean parameters
-! IKER: (0) waveform
-! (1) traveltime, cross-correlation, misfit
-! (2) amplitude, cross-correlation, misfit
-! (3) traveltime, multitaper
-! (4) amplitude, multitaper
-! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
-! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ !========================
+
+ ! boolean parameters
+ ! IKER: (0) waveform
+ ! (1) traveltime, cross-correlation, misfit
+ ! (2) amplitude, cross-correlation, misfit
+ ! (3) traveltime, multitaper
+ ! (4) amplitude, multitaper
+ ! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+ ! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
integer, parameter :: IKER = 1
integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
-! KEY: USE ONE LINE OR THE OTHER
+ ! KEY: USE ONE LINE OR THE OTHER
integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
-! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+ ! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
-! parameters controlling what to write to file
-! NOTE: for the tomography simulations, ALL of these can be .false.
+ ! parameters controlling what to write to file
+ ! NOTE: for the tomography simulations, ALL of these can be .false.
logical, parameter :: WRITE_STF_F = .false.
logical, parameter :: WRITE_SEISMO_F = .false.
-! logical, parameter :: WRITE_SPECTRA_F = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+ ! logical, parameter :: WRITE_SPECTRA_F = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
-! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+ ! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
logical, parameter :: WRITE_STF_A = .false.
logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
-! logical, parameter :: WRITE_SPECTRA_A = .false.
-! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+ ! logical, parameter :: WRITE_SPECTRA_A = .false.
+ ! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
-
-!--------------------------------------
-! INVERSION PARAMETERS
+ !--------------------------------------
+ ! INVERSION PARAMETERS
+
! whether you want to compute kernels, or simply the misfit function
logical, parameter :: COMPUTE_KERNELS = .true.
@@ -171,6 +181,8 @@
! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
integer, parameter :: NITERATION = 16
double precision, parameter :: VAR_RED_MIN = 8.0
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
! Gaussian errors containted in input file
! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
@@ -178,9 +190,6 @@
double precision, parameter :: SIGMA_DLNA = 1.0
double precision, parameter :: SIGMA_WAVEFORM = 1.0
logical, parameter :: ADD_DATA_ERRORS = .true.
-
- !double precision, parameter :: SIGMA_FAC = 2.0
- !double precision, parameter :: CONV_STOP = 1.0e-4
! order of interpolating polynomial in conjugate gradient algorithm
! using POLY_ORDER = 3 required computing the gradient of the test model
@@ -190,114 +199,115 @@
!integer, parameter :: GJI_PAPER = 1
! what to perturb, what to invert
- integer, parameter :: PERT_STRUCT = 1
+ ! (For the inverse tests, we only allow perturbations in beta.)
+ integer, parameter :: PERT_STRUCT_BETA = 1
integer, parameter :: PERT_SOURCE_T = 0
integer, parameter :: PERT_SOURCE_X = 0
- integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_STRUCT_BETA = 1
integer, parameter :: INV_SOURCE_T = 0
integer, parameter :: INV_SOURCE_X = 0
! whether to include the model norm term in the misfit function, which acts like damping
logical, parameter :: INCLUDE_MODEL_NORM = .true.
-
+
! log file showing source loactions
logical, parameter :: ISOURCE_LOG = .true.
-
+
! DO NOT CHANGE THESE
- integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ ! integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_STRUCT = 1 ! beta only (24-Jan-2010)
integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
-! parameterization for structure in the inversion
-! 1 : kappa-mu-rho
-! 2 : alpha-beta-rho
-! 3 : c-beta-rho
+ ! parameterization for structure in the inversion
+ ! 1 : kappa-mu-rho
+ ! 2 : alpha-beta-rho
+ ! 3 : c-beta-rho
integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
-! scalelength of smoothing Gaussian
-! GJI paper : 30,60,90
-! body waves : 10
- double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+ ! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 4.50e10 ! Pa
+ double precision, parameter :: RIGIDITY = 3.185e10 ! Pa
+ ! double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ ! double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
- ! parameters for smoothing
- logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
- logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+ ! compute additional parameters
+ double precision, parameter :: PWAVESPEED = &
+ sqrt( (INCOMPRESSIBILITY + (4.0/3.0)*RIGIDITY)/DENSITY )
+ double precision, parameter :: SWAVESPEED = sqrt( RIGIDITY/DENSITY )
+ double precision, parameter :: BWAVESPEED = sqrt( INCOMPRESSIBILITY/DENSITY )
-! homogeneous background model (S.I. units)
- double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
- double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
- double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+ !---------------------------------------------------------------
+ ! measurement windows
-!---------------------------------------------------------------
-! measurement windows
-
double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
-!---------------------------------------------------------------
-! CHT: do not change these
+ !---------------------------------------------------------------
+ ! CHT: do not change these
-! UTM zone for Southern California region
-! integer, parameter :: UTM_PROJECTION_ZONE = 11
+ ! UTM zone for Southern California region
+ ! integer, parameter :: UTM_PROJECTION_ZONE = 11
-! to suppress UTM projection for SCEC benchmarks
+ ! to suppress UTM projection for SCEC benchmarks
logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
-! flag for projection from latitude/longitude to UTM, and back
+ ! flag for projection from latitude/longitude to UTM, and back
integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
-! flag for projection from latitude/longitude to mesh-UTM, and back
+ ! flag for projection from latitude/longitude to mesh-UTM, and back
integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
-! max number of fake receivers
+ ! max number of fake receivers
integer, parameter :: MAX_SR_FAKE = 1000
-! max number of events, receivers, and phass
+ ! max number of events, receivers, and phass
integer, parameter :: MAX_EVENT = 50
integer, parameter :: MAX_SR = 1400
integer, parameter :: MAX_PHASE = 1
integer, parameter :: MAX_COMP = NCOMP
-!========================
-! GRID AND GLL POINTS
+ !========================
+ ! GRID AND GLL POINTS
integer, parameter :: NELE = MAX(NEX,NEZ)
integer, parameter :: NSPEC = NEX*NEZ
-! number of GLL points (polynomial degree plus one)
+ ! number of GLL points (polynomial degree plus one)
integer, parameter :: NGLLX = 5
integer, parameter :: NGLLZ = 5
integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
-! number of global points
+ ! number of global points
integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
-! number of nodes for 2D and 3D shape functions for hexahedra
-! we use 8-node mesh bricks, which are more stable than 27-node elements
+ ! number of nodes for 2D and 3D shape functions for hexahedra
+ ! we use 8-node mesh bricks, which are more stable than 27-node elements
integer, parameter :: NGNOD = 8, NGNOD2D = 4
-
-! number of iterations to solve the system for xi and eta
+
+ ! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 1
-! very large and very small values
+ ! very large and very small values
double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
-! for the Gauss-Lobatto-Legendre points and weights
+ ! for the Gauss-Lobatto-Legendre points and weights
double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
-!
-! CONSTANTS
-!
+ !
+ ! CONSTANTS
+ !
double precision, parameter :: PI = 3.141592653589793
double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
double precision, parameter :: ONE_THIRD = 1.0 / 3.0
double precision, parameter :: ONEOVERTWO = 0.5
double precision, parameter :: DEG = 180.0/PI
-! double precision, parameter :: EPS = 1.0e-35
+ ! double precision, parameter :: EPS = 1.0e-35
!!$! parameter for FFTW
!!$ integer, parameter :: NOUT = NSTEP/2 + 1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub.f90 2010-01-24 16:01:33 UTC (rev 16167)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_sub.f90 2010-01-24 23:32:07 UTC (rev 16168)
@@ -875,12 +875,13 @@
!----------------------------------------------------
! subroutine local2mvec(array1, array2, nmod_src, source_vec, nmod, mvec, array_fac0)
- subroutine local2mvec(array1, array2, nmod_src, source_vec, nmod, mvec)
+ !subroutine local2mvec(array1, array2, nmod_src, source_vec, nmod, mvec)
+ subroutine local2mvec(array1, nmod_src, source_vec, nmod, mvec)
! This inputs three arrays and creates a vector that can be used
! in the conjugate gradient algorithm.
- double precision, dimension(NGLLX,NGLLZ,NSPEC),intent(in) :: array1, array2
+ double precision, dimension(NGLLX,NGLLZ,NSPEC),intent(in) :: array1
integer, intent(in) :: nmod_src, nmod
double precision, dimension(nmod_src),intent(in) :: source_vec
double precision, dimension(nmod), intent(out) :: mvec
@@ -911,8 +912,8 @@
!if(ipar==1) mvec(k) = array1(i,j,ispec) * array_fac(i,j,ispec) ! alpha
!if(ipar==2) mvec(k) = array2(i,j,ispec) * array_fac(i,j,ispec) ! beta
- if(ipar==1) mvec(k) = array1(i,j,ispec) ! atype (kappa, alpha, c)
- if(ipar==2) mvec(k) = array2(i,j,ispec) ! btype (mu, beta, beta
+ if(ipar==1) mvec(k) = array1(i,j,ispec) ! btype (kappa, alpha, c)
+ !if(ipar==2) mvec(k) = array2(i,j,ispec) ! atype (mu, beta, beta
enddo
enddo
enddo
@@ -926,20 +927,20 @@
!----------------------------------------------------
- subroutine mvec2local(nmod, nmod_src, mvec, array1, array2, source_vec)
+ !subroutine mvec2local(nmod, nmod_src, mvec, array1, array2, source_vec)
+ subroutine mvec2local(nmod, nmod_src, mvec, array1, source_vec)
! This inputs a model vector and outputs the constituent arrays: alpha, beta, source.
integer, intent(in) :: nmod, nmod_src
double precision, dimension(nmod), intent(in) :: mvec
- double precision, dimension(NGLLX,NGLLZ,NSPEC), intent(out) :: array1, array2
+ double precision, dimension(NGLLX,NGLLZ,NSPEC), intent(out) :: array1
double precision, dimension(nmod_src), intent(out) :: source_vec
integer :: i,j,k,ipar,ispec
!----------
array1(:,:,:) = 0.0
- array2(:,:,:) = 0.0
source_vec(:) = 0.0
! structure arrays
@@ -949,8 +950,8 @@
do j = 1,NGLLZ
do i = 1,NGLLX
k = k+1
- if(ipar==1) array1(i,j,ispec) = mvec(k) ! atype (kappa, alpha, c)
- if(ipar==2) array2(i,j,ispec) = mvec(k) ! btype (mu, beta, beta)
+ if(ipar==1) array1(i,j,ispec) = mvec(k) ! btype (mu, beta, beta)
+ !if(ipar==2) array2(i,j,ispec) = mvec(k) ! atype (kappa, alpha, c)
enddo
enddo
enddo
@@ -986,8 +987,7 @@
norm_parts(:) = 0.0
! structure part of the norm -- BETA only
- norm_parts(1) = 0.0
- norm_parts(2) = sum( mvec(NLOCAL+1 : 2*NLOCAL)**2 / cov_model(NLOCAL+1 : 2*NLOCAL) )
+ norm_parts(1) = sum( mvec(1 : NLOCAL)**2 / cov_model(1 : NLOCAL) )
! source part of the norm -- only events that you are inverting for
do ievent = ievent_min, ievent_max
@@ -995,18 +995,18 @@
!itemp2 = 2*NLOCAL + (ievent-1)*3 + 2
!itemp3 = 2*NLOCAL + (ievent-1)*3 + 3
- itemp1 = 2*NLOCAL + index_source(1,ievent)
- itemp2 = 2*NLOCAL + index_source(2,ievent)
- itemp3 = 2*NLOCAL + index_source(3,ievent)
+ itemp1 = NLOCAL + index_source(1,ievent)
+ itemp2 = NLOCAL + index_source(2,ievent)
+ itemp3 = NLOCAL + index_source(3,ievent)
- norm_parts(3) = norm_parts(3) + mvec(itemp1)**2 / cov_model(itemp1)
- norm_parts(4) = norm_parts(4) + mvec(itemp2)**2 / cov_model(itemp2)
- norm_parts(5) = norm_parts(5) + mvec(itemp3)**2 / cov_model(itemp3)
+ norm_parts(2) = norm_parts(2) + mvec(itemp1)**2 / cov_model(itemp1)
+ norm_parts(3) = norm_parts(3) + mvec(itemp2)**2 / cov_model(itemp2)
+ norm_parts(4) = norm_parts(4) + mvec(itemp3)**2 / cov_model(itemp3)
enddo
- norm_struct = sum( norm_parts(1:2) )
- norm_source = sum( norm_parts(3:5) )
- norm_total = sum( norm_parts(1:5) )
+ norm_struct = norm_parts(1)
+ norm_source = sum( norm_parts(2:4) )
+ norm_total = sum( norm_parts(1:4) )
end subroutine compute_norm_sq
@@ -1018,13 +1018,14 @@
character(len=200), intent(in) :: filename
integer :: i
+ ! norm of each component
open(unit=19,file=filename,status='unknown')
do i = 1,NVAR
write(19,'(1i10,3e20.10)') i, norm_parts(i)
enddo
- write(19,'(1i10,1e20.10)') 1, sum( norm_parts(1:2) )
- write(19,'(1i10,1e20.10)') 2, sum( norm_parts(3:5) )
- write(19,'(1i10,1e20.10)') 1, sum( norm_parts(1:5) )
+ write(19,'(1i10,1e20.10)') 1, norm_parts(1) ! structure
+ write(19,'(1i10,1e20.10)') 2, sum( norm_parts(2:4) ) ! source
+ write(19,'(1i10,1e20.10)') 1, sum( norm_parts(1:4) ) ! total
end subroutine write_norm_sq
More information about the CIG-COMMITS
mailing list