[cig-commits] r16243 - in seismo/3D/ADJOINT_TOMO/iterate_adj/pangu: . smooth smooth_all sum_kernel sum_kernel/src topo_input/combine_vol_data
carltape at geodynamics.org
carltape at geodynamics.org
Tue Feb 9 15:44:07 PST 2010
Author: carltape
Date: 2010-02-09 15:44:05 -0800 (Tue, 09 Feb 2010)
New Revision: 16243
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/setup_sum.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/Makefile
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/README
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/combine_vol_data_mod.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/compile.bash
Removed:
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/gll_library.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/mpif.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/combine_vol_data.f90
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/README
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/go.bash
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/run.lsf
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth_all/run_dir.bash
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth_all/setup_smooth_all.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/go.bash
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/run.lsf
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/config.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/constants.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/exit_mpi.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/precision.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/sum_kernels
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/sum_kernels.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/values_from_mesher.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/config.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/constants.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/values_from_mesher.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/write_c_binary.c
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/xcombine_vol_data
Log:
Minor modifications to MPI scripts for operations within an adjoint-based tomographic inversion.
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/README 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/README 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,42 +1,67 @@
Carl Tape and Qinya Liu
-20-Oct-2008
+09-Nov-2009
From CIG SVN server:
/ADJOINT_TOMO/iterate_adj/pangu/README
This is designed to provide guidance for doing an inversion using event kernels
-derived from adjoint methods. At this time, THIS IS NOT PARTICULARLY USER FRIENDLY.
-Email me (carltape at gps.caltech.edu) with suggestions for improvements.
+derived from adjoint methods. Programs are contained for the basic operations needed for a tomographic inversion based on adjoint methods. However, several minor adjustments need to be made for each user.
+Email me (carltape at fas.harvard.edu) with suggestions for improvements.
+
Directions for user:
(1) Create a working directory
cd /ADJOINT_TOMO/iterate_adj
cp -r pangu pangu_work
+
(2) Copy slice file into topo_input/
-(3) Link your mesh topology files into topo/
-(4) Link your smoothing directory to smooth/kernel_smoothing
+
+(3) Compile combine_vol_data.f90 in /topo_input/combine_vol_data/
+
+(4) Link your mesh topology files into topo/
+ For example, for slice 162, we have these files:
+ proc000162_AVS_DXelements.txt
+ proc000162_AVS_DXpoints.txt
+ proc000162_ibool.bin
+ proc000162_jacobian2D.bin
+ proc000162_jacobian.bin
+ proc000162_x.bin
+ proc000162_y.bin
+ proc000162_z.bin
+ The Jacobian is used for taking the norm of the model vectors.
+
+(5) Link your smoothing directory to smooth/kernel_smoothing
cd smooth
- ln -s /ibrixfs1/scratch/carltape/kernel_smoothing_basin_new src
-(5) Link your kernels list
+ ln -s /---/ADJOINT_TOMO_EXTRA/kernel_smoothing_basin/smooth_sem_fun .
+
+(6) Link your kernels list
cd KERNELS_MODELS
- ln -s /ibrixfs1/scratch/carltape/KERNELS_MODELS/kernel_lists/kernels_m04 kernels_run
-(6) Link unsmoothed kernels into approriate directory
+ ln -s /---/KERNELS_MODELS/kernels_run .
+
+(7) Link unsmoothed kernels into approriate directory
cd KERNELS_MODELS/event_kernels/kernel_m04
ln -s /ibrixfs1/scratch/carltape/KERNELS_MODELS/event_kernels/kernel_m04/* .
-(7) Link present model into appropriate directory
+
+(8) Link present model into appropriate directory
cd KERNELS_MODELS/models/m04
ln -s /ibrixfs1/scratch/carltape/KERNELS_MODELS/models/m04/* .
-----------------
+NOTE: With mpich_gm, we needed to use the library file mpif.h.
+ With openmpi, this file is automatically present.
+ openmpi is our default and therefore no mpif.h file is present locally.
+-----------------
+
OPERATION DIRECTORIES:
-sum_kernel -- sum event kernels
-smooth -- smooth summed kernel or model update
-smooth_all -- setup directories for smoothing many event kernels at once
+sum_kernel -- sum event kernels
+smooth -- smooth summed kernel or model update
+smooth_all -- setup directories for smoothing many event kernels at once
subspace_hessian -- compute hessian for subspace method
-subspace_update -- compute model update, after computing mu vectors in Matlab
-model_vp_vs -- add a model update to make a new vp and vs model
+subspace_update -- compute model update, after computing mu vectors in Matlab
+model_vp_vs -- add a model update to make a new vp and vs model
+model_slice -- nearest neighbor algorithm for grid for plotting cross sections
cg_step -- compute the conjugate-gradient test model
@@ -52,15 +77,16 @@
-----------------
-Some programs may be linked to executables within these directories:
+See also Caltech executables in these Caltech directories:
/ibrixfs1/opt/seismo-util/source/kernel_MPI/
kernel_boun_int/
-kernel_integration/
kernel_resample/
kernel_rotate/
kernel_smoothing/
kernel_smoothing_basin/
+kernel_smooth_V4.0/
+kernel_vol_int/
sum_kernels/
-----------------
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/go.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/go.bash 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/go.bash 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,6 +1,6 @@
#!/bin/bash -v
#BSUB -o OUTPUT_FILES/%J.o
-#BSUB -a mpich_gm
+#BSUB -a openmpi
#BSUB -J smooth_model
echo "$LSB_MCPU_HOSTS" > OUTPUT_FILES/lsf_machines
@@ -19,12 +19,14 @@
# the last two entries are the directories for (1) the unsmoothed and smoothed files and (2) topology
-mpirun.lsf --gm-no-shmem --gm-copy-env smooth_sem_fun 14 12 0 2 6 1 mu_kernel inout_smooth topo
-sleep 5s
+#--------------------------------------------------
-mpirun.lsf --gm-no-shmem --gm-copy-env smooth_sem_fun 14 12 0 2 6 1 kappa_kernel inout_smooth topo
-sleep 5s
+#mpirun.lsf --gm-no-shmem --gm-copy-env smooth_sem_fun 14 12 0 2 6 1 mu_kernel inout_smooth topo
+#sleep 5s
+#mpirun.lsf --gm-no-shmem --gm-copy-env smooth_sem_fun 14 12 0 2 6 1 kappa_kernel inout_smooth topo
+#sleep 5s
+
#xcombine_vol_data slice_file mu_kernel_smooth topo inout_smooth . 0
#sleep 5s
#mv mu_kernel_smooth.mesh mu_kernel_smooth_h006km_v001km.mesh
@@ -34,3 +36,32 @@
#sleep 5s
#mv kappa_kernel_smooth.mesh kappa_kernel_smooth_h006km_v001km.mesh
#sleep 5s
+
+#--------------------------------------------
+
+# new command for Harvard cluster
+#mpirun smooth_sem_fun 14 12 0 2 12 2 vs_m16 inout_smooth topo
+#sleep 5s
+#xcombine_vol_data slice_file vs_m16_smooth topo inout_smooth . 0
+#sleep 5s
+#mv vs_m16_smooth.mesh vs_m16_smooth_h012km_v002km.mesh
+#sleep 5s
+
+# new command for Harvard cluster
+#mpirun smooth_sem_fun 14 12 0 2 12 2 vp_m16 inout_smooth topo
+#sleep 5s
+#xcombine_vol_data slice_file vp_m16_smooth topo inout_smooth . 0
+#sleep 5s
+#mv vp_m16_smooth.mesh vp_m16_smooth_h012km_v002km.mesh
+#sleep 5s
+
+#--------------------------------------------
+
+# new command for Harvard cluster
+mpirun smooth_sem_fun 14 12 0 2 4 1 rho_cbr_kernel inout_smooth topo
+sleep 5s
+xcombine_vol_data slice_file rho_cbr_kernel_smooth topo inout_smooth . 0
+sleep 5s
+mv rho_cbr_kernel_smooth.mesh rho_cbr_kernel_smooth_h004km_v001km.mesh
+sleep 5s
+
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/run.lsf
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/run.lsf 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/run.lsf 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,8 +1,8 @@
#!/bin/bash
-#compile.bash
+
mkdir -p OUTPUT_FILES
rm OUTPUT_FILES/*
date
-#bsub -q test -W 30 -n 168 < go.bash
-bsub -q normal -W 90 -n 168 < go.bash
+bsub -q short_parallel -W 30 -n 168 < go.bash
+#bsub -q sophrosyne -W 30 -n 168 -R span[ptile=8] < go.bash
date
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth_all/run_dir.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth_all/run_dir.bash 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth_all/run_dir.bash 2010-02-09 23:44:05 UTC (rev 16243)
@@ -11,7 +11,9 @@
cd $eid
echo "==== $eid ====="
- run.lsf
+ #run.lsf
+ run.lsf_hero
+ #sleep 15m
sleep 5s
done
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth_all/setup_smooth_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth_all/setup_smooth_all.pl 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth_all/setup_smooth_all.pl 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,32 +1,34 @@
#!/usr/bin/perl -w
#-----------------------------------
-# Carl Tape, 29-Jan-2009
+# Carl Tape, 10-Nov-2009
# setup_smooth_all.pl
#
-# This script copies directories with binary files of a specified parameter
-# from one place to a local directory on pangu.
-#
+# This script takes a template smoothing directory and sets up many copies.
#
# EXAMPLE:
-# setup_smooth_all.pl m14 168 mu_kernel ; setup_smooth_all.pl m14 168 kappa_kernel
+# setup_smooth_all.pl m16 168 rho_cbr_kernel
#
#-----------------------------------
if (@ARGV < 3) {die("Usage: setup_smooth_all.pl smodel nproc ftag\n")}
($smodel,$nproc,$ftag) = @ARGV;
-$basedir = "/ibrixfs1/home/carltape/ADJOINT_TOMO/iterate_adj/pangu_work";
-$edir = "$basedir/KERNELS_MODELS/event_kernels/kernel_${smodel}";
-$edir_smooth = "$basedir/KERNELS_MODELS/event_kernels_smooth/kernel_${smodel}";
-$dir_done = "/net/sierra/raid1/carltape/socal/socal_3D/RUNS/SMOOTH_EVENT_KERNELS/$smodel";
+#==========================================================
+# USER INPUT: will depend on how files are organized
+$basedir = "/n/home/ctape/ADJOINT_TOMO/iterate_adj/pangu_work";
+$edirall = "/n/nobackup1/shaw_lab/ctape/SEM_RUNS/kernels_${smodel}_imaging";
+$edirall = "${edirall}_target1";
+
+$dtag1 = "SEM/MESH_T003_T010";
+$dtag2 = "collect";
+
+#==========================================================
+
$masterdir = "$basedir/smooth";
-
-if(not -e $edir) {die("edir $edir does not exist");}
-if(not -e $edir_smooth) {`mkdir $edir_smooth`;}
+if(not -e $edirall) {die("edirall $edirall does not exist");}
if(not -e $masterdir) {die("masterdir $masterdir does not exist");}
-if(not -e ${dir_done}) {die("dir_done $dir_done does not exist");}
# read in the list of kernels that you want to sum
$kernel_list = "$basedir/KERNELS_MODELS/kernels_run";
@@ -37,69 +39,82 @@
print "number of event kernels : $neid\n";
$imin = 1; $imax = $neid; # default
-#$imin = 1; $imax = 20;
+#$imin = 1; $imax = 2;
#$imin = 14; $imax = $imin;
for ($i = $imin; $i <= $imax; $i++) {
$eid = $eids[$i-1]; chomp($eid);
print "======== $eid -- $i out of $imax ========\n";
- $edir_in = "$edir/$eid";
- ${efiles_in} = "${edir_in}/*${ftag}.bin";
- $nfiles_in = `ls -1 $efiles_in | wc | awk '{print \$1}'`; chomp($nfiles_in);
+
+ # unsmoothed files
+ $edir0 = "$edirall/$eid/$dtag1";
+ $edir = "$edir0/$dtag2";
+ $efiles_in = "$edir/*${ftag}.bin";
+ $nfiles_in = `ls -1 ${efiles_in} | wc | awk '{print \$1}'`; chomp($nfiles_in);
+ if ($nfiles_in == $nproc) {$boolin = 1;} else {$boolin = 0;}
print "nfiles_in = ${nfiles_in}\n";
- if (${nfiles_in} != $nproc ) {
- #die("no $nproc files in ${efiles_in}");
- print "no $nproc files in ${efiles_in}\n";
+ # smoothed files
+ $efiles_out = "$edir/*${ftag}_smooth.bin";
+ $nfiles_out = `ls -1 ${efiles_out} | wc | awk '{print \$1}'`; chomp($nfiles_out);
+ if ($nfiles_out == $nproc) {$boolout = 1;} else {$boolout = 0;}
+ print "nfiles_out = ${nfiles_out}\n";
- } else {
+ # local run directory
+ $edir_local = "$eid";
- $edir_done = "${dir_done}/$eid";
- if ( (-e "${dir_done}/$eid") || ( (-e "TO_COMBINE/$eid") || (-e "TO_COMBINE/DONE/$eid") ) ) {
- print "--> kernel is already smoothed\n";
+ # Four possible cases:
+ # 1. both smoothed and unsmoothed are done
+ # 2. smoothed done, but missing unsmoothed
+ # 3. neither unsmoothed nor smoothed done
+ # 4. unsmoothed done, but not smoothed ==> SETUP
- } else {
+ if ($boolin==1 && $boolout==1) {
- $edir_copy = "$eid";
- $edir_out = "${edir_copy}/inout_smooth";
- $efiles_out = "${edir_out}/*${ftag}.bin";
- $nfiles_out = `ls -1 ${efiles_out} | wc | awk '{print \$1}'`; chomp($nfiles_out);
- print "nfiles_out = ${nfiles_out}\n";
+ # check if mesh files have been made
+ $mesh_local = "${edir_local}/${ftag}_smooth*";
+ $mesh_out = "$edir0/${ftag}_smooth*";
+ $nmesh_local = `ls -1 ${mesh_local} | wc | awk '{print \$1}'`; chomp($nmesh_local);
+ $nmesh_out = `ls -1 ${mesh_out} | wc | awk '{print \$1}'`; chomp($nmesh_out);
+ if ($nmesh_out == 1) {
+ print "--> DONE\n";
+ } elsif ( $nmesh_out==0 && $nmesh_local==1 ) {
+ `cp $mesh_local $edir0`; `sleep 1s`;
+ } elsif ( $nmesh_out==0 && $nmesh_local==0 ) {
+ print "--> STILL NEED TO COMBINE BIN FILES\n";
+ }
- if ( ${nfiles_out} == $nproc ) {
- print "--> $nfiles_out bin files exist -- ready to go\n";
+ } elsif ($boolin==0 && $boolout==1) {
+ print "missing unsmoothed files, but smoothed files are present\n";
- } else {
+ } elsif ($boolin==0 && $boolout==0) {
+ print "no unsmoothed files or smoothed files\n";
- `mkdir -p $edir_smooth/$eid`;
- `mkdir -p $edir_copy`;
- `mkdir -p $eid/OUTPUT_FILES`;
- `mkdir -p $eid/inout_smooth`;
+ } else {
+ print "unsmoothed files present but no smoothed files\n";
- $sfiles = "${edir_out}/*${ftag}_smooth.bin";
- $nsfiles = `ls -1 $sfiles | wc | awk '{print \$1}'`; chomp($nsfiles);
- if ( $nsfiles == $nproc ) {
- print "--> $nsfiles bin files already exist -- no need to run\n";
+ if (-e $edir_local) {
+ print "--> READY TO EXECUTE SMOOTHING\n";
+ } else {
- } else {
- # link input event kernels
- `ln -s $efiles_in $edir_out`;
+ `mkdir -p ${edir_local}`;
+ `mkdir -p ${edir_local}/OUTPUT_FILES`;
+ `ln -s $edir ${edir_local}/inout_smooth`; # KEY LINK
+ # NOTE: the link leaves the danger of over-writing input files
- # copy run scripts and link executable file
- `cp $masterdir/run.lsf $edir_copy`;
- `cp $masterdir/go.bash $edir_copy`;
- `ln -s $masterdir/topo $edir_copy`;
- `ln -s $masterdir/smooth_sem_fun $edir_copy`;
- `ln -s $masterdir/xcombine_vol_data $edir_copy`;
- `ln -s $masterdir/slice_file $edir_copy`;
- `ln -s $masterdir/combine.bash $edir_copy`;
+ # copy run scripts and link executable file
+ `cp $masterdir/run* $edir_local`;
+ `cp $masterdir/go.bash $edir_local`;
+ `ln -s $masterdir/topo $edir_local`;
+ `ln -s $masterdir/smooth_sem_fun $edir_local`;
+ `ln -s $masterdir/xcombine_vol_data $edir_local`;
+ `ln -s $masterdir/slice_file $edir_local`;
+ `ln -s $masterdir/combine.bash $edir_local`;
- # change the label for the run
- `sed "/smooth_model/s/smooth_model/${eid}_smooth/" ${edir_copy}/go.bash > ${edir_copy}/run.tmp`;
- `mv -f ${edir_copy}/run.tmp ${edir_copy}/go.bash`;
- }
- }
+ # change the label for the run
+ `sed "/smooth_model/s/smooth_model/${eid}S/" ${edir_local}/go.bash > ${edir_local}/run.tmp`;
+ `mv -f ${edir_local}/run.tmp ${edir_local}/go.bash`;
}
}
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/go.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/go.bash 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/go.bash 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,11 +1,15 @@
#!/bin/bash -v
#BSUB -o OUTPUT_FILES/%J.o
-#BSUB -a mpich_gm
+#BSUB -a openmpi
#BSUB -J sum_kernels
echo "$LSB_MCPU_HOSTS" > OUTPUT_FILES/lsf_machines
echo "$LSB_JOBID" > OUTPUT_FILES/jobid
remap_lsf_machines.pl OUTPUT_FILES/lsf_machines >OUTPUT_FILES/machines
-mpirun.lsf --gm-no-shmem --gm-copy-env src/sum_kernels
+#mpirun src/sum_kernels
+mpirun.lsf src/sum_kernels
+# Caltech command, used with mpich_gm
+#mpirun.lsf --gm-no-shmem --gm-copy-env src/sum_kernels
+
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/run.lsf
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/run.lsf 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/run.lsf 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,8 +1,7 @@
#!/bin/bash
-#compile.bash
mkdir -p OUTPUT_FILES
rm OUTPUT_FILES/*
date
-#bsub -q normal -W 30 -n 168 < go.bash
-bsub -q debug -W 30 -n 168 < go.bash
+#bsub -q short_parallel -W 30 -n 168 < go.bash
+bsub -q sophrosyne -W 30 -n 168 -R span[ptile=8] < go.bash
date
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/setup_sum.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/setup_sum.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/setup_sum.pl 2010-02-09 23:44:05 UTC (rev 16243)
@@ -0,0 +1,79 @@
+#!/usr/bin/perl -w
+
+#-----------------------------------
+# Carl Tape, 10-Nov-2009
+# setup_sum.pl
+#
+# Example script for setting up this directory; this will vary for each user.
+#
+# EXAMPLE:
+# setup_sum.pl 168 rho_cbr_kernel
+#
+#-----------------------------------
+
+if (@ARGV < 2) {die("Usage: setup_sum.pl xxx\n")}
+($nproc,$ktag) = @ARGV;
+
+#-------------------------------------
+# USER INPUT
+
+# directory containing the SEM runs
+$basedir = "/n/nobackup1/shaw_lab/ctape/SEM_RUNS/kernels_m16_imaging";
+if(not -e $basedir) {die("basedir $basedir does not exist");}
+
+$dtag = "SEM/MESH_T003_T010/collect";
+
+# defaults
+$idir = "INPUT_KERNELS";
+$kernel_list = "kernels_run";
+$kernel_list_sub = "${kernel_list}_sub";
+$otag = "ftags";
+
+#-------------------------------------
+
+# write kernel tag to be read in
+open(OUT,">$otag"); print OUT "$ktag"; close(OUT);
+
+# delete all symbolic links in input directory
+`rm -rf $idir/*`; `sleep 1s`;
+
+# read in the list of kernels that you want to sum
+if (not -e ${kernel_list}) {die("check if kernel_list $kernel_list exist or not\n")}
+open(IN,"${kernel_list}"); @eids = <IN>; close(IN);
+$neid = @eids;
+
+# open list for subset of EIDs
+open(KEIDS,">${kernel_list_sub}");
+
+$imin = 1; $imax = $neid; # default
+#$imin = 1; $imax = 2;
+#$imin = 18; $imax = $imin;
+
+for ($i = $imin; $i <= $imax; $i++) {
+
+ $eid = $eids[$i-1]; chomp($eid);
+ print "======== $eid -- $i out of $imax ========\n";
+ $edir_in = "$basedir/$eid/$dtag";
+ if(not -e $edir_in) {die("edir_in $edir_in does not exist");}
+
+ # print EID to file
+ print KEIDS "$eid\n";
+
+ # link directory containing .bin files to the local directory
+ `ln -s $edir_in $idir/$eid`;
+
+ # check that all the linked files exist
+ $efiles = "$idir/$eid/*${ktag}*";
+ $nfiles = `ls -1 $efiles | wc | awk '{print \$1}'`; chomp($nfiles);
+ print "nfiles = ${nfiles} for $ktag\n";
+ if ($nfiles != $nproc) {die("you do not have $nproc $ktag files setup");}
+
+}
+
+close(KEIDS);
+
+# instructions to make mesh file
+print "to combine into a low-res mesh, use this command:\n";
+print " xcombine_vol_data ../topo_input/slice_file $ktag topo OUTPUT_SUM . 0\n";
+
+#=================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/setup_sum.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/Makefile
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/Makefile (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/Makefile 2010-02-09 23:44:05 UTC (rev 16243)
@@ -0,0 +1,13 @@
+FC=mpif90
+FFLAGS=-O3
+
+OBJS = sum_kernels.o exit_mpi.o
+sum_kernels: $(OBJS)
+ $(FC) $(FFLAGS) -o $@ $(OBJS)
+
+.SUFFIXES: .o .f90
+.f90.o:
+ $(FC) $(FFLAGS) -c $<
+
+clean:
+ rm -rf *.o *~ sum_kernels
\ No newline at end of file
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/config.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/config.h 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/config.h 2010-02-09 23:44:05 UTC (rev 16243)
@@ -23,13 +23,13 @@
#define PACKAGE_BUGREPORT "jtromp AT caltech.edu"
/* Define to the full name of this package. */
-#define PACKAGE_NAME "Specfem 3D Basin"
+#define PACKAGE_NAME "Specfem 3D"
/* Define to the full name and version of this package. */
-#define PACKAGE_STRING "Specfem 3D Basin 1.4.0"
+#define PACKAGE_STRING "Specfem 3D 1.4.0"
/* Define to the one symbol short name of this package. */
-#define PACKAGE_TARNAME "Specfem3DBasin"
+#define PACKAGE_TARNAME "Specfem3D"
/* Define to the version of this package. */
#define PACKAGE_VERSION "1.4.0"
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/constants.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/constants.h 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/constants.h 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,18 +1,26 @@
!=====================================================================
!
-! S p e c f e m 3 D B a s i n V e r s i o n 1 . 4
-! --------------------------------------------------
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
!
! Dimitri Komatitsch and Jeroen Tromp
! Seismological Laboratory - California Institute of Technology
! (c) California Institute of Technology July 2005
!
-! A signed non-commercial agreement is required to use this program.
-! Please check http://www.gps.caltech.edu/research/jtromp for details.
-! Free for non-commercial academic research ONLY.
-! This program is distributed WITHOUT ANY WARRANTY whatsoever.
-! Do not redistribute this program without written permission.
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
!=====================================================================
! constants.h. Generated from constants.h.in by configure.
@@ -65,7 +73,7 @@
! density of sea water
real(kind=CUSTOM_REAL), parameter :: RHO_OCEANS = 1020.0
-! extend basin model below threshold and above topography to make sure
+! extend model below threshold and above topography to make sure
! there is no small gap between interpolated maps and sediments
logical, parameter :: EXTEND_VOXET_BELOW_BASEMENT = .true.
logical, parameter :: EXTEND_VOXET_ABOVE_TOPO = .true.
@@ -85,24 +93,27 @@
! number of sources to be gathered by MPI_Gather
integer, parameter :: NGATHER_SOURCES = 10000
-! source decay rate
- double precision, parameter :: SOURCE_DECAY_RATE = 1.628d0
+! we mimic a triangle of half duration equal to half_duration_triangle
+! using a Gaussian having a very close shape, as explained in Figure 4.2
+! of the manual. This source decay rate to mimic an equivalent triangle
+! was found by trial and error
+ double precision, parameter :: SOURCE_DECAY_MIMIC_TRIANGLE = 1.628d0
! ---------------------------------------------------------------------------------------
-! LQY -- Following 3 variables stays here temporarily,
+! LQY -- Following 3 variables stays here temporarily,
! we need to move them to Par_file at a proper time
! ---------------------------------------------------------------------------------------
! save moho mesh and compute Moho boundary kernels
- logical, parameter :: SAVE_MOHO_MESH = .true.
+ logical, parameter :: SAVE_MOHO_MESH = .false.
! number of steps to save the state variables in the forward simulation,
! to be used in the backward reconstruction in the presence of attenuation
integer, parameter :: NSTEP_Q_SAVE = 200
-! the scratch disk to save the state variables saved in the forward
+! the scratch disk to save the state variables saved in the forward
! simulation, this can be a global scratch disk in case you run out of
! space on the local scratch disk
- character(len=150), parameter :: LOCAL_PATH_Q = '/ibrixfs1/scratch/lqy'
+ character(len=150), parameter :: LOCAL_PATH_Q = '/ibrixfs1/scratch/lqy/DATABASES_MPI_Q/'
!------------------------------------------------------
!----------- do not modify anything below -------------
@@ -223,13 +234,6 @@
double precision, parameter :: DEGREES_PER_CELL_TOPO_SOCAL = 5.d0 / 1000.d0
character(len=100), parameter :: TOPO_FILE_SOCAL = 'DATA/la_topography/topo_bathy_final.dat'
-! size of topography and bathymetry file for Lacq (France) gas field
- integer, parameter :: NX_TOPO_LACQ = 499,NY_TOPO_LACQ = 401
- double precision, parameter :: ORIG_LAT_TOPO_LACQ = 100000.d0
- double precision, parameter :: ORIG_LONG_TOPO_LACQ = 340400.d0
- double precision, parameter :: DEGREES_PER_CELL_TOPO_LACQ = 100.d0
- character(len=100), parameter :: TOPO_FILE_LACQ = 'DATA/lacq_thomas/mnt_Lacq_Lambert_final_dimitri.dat'
-
! size of Lupei Zhu's Moho map file for Southern California
integer, parameter :: NX_MOHO = 71,NY_MOHO = 51
double precision, parameter :: ORIG_LAT_MOHO = 32.d0
@@ -309,7 +313,7 @@
!
! integer, parameter :: NGRID_NEW_HAUKSSON = 201
!
-!! corners of new Hauksson's interpolated grid
+!! corners of new Hauksson interpolated grid
! double precision, parameter :: UTM_X_ORIG_HAUKSSON = 122035.012d0
! double precision, parameter :: UTM_X_END_HAUKSSON = 766968.628d0
! double precision, parameter :: UTM_Y_ORIG_HAUKSSON = 3547232.986d0
@@ -330,7 +334,7 @@
integer, parameter :: NGRID_NEW_HAUKSSON = 241
-! corners of new Hauksson's interpolated grid
+! corners of new Hauksson interpolated grid
double precision, parameter :: UTM_X_ORIG_HAUKSSON = 88021.568d0
double precision, parameter :: UTM_X_END_HAUKSSON = 861517.886d0
double precision, parameter :: UTM_Y_ORIG_HAUKSSON = 3404059.875d0
@@ -347,13 +351,7 @@
double precision, parameter :: DEPTH_5p5km_SOCAL = -5500.d0
double precision, parameter :: DEPTH_16km_SOCAL = -16000.d0
double precision, parameter :: DEPTH_MOHO_SOCAL = -32000.d0
-! double precision, parameter :: DEPTH_5p5km_SOCAL = -7700.d0
-! double precision, parameter :: DEPTH_16km_SOCAL = -18200.d0
-! double precision, parameter :: DEPTH_MOHO_SOCAL = -34200.d0
-! layer for Lacq (France) gas field
- double precision, parameter :: DEPTH_INTERFACE_LACQ = -7000.d0
-
! reference surface of the model before adding topography
double precision, parameter :: Z_SURFACE = 0.d0
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/exit_mpi.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/exit_mpi.f90 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/exit_mpi.f90 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,30 +1,34 @@
!=====================================================================
!
-! S p e c f e m 3 D G l o b e V e r s i o n 3 . 6
-! --------------------------------------------------
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
!
! Dimitri Komatitsch and Jeroen Tromp
! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
+! (c) California Institute of Technology September 2006
!
-! A signed non-commercial agreement is required to use this program.
-! Please check http://www.gps.caltech.edu/research/jtromp for details.
-! Free for non-commercial academic research ONLY.
-! This program is distributed WITHOUT ANY WARRANTY whatsoever.
-! Do not redistribute this program without written permission.
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
!=====================================================================
! end the simulation and exit MPI
-! version with rank number printed in the error message
subroutine exit_MPI(myrank,error_msg)
implicit none
-! standard include of the MPI library
- include 'mpif.h'
-
include "constants.h"
! identifier for error message file
@@ -35,21 +39,27 @@
integer ier
character(len=80) outputname
+ !character(len=150) OUTPUT_FILES
! write error message to screen
write(*,*) error_msg(1:len(error_msg))
write(*,*) 'Error detected, aborting MPI... proc ',myrank
! write error message to file
+
+ !call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
+ !write(outputname,"('/error_message',i6.6,'.txt')") myrank
+ !open(unit=IERROR,file=trim(OUTPUT_FILES)//outputname,status='unknown')
write(outputname,"('OUTPUT_FILES/error_message',i6.6,'.txt')") myrank
open(unit=IERROR,file=trim(outputname),status='unknown')
+
write(IERROR,*) error_msg(1:len(error_msg))
write(IERROR,*) 'Error detected, aborting MPI... proc ',myrank
close(IERROR)
-
-! stop all the MPI processes, and exit
-! on some machines, MPI_FINALIZE needs to be called before MPI_ABORT
+! close output file
+! if(myrank == 0 .and. IMAIN /= ISTANDARD_OUTPUT) close(IMAIN)
+! call stop_all()
call MPI_FINALIZE(ier)
call MPI_ABORT(ier)
stop 'error, program ended in exit_MPI'
@@ -59,3 +69,25 @@
!
!----
!
+
+! version without rank number printed in the error message
+ subroutine exit_MPI_without_rank(error_msg)
+
+ implicit none
+
+ include "constants.h"
+
+ integer ier
+ character(len=*) error_msg
+
+! write error message to screen
+ write(*,*) error_msg(1:len(error_msg))
+ write(*,*) 'Error detected, aborting MPI...'
+
+ !call stop_all()
+ call MPI_FINALIZE(ier)
+ call MPI_ABORT(ier)
+ stop 'error, program ended in exit_MPI'
+
+ end subroutine exit_MPI_without_rank
+
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/gll_library.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/gll_library.f90 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/gll_library.f90 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,529 +0,0 @@
-
-!=======================================================================
-!
-! Library to compute the Gauss-Lobatto-Legendre points and weights
-! Based on Gauss-Lobatto routines from M.I.T.
-! Department of Mechanical Engineering
-!
-!=======================================================================
-
- double precision function endw1(n,alpha,beta)
-
- implicit none
-
- integer n
- double precision alpha,beta
-
- double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0,three=3.d0,four=4.d0
- double precision apb,f1,fint1,fint2,f2,di,abn,abnn,a1,a2,a3,f3
- double precision, external :: gammaf
- integer i
-
- f3 = zero
- apb = alpha+beta
- if (n == 0) then
- endw1 = zero
- return
- endif
- f1 = gammaf(alpha+two)*gammaf(beta+one)/gammaf(apb+three)
- f1 = f1*(apb+two)*two**(apb+two)/two
- if (n == 1) then
- endw1 = f1
- return
- endif
- fint1 = gammaf(alpha+two)*gammaf(beta+one)/gammaf(apb+three)
- fint1 = fint1*two**(apb+two)
- fint2 = gammaf(alpha+two)*gammaf(beta+two)/gammaf(apb+four)
- fint2 = fint2*two**(apb+three)
- f2 = (-two*(beta+two)*fint1 + (apb+four)*fint2) * (apb+three)/four
- if (n == 2) then
- endw1 = f2
- return
- endif
- do i=3,n
- di = dble(i-1)
- abn = alpha+beta+di
- abnn = abn+di
- a1 = -(two*(di+alpha)*(di+beta))/(abn*abnn*(abnn+one))
- a2 = (two*(alpha-beta))/(abnn*(abnn+two))
- a3 = (two*(abn+one))/((abnn+two)*(abnn+one))
- f3 = -(a2*f2+a1*f1)/a3
- f1 = f2
- f2 = f3
- enddo
- endw1 = f3
-
- end function endw1
-
-!
-!=======================================================================
-!
-
- double precision function endw2(n,alpha,beta)
-
- implicit none
-
- integer n
- double precision alpha,beta
-
- double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0,three=3.d0,four=4.d0
- double precision apb,f1,fint1,fint2,f2,di,abn,abnn,a1,a2,a3,f3
- double precision, external :: gammaf
- integer i
-
- apb = alpha+beta
- f3 = zero
- if (n == 0) then
- endw2 = zero
- return
- endif
- f1 = gammaf(alpha+one)*gammaf(beta+two)/gammaf(apb+three)
- f1 = f1*(apb+two)*two**(apb+two)/two
- if (n == 1) then
- endw2 = f1
- return
- endif
- fint1 = gammaf(alpha+one)*gammaf(beta+two)/gammaf(apb+three)
- fint1 = fint1*two**(apb+two)
- fint2 = gammaf(alpha+two)*gammaf(beta+two)/gammaf(apb+four)
- fint2 = fint2*two**(apb+three)
- f2 = (two*(alpha+two)*fint1 - (apb+four)*fint2) * (apb+three)/four
- if (n == 2) then
- endw2 = f2
- return
- endif
- do i=3,n
- di = dble(i-1)
- abn = alpha+beta+di
- abnn = abn+di
- a1 = -(two*(di+alpha)*(di+beta))/(abn*abnn*(abnn+one))
- a2 = (two*(alpha-beta))/(abnn*(abnn+two))
- a3 = (two*(abn+one))/((abnn+two)*(abnn+one))
- f3 = -(a2*f2+a1*f1)/a3
- f1 = f2
- f2 = f3
- enddo
- endw2 = f3
-
- end function endw2
-
-!
-!=======================================================================
-!
-
- double precision function gammaf (x)
-
- implicit none
-
- double precision, parameter :: pi = 3.141592653589793d0
-
- double precision x
-
- double precision, parameter :: half=0.5d0,one=1.d0,two=2.d0
-
- gammaf = one
-
- if (x == -half) gammaf = -two*dsqrt(pi)
- if (x == half) gammaf = dsqrt(pi)
- if (x == one ) gammaf = one
- if (x == two ) gammaf = one
- if (x == 1.5d0) gammaf = dsqrt(pi)/2.d0
- if (x == 2.5d0) gammaf = 1.5d0*dsqrt(pi)/2.d0
- if (x == 3.5d0) gammaf = 2.5d0*1.5d0*dsqrt(pi)/2.d0
- if (x == 3.d0 ) gammaf = 2.d0
- if (x == 4.d0 ) gammaf = 6.d0
- if (x == 5.d0 ) gammaf = 24.d0
- if (x == 6.d0 ) gammaf = 120.d0
-
- end function gammaf
-
-!
-!=====================================================================
-!
-
- subroutine jacg (xjac,np,alpha,beta)
-
-!=======================================================================
-!
-! computes np Gauss points, which are the zeros of the
-! Jacobi polynomial with parameters alpha and beta
-!
-! .alpha = beta = 0.0 -> Legendre points
-! .alpha = beta = -0.5 -> Chebyshev points
-!
-!=======================================================================
-
- implicit none
-
- integer np
- double precision alpha,beta
- double precision xjac(np)
-
- integer k,j,i,jmin,jm,n
- double precision xlast,dth,x,x1,x2,recsum,delx,xmin,swap
- double precision p,pd,pm1,pdm1,pm2,pdm2
-
- integer, parameter :: K_MAX_ITER = 10
- double precision, parameter :: zero = 0.d0, eps = 1.0d-12
-
- pm1 = zero
- pm2 = zero
- pdm1 = zero
- pdm2 = zero
-
- xlast = 0.d0
- n = np-1
- dth = 4.d0*datan(1.d0)/(2.d0*dble(n)+2.d0)
- p = 0.d0
- pd = 0.d0
- jmin = 0
- do j=1,np
- if(j == 1) then
- x = dcos((2.d0*(dble(j)-1.d0)+1.d0)*dth)
- else
- x1 = dcos((2.d0*(dble(j)-1.d0)+1.d0)*dth)
- x2 = xlast
- x = (x1+x2)/2.d0
- endif
- do k=1,K_MAX_ITER
- call jacobf (p,pd,pm1,pdm1,pm2,pdm2,np,alpha,beta,x)
- recsum = 0.d0
- jm = j-1
- do i=1,jm
- recsum = recsum+1.d0/(x-xjac(np-i+1))
- enddo
- delx = -p/(pd-recsum*p)
- x = x+delx
- if(abs(delx) < eps) goto 31
- enddo
- 31 continue
- xjac(np-j+1) = x
- xlast = x
- enddo
- do i=1,np
- xmin = 2.d0
- do j=i,np
- if(xjac(j) < xmin) then
- xmin = xjac(j)
- jmin = j
- endif
- enddo
- if(jmin /= i) then
- swap = xjac(i)
- xjac(i) = xjac(jmin)
- xjac(jmin) = swap
- endif
- enddo
-
- end subroutine jacg
-
-!
-!=====================================================================
-!
-
- subroutine jacobf (poly,pder,polym1,pderm1,polym2,pderm2,n,alp,bet,x)
-
-!=======================================================================
-!
-! Computes the Jacobi polynomial of degree n and its derivative at x
-!
-!=======================================================================
-
- implicit none
-
- double precision poly,pder,polym1,pderm1,polym2,pderm2,alp,bet,x
- integer n
-
- double precision apb,polyl,pderl,dk,a1,a2,b3,a3,a4,polyn,pdern,psave,pdsave
- integer k
-
- apb = alp+bet
- poly = 1.d0
- pder = 0.d0
- psave = 0.d0
- pdsave = 0.d0
-
- if (n == 0) return
-
- polyl = poly
- pderl = pder
- poly = (alp-bet+(apb+2.d0)*x)/2.d0
- pder = (apb+2.d0)/2.d0
- if (n == 1) return
-
- do k=2,n
- dk = dble(k)
- a1 = 2.d0*dk*(dk+apb)*(2.d0*dk+apb-2.d0)
- a2 = (2.d0*dk+apb-1.d0)*(alp**2-bet**2)
- b3 = (2.d0*dk+apb-2.d0)
- a3 = b3*(b3+1.d0)*(b3+2.d0)
- a4 = 2.d0*(dk+alp-1.d0)*(dk+bet-1.d0)*(2.d0*dk+apb)
- polyn = ((a2+a3*x)*poly-a4*polyl)/a1
- pdern = ((a2+a3*x)*pder-a4*pderl+a3*poly)/a1
- psave = polyl
- pdsave = pderl
- polyl = poly
- poly = polyn
- pderl = pder
- pder = pdern
- enddo
-
- polym1 = polyl
- pderm1 = pderl
- polym2 = psave
- pderm2 = pdsave
-
- end subroutine jacobf
-
-!
-!------------------------------------------------------------------------
-!
-
- double precision FUNCTION PNDLEG (Z,N)
-
-!------------------------------------------------------------------------
-!
-! Compute the derivative of the Nth order Legendre polynomial at Z.
-! Based on the recursion formula for the Legendre polynomials.
-!
-!------------------------------------------------------------------------
- implicit none
-
- double precision z
- integer n
-
- double precision P1,P2,P1D,P2D,P3D,FK,P3
- integer k
-
- P1 = 1.d0
- P2 = Z
- P1D = 0.d0
- P2D = 1.d0
- P3D = 1.d0
-
- do K = 1, N-1
- FK = dble(K)
- P3 = ((2.d0*FK+1.d0)*Z*P2 - FK*P1)/(FK+1.d0)
- P3D = ((2.d0*FK+1.d0)*P2 + (2.d0*FK+1.d0)*Z*P2D - FK*P1D) / (FK+1.d0)
- P1 = P2
- P2 = P3
- P1D = P2D
- P2D = P3D
- enddo
-
- PNDLEG = P3D
-
- end function pndleg
-
-!
-!------------------------------------------------------------------------
-!
-
- double precision FUNCTION PNLEG (Z,N)
-
-!------------------------------------------------------------------------
-!
-! Compute the value of the Nth order Legendre polynomial at Z.
-! Based on the recursion formula for the Legendre polynomials.
-!
-!------------------------------------------------------------------------
- implicit none
-
- double precision z
- integer n
-
- double precision P1,P2,P3,FK
- integer k
-
- P1 = 1.d0
- P2 = Z
- P3 = P2
-
- do K = 1, N-1
- FK = dble(K)
- P3 = ((2.d0*FK+1.d0)*Z*P2 - FK*P1)/(FK+1.d0)
- P1 = P2
- P2 = P3
- enddo
-
- PNLEG = P3
-
- end function pnleg
-
-!
-!------------------------------------------------------------------------
-!
-
- double precision function pnormj (n,alpha,beta)
-
- implicit none
-
- double precision alpha,beta
- integer n
-
- double precision one,two,dn,const,prod,dindx,frac
- double precision, external :: gammaf
- integer i
-
- one = 1.d0
- two = 2.d0
- dn = dble(n)
- const = alpha+beta+one
-
- if (n <= 1) then
- prod = gammaf(dn+alpha)*gammaf(dn+beta)
- prod = prod/(gammaf(dn)*gammaf(dn+alpha+beta))
- pnormj = prod * two**const/(two*dn+const)
- return
- endif
-
- prod = gammaf(alpha+one)*gammaf(beta+one)
- prod = prod/(two*(one+const)*gammaf(const+one))
- prod = prod*(one+alpha)*(two+alpha)
- prod = prod*(one+beta)*(two+beta)
-
- do i=3,n
- dindx = dble(i)
- frac = (dindx+alpha)*(dindx+beta)/(dindx*(dindx+alpha+beta))
- prod = prod*frac
- enddo
-
- pnormj = prod * two**const/(two*dn+const)
-
- end function pnormj
-
-!
-!------------------------------------------------------------------------
-!
-
- subroutine zwgjd(z,w,np,alpha,beta)
-
-!=======================================================================
-!
-! Z w g j d : Generate np Gauss-Jacobi points and weights
-! associated with Jacobi polynomial of degree n = np-1
-!
-! Note : Coefficients alpha and beta must be greater than -1.
-! ----
-!=======================================================================
-
- implicit none
-
- double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0
-
- integer np
- double precision z(np),w(np)
- double precision alpha,beta
-
- integer n,np1,np2,i
- double precision p,pd,pm1,pdm1,pm2,pdm2
- double precision apb,dnp1,dnp2,fac1,fac2,fac3,fnorm,rcoef
- double precision, external :: gammaf,pnormj
-
- pd = zero
- pm1 = zero
- pm2 = zero
- pdm1 = zero
- pdm2 = zero
-
- n = np-1
- apb = alpha+beta
- p = zero
- pdm1 = zero
-
- if (np <= 0) stop 'minimum number of Gauss points is 1'
-
- if ((alpha <= -one) .or. (beta <= -one)) stop 'alpha and beta must be greater than -1'
-
- if (np == 1) then
- z(1) = (beta-alpha)/(apb+two)
- w(1) = gammaf(alpha+one)*gammaf(beta+one)/gammaf(apb+two) * two**(apb+one)
- return
- endif
-
- call jacg(z,np,alpha,beta)
-
- np1 = n+1
- np2 = n+2
- dnp1 = dble(np1)
- dnp2 = dble(np2)
- fac1 = dnp1+alpha+beta+one
- fac2 = fac1+dnp1
- fac3 = fac2+one
- fnorm = pnormj(np1,alpha,beta)
- rcoef = (fnorm*fac2*fac3)/(two*fac1*dnp2)
- do i=1,np
- call jacobf(p,pd,pm1,pdm1,pm2,pdm2,np2,alpha,beta,z(i))
- w(i) = -rcoef/(p*pdm1)
- enddo
-
- end subroutine zwgjd
-
-!
-!------------------------------------------------------------------------
-!
-
- subroutine zwgljd(z,w,np,alpha,beta)
-
-!=======================================================================
-!
-! Z w g l j d : Generate np Gauss-Lobatto-Jacobi points and the
-! ----------- weights associated with Jacobi polynomials of degree
-! n = np-1.
-!
-! Note : alpha and beta coefficients must be greater than -1.
-! Legendre polynomials are special case of Jacobi polynomials
-! just by setting alpha and beta to 0.
-!
-!=======================================================================
-
- implicit none
-
- double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0
-
- integer np
- double precision alpha,beta
- double precision z(np), w(np)
-
- integer n,nm1,i
- double precision p,pd,pm1,pdm1,pm2,pdm2
- double precision alpg,betg
- double precision, external :: endw1,endw2
-
- p = zero
- pm1 = zero
- pm2 = zero
- pdm1 = zero
- pdm2 = zero
-
- n = np-1
- nm1 = n-1
- pd = zero
-
- if (np <= 1) stop 'minimum number of Gauss-Lobatto points is 2'
-
-! with spectral elements, use at least 3 points
- if (np <= 2) stop 'minimum number of Gauss-Lobatto points for the SEM is 3'
-
- if ((alpha <= -one) .or. (beta <= -one)) stop 'alpha and beta must be greater than -1'
-
- if (nm1 > 0) then
- alpg = alpha+one
- betg = beta+one
- call zwgjd(z(2),w(2),nm1,alpg,betg)
- endif
-
- z(1) = - one
- z(np) = one
-
- do i=2,np-1
- w(i) = w(i)/(one-z(i)**2)
- enddo
-
- call jacobf(p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(1))
- w(1) = endw1(n,alpha,beta)/(two*pd)
- call jacobf(p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(np))
- w(np) = endw2(n,alpha,beta)/(two*pd)
-
- end subroutine zwgljd
-
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/mpif.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/mpif.h 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/mpif.h 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,265 +0,0 @@
-! mpif.h. Generated from /opt/mpich/myrinet/intel//include/mpif.h by configure.
-
-!
-!
-! (C) 1993 by Argonne National Laboratory and Mississipi State University.
-! All rights reserved. See COPYRIGHT in top-level directory.
-!
-!
-! user include file for MPI programs, with no dependencies
-!
-! It really is not possible to make a perfect include file that can
-! be used by both F77 and F90 compilers, but this is close. We have removed
-! continuation lines (allows free form input in F90); systems whose
-! Fortran compilers support ! instead of just C or * for comments can
-! globally replace a C in the first column with !; the resulting file
-! should work for both Fortran 77 and Fortran 90.
-!
-! If your Fortran compiler supports ! for comments, you can run this
-! through sed with
-! sed -e 's/^C/\!/g'
-!
-! We have also removed the use of contractions (involving the single quote)
-! character because some users use .F instead of .f files (to invoke the
-! cpp preprocessor) and further, their preprocessor is determined to find
-! matching single quote pairs (and probably double quotes; given the
-! different rules in C and Fortran, this sounds like a disaster). Rather than
-! take the position that the poor users should get a better system, we
-! have removed the text that caused problems. Of course, the users SHOULD
-! get a better system...
-!
-! return codes
- INTEGER MPI_SUCCESS,MPI_ERR_BUFFER,MPI_ERR_COUNT,MPI_ERR_TYPE
- INTEGER MPI_ERR_TAG,MPI_ERR_COMM,MPI_ERR_RANK,MPI_ERR_ROOT
- INTEGER MPI_ERR_GROUP
- INTEGER MPI_ERR_OP,MPI_ERR_TOPOLOGY,MPI_ERR_DIMS,MPI_ERR_ARG
- INTEGER MPI_ERR_UNKNOWN,MPI_ERR_TRUNCATE,MPI_ERR_OTHER
- INTEGER MPI_ERR_INTERN,MPI_ERR_IN_STATUS,MPI_ERR_PENDING
- INTEGER MPI_ERR_REQUEST, MPI_ERR_LASTCODE
- PARAMETER (MPI_SUCCESS=0,MPI_ERR_BUFFER=1,MPI_ERR_COUNT=2)
- PARAMETER (MPI_ERR_TYPE=3,MPI_ERR_TAG=4,MPI_ERR_COMM=5)
- PARAMETER (MPI_ERR_RANK=6,MPI_ERR_ROOT=7,MPI_ERR_GROUP=8)
- PARAMETER (MPI_ERR_OP=9,MPI_ERR_TOPOLOGY=10,MPI_ERR_DIMS=11)
- PARAMETER (MPI_ERR_ARG=12,MPI_ERR_UNKNOWN=13)
- PARAMETER (MPI_ERR_TRUNCATE=14,MPI_ERR_OTHER=15)
- PARAMETER (MPI_ERR_INTERN=16,MPI_ERR_IN_STATUS=17)
- PARAMETER (MPI_ERR_PENDING=18,MPI_ERR_REQUEST=19)
- PARAMETER (MPI_ERR_LASTCODE=1073741823)
-!
- INTEGER MPI_UNDEFINED
- parameter (MPI_UNDEFINED = (-32766))
-!
- INTEGER MPI_GRAPH, MPI_CART
- PARAMETER (MPI_GRAPH = 1, MPI_CART = 2)
- INTEGER MPI_PROC_NULL
- PARAMETER ( MPI_PROC_NULL = (-1) )
-!
- INTEGER MPI_BSEND_OVERHEAD
- PARAMETER ( MPI_BSEND_OVERHEAD = 512 )
-
- INTEGER MPI_SOURCE, MPI_TAG, MPI_ERROR
- PARAMETER(MPI_SOURCE=2, MPI_TAG=3, MPI_ERROR=4)
- INTEGER MPI_STATUS_SIZE
- PARAMETER (MPI_STATUS_SIZE=4)
- INTEGER MPI_MAX_PROCESSOR_NAME, MPI_MAX_ERROR_STRING
- PARAMETER (MPI_MAX_PROCESSOR_NAME=256)
- PARAMETER (MPI_MAX_ERROR_STRING=512)
- INTEGER MPI_MAX_NAME_STRING
- PARAMETER (MPI_MAX_NAME_STRING=63)
- INTEGER MPI_MAX_PORT_NAME
- PARAMETER (MPI_MAX_PORT_NAME=256)
-!
- INTEGER MPI_COMM_NULL
- PARAMETER (MPI_COMM_NULL=0)
-!
- INTEGER MPI_DATATYPE_NULL
- PARAMETER (MPI_DATATYPE_NULL = 0)
-
- INTEGER MPI_ERRHANDLER_NULL
- PARAMETER (MPI_ERRHANDLER_NULL = 0)
-
- INTEGER MPI_GROUP_NULL
- PARAMETER (MPI_GROUP_NULL = 0)
-
- INTEGER MPI_KEYVAL_INVALID
- PARAMETER (MPI_KEYVAL_INVALID = 0)
-
- INTEGER MPI_REQUEST_NULL
- PARAMETER (MPI_REQUEST_NULL = 0)
-!
- INTEGER MPI_IDENT, MPI_CONGRUENT, MPI_SIMILAR, MPI_UNEQUAL
- PARAMETER (MPI_IDENT=0, MPI_CONGRUENT=1, MPI_SIMILAR=2)
- PARAMETER (MPI_UNEQUAL=3)
-!
-! MPI_BOTTOM needs to be a known address; here we put it at the
-! beginning of the common block. The point-to-point and collective
-! routines know about MPI_BOTTOM, but MPI_TYPE_STRUCT as yet does not.
-!
-! MPI_STATUS_IGNORE and MPI_STATUSES_IGNORE are similar objects
-! Until the underlying MPI library implements the C version of these
-! (a null pointer), these are declared as arrays of MPI_STATUS_SIZE
-!
-! The types MPI_INTEGER1,2,4 and MPI_REAL4,8 are OPTIONAL.
-! Their values are zero if they are not available. Note that
-! using these reduces the portability of code (though may enhance
-! portability between Crays and other systems)
-!
- INTEGER MPI_TAG_UB, MPI_HOST, MPI_IO
- INTEGER MPI_BOTTOM
- INTEGER MPI_STATUS_IGNORE(MPI_STATUS_SIZE)
- INTEGER MPI_STATUSES_IGNORE(MPI_STATUS_SIZE)
- INTEGER MPI_INTEGER, MPI_REAL, MPI_DOUBLE_PRECISION
- INTEGER MPI_COMPLEX, MPI_DOUBLE_COMPLEX,MPI_LOGICAL
- INTEGER MPI_CHARACTER, MPI_BYTE, MPI_2INTEGER, MPI_2REAL
- INTEGER MPI_2DOUBLE_PRECISION, MPI_2COMPLEX, MPI_2DOUBLE_COMPLEX
- INTEGER MPI_UB, MPI_LB
- INTEGER MPI_PACKED, MPI_WTIME_IS_GLOBAL
- INTEGER MPI_COMM_WORLD, MPI_COMM_SELF, MPI_GROUP_EMPTY
- INTEGER MPI_SUM, MPI_MAX, MPI_MIN, MPI_PROD, MPI_LAND, MPI_BAND
- INTEGER MPI_LOR, MPI_BOR, MPI_LXOR, MPI_BXOR, MPI_MINLOC
- INTEGER MPI_MAXLOC
- INTEGER MPI_OP_NULL
- INTEGER MPI_ERRORS_ARE_FATAL, MPI_ERRORS_RETURN
-!
- PARAMETER (MPI_ERRORS_ARE_FATAL=119)
- PARAMETER (MPI_ERRORS_RETURN=120)
-!
- PARAMETER (MPI_COMPLEX=23,MPI_DOUBLE_COMPLEX=24,MPI_LOGICAL=25)
- PARAMETER (MPI_REAL=26,MPI_DOUBLE_PRECISION=27,MPI_INTEGER=28)
- PARAMETER (MPI_2INTEGER=29,MPI_2COMPLEX=30,MPI_2DOUBLE_COMPLEX=31)
- PARAMETER (MPI_2REAL=32,MPI_2DOUBLE_PRECISION=33,MPI_CHARACTER=1)
- PARAMETER (MPI_BYTE=3,MPI_UB=16,MPI_LB=15,MPI_PACKED=14)
-
- INTEGER MPI_ORDER_C, MPI_ORDER_FORTRAN
- PARAMETER (MPI_ORDER_C=56, MPI_ORDER_FORTRAN=57)
- INTEGER MPI_DISTRIBUTE_BLOCK, MPI_DISTRIBUTE_CYCLIC
- INTEGER MPI_DISTRIBUTE_NONE, MPI_DISTRIBUTE_DFLT_DARG
- PARAMETER (MPI_DISTRIBUTE_BLOCK=121, MPI_DISTRIBUTE_CYCLIC=122)
- PARAMETER (MPI_DISTRIBUTE_NONE=123)
- PARAMETER (MPI_DISTRIBUTE_DFLT_DARG=-49767)
- INTEGER MPI_MAX_INFO_KEY, MPI_MAX_INFO_VAL
- PARAMETER (MPI_MAX_INFO_KEY=255, MPI_MAX_INFO_VAL=1024)
- INTEGER MPI_INFO_NULL
- PARAMETER (MPI_INFO_NULL=0)
-
-!
-! Optional Fortran Types. Configure attempts to determine these.
-!
- INTEGER MPI_INTEGER1, MPI_INTEGER2, MPI_INTEGER4, MPI_INTEGER8
- INTEGER MPI_INTEGER16
- INTEGER MPI_REAL4, MPI_REAL8, MPI_REAL16
- INTEGER MPI_COMPLEX8, MPI_COMPLEX16, MPI_COMPLEX32
- PARAMETER (MPI_INTEGER1=1,MPI_INTEGER2=4)
- PARAMETER (MPI_INTEGER4=6)
- PARAMETER (MPI_INTEGER8=8)
- PARAMETER (MPI_INTEGER16=0)
- PARAMETER (MPI_REAL4=10)
- PARAMETER (MPI_REAL8=11)
- PARAMETER (MPI_REAL16=12)
- PARAMETER (MPI_COMPLEX8=23)
- PARAMETER (MPI_COMPLEX16=24)
- PARAMETER (MPI_COMPLEX32=0)
-!
-! This is now handled with either the "pointer" extension or this same
-! code, appended at the end.
-! COMMON /MPIPRIV/ MPI_BOTTOM,MPI_STATUS_IGNORE,MPI_STATUSES_IGNORE
-!C
-!C Without this save, some Fortran implementations may make the common
-!C dynamic!
-!C
-!C For a Fortran90 module, we might replace /MPIPRIV/ with a simple
-!C SAVE MPI_BOTTOM
-!C
-! SAVE /MPIPRIV/
-!
- PARAMETER (MPI_MAX=100,MPI_MIN=101,MPI_SUM=102,MPI_PROD=103)
- PARAMETER (MPI_LAND=104,MPI_BAND=105,MPI_LOR=106,MPI_BOR=107)
- PARAMETER (MPI_LXOR=108,MPI_BXOR=109,MPI_MINLOC=110)
- PARAMETER (MPI_MAXLOC=111, MPI_OP_NULL=0)
-!
- PARAMETER (MPI_GROUP_EMPTY=90,MPI_COMM_WORLD=91,MPI_COMM_SELF=92)
- PARAMETER (MPI_TAG_UB=80,MPI_HOST=82,MPI_IO=84)
- PARAMETER (MPI_WTIME_IS_GLOBAL=86)
-!
- INTEGER MPI_ANY_SOURCE
- PARAMETER (MPI_ANY_SOURCE = (-2))
- INTEGER MPI_ANY_TAG
- PARAMETER (MPI_ANY_TAG = (-1))
-!
- INTEGER MPI_VERSION, MPI_SUBVERSION
- PARAMETER (MPI_VERSION = 1, MPI_SUBVERSION = 2)
-!
-! There are additional MPI-2 constants
- INTEGER MPI_ADDRESS_KIND, MPI_OFFSET_KIND
- PARAMETER (MPI_ADDRESS_KIND=8)
- PARAMETER (MPI_OFFSET_KIND=8)
-!
-! All other MPI routines are subroutines
-! This may cause some Fortran compilers to complain about defined and
-! not used. Such compilers should be improved.
-!
-! Some Fortran compilers will not link programs that contain
-! external statements to routines that are not provided, even if
-! the routine is never called. Remove PMPI_WTIME and PMPI_WTICK
-! if you have trouble with them.
-!
- DOUBLE PRECISION MPI_WTIME, MPI_WTICK,PMPI_WTIME,PMPI_WTICK
- EXTERNAL MPI_WTIME, MPI_WTICK,PMPI_WTIME,PMPI_WTICK
-!
-! The attribute copy/delete subroutines are symbols that can be passed
-! to MPI routines
-!
- EXTERNAL MPI_NULL_COPY_FN, MPI_NULL_DELETE_FN, MPI_DUP_FN
- COMMON /MPIPRIV/ MPI_BOTTOM,MPI_STATUS_IGNORE,MPI_STATUSES_IGNORE
-!
-! Without this save, some Fortran implementations may make the common
-! dynamic!
-!
-! For a Fortran90 module, we might replace /MPIPRIV/ with a simple
-! SAVE MPI_BOTTOM
-!
- SAVE /MPIPRIV/
-!
-! $Id: mpiof.h.in,v 1.3 1999/08/06 18:33:09 thakur Exp $
-!
-! Copyright (C) 1997 University of Chicago.
-! See COPYRIGHT notice in top-level directory.
-!
-!
-! user include file for Fortran MPI-IO programs
-!
- INTEGER MPI_MODE_RDONLY, MPI_MODE_RDWR, MPI_MODE_WRONLY
- INTEGER MPI_MODE_DELETE_ON_CLOSE, MPI_MODE_UNIQUE_OPEN
- INTEGER MPI_MODE_CREATE, MPI_MODE_EXCL
- INTEGER MPI_MODE_APPEND, MPI_MODE_SEQUENTIAL
- PARAMETER (MPI_MODE_RDONLY=2, MPI_MODE_RDWR=8, MPI_MODE_WRONLY=4)
- PARAMETER (MPI_MODE_CREATE=1, MPI_MODE_DELETE_ON_CLOSE=16)
- PARAMETER (MPI_MODE_UNIQUE_OPEN=32, MPI_MODE_EXCL=64)
- PARAMETER (MPI_MODE_APPEND=128, MPI_MODE_SEQUENTIAL=256)
-!
- INTEGER MPI_FILE_NULL
- PARAMETER (MPI_FILE_NULL=0)
-!
- INTEGER MPI_MAX_DATAREP_STRING
- PARAMETER (MPI_MAX_DATAREP_STRING=128)
-!
- INTEGER MPI_SEEK_SET, MPI_SEEK_CUR, MPI_SEEK_END
- PARAMETER (MPI_SEEK_SET=600, MPI_SEEK_CUR=602, MPI_SEEK_END=604)
-!
- INTEGER MPIO_REQUEST_NULL
- PARAMETER (MPIO_REQUEST_NULL=0)
-!
-!
-!
-
-
-
-
-
-
-
-!
-!
-!
-!
-!
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/precision.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/precision.h 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/precision.h 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,18 +1,26 @@
!=====================================================================
!
-! S p e c f e m 3 D B a s i n V e r s i o n 1 . 4
-! --------------------------------------------------
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
!
! Dimitri Komatitsch and Jeroen Tromp
! Seismological Laboratory - California Institute of Technology
! (c) California Institute of Technology July 2005
!
-! A signed non-commercial agreement is required to use this program.
-! Please check http://www.gps.caltech.edu/research/jtromp for details.
-! Free for non-commercial academic research ONLY.
-! This program is distributed WITHOUT ANY WARRANTY whatsoever.
-! Do not redistribute this program without written permission.
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
!=====================================================================
! precision.h. Generated from precision.h.in by configure.
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/sum_kernels
===================================================================
(Binary files differ)
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/sum_kernels.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/sum_kernels.f90 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/sum_kernels.f90 2010-02-09 23:44:05 UTC (rev 16243)
@@ -8,9 +8,9 @@
! ======================================================
- integer, parameter :: NSPEC=NSPEC_AB
+ integer, parameter :: NSPEC=NSPEC_AB_VAL
- character(len=150) :: kernel_file_list, kernel_list(1000), sline, k_file, kernel_name
+ character(len=150) :: kernel_file_list, kernel_list(1000), sline, k_file, kernel_name, ftagfile, ftag
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: kernel, total_kernel
integer :: iker, nker, myrank, sizeprocs, ier
integer :: i, j, k,ispec, iglob, ishell, n, it, j1, ib, npts_sem, ios
@@ -21,8 +21,10 @@
call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
- kernel_file_list='kernels_run'
+ kernel_file_list='kernels_run_sub'
+ ftagfile = 'ftags'
+ ! read in the list of event IDs
nker=0
open(unit = 20, file = trim(kernel_file_list), status = 'old',iostat = ios)
if (ios /= 0) then
@@ -37,8 +39,22 @@
enddo
close(20)
+ ! read in the name of the kernel
+ open(unit = 20, file = trim(ftagfile), status = 'old',iostat = ios)
+ if (ios /= 0) then
+ print *,'Error opening ',trim(ftagfile)
+ stop
+ endif
+ read(20,'(a)',iostat=ios) ftag
+ if (ios /= 0) then
+ print *,'Error reading ',trim(ftagfile)
+ stop
+ endif
+ close(20)
+
!kernel_name = 'mu_kernel'
- kernel_name = 'mu_kernel_smooth'
+ !kernel_name = 'mu_kernel_smooth'
+ kernel_name = ftag
total_kernel=0.
do iker = 1, nker
@@ -59,28 +75,28 @@
write(12) total_kernel(:,:,:,1:nspec)
close(12)
- !kernel_name = 'kappa_kernel'
- kernel_name = 'kappa_kernel_smooth'
+!!$ !kernel_name = 'kappa_kernel'
+!!$ kernel_name = 'kappa_kernel_smooth'
+!!$
+!!$ total_kernel=0.
+!!$ do iker = 1, nker
+!!$ if(myrank==1) write(*,*) 'reading in event kernel for kappa: ', iker, ' out of ', nker
+!!$ write(k_file,'(a,i6.6,a)') 'INPUT_KERNELS/'//trim(kernel_list(iker))//'/proc',myrank,'_'//trim(kernel_name)//'.bin'
+!!$
+!!$ open(12,file=trim(k_file),status='old',form='unformatted')
+!!$ read(12) kernel(:,:,:,1:nspec)
+!!$ close(12)
+!!$
+!!$ !total_kernel(:,:,:,1:nspec) = total_kernel(:,:,:,1:nspec) + abs( kernel(:,:,:,1:nspec) )
+!!$ total_kernel(:,:,:,1:nspec) = total_kernel(:,:,:,1:nspec) + kernel(:,:,:,1:nspec)
+!!$
+!!$ enddo
+!!$ if(myrank==1) write(*,*) 'writing out summed kernel for kappa'
+!!$ write(k_file,'(a,i6.6,a)') 'OUTPUT_SUM/proc',myrank,'_'//trim(kernel_name)//'.bin'
+!!$ open(12,file=trim(k_file),form='unformatted')
+!!$ write(12) total_kernel(:,:,:,1:nspec)
+!!$ close(12)
- total_kernel=0.
- do iker = 1, nker
- if(myrank==1) write(*,*) 'reading in event kernel for kappa: ', iker, ' out of ', nker
- write(k_file,'(a,i6.6,a)') 'INPUT_KERNELS/'//trim(kernel_list(iker))//'/proc',myrank,'_'//trim(kernel_name)//'.bin'
-
- open(12,file=trim(k_file),status='old',form='unformatted')
- read(12) kernel(:,:,:,1:nspec)
- close(12)
-
- !total_kernel(:,:,:,1:nspec) = total_kernel(:,:,:,1:nspec) + abs( kernel(:,:,:,1:nspec) )
- total_kernel(:,:,:,1:nspec) = total_kernel(:,:,:,1:nspec) + kernel(:,:,:,1:nspec)
-
- enddo
- if(myrank==1) write(*,*) 'writing out summed kernel for kappa'
- write(k_file,'(a,i6.6,a)') 'OUTPUT_SUM/proc',myrank,'_'//trim(kernel_name)//'.bin'
- open(12,file=trim(k_file),form='unformatted')
- write(12) total_kernel(:,:,:,1:nspec)
- close(12)
-
if(myrank==1) write(*,*) 'done writing all kernels, now finishing...'
! stop all the MPI processes, and exit
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/values_from_mesher.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/values_from_mesher.h 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/values_from_mesher.h 2010-02-09 23:44:05 UTC (rev 16243)
@@ -10,6 +10,9 @@
! number of ES nodes = 21.00000
! percentage of total 640 ES nodes = 3.281250 %
! total memory available on these ES nodes (Gb) = 336.0000
+ integer, parameter :: NPROC_VAL = 168
+ integer, parameter :: NPROC_XI_VAL = 14
+ integer, parameter :: NPROC_ETA_VAL = 12
!
! max points per processor = max vector length = 163941
! min vector length = 25
@@ -42,13 +45,29 @@
! average distance between points along Y in m = 436.8476
!
- integer, parameter :: NSPEC_AB = 2412
- integer, parameter :: NGLOB_AB = 163941
+ integer, parameter :: NSPEC_AB_VAL = 2412
+ integer, parameter :: NGLOB_AB_VAL = 163941
!
- ! number of time steps = 27300
+ ! number of time steps = 18200
!
- integer, parameter :: NSPEC_ATTENUATION = 1
+ integer, parameter :: NSPEC_ATTENUATION = 1
logical, parameter :: ATTENUATION_VAL = .false.
+
+ integer, parameter :: NSPEC_ANISO = 1
logical, parameter :: ANISOTROPY_VAL = .false.
+ integer, parameter :: NSPEC_ATT_AND_KERNEL = 1
+ integer, parameter :: NSPEC_ADJOINT = 2412
+ integer, parameter :: NGLOB_ADJOINT = 163941
+
+ integer, parameter :: NSPEC2DMAX_XMIN_XMAX_VAL = 150
+ integer, parameter :: NSPEC2DMAX_YMIN_YMAX_VAL = 150
+ integer, parameter :: NSPEC2D_BOTTOM_VAL = 36
+ integer, parameter :: NSPEC2D_TOP_VAL = 576
+ integer, parameter :: NPOIN2DMAX_XMIN_XMAX_VAL = 3751
+ integer, parameter :: NPOIN2DMAX_YMIN_YMAX_VAL = 3751
+ integer, parameter :: NPOIN2DMAX_XY_VAL = 3751
+
+ integer, parameter :: NSPEC2D_MOHO_BOUN = 1
+ integer, parameter :: NSPEC_BOUN = 1
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/README (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/README 2010-02-09 23:44:05 UTC (rev 16243)
@@ -0,0 +1,20 @@
+Carl Tape
+09-Nov-2009
+
+This version of combine_vol_data.f90 is modified by Qinya Liu to allow the topology files to be read in from a base directory. Probably this version should update the one residing in SPECFEM3D-1.4.3.
+
+To avoid confusion, the file has been renamed combine_vol_data_mod.f90
+
+Compiling with compile.bash generated the following warning messages:
+
+write_c_binary.c: In function 'open_file_':
+write_c_binary.c:45: warning: incompatible implicit declaration of built-in function 'exit'
+/n/sw/intel/ifort/11.1/046/lib/intel64/libimf.so: warning: warning: feupdateenv is not implemented and will always fail
+
+These do not appear to cause any problems.
+
+Note that the following files will differ for each user:
+ values_from_mesher.h
+ constants.h
+
+------------------------------------------
\ No newline at end of file
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/combine_vol_data.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/combine_vol_data.f90 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/combine_vol_data.f90 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,427 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D B a s i n V e r s i o n 1 . 4
-! --------------------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! A signed non-commercial agreement is required to use this program.
-! Please check http://www.gps.caltech.edu/research/jtromp for details.
-! Free for non-commercial academic research ONLY.
-! This program is distributed WITHOUT ANY WARRANTY whatsoever.
-! Do not redistribute this program without written permission.
-!
-!=====================================================================
-
- program combine_paraview_data
-
-! puts the output of SPECFEM3D in ParaView format.
-! see http://www.paraview.org for details
-
-! combines the database files on several slices.
-! the local database file needs to have been collected onto the frontend (copy_local_database.pl)
-
- implicit none
-
- include 'constants.h'
- include 'values_from_mesher.h'
-
- integer i,j,k,ispec, ios, it
- integer iproc, proc1, proc2, num_node, node_list(300), nspec, nglob
- integer np, ne, npp, nee, npoint, nelement, njunk, njunk2, n1, n2, n3, n4, n5, n6, n7, n8
- integer ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB)
- integer numpoin, iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8, iglob
- logical mask_ibool(NGLOB_AB)
- real(kind=CUSTOM_REAL) data(NGLLX,NGLLY,NGLLZ,NSPEC_AB)
- real(kind=CUSTOM_REAL),dimension(NGLOB_AB) :: xstore, ystore, zstore
- real x, y, z, dat(NGLLX,NGLLY,NGLLZ,NSPEC_AB)
- character(len=150) :: sline, arg(7), filename, in_file_dir, in_topo_dir, outdir, prname_topo,prname_file
- character(len=150) :: mesh_file, local_point_file, local_element_file, local_file, local_data_file, local_ibool_file
- integer :: num_ibool(NGLOB_AB)
- logical :: HIGH_RESOLUTION_MESH
- integer :: ires
-
- print *
- print *,'Recombining ParaView data for slices'
- print *
-
- do i = 1, 6
- call getarg(i,arg(i))
- if (trim(arg(i)) == '') then
- print *, 'Usage: xcombine_data slice_list filename input_topo_dir input_kernel_dir output_dir high/low-resolution'
- print *, ' possible filenames are '
- print *, ' rho_vp, rho_vs, kappastore, mustore etc'
- print *, ' that are stored in the local directory as real(kind=CUSTOM_REAL) filename(NGLLX,NGLLY,NGLLZ,nspec) '
- print *, ' in filename.bin'
- print *, ' topology files have been collected in input_topo_dir'
- print *, ' files have been collected in input_file_dir, output mesh file goes to output_dir '
- print *, ' give 0 for low resolution and 1 for high resolution'
- stop ' Reenter command line options'
- endif
- enddo
-
-! get slice list
- num_node = 0
- open(unit = 20, file = trim(arg(1)), status = 'unknown',iostat = ios)
- if (ios /= 0) then
- print *,'Error opening ',trim(arg(1))
- stop
- endif
- do while ( 1 == 1)
- read(20,'(a)',iostat=ios) sline
- if (ios /= 0) exit
- read(sline,*,iostat=ios) njunk
- if (ios /= 0) exit
- num_node = num_node + 1
- node_list(num_node) = njunk
- enddo
- close(20)
- filename = arg(2)
- in_topo_dir=arg(3)
- in_file_dir= arg(4)
- outdir = arg(5)
- read(arg(6),*) ires
-
- if (ires == 0) then
- HIGH_RESOLUTION_MESH = .false.
- else
- HIGH_RESOLUTION_MESH = .true.
- endif
-
- print *, 'Slice list: '
- print *, node_list(1:num_node)
-
- ! open paraview output mesh file
- mesh_file = trim(outdir) // '/' // trim(filename)//'.mesh'
- call open_file(trim(mesh_file)//char(0))
-
- nspec = NSPEC_AB
- nglob = NGLOB_AB
-
- np = 0
-
- ! write point and scalar information
-
- do it = 1, num_node
-
- iproc = node_list(it)
-
- print *, ' '
- print *, 'Reading slice ', iproc
- write(prname_topo,'(a,i6.6,a)') trim(in_topo_dir)//'/proc',iproc,'_'
- write(prname_file,'(a,i6.6,a)') trim(in_file_dir)//'/proc',iproc,'_'
-
- ! data file
- local_data_file = trim(prname_file) // trim(filename) // '.bin'
- open(unit = 27,file = trim(local_data_file),status='old',action='read', iostat = ios,form ='unformatted')
- if (ios /= 0) then
- print *,'Error opening ',trim(local_data_file)
- stop
- endif
- read(27) data
- close(27)
- print *, trim(local_data_file)
-
- dat = data
-
- ! ibool file
- local_ibool_file = trim(prname_topo) // 'ibool' // '.bin'
- open(unit = 28,file = trim(local_ibool_file),status='old',action='read', iostat = ios, form='unformatted')
- if (ios /= 0) then
- print *,'Error opening ',trim(local_data_file)
- stop
- endif
- read(28) ibool
- close(28)
- print *, trim(local_ibool_file)
-
- mask_ibool(:) = .false.
- numpoin = 0
-
- if (.not. HIGH_RESOLUTION_MESH) then
-
- local_point_file = trim(prname_topo) // 'AVS_DXpoints.txt'
- open(unit = 25, file = trim(local_point_file), status = 'old', iostat = ios)
- if (ios /= 0) then
- print *,'Error opening ',trim(local_point_file)
- stop
- endif
- read(25,*) npoint
-
- if (it == 1) then
- npp = npoint * num_node
- call write_integer(npp)
- endif
-
- do ispec=1,nspec
- iglob1=ibool(1,1,1,ispec)
- iglob2=ibool(NGLLX,1,1,ispec)
- iglob3=ibool(NGLLX,NGLLY,1,ispec)
- iglob4=ibool(1,NGLLY,1,ispec)
- iglob5=ibool(1,1,NGLLZ,ispec)
- iglob6=ibool(NGLLX,1,NGLLZ,ispec)
- iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
- iglob8=ibool(1,NGLLY,NGLLZ,ispec)
-
- if(.not. mask_ibool(iglob1)) then
- numpoin = numpoin + 1
- read(25,*) njunk, x, y, z
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(1,1,1,ispec))
- mask_ibool(iglob1) = .true.
- endif
- if(.not. mask_ibool(iglob2)) then
- numpoin = numpoin + 1
- read(25,*) njunk, x, y, z
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(NGLLX,1,1,ispec))
- mask_ibool(iglob2) = .true.
- endif
- if(.not. mask_ibool(iglob3)) then
- numpoin = numpoin + 1
- read(25,*) njunk, x, y, z
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(NGLLX,NGLLY,1,ispec))
- mask_ibool(iglob3) = .true.
- endif
- if(.not. mask_ibool(iglob4)) then
- numpoin = numpoin + 1
- read(25,*) njunk, x, y, z
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(1,NGLLY,1,ispec))
- mask_ibool(iglob4) = .true.
- endif
- if(.not. mask_ibool(iglob5)) then
- numpoin = numpoin + 1
- read(25,*) njunk, x, y, z
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(1,1,NGLLZ,ispec))
- mask_ibool(iglob5) = .true.
- endif
- if(.not. mask_ibool(iglob6)) then
- numpoin = numpoin + 1
- read(25,*) njunk, x, y, z
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(NGLLX,1,NGLLZ,ispec))
- mask_ibool(iglob6) = .true.
- endif
- if(.not. mask_ibool(iglob7)) then
- numpoin = numpoin + 1
- read(25,*) njunk, x, y, z
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(NGLLX,NGLLY,NGLLZ,ispec))
- mask_ibool(iglob7) = .true.
- endif
- if(.not. mask_ibool(iglob8)) then
- numpoin = numpoin + 1
- read(25,*) njunk, x, y, z
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(1,NGLLY,NGLLZ,ispec))
- mask_ibool(iglob8) = .true.
- endif
- enddo ! ispec
- close(25)
-
- else ! high resolution
-
- if (it == 1) then
- npp = nglob * num_node
- npoint = nglob
- call write_integer(npp)
- endif
-
- local_file = trim(prname_topo)//'x.bin'
- open(unit = 27,file = trim(prname_topo)//'x.bin',status='old',action='read', iostat = ios,form ='unformatted')
- if (ios /= 0) then
- print *,'Error opening ',trim(local_file)
- stop
- endif
- read(27) xstore
- close(27)
- local_file = trim(prname_topo)//'y.bin'
- open(unit = 27,file = trim(prname_topo)//'y.bin',status='old',action='read', iostat = ios,form ='unformatted')
- if (ios /= 0) then
- print *,'Error opening ',trim(local_file)
- stop
- endif
- read(27) ystore
- close(27)
- local_file = trim(prname_topo)//'z.bin'
- open(unit = 27,file = trim(prname_topo)//'z.bin',status='old',action='read', iostat = ios,form ='unformatted')
- if (ios /= 0) then
- print *,'Error opening ',trim(local_file)
- stop
- endif
- read(27) zstore
- close(27)
-
- do ispec=1,nspec
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob = ibool(i,j,k,ispec)
- if(.not. mask_ibool(iglob)) then
- numpoin = numpoin + 1
- x = xstore(iglob)
- y = ystore(iglob)
- z = zstore(iglob)
- call write_real(x)
- call write_real(y)
- call write_real(z)
- call write_real(dat(i,j,k,ispec))
- mask_ibool(iglob) = .true.
- endif
- enddo ! i
- enddo ! j
- enddo ! k
- enddo !ispec
- endif
-
- if (numpoin /= npoint) stop 'Error: number of points are not consistent'
- np = np + npoint
-
- enddo ! all slices for points
-
- if (np /= npp) stop 'Error: Number of total points are not consistent'
- print *, 'Total number of points: ', np
- print *, ' '
-
-
- ne = 0
-! write element information
- do it = 1, num_node
-
- iproc = node_list(it)
-
- print *, 'Reading slice ', iproc
- write(prname_topo,'(a,i6.6,a)') trim(in_topo_dir)//'/proc',iproc,'_'
- write(prname_file,'(a,i6.6,a)') trim(in_file_dir)//'/proc',iproc,'_'
-
- np = npoint * (it-1)
-
- if (.not. HIGH_RESOLUTION_MESH) then
-
- local_element_file = trim(prname_topo) // 'AVS_DXelements.txt'
- open(unit = 26, file = trim(local_element_file), status = 'old', iostat = ios)
- if (ios /= 0) then
- print *,'Error opening ',trim(local_element_file)
- stop
- endif
- print *, trim(local_element_file)
-
- read(26, *) nelement
- if (it == 1) then
- nee = nelement * num_node
- call write_integer(nee)
- end if
-
- do i = 1, nelement
- read(26,*) njunk, njunk2, n1, n2, n3, n4, n5, n6, n7, n8
- n1 = n1+np-1
- n2 = n2+np-1
- n3 = n3+np-1
- n4 = n4+np-1
- n5 = n5+np-1
- n6 = n6+np-1
- n7 = n7+np-1
- n8 = n8+np-1
- call write_integer(n1)
- call write_integer(n2)
- call write_integer(n3)
- call write_integer(n4)
- call write_integer(n5)
- call write_integer(n6)
- call write_integer(n7)
- call write_integer(n8)
- enddo
- close(26)
-
- else ! high resolution mesh
-
- if (it == 1) then
- nelement = nspec * (NGLLX-1) * (NGLLY-1) * (NGLLZ-1)
- nee = nelement * num_node
- call write_integer(nee)
- endif
-
- numpoin = 0
- mask_ibool = .false.
- do ispec=1,nspec
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- iglob = ibool(i,j,k,ispec)
- if(.not. mask_ibool(iglob)) then
- numpoin = numpoin + 1
- num_ibool(iglob) = numpoin
- mask_ibool(iglob) = .true.
- endif
- enddo ! i
- enddo ! j
- enddo ! k
- enddo !ispec
-
- do ispec = 1, nspec
- do k = 1, NGLLZ-1
- do j = 1, NGLLY-1
- do i = 1, NGLLX-1
- iglob1 = ibool(i,j,k,ispec)
- iglob2 = ibool(i+1,j,k,ispec)
- iglob3 = ibool(i+1,j+1,k,ispec)
- iglob4 = ibool(i,j+1,k,ispec)
- iglob5 = ibool(i,j,k+1,ispec)
- iglob6 = ibool(i+1,j,k+1,ispec)
- iglob7 = ibool(i+1,j+1,k+1,ispec)
- iglob8 = ibool(i,j+1,k+1,ispec)
- n1 = num_ibool(iglob1)+np-1
- n2 = num_ibool(iglob2)+np-1
- n3 = num_ibool(iglob3)+np-1
- n4 = num_ibool(iglob4)+np-1
- n5 = num_ibool(iglob5)+np-1
- n6 = num_ibool(iglob6)+np-1
- n7 = num_ibool(iglob7)+np-1
- n8 = num_ibool(iglob8)+np-1
- call write_integer(n1)
- call write_integer(n2)
- call write_integer(n3)
- call write_integer(n4)
- call write_integer(n5)
- call write_integer(n6)
- call write_integer(n7)
- call write_integer(n8)
- enddo
- enddo
- enddo
- enddo
-
- endif
- ne = ne + nelement
-
- enddo ! num_node
- if (ne /= nee) stop 'Number of total elements are not consistent'
- print *, 'Total number of elements: ', ne
-
- call close_file()
-
- print *, 'Done writing '//trim(mesh_file)
-
- end program combine_paraview_data
-
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/combine_vol_data_mod.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/combine_vol_data_mod.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/combine_vol_data_mod.f90 2010-02-09 23:44:05 UTC (rev 16243)
@@ -0,0 +1,434 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ program combine_paraview_data
+
+! puts the output of SPECFEM3D in ParaView format.
+! see http://www.paraview.org for details
+
+! combines the database files on several slices.
+! the local database file needs to have been collected onto the frontend (copy_local_database.pl)
+
+ implicit none
+
+ include 'constants.h'
+ include 'values_from_mesher.h'
+
+ integer i,j,k,ispec, ios, it
+ integer iproc, proc1, proc2, num_node, node_list(300), nspec, nglob
+ integer np, ne, npp, nee, npoint, nelement, njunk, njunk2, n1, n2, n3, n4, n5, n6, n7, n8
+ integer ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL)
+ integer numpoin, iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8, iglob
+ logical mask_ibool(NGLOB_AB_VAL)
+ real(kind=CUSTOM_REAL) data(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL)
+ real(kind=CUSTOM_REAL),dimension(NGLOB_AB_VAL) :: xstore, ystore, zstore
+ real x, y, z, dat(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL)
+ character(len=150) :: sline, arg(7), filename, in_file_dir, in_topo_dir, outdir, prname_topo,prname_file
+ character(len=150) :: mesh_file, local_point_file, local_element_file, local_file, local_data_file, local_ibool_file
+ integer :: num_ibool(NGLOB_AB_VAL)
+ logical :: HIGH_RESOLUTION_MESH
+ integer :: ires
+
+ print *
+ print *,'Recombining ParaView data for slices'
+ print *
+
+ do i = 1, 6
+ call getarg(i,arg(i))
+ if (trim(arg(i)) == '') then
+ print *, 'Usage: xcombine_data slice_list filename input_topo_dir input_kernel_dir output_dir high/low-resolution'
+ print *, ' possible filenames are '
+ print *, ' rho_vp, rho_vs, kappastore, mustore etc'
+ print *, ' that are stored in the local directory as real(kind=CUSTOM_REAL) filename(NGLLX,NGLLY,NGLLZ,nspec) '
+ print *, ' in filename.bin'
+ print *, ' topology files have been collected in input_topo_dir'
+ print *, ' files have been collected in input_file_dir, output mesh file goes to output_dir '
+ print *, ' give 0 for low resolution and 1 for high resolution'
+ stop ' Reenter command line options'
+ endif
+ enddo
+
+! get slice list
+ num_node = 0
+ open(unit = 20, file = trim(arg(1)), status = 'unknown',iostat = ios)
+ if (ios /= 0) then
+ print *,'Error opening ',trim(arg(1))
+ stop
+ endif
+ do while ( 1 == 1)
+ read(20,'(a)',iostat=ios) sline
+ if (ios /= 0) exit
+ read(sline,*,iostat=ios) njunk
+ if (ios /= 0) exit
+ num_node = num_node + 1
+ node_list(num_node) = njunk
+ enddo
+ close(20)
+ filename = arg(2)
+ in_topo_dir=arg(3)
+ in_file_dir= arg(4)
+ outdir = arg(5)
+ read(arg(6),*) ires
+
+ if (ires == 0) then
+ HIGH_RESOLUTION_MESH = .false.
+ else
+ HIGH_RESOLUTION_MESH = .true.
+ endif
+
+ print *, 'Slice list: '
+ print *, node_list(1:num_node)
+
+ ! open paraview output mesh file
+ mesh_file = trim(outdir) // '/' // trim(filename)//'.mesh'
+ call open_file(trim(mesh_file)//char(0))
+
+ nspec = NSPEC_AB_VAL
+ nglob = NGLOB_AB_VAL
+
+ np = 0
+
+ ! write point and scalar information
+
+ do it = 1, num_node
+
+ iproc = node_list(it)
+
+ print *, ' '
+ print *, 'Reading slice ', iproc
+ write(prname_topo,'(a,i6.6,a)') trim(in_topo_dir)//'/proc',iproc,'_'
+ write(prname_file,'(a,i6.6,a)') trim(in_file_dir)//'/proc',iproc,'_'
+
+ ! data file
+ local_data_file = trim(prname_file) // trim(filename) // '.bin'
+ open(unit = 27,file = trim(local_data_file),status='old',action='read', iostat = ios,form ='unformatted')
+ if (ios /= 0) then
+ print *,'Error opening ',trim(local_data_file)
+ stop
+ endif
+ read(27) data
+ close(27)
+ print *, trim(local_data_file)
+
+ dat = data
+
+ ! ibool file
+ local_ibool_file = trim(prname_topo) // 'ibool' // '.bin'
+ open(unit = 28,file = trim(local_ibool_file),status='old',action='read', iostat = ios, form='unformatted')
+ if (ios /= 0) then
+ print *,'Error opening ',trim(local_data_file)
+ stop
+ endif
+ read(28) ibool
+ close(28)
+ print *, trim(local_ibool_file)
+
+ mask_ibool(:) = .false.
+ numpoin = 0
+
+ if (.not. HIGH_RESOLUTION_MESH) then
+
+ local_point_file = trim(prname_topo) // 'AVS_DXpoints.txt'
+ open(unit = 25, file = trim(local_point_file), status = 'old', iostat = ios)
+ if (ios /= 0) then
+ print *,'Error opening ',trim(local_point_file)
+ stop
+ endif
+ read(25,*) npoint
+
+ if (it == 1) then
+ npp = npoint * num_node
+ call write_integer(npp)
+ endif
+
+ do ispec=1,nspec
+ iglob1=ibool(1,1,1,ispec)
+ iglob2=ibool(NGLLX,1,1,ispec)
+ iglob3=ibool(NGLLX,NGLLY,1,ispec)
+ iglob4=ibool(1,NGLLY,1,ispec)
+ iglob5=ibool(1,1,NGLLZ,ispec)
+ iglob6=ibool(NGLLX,1,NGLLZ,ispec)
+ iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglob8=ibool(1,NGLLY,NGLLZ,ispec)
+
+ if(.not. mask_ibool(iglob1)) then
+ numpoin = numpoin + 1
+ read(25,*) njunk, x, y, z
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(1,1,1,ispec))
+ mask_ibool(iglob1) = .true.
+ endif
+ if(.not. mask_ibool(iglob2)) then
+ numpoin = numpoin + 1
+ read(25,*) njunk, x, y, z
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(NGLLX,1,1,ispec))
+ mask_ibool(iglob2) = .true.
+ endif
+ if(.not. mask_ibool(iglob3)) then
+ numpoin = numpoin + 1
+ read(25,*) njunk, x, y, z
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(NGLLX,NGLLY,1,ispec))
+ mask_ibool(iglob3) = .true.
+ endif
+ if(.not. mask_ibool(iglob4)) then
+ numpoin = numpoin + 1
+ read(25,*) njunk, x, y, z
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(1,NGLLY,1,ispec))
+ mask_ibool(iglob4) = .true.
+ endif
+ if(.not. mask_ibool(iglob5)) then
+ numpoin = numpoin + 1
+ read(25,*) njunk, x, y, z
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(1,1,NGLLZ,ispec))
+ mask_ibool(iglob5) = .true.
+ endif
+ if(.not. mask_ibool(iglob6)) then
+ numpoin = numpoin + 1
+ read(25,*) njunk, x, y, z
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(NGLLX,1,NGLLZ,ispec))
+ mask_ibool(iglob6) = .true.
+ endif
+ if(.not. mask_ibool(iglob7)) then
+ numpoin = numpoin + 1
+ read(25,*) njunk, x, y, z
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(NGLLX,NGLLY,NGLLZ,ispec))
+ mask_ibool(iglob7) = .true.
+ endif
+ if(.not. mask_ibool(iglob8)) then
+ numpoin = numpoin + 1
+ read(25,*) njunk, x, y, z
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(1,NGLLY,NGLLZ,ispec))
+ mask_ibool(iglob8) = .true.
+ endif
+ enddo ! ispec
+ close(25)
+
+ else ! high resolution
+
+ if (it == 1) then
+ npp = nglob * num_node
+ npoint = nglob
+ call write_integer(npp)
+ endif
+
+ local_file = trim(prname_topo)//'x.bin'
+ open(unit = 27,file = trim(prname_topo)//'x.bin',status='old',action='read', iostat = ios,form ='unformatted')
+ if (ios /= 0) then
+ print *,'Error opening ',trim(local_file)
+ stop
+ endif
+ read(27) xstore
+ close(27)
+ local_file = trim(prname_topo)//'y.bin'
+ open(unit = 27,file = trim(prname_topo)//'y.bin',status='old',action='read', iostat = ios,form ='unformatted')
+ if (ios /= 0) then
+ print *,'Error opening ',trim(local_file)
+ stop
+ endif
+ read(27) ystore
+ close(27)
+ local_file = trim(prname_topo)//'z.bin'
+ open(unit = 27,file = trim(prname_topo)//'z.bin',status='old',action='read', iostat = ios,form ='unformatted')
+ if (ios /= 0) then
+ print *,'Error opening ',trim(local_file)
+ stop
+ endif
+ read(27) zstore
+ close(27)
+
+ do ispec=1,nspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+ if(.not. mask_ibool(iglob)) then
+ numpoin = numpoin + 1
+ x = xstore(iglob)
+ y = ystore(iglob)
+ z = zstore(iglob)
+ call write_real(x)
+ call write_real(y)
+ call write_real(z)
+ call write_real(dat(i,j,k,ispec))
+ mask_ibool(iglob) = .true.
+ endif
+ enddo ! i
+ enddo ! j
+ enddo ! k
+ enddo !ispec
+ endif
+
+ if (numpoin /= npoint) stop 'Error: number of points are not consistent'
+ np = np + npoint
+
+ enddo ! all slices for points
+
+ if (np /= npp) stop 'Error: Number of total points are not consistent'
+ print *, 'Total number of points: ', np
+ print *, ' '
+
+
+ ne = 0
+! write element information
+ do it = 1, num_node
+
+ iproc = node_list(it)
+
+ print *, 'Reading slice ', iproc
+ write(prname_topo,'(a,i6.6,a)') trim(in_topo_dir)//'/proc',iproc,'_'
+ write(prname_file,'(a,i6.6,a)') trim(in_file_dir)//'/proc',iproc,'_'
+
+ np = npoint * (it-1)
+
+ if (.not. HIGH_RESOLUTION_MESH) then
+
+ local_element_file = trim(prname_topo) // 'AVS_DXelements.txt'
+ open(unit = 26, file = trim(local_element_file), status = 'old', iostat = ios)
+ if (ios /= 0) then
+ print *,'Error opening ',trim(local_element_file)
+ stop
+ endif
+ print *, trim(local_element_file)
+
+ read(26, *) nelement
+ if (it == 1) then
+ nee = nelement * num_node
+ call write_integer(nee)
+ end if
+
+ do i = 1, nelement
+ read(26,*) njunk, njunk2, n1, n2, n3, n4, n5, n6, n7, n8
+ n1 = n1+np-1
+ n2 = n2+np-1
+ n3 = n3+np-1
+ n4 = n4+np-1
+ n5 = n5+np-1
+ n6 = n6+np-1
+ n7 = n7+np-1
+ n8 = n8+np-1
+ call write_integer(n1)
+ call write_integer(n2)
+ call write_integer(n3)
+ call write_integer(n4)
+ call write_integer(n5)
+ call write_integer(n6)
+ call write_integer(n7)
+ call write_integer(n8)
+ enddo
+ close(26)
+
+ else ! high resolution mesh
+
+ if (it == 1) then
+ nelement = nspec * (NGLLX-1) * (NGLLY-1) * (NGLLZ-1)
+ nee = nelement * num_node
+ call write_integer(nee)
+ endif
+
+ numpoin = 0
+ mask_ibool = .false.
+ do ispec=1,nspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+ if(.not. mask_ibool(iglob)) then
+ numpoin = numpoin + 1
+ num_ibool(iglob) = numpoin
+ mask_ibool(iglob) = .true.
+ endif
+ enddo ! i
+ enddo ! j
+ enddo ! k
+ enddo !ispec
+
+ do ispec = 1, nspec
+ do k = 1, NGLLZ-1
+ do j = 1, NGLLY-1
+ do i = 1, NGLLX-1
+ iglob1 = ibool(i,j,k,ispec)
+ iglob2 = ibool(i+1,j,k,ispec)
+ iglob3 = ibool(i+1,j+1,k,ispec)
+ iglob4 = ibool(i,j+1,k,ispec)
+ iglob5 = ibool(i,j,k+1,ispec)
+ iglob6 = ibool(i+1,j,k+1,ispec)
+ iglob7 = ibool(i+1,j+1,k+1,ispec)
+ iglob8 = ibool(i,j+1,k+1,ispec)
+ n1 = num_ibool(iglob1)+np-1
+ n2 = num_ibool(iglob2)+np-1
+ n3 = num_ibool(iglob3)+np-1
+ n4 = num_ibool(iglob4)+np-1
+ n5 = num_ibool(iglob5)+np-1
+ n6 = num_ibool(iglob6)+np-1
+ n7 = num_ibool(iglob7)+np-1
+ n8 = num_ibool(iglob8)+np-1
+ call write_integer(n1)
+ call write_integer(n2)
+ call write_integer(n3)
+ call write_integer(n4)
+ call write_integer(n5)
+ call write_integer(n6)
+ call write_integer(n7)
+ call write_integer(n8)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ endif
+ ne = ne + nelement
+
+ enddo ! num_node
+ if (ne /= nee) stop 'Number of total elements are not consistent'
+ print *, 'Total number of elements: ', ne
+
+ call close_file()
+
+ print *, 'Done writing '//trim(mesh_file)
+
+ end program combine_paraview_data
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/compile.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/compile.bash (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/compile.bash 2010-02-09 23:44:05 UTC (rev 16243)
@@ -0,0 +1,4 @@
+rm -rf xcombine_vol_data *.o
+mpif90 -O3 -c -o combine_vol_data_mod.o combine_vol_data_mod.f90
+cc -c -o write_c_binary.o write_c_binary.c
+mpif90 -O3 -o xcombine_vol_data combine_vol_data_mod.o write_c_binary.o
\ No newline at end of file
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/compile.bash
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/config.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/config.h 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/config.h 2010-02-09 23:44:05 UTC (rev 16243)
@@ -23,13 +23,13 @@
#define PACKAGE_BUGREPORT "jtromp AT caltech.edu"
/* Define to the full name of this package. */
-#define PACKAGE_NAME "Specfem 3D Basin"
+#define PACKAGE_NAME "Specfem 3D"
/* Define to the full name and version of this package. */
-#define PACKAGE_STRING "Specfem 3D Basin 1.4.0"
+#define PACKAGE_STRING "Specfem 3D 1.4.0"
/* Define to the one symbol short name of this package. */
-#define PACKAGE_TARNAME "Specfem3DBasin"
+#define PACKAGE_TARNAME "Specfem3D"
/* Define to the version of this package. */
#define PACKAGE_VERSION "1.4.0"
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/constants.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/constants.h 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/constants.h 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,18 +1,26 @@
!=====================================================================
!
-! S p e c f e m 3 D B a s i n V e r s i o n 1 . 4
-! --------------------------------------------------
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
!
! Dimitri Komatitsch and Jeroen Tromp
! Seismological Laboratory - California Institute of Technology
! (c) California Institute of Technology July 2005
!
-! A signed non-commercial agreement is required to use this program.
-! Please check http://www.gps.caltech.edu/research/jtromp for details.
-! Free for non-commercial academic research ONLY.
-! This program is distributed WITHOUT ANY WARRANTY whatsoever.
-! Do not redistribute this program without written permission.
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
!=====================================================================
! constants.h. Generated from constants.h.in by configure.
@@ -65,7 +73,7 @@
! density of sea water
real(kind=CUSTOM_REAL), parameter :: RHO_OCEANS = 1020.0
-! extend basin model below threshold and above topography to make sure
+! extend model below threshold and above topography to make sure
! there is no small gap between interpolated maps and sediments
logical, parameter :: EXTEND_VOXET_BELOW_BASEMENT = .true.
logical, parameter :: EXTEND_VOXET_ABOVE_TOPO = .true.
@@ -85,24 +93,27 @@
! number of sources to be gathered by MPI_Gather
integer, parameter :: NGATHER_SOURCES = 10000
-! source decay rate
- double precision, parameter :: SOURCE_DECAY_RATE = 1.628d0
+! we mimic a triangle of half duration equal to half_duration_triangle
+! using a Gaussian having a very close shape, as explained in Figure 4.2
+! of the manual. This source decay rate to mimic an equivalent triangle
+! was found by trial and error
+ double precision, parameter :: SOURCE_DECAY_MIMIC_TRIANGLE = 1.628d0
! ---------------------------------------------------------------------------------------
-! LQY -- Following 3 variables stays here temporarily,
+! LQY -- Following 3 variables stays here temporarily,
! we need to move them to Par_file at a proper time
! ---------------------------------------------------------------------------------------
! save moho mesh and compute Moho boundary kernels
- logical, parameter :: SAVE_MOHO_MESH = .true.
+ logical, parameter :: SAVE_MOHO_MESH = .false.
! number of steps to save the state variables in the forward simulation,
! to be used in the backward reconstruction in the presence of attenuation
integer, parameter :: NSTEP_Q_SAVE = 200
-! the scratch disk to save the state variables saved in the forward
+! the scratch disk to save the state variables saved in the forward
! simulation, this can be a global scratch disk in case you run out of
! space on the local scratch disk
- character(len=150), parameter :: LOCAL_PATH_Q = '/ibrixfs1/scratch/lqy'
+ character(len=150), parameter :: LOCAL_PATH_Q = '/ibrixfs1/scratch/lqy/DATABASES_MPI_Q/'
!------------------------------------------------------
!----------- do not modify anything below -------------
@@ -223,13 +234,6 @@
double precision, parameter :: DEGREES_PER_CELL_TOPO_SOCAL = 5.d0 / 1000.d0
character(len=100), parameter :: TOPO_FILE_SOCAL = 'DATA/la_topography/topo_bathy_final.dat'
-! size of topography and bathymetry file for Lacq (France) gas field
- integer, parameter :: NX_TOPO_LACQ = 499,NY_TOPO_LACQ = 401
- double precision, parameter :: ORIG_LAT_TOPO_LACQ = 100000.d0
- double precision, parameter :: ORIG_LONG_TOPO_LACQ = 340400.d0
- double precision, parameter :: DEGREES_PER_CELL_TOPO_LACQ = 100.d0
- character(len=100), parameter :: TOPO_FILE_LACQ = 'DATA/lacq_thomas/mnt_Lacq_Lambert_final_dimitri.dat'
-
! size of Lupei Zhu's Moho map file for Southern California
integer, parameter :: NX_MOHO = 71,NY_MOHO = 51
double precision, parameter :: ORIG_LAT_MOHO = 32.d0
@@ -309,7 +313,7 @@
!
! integer, parameter :: NGRID_NEW_HAUKSSON = 201
!
-!! corners of new Hauksson's interpolated grid
+!! corners of new Hauksson interpolated grid
! double precision, parameter :: UTM_X_ORIG_HAUKSSON = 122035.012d0
! double precision, parameter :: UTM_X_END_HAUKSSON = 766968.628d0
! double precision, parameter :: UTM_Y_ORIG_HAUKSSON = 3547232.986d0
@@ -330,7 +334,7 @@
integer, parameter :: NGRID_NEW_HAUKSSON = 241
-! corners of new Hauksson's interpolated grid
+! corners of new Hauksson interpolated grid
double precision, parameter :: UTM_X_ORIG_HAUKSSON = 88021.568d0
double precision, parameter :: UTM_X_END_HAUKSSON = 861517.886d0
double precision, parameter :: UTM_Y_ORIG_HAUKSSON = 3404059.875d0
@@ -347,13 +351,7 @@
double precision, parameter :: DEPTH_5p5km_SOCAL = -5500.d0
double precision, parameter :: DEPTH_16km_SOCAL = -16000.d0
double precision, parameter :: DEPTH_MOHO_SOCAL = -32000.d0
-! double precision, parameter :: DEPTH_5p5km_SOCAL = -7700.d0
-! double precision, parameter :: DEPTH_16km_SOCAL = -18200.d0
-! double precision, parameter :: DEPTH_MOHO_SOCAL = -34200.d0
-! layer for Lacq (France) gas field
- double precision, parameter :: DEPTH_INTERFACE_LACQ = -7000.d0
-
! reference surface of the model before adding topography
double precision, parameter :: Z_SURFACE = 0.d0
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/values_from_mesher.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/values_from_mesher.h 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/values_from_mesher.h 2010-02-09 23:44:05 UTC (rev 16243)
@@ -10,6 +10,9 @@
! number of ES nodes = 21.00000
! percentage of total 640 ES nodes = 3.281250 %
! total memory available on these ES nodes (Gb) = 336.0000
+ integer, parameter :: NPROC_VAL = 168
+ integer, parameter :: NPROC_XI_VAL = 14
+ integer, parameter :: NPROC_ETA_VAL = 12
!
! max points per processor = max vector length = 163941
! min vector length = 25
@@ -42,13 +45,29 @@
! average distance between points along Y in m = 436.8476
!
- integer, parameter :: NSPEC_AB = 2412
- integer, parameter :: NGLOB_AB = 163941
+ integer, parameter :: NSPEC_AB_VAL = 2412
+ integer, parameter :: NGLOB_AB_VAL = 163941
!
! number of time steps = 27300
!
- integer, parameter :: NSPEC_ATTENUATION = 1
+ integer, parameter :: NSPEC_ATTENUATION = 1
logical, parameter :: ATTENUATION_VAL = .false.
+
+ integer, parameter :: NSPEC_ANISO = 1
logical, parameter :: ANISOTROPY_VAL = .false.
+ integer, parameter :: NSPEC_ATT_AND_KERNEL = 1
+ integer, parameter :: NSPEC_ADJOINT = 1
+ integer, parameter :: NGLOB_ADJOINT = 1
+
+ integer, parameter :: NSPEC2DMAX_XMIN_XMAX_VAL = 150
+ integer, parameter :: NSPEC2DMAX_YMIN_YMAX_VAL = 150
+ integer, parameter :: NSPEC2D_BOTTOM_VAL = 36
+ integer, parameter :: NSPEC2D_TOP_VAL = 576
+ integer, parameter :: NPOIN2DMAX_XMIN_XMAX_VAL = 3751
+ integer, parameter :: NPOIN2DMAX_YMIN_YMAX_VAL = 3751
+ integer, parameter :: NPOIN2DMAX_XY_VAL = 3751
+
+ integer, parameter :: NSPEC2D_MOHO_BOUN = 1
+ integer, parameter :: NSPEC_BOUN = 1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/write_c_binary.c
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/write_c_binary.c 2010-02-09 22:51:38 UTC (rev 16242)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/write_c_binary.c 2010-02-09 23:44:05 UTC (rev 16243)
@@ -1,19 +1,27 @@
/*
!=====================================================================
!
-! S p e c f e m 3 D B a s i n V e r s i o n 1 . 4
-! --------------------------------------------------
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
!
! Dimitri Komatitsch and Jeroen Tromp
! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology July 2005
+! (c) California Institute of Technology September 2006
!
-! A signed non-commercial agreement is required to use this program.
-! Please check http://www.gps.caltech.edu/research/jtromp for details.
-! Free for non-commercial academic research ONLY.
-! This program is distributed WITHOUT ANY WARRANTY whatsoever.
-! Do not redistribute this program without written permission.
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
!=====================================================================
*/
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/topo_input/combine_vol_data/xcombine_vol_data
===================================================================
(Binary files differ)
More information about the CIG-COMMITS
mailing list