[cig-commits] r16343 - in seismo/3D/SPECFEM3D_SESAME/trunk: . DATA EXAMPLES/homogeneous_halfspace EXAMPLES/layered_halfspace EXAMPLES/waterlayered_halfspace decompose_mesh_SCOTCH
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Thu Feb 25 15:07:31 PST 2010
Author: danielpeter
Date: 2010-02-25 15:07:29 -0800 (Thu, 25 Feb 2010)
New Revision: 16343
Added:
seismo/3D/SPECFEM3D_SESAME/trunk/DATA/CMTSOLUTION.in
seismo/3D/SPECFEM3D_SESAME/trunk/DATA/Par_file.in
seismo/3D/SPECFEM3D_SESAME/trunk/DATA/STATIONS.in
seismo/3D/SPECFEM3D_SESAME/trunk/DATABASES_MPI/
seismo/3D/SPECFEM3D_SESAME/trunk/create_mass_matrices.f90
seismo/3D/SPECFEM3D_SESAME/trunk/get_MPI.f90
seismo/3D/SPECFEM3D_SESAME/trunk/get_absorbing_boundary.f90
seismo/3D/SPECFEM3D_SESAME/trunk/get_coupling_surfaces.f90
seismo/3D/SPECFEM3D_SESAME/trunk/get_model.f90
seismo/3D/SPECFEM3D_SESAME/trunk/model_aniso.f90
seismo/3D/SPECFEM3D_SESAME/trunk/model_external_values.f90
seismo/3D/SPECFEM3D_SESAME/trunk/save_moho_arrays.f90
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/README
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/README
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/README
seismo/3D/SPECFEM3D_SESAME/trunk/INSTALL
seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in
seismo/3D/SPECFEM3D_SESAME/trunk/configure
seismo/3D/SPECFEM3D_SESAME/trunk/configure.ac
seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in
seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90
Log:
added a generic file model_external_values.f90 to superimpose external velocities on the mesh points when flag is set for USE_MODEL_EXTERNAL_VALUES in constants.h; updated READMEs in EXAMPLES/ and configure script
Added: seismo/3D/SPECFEM3D_SESAME/trunk/DATA/CMTSOLUTION.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/DATA/CMTSOLUTION.in (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/DATA/CMTSOLUTION.in 2010-02-25 23:07:29 UTC (rev 16343)
@@ -0,0 +1,13 @@
+PDE 1999 01 01 00 00 00.00 67000 67000 -25000 4.2 4.2 hom_explosion
+event name: hom_explosion
+time shift: 0.0000
+half duration: 5.0
+latitude: 67000.0
+longitude: 67000.0
+depth: -25000.0
+Mrr: 1.000000e+23
+Mtt: 1.000000e+23
+Mpp: 1.000000e+23
+Mrt: 0.000000
+Mrp: 0.000000
+Mtp: 0.000000
Added: seismo/3D/SPECFEM3D_SESAME/trunk/DATA/Par_file.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/DATA/Par_file.in (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/DATA/Par_file.in 2010-02-25 23:07:29 UTC (rev 16343)
@@ -0,0 +1,51 @@
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint, 3 = both simultaneously
+SAVE_FORWARD = .false.
+
+# UTM projection parameters
+UTM_PROJECTION_ZONE = 11
+SUPPRESS_UTM_PROJECTION = .true.
+
+# number of MPI processors
+NPROC = 4
+
+# time step parameters
+NSTEP = 1000
+DT = 0.05d0
+
+# parameters describing the model
+OCEANS = .false.
+TOPOGRAPHY = .false.
+ATTENUATION = .false.
+USE_OLSEN_ATTENUATION = .false.
+ANISOTROPY = .false.
+
+# absorbing boundary conditions for a regional simulation
+ABSORBING_CONDITIONS = .false.
+
+# save AVS or OpenDX movies
+MOVIE_SURFACE = .false.
+MOVIE_VOLUME = .false.
+NTSTEP_BETWEEN_FRAMES = 200
+CREATE_SHAKEMAP = .false.
+SAVE_DISPLACEMENT = .false.
+USE_HIGHRES_FOR_MOVIES = .false.
+HDUR_MOVIE = 0.0
+
+# save AVS or OpenDX mesh files to check the mesh
+SAVE_MESH_FILES = .true.
+
+# path to store the local database file on each node
+LOCAL_PATH = DATABASES_MPI
+
+# interval at which we output time step info and max of norm of displacement
+NTSTEP_BETWEEN_OUTPUT_INFO = 500
+
+# interval in time steps for writing of seismograms
+NTSTEP_BETWEEN_OUTPUT_SEISMOS = 10000
+
+# print source time function
+PRINT_SOURCE_TIME_FUNCTION = .false.
+
+
Added: seismo/3D/SPECFEM3D_SESAME/trunk/DATA/STATIONS.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/DATA/STATIONS.in (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/DATA/STATIONS.in 2010-02-25 23:07:29 UTC (rev 16343)
@@ -0,0 +1,224 @@
+X1 DB 67000.00 0.000000 0.0 0.0
+X2 DB 67000.00 1196.429 0.0 0.0
+X3 DB 67000.00 2392.857 0.0 0.0
+X4 DB 67000.00 3589.286 0.0 0.0
+X5 DB 67000.00 4785.714 0.0 0.0
+X6 DB 67000.00 5982.143 0.0 0.0
+X7 DB 67000.00 7178.571 0.0 0.0
+X8 DB 67000.00 8375.000 0.0 0.0
+X9 DB 67000.00 9571.429 0.0 0.0
+X10 DB 67000.00 10767.86 0.0 0.0
+X11 DB 67000.00 11964.29 0.0 0.0
+X12 DB 67000.00 13160.71 0.0 0.0
+X13 DB 67000.00 14357.14 0.0 0.0
+X14 DB 67000.00 15553.57 0.0 0.0
+X15 DB 67000.00 16750.00 0.0 0.0
+X16 DB 67000.00 17946.43 0.0 0.0
+X17 DB 67000.00 19142.86 0.0 0.0
+X18 DB 67000.00 20339.29 0.0 0.0
+X19 DB 67000.00 21535.71 0.0 0.0
+X20 DB 67000.00 22732.14 0.0 0.0
+X21 DB 67000.00 23928.57 0.0 0.0
+X22 DB 67000.00 25125.00 0.0 0.0
+X23 DB 67000.00 26321.43 0.0 0.0
+X24 DB 67000.00 27517.86 0.0 0.0
+X25 DB 67000.00 28714.29 0.0 0.0
+X26 DB 67000.00 29910.71 0.0 0.0
+X27 DB 67000.00 31107.14 0.0 0.0
+X28 DB 67000.00 32303.57 0.0 0.0
+X29 DB 67000.00 33500.00 0.0 0.0
+X30 DB 67000.00 34696.43 0.0 0.0
+X31 DB 67000.00 35892.86 0.0 0.0
+X32 DB 67000.00 37089.29 0.0 0.0
+X33 DB 67000.00 38285.71 0.0 0.0
+X34 DB 67000.00 39482.14 0.0 0.0
+X35 DB 67000.00 40678.57 0.0 0.0
+X36 DB 67000.00 41875.00 0.0 0.0
+X37 DB 67000.00 43071.43 0.0 0.0
+X38 DB 67000.00 44267.86 0.0 0.0
+X39 DB 67000.00 45464.29 0.0 0.0
+X40 DB 67000.00 46660.71 0.0 0.0
+X41 DB 67000.00 47857.14 0.0 0.0
+X42 DB 67000.00 49053.57 0.0 0.0
+X43 DB 67000.00 50250.00 0.0 0.0
+X44 DB 67000.00 51446.43 0.0 0.0
+X45 DB 67000.00 52642.86 0.0 0.0
+X46 DB 67000.00 53839.29 0.0 0.0
+X47 DB 67000.00 55035.71 0.0 0.0
+X48 DB 67000.00 56232.14 0.0 0.0
+X49 DB 67000.00 57428.57 0.0 0.0
+X50 DB 67000.00 58625.00 0.0 0.0
+X51 DB 67000.00 59821.43 0.0 0.0
+X52 DB 67000.00 61017.86 0.0 0.0
+X53 DB 67000.00 62214.29 0.0 0.0
+X54 DB 67000.00 63410.71 0.0 0.0
+X55 DB 67000.00 64607.14 0.0 0.0
+X56 DB 67000.00 65803.57 0.0 0.0
+X57 DB 67000.00 67000.00 0.0 0.0
+X58 DB 67000.00 68196.43 0.0 0.0
+X59 DB 67000.00 69392.86 0.0 0.0
+X60 DB 67000.00 70589.29 0.0 0.0
+X61 DB 67000.00 71785.71 0.0 0.0
+X62 DB 67000.00 72982.14 0.0 0.0
+X63 DB 67000.00 74178.57 0.0 0.0
+X64 DB 67000.00 75375.00 0.0 0.0
+X65 DB 67000.00 76571.43 0.0 0.0
+X66 DB 67000.00 77767.86 0.0 0.0
+X67 DB 67000.00 78964.29 0.0 0.0
+X68 DB 67000.00 80160.71 0.0 0.0
+X69 DB 67000.00 81357.14 0.0 0.0
+X70 DB 67000.00 82553.57 0.0 0.0
+X71 DB 67000.00 83750.00 0.0 0.0
+X72 DB 67000.00 84946.43 0.0 0.0
+X73 DB 67000.00 86142.86 0.0 0.0
+X74 DB 67000.00 87339.29 0.0 0.0
+X75 DB 67000.00 88535.71 0.0 0.0
+X76 DB 67000.00 89732.14 0.0 0.0
+X77 DB 67000.00 90928.57 0.0 0.0
+X78 DB 67000.00 92125.00 0.0 0.0
+X79 DB 67000.00 93321.43 0.0 0.0
+X80 DB 67000.00 94517.86 0.0 0.0
+X81 DB 67000.00 95714.29 0.0 0.0
+X82 DB 67000.00 96910.71 0.0 0.0
+X83 DB 67000.00 98107.14 0.0 0.0
+X84 DB 67000.00 99303.57 0.0 0.0
+X85 DB 67000.00 100500.0 0.0 0.0
+X86 DB 67000.00 101696.4 0.0 0.0
+X87 DB 67000.00 102892.9 0.0 0.0
+X88 DB 67000.00 104089.3 0.0 0.0
+X89 DB 67000.00 105285.7 0.0 0.0
+X90 DB 67000.00 106482.1 0.0 0.0
+X91 DB 67000.00 107678.6 0.0 0.0
+X92 DB 67000.00 108875.0 0.0 0.0
+X93 DB 67000.00 110071.4 0.0 0.0
+X94 DB 67000.00 111267.9 0.0 0.0
+X95 DB 67000.00 112464.3 0.0 0.0
+X96 DB 67000.00 113660.7 0.0 0.0
+X97 DB 67000.00 114857.1 0.0 0.0
+X98 DB 67000.00 116053.6 0.0 0.0
+X99 DB 67000.00 117250.0 0.0 0.0
+X100 DB 67000.00 118446.4 0.0 0.0
+X101 DB 67000.00 119642.9 0.0 0.0
+X102 DB 67000.00 120839.3 0.0 0.0
+X103 DB 67000.00 122035.7 0.0 0.0
+X104 DB 67000.00 123232.1 0.0 0.0
+X105 DB 67000.00 124428.6 0.0 0.0
+X106 DB 67000.00 125625.0 0.0 0.0
+X107 DB 67000.00 126821.4 0.0 0.0
+X108 DB 67000.00 128017.9 0.0 0.0
+X109 DB 67000.00 129214.3 0.0 0.0
+X110 DB 67000.00 130410.7 0.0 0.0
+X111 DB 67000.00 131607.1 0.0 0.0
+X112 DB 67000.00 132803.6 0.0 0.0
+Y1 DB 0.000000 67000.00 0.0 0.0
+Y2 DB 1196.429 67000.00 0.0 0.0
+Y3 DB 2392.857 67000.00 0.0 0.0
+Y4 DB 3589.286 67000.00 0.0 0.0
+Y5 DB 4785.714 67000.00 0.0 0.0
+Y6 DB 5982.143 67000.00 0.0 0.0
+Y7 DB 7178.571 67000.00 0.0 0.0
+Y8 DB 8375.000 67000.00 0.0 0.0
+Y9 DB 9571.429 67000.00 0.0 0.0
+Y10 DB 10767.86 67000.00 0.0 0.0
+Y11 DB 11964.29 67000.00 0.0 0.0
+Y12 DB 13160.71 67000.00 0.0 0.0
+Y13 DB 14357.14 67000.00 0.0 0.0
+Y14 DB 15553.57 67000.00 0.0 0.0
+Y15 DB 16750.00 67000.00 0.0 0.0
+Y16 DB 17946.43 67000.00 0.0 0.0
+Y17 DB 19142.86 67000.00 0.0 0.0
+Y18 DB 20339.29 67000.00 0.0 0.0
+Y19 DB 21535.71 67000.00 0.0 0.0
+Y20 DB 22732.14 67000.00 0.0 0.0
+Y21 DB 23928.57 67000.00 0.0 0.0
+Y22 DB 25125.00 67000.00 0.0 0.0
+Y23 DB 26321.43 67000.00 0.0 0.0
+Y24 DB 27517.86 67000.00 0.0 0.0
+Y25 DB 28714.29 67000.00 0.0 0.0
+Y26 DB 29910.71 67000.00 0.0 0.0
+Y27 DB 31107.14 67000.00 0.0 0.0
+Y28 DB 32303.57 67000.00 0.0 0.0
+Y29 DB 33500.00 67000.00 0.0 0.0
+Y30 DB 34696.43 67000.00 0.0 0.0
+Y31 DB 35892.86 67000.00 0.0 0.0
+Y32 DB 37089.29 67000.00 0.0 0.0
+Y33 DB 38285.71 67000.00 0.0 0.0
+Y34 DB 39482.14 67000.00 0.0 0.0
+Y35 DB 40678.57 67000.00 0.0 0.0
+Y36 DB 41875.00 67000.00 0.0 0.0
+Y37 DB 43071.43 67000.00 0.0 0.0
+Y38 DB 44267.86 67000.00 0.0 0.0
+Y39 DB 45464.29 67000.00 0.0 0.0
+Y40 DB 46660.71 67000.00 0.0 0.0
+Y41 DB 47857.14 67000.00 0.0 0.0
+Y42 DB 49053.57 67000.00 0.0 0.0
+Y43 DB 50250.00 67000.00 0.0 0.0
+Y44 DB 51446.43 67000.00 0.0 0.0
+Y45 DB 52642.86 67000.00 0.0 0.0
+Y46 DB 53839.29 67000.00 0.0 0.0
+Y47 DB 55035.71 67000.00 0.0 0.0
+Y48 DB 56232.14 67000.00 0.0 0.0
+Y49 DB 57428.57 67000.00 0.0 0.0
+Y50 DB 58625.00 67000.00 0.0 0.0
+Y51 DB 59821.43 67000.00 0.0 0.0
+Y52 DB 61017.86 67000.00 0.0 0.0
+Y53 DB 62214.29 67000.00 0.0 0.0
+Y54 DB 63410.71 67000.00 0.0 0.0
+Y55 DB 64607.14 67000.00 0.0 0.0
+Y56 DB 65803.57 67000.00 0.0 0.0
+Y57 DB 67000.00 67000.00 0.0 0.0
+Y58 DB 68196.43 67000.00 0.0 0.0
+Y59 DB 69392.86 67000.00 0.0 0.0
+Y60 DB 70589.29 67000.00 0.0 0.0
+Y61 DB 71785.71 67000.00 0.0 0.0
+Y62 DB 72982.14 67000.00 0.0 0.0
+Y63 DB 74178.57 67000.00 0.0 0.0
+Y64 DB 75375.00 67000.00 0.0 0.0
+Y65 DB 76571.43 67000.00 0.0 0.0
+Y66 DB 77767.86 67000.00 0.0 0.0
+Y67 DB 78964.29 67000.00 0.0 0.0
+Y68 DB 80160.71 67000.00 0.0 0.0
+Y69 DB 81357.14 67000.00 0.0 0.0
+Y70 DB 82553.57 67000.00 0.0 0.0
+Y71 DB 83750.00 67000.00 0.0 0.0
+Y72 DB 84946.43 67000.00 0.0 0.0
+Y73 DB 86142.86 67000.00 0.0 0.0
+Y74 DB 87339.29 67000.00 0.0 0.0
+Y75 DB 88535.71 67000.00 0.0 0.0
+Y76 DB 89732.14 67000.00 0.0 0.0
+Y77 DB 90928.57 67000.00 0.0 0.0
+Y78 DB 92125.00 67000.00 0.0 0.0
+Y79 DB 93321.43 67000.00 0.0 0.0
+Y80 DB 94517.86 67000.00 0.0 0.0
+Y81 DB 95714.29 67000.00 0.0 0.0
+Y82 DB 96910.71 67000.00 0.0 0.0
+Y83 DB 98107.14 67000.00 0.0 0.0
+Y84 DB 99303.57 67000.00 0.0 0.0
+Y85 DB 100500.0 67000.00 0.0 0.0
+Y86 DB 101696.4 67000.00 0.0 0.0
+Y87 DB 102892.9 67000.00 0.0 0.0
+Y88 DB 104089.3 67000.00 0.0 0.0
+Y89 DB 105285.7 67000.00 0.0 0.0
+Y90 DB 106482.1 67000.00 0.0 0.0
+Y91 DB 107678.6 67000.00 0.0 0.0
+Y92 DB 108875.0 67000.00 0.0 0.0
+Y93 DB 110071.4 67000.00 0.0 0.0
+Y94 DB 111267.9 67000.00 0.0 0.0
+Y95 DB 112464.3 67000.00 0.0 0.0
+Y96 DB 113660.7 67000.00 0.0 0.0
+Y97 DB 114857.1 67000.00 0.0 0.0
+Y98 DB 116053.6 67000.00 0.0 0.0
+Y99 DB 117250.0 67000.00 0.0 0.0
+Y100 DB 118446.4 67000.00 0.0 0.0
+Y101 DB 119642.9 67000.00 0.0 0.0
+Y102 DB 120839.3 67000.00 0.0 0.0
+Y103 DB 122035.7 67000.00 0.0 0.0
+Y104 DB 123232.1 67000.00 0.0 0.0
+Y105 DB 124428.6 67000.00 0.0 0.0
+Y106 DB 125625.0 67000.00 0.0 0.0
+Y107 DB 126821.4 67000.00 0.0 0.0
+Y108 DB 128017.9 67000.00 0.0 0.0
+Y109 DB 129214.3 67000.00 0.0 0.0
+Y110 DB 130410.7 67000.00 0.0 0.0
+Y111 DB 131607.1 67000.00 0.0 0.0
+Y112 DB 132803.6 67000.00 0.0 0.0
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/README
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/README 2010-02-25 21:58:34 UTC (rev 16342)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/README 2010-02-25 23:07:29 UTC (rev 16343)
@@ -5,14 +5,14 @@
step-by-step tutorial:
0. check that all software is available (or that modules are loaded):
- intel/openmpi, cubit, scotch, gnuplot
+ intel/openmpi, cubit, scotch, python, gnuplot
1. configure package:
- - From the SPECFEM3D_SESAME root directory /SPECFEM3D_SESAME/
+ - From the SPECFEM3D_SESAME root directory SPECFEM3D_SESAME/
configure the package, e.g. using intel's ifort compiler:
- > cd /SPECFEM3D_SESAME
+ > cd SPECFEM3D_SESAME
> ./configure F90=ifort
If successful, this will generate the files:
@@ -20,11 +20,16 @@
- make a directory for the databases
> mkdir DATABASES_MPI
+
+ - copy two run scripts from SPECFEM3D_SESAME/UTILS/Cluster/
+ into SPECFEM3D_SESAME/, e.g.,
+ pbs/go_generate_databases_pbs.bash
+ pbs/go_solver_pbs.bash
2. create mesh:
- - change to the examples directory /SPECFEM3D_SESAME/EXAMPLES/homogeneous_halfspace:
+ - change to the examples directory SPECFEM3D_SESAME/EXAMPLES/homogeneous_halfspace:
> cd EXAMPLES/homogeneous_halfspace
- open the cubit GUI:
@@ -48,7 +53,7 @@
3. decompose mesh files:
- - compile decomposer in directory /SPECFEM3D_SESAME/decompose_mesh_SCOTCH/:
+ - compile decomposer in directory SPECFEM3D_SESAME/decompose_mesh_SCOTCH/:
> make
NOTE 1: check that the two scotch libraries are properly specified in Makefile
@@ -66,17 +71,17 @@
4. generate databases:
- copy three files -- Par_file CMTSOLUTION STATIONS -- from
- /SPECFEM3D_SESAME/EXAMPLES/homogeneous_halfspace/ to /SPECFEM3D_SESAME/DATA/
+ SPECFEM3D_SESAME/EXAMPLES/homogeneous_halfspace/ to SPECFEM3D_SESAME/DATA/
- - compile generate_databases from /SPECFEM3D_SESAME/ :
- > cd /SPECFEM3D_SESAME
+ - compile generate_databases from SPECFEM3D_SESAME/ :
+ > cd SPECFEM3D_SESAME
> make xgenerate_databases
- submit job script
- > qsub go_generate_database_pbs.sesame.bash
+ > qsub go_generate_databases_pbs.bash
NOTE: this script will need to be tailored to your cluster, e.g.,
- > bsub < go_generate_database_lsf.sesame.bash
+ > bsub < go_generate_databases_lsf.bash
this will create binary mesh files, e.g. "proc000***_external_mesh.bin"
in directory DATABASES_MPI/.
@@ -88,20 +93,23 @@
> make xspecfem3D
- submit job script:
- > qsub go_solver_pbs.sesame.bash
+ > qsub go_solver_pbs.bash
NOTE 1: this script will need to be tailored to your cluster, e.g.,
- > bsub < go_solver_lsf.sesame.bash
+ > bsub < go_solver_lsf.bash
NOTE 2: the simulation runs on 4 cores and should take about 15 minutes,
and you can track the progress with the timestamp files
- NOTE 3: you should have 672 x 3 (semd,semv,sema) seismograms in OUTPUT_FILES
+ - when the job is complete, you should have 3 sets (semd,semv,sema)
+ of 672 (ls -1 *semd | wc) in the directory OUTPUT_FILES,
+ as well as 11 timestamp****** files
+
6. check with 6 reference seismograms in
- /SPECFEM3D_SESAME/EXAMPLES/homogeneous_halfspace/REF_SEIS/
+ SPECFEM3D_SESAME/EXAMPLES/homogeneous_halfspace/REF_SEIS/
- execute gnuplot script
- > cd /SPECFEM3D_SESAME/EXAMPLES/homogeneous_halfspace/REF_SEIS
+ > cd SPECFEM3D_SESAME/EXAMPLES/homogeneous_halfspace/REF_SEIS
> gnuplot
gnuplot> load "X1-50.BHZ.gnuplot"
@@ -122,7 +130,7 @@
gnuplot> quit
-7. plot your output seismograms in /SPECFEM3D_SESAME/OUTPUT_FILES/
+7. plot your output seismograms in SPECFEM3D_SESAME/OUTPUT_FILES/
- copy gnuplot script X1-50.BHZ.gnuplot to OUTPUT_FILES/
- execute same commands as above
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/README
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/README 2010-02-25 21:58:34 UTC (rev 16342)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/README 2010-02-25 23:07:29 UTC (rev 16343)
@@ -62,7 +62,7 @@
> make xgenerate_databases
- submit job script:
- > qsub go_generate_database_pbs.sesame.bash
+ > qsub go_generate_databases_pbs.bash
(note: if execution fails due to memory shortage - most likely fails when calling routine to
create regional mesh - then try to increase process memory stack size: ulimit -s 2000000 (2GB) )
@@ -74,7 +74,7 @@
> make xspecfem3D
- submit job script:
- > qsub go_solver_pbs.sesame.bash
+ > qsub go_solver_pbs.bash
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/README
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/README 2010-02-25 21:58:34 UTC (rev 16342)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/README 2010-02-25 23:07:29 UTC (rev 16343)
@@ -39,7 +39,7 @@
> make xgenerate_databases
- submit job script:
- > qsub go_generate_database_pbs.sesame.bash
+ > qsub go_generate_databases_pbs.bash
(note: if execution fails due to memory shortage - most likely fails when calling routine to
create regional mesh - then try to increase process memory stack size: ulimit -s 2000000 (2GB) )
@@ -51,7 +51,7 @@
> make xspecfem3D
- submit job script:
- > qsub go_solver_pbs.sesame.bash
+ > qsub go_solver_pbs.bash
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/INSTALL
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/INSTALL 2010-02-25 21:58:34 UTC (rev 16342)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/INSTALL 2010-02-25 23:07:29 UTC (rev 16343)
@@ -61,6 +61,9 @@
./Makefile
./constants.h
./precision.h
+ ./DATA/Par_file
+ ./DATA/CMTSOLUTION
+ ./DATA/STATIONS
please make sure, your installation is working and that
the created files 'constant.h' and 'Makefile' satisfy
@@ -76,10 +79,7 @@
./DATA/STATIONS
default simulation files are provided in the directory DATA/ which you can
- use as the ones needed above:
- > cp DATA/Par_file.default DATA/Par_file
- > cp DATA/CMTSOLUTION.default DATA/CMTSOLUTION
- > cp DATA/STATIONS.default DATA/STATIONS
+ use as the ones needed above.
then type
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in 2010-02-25 21:58:34 UTC (rev 16342)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in 2010-02-25 23:07:29 UTC (rev 16343)
@@ -57,7 +57,6 @@
O = obj
libspecfem_a_OBJECTS = \
- $O/aniso_model.o \
$O/assemble_MPI_scalar.o \
$O/calc_jacobian.o \
$O/check_mesh_resolution.o \
@@ -67,10 +66,15 @@
$O/create_header_file.o \
$O/create_name_database.o \
$O/create_regions_mesh.o \
+ $O/create_mass_matrices.o \
$O/create_serial_name_database.o \
$O/define_derivation_matrices.o \
$O/detect_surface.o \
$O/exit_mpi.o \
+ $O/get_absorbing_boundary.o \
+ $O/get_coupling_surfaces.o \
+ $O/get_model.o \
+ $O/get_MPI.o \
$O/get_MPI_cutplanes_eta.o \
$O/get_MPI_cutplanes_xi.o \
$O/get_attenuation_model.o \
@@ -87,6 +91,8 @@
$O/locate_receivers.o \
$O/locate_source.o \
$O/generate_databases.o \
+ $O/model_external_values.o \
+ $O/model_aniso.o \
$O/netlib_specfun_erf.o \
$O/param_reader.o \
$O/prepare_assemble_MPI.o \
@@ -96,6 +102,7 @@
$O/recompute_jacobian.o \
$O/save_arrays_solver.o \
$O/save_header_file.o \
+ $O/save_moho_arrays.o \
$O/sort_array_coordinates.o \
$O/utm_geo.o \
$O/write_AVS_DX_global_data.o \
@@ -197,7 +204,7 @@
all: clean default
backup:
- cp *f90 *h README_SPECFEM3D DATA/Par_file* Makefile go_mesher go_solver mymachines bak
+ cp *f90 *h README_SPECFEM3D DATA/Par_file* Makefile go_generate_databases* go_mesher* go_solver* mymachines bak
bak: backup
@@ -382,12 +389,16 @@
$O/parallel.o: constants.h parallel.f90
${MPIFCCOMPILE_CHECK} -c -o $O/parallel.o parallel.f90
+$O/model_external_values.o: constants.h model_external_values.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/model_external_values.o model_external_values.f90
+
+
###
### serial compilation without optimization
###
-$O/aniso_model.o: constants.h aniso_model.f90
- ${FCCOMPILE_CHECK} -c -o $O/aniso_model.o aniso_model.f90
+$O/model_aniso.o: constants.h model_aniso.f90
+ ${FCCOMPILE_CHECK} -c -o $O/model_aniso.o model_aniso.f90
$O/serial.o: constants.h exit_mpi.f90
${FCCOMPILE_CHECK} -c -o $O/serial.o serial.f90
@@ -524,9 +535,29 @@
$O/recompute_jacobian.o: constants.h recompute_jacobian.f90
${FCCOMPILE_CHECK} -c -o $O/recompute_jacobian.o recompute_jacobian.f90
+
$O/create_regions_mesh.o: constants.h create_regions_mesh.f90
${FCCOMPILE_CHECK} -c -o $O/create_regions_mesh.o create_regions_mesh.f90
+$O/create_mass_matrices.o: constants.h create_mass_matrices.f90
+ ${FCCOMPILE_CHECK} -c -o $O/create_mass_matrices.o create_mass_matrices.f90
+
+$O/get_absorbing_boundary.o: constants.h get_absorbing_boundary.f90
+ ${FCCOMPILE_CHECK} -c -o $O/get_absorbing_boundary.o get_absorbing_boundary.f90
+
+$O/get_coupling_surfaces.o: constants.h get_coupling_surfaces.f90
+ ${FCCOMPILE_CHECK} -c -o $O/get_coupling_surfaces.o get_coupling_surfaces.f90
+
+$O/get_model.o: constants.h get_model.f90
+ ${FCCOMPILE_CHECK} -c -o $O/get_model.o get_model.f90
+
+$O/get_MPI.o: constants.h get_MPI.f90
+ ${FCCOMPILE_CHECK} -c -o $O/get_MPI.o get_MPI.f90
+
+$O/save_moho_arrays.o: constants.h save_moho_arrays.f90
+ ${FCCOMPILE_CHECK} -c -o $O/save_moho_arrays.o save_moho_arrays.f90
+
+
$O/create_name_database.o: constants.h create_name_database.f90
${FCCOMPILE_CHECK} -c -o $O/create_name_database.o create_name_database.f90
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/configure
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/configure 2010-02-25 21:58:34 UTC (rev 16342)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/configure 2010-02-25 23:07:29 UTC (rev 16343)
@@ -5397,7 +5397,7 @@
# Checks for library functions.
-ac_config_files="$ac_config_files Makefile constants.h precision.h"
+ac_config_files="$ac_config_files Makefile constants.h precision.h DATA/Par_file DATA/CMTSOLUTION DATA/STATIONS"
cat >confcache <<\_ACEOF
# This file is a shell script that caches the results of configure
@@ -5969,6 +5969,9 @@
"Makefile") CONFIG_FILES="$CONFIG_FILES Makefile" ;;
"constants.h") CONFIG_FILES="$CONFIG_FILES constants.h" ;;
"precision.h") CONFIG_FILES="$CONFIG_FILES precision.h" ;;
+ "DATA/Par_file") CONFIG_FILES="$CONFIG_FILES DATA/Par_file" ;;
+ "DATA/CMTSOLUTION") CONFIG_FILES="$CONFIG_FILES DATA/CMTSOLUTION" ;;
+ "DATA/STATIONS") CONFIG_FILES="$CONFIG_FILES DATA/STATIONS" ;;
*) { { echo "$as_me:$LINENO: error: invalid argument: $ac_config_target" >&5
echo "$as_me: error: invalid argument: $ac_config_target" >&2;}
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/configure.ac
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/configure.ac 2010-02-25 21:58:34 UTC (rev 16342)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/configure.ac 2010-02-25 23:07:29 UTC (rev 16343)
@@ -173,7 +173,7 @@
# Checks for library functions.
-AC_CONFIG_FILES([Makefile constants.h precision.h])
+AC_CONFIG_FILES([Makefile constants.h precision.h DATA/Par_file DATA/CMTSOLUTION DATA/STATIONS])
AC_OUTPUT
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in 2010-02-25 21:58:34 UTC (rev 16342)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in 2010-02-25 23:07:29 UTC (rev 16343)
@@ -48,6 +48,9 @@
! if running on a Beowulf, also modify name of nodes in filter_machine_file.f90
logical, parameter :: LOCAL_PATH_IS_ALSO_GLOBAL = . at LOCAL_PATH_IS_ALSO_GLOBAL@.
+! adds/superimposes velocity model values from 'model_external_values.f90'
+ logical, parameter :: USE_MODEL_EXTERNAL_VALUES = .false.
+
! use inlined products of Deville et al. (2002) to speedup the calculations to compute internal forces
logical, parameter :: USE_DEVILLE_PRODUCTS = .true.
Added: seismo/3D/SPECFEM3D_SESAME/trunk/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_mass_matrices.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_mass_matrices.f90 2010-02-25 23:07:29 UTC (rev 16343)
@@ -0,0 +1,265 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+ subroutine create_mass_matrices(nglob,nspec,ibool)
+
+! returns precomputed mass matrix in rmass array
+
+ use create_regions_mesh_ext_par
+ implicit none
+
+! number of spectral elements in each block
+ integer :: nspec
+ integer :: nglob
+
+! arrays with the mesh global indices
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+! local parameters
+ double precision :: weight
+ real(kind=CUSTOM_REAL) :: jacobianl
+ integer :: ispec,i,j,k,iglob,ier
+
+! allocates memory
+ allocate(rmass(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
+ allocate(rmass_acoustic(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
+ allocate(rmass_solid_poroelastic(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
+ allocate(rmass_fluid_poroelastic(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
+
+! creates mass matrix
+ rmass(:) = 0._CUSTOM_REAL
+ rmass_acoustic(:) = 0._CUSTOM_REAL
+ rmass_solid_poroelastic(:) = 0._CUSTOM_REAL
+ rmass_fluid_poroelastic(:) = 0._CUSTOM_REAL
+
+ do ispec=1,nspec
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+ jacobianl = jacobianstore(i,j,k,ispec)
+
+! acoustic mass matrix
+ if( ispec_is_acoustic(ispec) ) then
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ sngl( dble(jacobianl) * weight / dble(kappastore(i,j,k,ispec)) )
+ else
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ jacobianl * weight / kappastore(i,j,k,ispec)
+ endif
+ endif
+
+! elastic mass matrix
+ if( ispec_is_elastic(ispec) ) then
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass(iglob) = rmass(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(rhostore(i,j,k,ispec)) )
+ else
+ rmass(iglob) = rmass(iglob) + &
+ jacobianl * weight * rhostore(i,j,k,ispec)
+ endif
+ endif
+
+! poroelastic mass matrices
+ if( ispec_is_poroelastic(ispec) ) then
+
+ stop 'poroelastic mass matrices not implemented yet'
+
+ !rho_solid = density(1,kmato(ispec))
+ !rho_fluid = density(2,kmato(ispec))
+ !phi = porosity(kmato(ispec))
+ !tort = tortuosity(kmato(ispec))
+ !rho_bar = (1._CUSTOM_REAL-phil)*rhol_s + phil*rhol_f
+ !
+ !if(CUSTOM_REAL == SIZE_REAL) then
+ ! ! for the solid mass matrix
+ ! rmass_solid_poroelastic(iglob) = rmass_solid_poroelastic(iglob) + &
+ ! sngl( dble(jacobianl) * weight * dble(rho_bar - phi*rho_fluid/tort) )
+ !
+ ! ! for the fluid mass matrix
+ ! rmass_fluid_poroelastic(iglob) = rmass_fluid_poroelastic(iglob) + &
+ ! sngl( dble(jacobianl) * weight * dble(rho_bar*rho_fluid*tort - &
+ ! phi*rho_fluid*rho_fluid)/dble(rho_bar*phi) )
+ !else
+ ! rmass_solid_poroelastic(iglob) = rmass_solid_poroelastic(iglob) + &
+ ! jacobianl * weight * (rho_bar - phi*rho_fluid/tort)
+ !
+ ! rmass_fluid_poroelastic(iglob) = rmass_fluid_poroelastic(iglob) + &
+ ! jacobianl * weight * (rho_bar*rho_fluid*tort - &
+ ! phi*rho_fluid*rho_fluid) / (rho_bar*phi)
+ !endif
+ endif
+
+ enddo
+ enddo
+ enddo
+ enddo ! nspec
+
+ end subroutine create_mass_matrices
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine create_mass_matrices_ocean_load(nglob,nspec,ibool,OCEANS,&
+ UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION,NX_TOPO,NY_TOPO, &
+ ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO, &
+ itopo_bathy)
+
+! returns precomputed mass matrix in rmass array
+
+ use create_regions_mesh_ext_par
+ implicit none
+
+! number of spectral elements in each block
+ integer :: nspec
+ integer :: nglob
+
+! arrays with the mesh global indices
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ logical :: OCEANS
+
+! use integer array to store topography values
+ integer :: UTM_PROJECTION_ZONE
+ logical :: SUPPRESS_UTM_PROJECTION
+ integer :: NX_TOPO,NY_TOPO
+ double precision :: ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO
+ integer, dimension(NX_TOPO,NY_TOPO) :: itopo_bathy
+
+
+! local parameters
+ double precision :: weight
+ double precision :: xval,yval,long,lat,elevation
+ double precision :: height_oceans
+ double precision :: long_corner,lat_corner,ratio_xi,ratio_eta
+ integer :: ix_oceans,iy_oceans,iz_oceans,ispec_oceans,ispec2D,igll,iglobnum
+ integer :: icornerlong,icornerlat
+
+! creates ocean load mass matrix
+ if(OCEANS) then
+
+ ! adding ocean load mass matrix at ocean bottom
+ NGLOB_OCEAN = nglob
+ allocate(rmass_ocean_load(NGLOB_OCEAN))
+
+ ! create ocean load mass matrix for degrees of freedom at ocean bottom
+ rmass_ocean_load(:) = 0._CUSTOM_REAL
+
+ ! add contribution of the oceans for surface elements exactly at ocean bottom
+ do ispec2D = 1,num_free_surface_faces
+
+ ispec_oceans = free_surface_ispec(ispec2D)
+
+ ! only adds contribution if top boundary is elastic, no need to add this approximate calculation
+ ! if top is already acoustic/poroelastic
+ if( ispec_is_elastic(ispec_oceans) ) then
+
+ do igll=1,NGLLSQUARE
+ ix_oceans = free_surface_ijk(1,igll,ispec2D)
+ iy_oceans = free_surface_ijk(1,igll,ispec2D)
+ iz_oceans = free_surface_ijk(1,igll,ispec2D)
+
+ iglobnum=ibool(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
+
+ ! compute local height of oceans
+
+ ! get coordinates of current point
+ xval = xstore_dummy(iglobnum)
+ yval = ystore_dummy(iglobnum)
+
+ ! project x and y in UTM back to long/lat since topo file is in long/lat
+ call utm_geo(long,lat,xval,yval,UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
+
+ ! get coordinate of corner in bathy/topo model
+ icornerlong = int((long - ORIG_LONG_TOPO) / DEGREES_PER_CELL_TOPO) + 1
+ icornerlat = int((lat - ORIG_LAT_TOPO) / DEGREES_PER_CELL_TOPO) + 1
+
+ ! avoid edge effects and extend with identical point if outside model
+ if(icornerlong < 1) icornerlong = 1
+ if(icornerlong > NX_TOPO-1) icornerlong = NX_TOPO-1
+ if(icornerlat < 1) icornerlat = 1
+ if(icornerlat > NY_TOPO-1) icornerlat = NY_TOPO-1
+
+ ! compute coordinates of corner
+ long_corner = ORIG_LONG_TOPO + (icornerlong-1)*DEGREES_PER_CELL_TOPO
+ lat_corner = ORIG_LAT_TOPO + (icornerlat-1)*DEGREES_PER_CELL_TOPO
+
+ ! compute ratio for interpolation
+ ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO
+ ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO
+
+ ! avoid edge effects
+ if(ratio_xi < 0.) ratio_xi = 0.
+ if(ratio_xi > 1.) ratio_xi = 1.
+ if(ratio_eta < 0.) ratio_eta = 0.
+ if(ratio_eta > 1.) ratio_eta = 1.
+
+ ! interpolate elevation at current point
+ elevation = &
+ itopo_bathy(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
+ itopo_bathy(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
+ itopo_bathy(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
+ itopo_bathy(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
+
+ ! suppress positive elevation, which means no oceans
+ if(elevation >= - MINIMUM_THICKNESS_3D_OCEANS) then
+ height_oceans = 0.d0
+ else
+ height_oceans = dabs(elevation)
+ endif
+
+ ! take into account inertia of water column
+ weight = dble( free_surface_jacobian2Dw(igll,ispec2D)) &
+ * dble(RHO_OCEANS) * height_oceans
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_ocean_load(iglobnum) = rmass_ocean_load(iglobnum) + sngl(weight)
+ else
+ rmass_ocean_load(iglobnum) = rmass_ocean_load(iglobnum) + weight
+ endif
+
+ enddo ! igll
+ endif ! ispec_is_elastic
+ enddo ! num_free_surface_faces
+
+ ! add regular mass matrix to ocean load contribution
+ rmass_ocean_load(:) = rmass_ocean_load(:) + rmass(:)
+
+ else
+
+ ! allocate dummy array if no oceans
+ NGLOB_OCEAN = 1
+ allocate(rmass_ocean_load(NGLOB_OCEAN))
+
+ endif
+
+ end subroutine create_mass_matrices_ocean_load
+
+
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2010-02-25 21:58:34 UTC (rev 16342)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2010-02-25 23:07:29 UTC (rev 16343)
@@ -23,7 +23,6 @@
!
!=====================================================================
-
module create_regions_mesh_ext_par
include 'constants.h'
@@ -259,18 +258,7 @@
xstore,ystore,zstore,nspec, &
nodes_coords_ext_mesh,nnodes_ext_mesh,&
elmnts_ext_mesh,nelmnts_ext_mesh)
-
-! sets material velocities
- call sync_all()
- if( myrank == 0) then
- write(IMAIN,*) ' ...determining velocity model'
- endif
- call crm_ext_determine_velocity(nspec,&
- mat_ext_mesh,nelmnts_ext_mesh, &
- materials_ext_mesh,nmat_ext_mesh, &
- undef_mat_prop,nundefMat_ext_mesh, &
- ANISOTROPY)
-
+
! creates ibool index array for projection from local to global points
call sync_all()
if( myrank == 0) then
@@ -285,7 +273,7 @@
if( myrank == 0) then
write(IMAIN,*) ' ...preparing MPI interfaces '
endif
- call crm_ext_prepare_MPI(myrank,nglob,nspec,ibool, &
+ call get_MPI(myrank,nglob,nspec,ibool, &
nelmnts_ext_mesh,elmnts_ext_mesh, &
my_nelmnts_neighbours_ext_mesh, my_interfaces_ext_mesh, &
ibool_interfaces_ext_mesh, &
@@ -293,19 +281,22 @@
num_interfaces_ext_mesh,max_interface_size_ext_mesh,&
my_neighbours_ext_mesh,NPROC)
-! creates mass matrix
+! sets material velocities
call sync_all()
if( myrank == 0) then
- write(IMAIN,*) ' ...creating mass matrix '
+ write(IMAIN,*) ' ...determining velocity model'
endif
- call crm_ext_create_mass_matrix(nglob,nspec,ibool)
+ call get_model(nspec,ibool,mat_ext_mesh,nelmnts_ext_mesh, &
+ materials_ext_mesh,nmat_ext_mesh, &
+ undef_mat_prop,nundefMat_ext_mesh, &
+ ANISOTROPY)
! sets up absorbing/free surface boundaries
call sync_all()
if( myrank == 0) then
write(IMAIN,*) ' ...setting up absorbing boundaries '
endif
- call crm_ext_setup_abs_boundary(myrank,nspec,nglob,ibool, &
+ call get_absorbing_boundary(myrank,nspec,nglob,ibool, &
nodes_coords_ext_mesh,nnodes_ext_mesh, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
nodes_ibelm_xmin,nodes_ibelm_xmax,nodes_ibelm_ymin,nodes_ibelm_ymax, &
@@ -318,23 +309,29 @@
if( myrank == 0) then
write(IMAIN,*) ' ...detecting acoustic-elastic surfaces '
endif
- call crm_ext_detect_ac_el_surface(myrank, &
+ call get_coupling_surfaces(myrank, &
nspec,nglob,ibool,NPROC, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
num_interfaces_ext_mesh,max_interface_size_ext_mesh, &
my_neighbours_ext_mesh)
+! creates mass matrix
+ call sync_all()
+ if( myrank == 0) then
+ write(IMAIN,*) ' ...creating mass matrix '
+ endif
+ call create_mass_matrices(nglob,nspec,ibool)
+
! creates ocean load mass matrix
call sync_all()
if( myrank == 0) then
write(IMAIN,*) ' ...creating ocean load mass matrix '
endif
- call crm_ext_create_ocean_load_mass(nglob,nspec,ibool,OCEANS,&
+ call create_mass_matrices_ocean_load(nglob,nspec,ibool,OCEANS,&
UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION,NX_TOPO,NY_TOPO, &
ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO, &
itopo_bathy)
-
! saves the binary mesh files
call sync_all()
if( myrank == 0) then
@@ -681,368 +678,6 @@
!-------------------------------------------------------------------------------------------------
!
-subroutine crm_ext_determine_velocity(nspec,&
- mat_ext_mesh,nelmnts_ext_mesh, &
- materials_ext_mesh,nmat_ext_mesh, &
- undef_mat_prop,nundefMat_ext_mesh, &
- ANISOTROPY)
-
- use create_regions_mesh_ext_par
- implicit none
-
-! number of spectral elements in each block
- integer :: nspec
-
-! external mesh
- integer :: nelmnts_ext_mesh
- integer :: nmat_ext_mesh,nundefMat_ext_mesh
-
- integer, dimension(2,nelmnts_ext_mesh) :: mat_ext_mesh
- double precision, dimension(6,nmat_ext_mesh) :: materials_ext_mesh
- character (len=30), dimension(6,nundefMat_ext_mesh):: undef_mat_prop
-
-! anisotropy
- logical :: ANISOTROPY
-
-! local parameters
- real(kind=CUSTOM_REAL) :: vp,vs,rho
- real(kind=CUSTOM_REAL) :: c11,c12,c13,c14,c15,c16,c22,c23,c24,c25, &
- c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
-
- integer :: ispec,i,j,k,iundef,iflag_atten
- integer :: iflag,flag_below,flag_above
- integer :: iflag_aniso,idomain_id
-
-! ! Piero, read bedrock file
-! allocate(ibedrock(NX_TOPO_ANT,NY_TOPO_ANT))
-! if(myrank == 0) then
-! call read_bedrock_file(ibedrock)
-! ! write(IMAIN,*)
-! ! write(IMAIN,*) 'regional bedrock file read ranges in m from ',minval(ibedrock),' to ',maxval(ibedrock)
-! ! write(IMAIN,*)
-! endif
-! ! broadcast the information read on the master to the nodes
-! ! call MPI_BCAST(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT,MPI_REAL,0,MPI_COMM_WORLD,ier)
-! call bcast_all_cr(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT)
-
- ispec_is_acoustic(:) = .false.
- ispec_is_elastic(:) = .false.
- ispec_is_poroelastic(:) = .false.
-
-! material properties on all GLL points: taken from material values defined for
-! each spectral element in input mesh
- do ispec = 1, nspec
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
-! check if the material is known or unknown
- if (mat_ext_mesh(1,ispec) > 0) then
- ! density
- ! materials_ext_mesh format: #index1 = rho #index2 = vp #index3 = vs #index4 = Q_flag #index5 = 0
- !rhostore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(1,ispec))
- rho = materials_ext_mesh(1,mat_ext_mesh(1,ispec))
-
- ! isotropic values: vp, vs
- !vpstore(i,j,k,ispec) = materials_ext_mesh(2,mat_ext_mesh(1,ispec))
- !vsstore(i,j,k,ispec) = materials_ext_mesh(3,mat_ext_mesh(1,ispec))
- vp = materials_ext_mesh(2,mat_ext_mesh(1,ispec))
- vs = materials_ext_mesh(3,mat_ext_mesh(1,ispec))
-
- ! attenuation
- iflag_atten = materials_ext_mesh(4,mat_ext_mesh(1,ispec))
- !change for piero :
- !if(mat_ext_mesh(1,ispec) == 1) then
- ! iflag_attenuation_store(i,j,k,ispec) = 1
- !else
- ! iflag_attenuation_store(i,j,k,ispec) = 2
- !endif
-
- ! anisotropy
- iflag_aniso = materials_ext_mesh(5,mat_ext_mesh(1,ispec))
-
- ! material domain_id
- idomain_id = materials_ext_mesh(6,mat_ext_mesh(1,ispec))
-
- else if (mat_ext_mesh(2,ispec) == 1) then
- stop 'material: interface not implemented yet'
-
- do iundef = 1,nundefMat_ext_mesh
- if(trim(undef_mat_prop(2,iundef)) == 'interface') then
- read(undef_mat_prop(4,iundef),'(1i3)') flag_below
- read(undef_mat_prop(5,iundef),'(1i3)') flag_above
- endif
- enddo
- !call interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
- iflag = 1
- !rhostore(i,j,k,ispec) = materials_ext_mesh(1,iflag)
- !vpstore(i,j,k,ispec) = materials_ext_mesh(2,iflag)
- !vsstore(i,j,k,ispec) = materials_ext_mesh(3,iflag)
- rho = materials_ext_mesh(1,iflag)
- vp = materials_ext_mesh(2,iflag)
- vs = materials_ext_mesh(3,iflag)
- iflag_atten = materials_ext_mesh(4,iflag)
- !change for piero :
- ! if(iflag == 1) then
- ! iflag_attenuation_store(i,j,k,ispec) = 1
- ! else
- ! iflag_attenuation_store(i,j,k,ispec) = 2
- ! endif
- iflag_aniso = materials_ext_mesh(5,iflag)
- idomain_id = materials_ext_mesh(6,iflag)
- else
- stop 'material: tomography not implemented yet'
- ! call tomography()
- end if
-
-! adds anisotropic perturbation to vp, vs
- if( ANISOTROPY ) then
- call aniso_model(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
- c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
-
- c11store(i,j,k,ispec) = c11
- c12store(i,j,k,ispec) = c12
- c13store(i,j,k,ispec) = c13
- c14store(i,j,k,ispec) = c14
- c15store(i,j,k,ispec) = c15
- c16store(i,j,k,ispec) = c16
- c22store(i,j,k,ispec) = c22
- c23store(i,j,k,ispec) = c23
- c24store(i,j,k,ispec) = c24
- c25store(i,j,k,ispec) = c25
- c26store(i,j,k,ispec) = c26
- c33store(i,j,k,ispec) = c33
- c34store(i,j,k,ispec) = c34
- c35store(i,j,k,ispec) = c35
- c36store(i,j,k,ispec) = c36
- c44store(i,j,k,ispec) = c44
- c45store(i,j,k,ispec) = c45
- c46store(i,j,k,ispec) = c46
- c55store(i,j,k,ispec) = c55
- c56store(i,j,k,ispec) = c56
- c66store(i,j,k,ispec) = c66
-
- endif
-! density
- rhostore(i,j,k,ispec) = rho
-
-! kappa, mu
- kappastore(i,j,k,ispec) = rho*( vp*vp - FOUR_THIRDS*vs*vs )
- mustore(i,j,k,ispec) = rho*vs*vs
-
-! attenuation
- iflag_attenuation_store(i,j,k,ispec) = iflag_atten
- ! Stacey, a completer par la suite
- rho_vp(i,j,k,ispec) = rho*vp
- rho_vs(i,j,k,ispec) = rho*vs
- !end pll
-
-! material domain
- !print*,'velocity model:',ispec,idomain_id
- if( idomain_id == IDOMAIN_ACOUSTIC ) then
- ispec_is_acoustic(ispec) = .true.
- else if( idomain_id == IDOMAIN_ELASTIC ) then
- ispec_is_elastic(ispec) = .true.
- else if( idomain_id == IDOMAIN_POROELASTIC ) then
- stop 'poroelastic material domain not implemented yet'
- ispec_is_poroelastic(ispec) = .true.
- else
- stop 'error material domain index'
- endif
- enddo
- enddo
- enddo
- !print*,myrank,'ispec:',ispec,'rho:',rhostore(1,1,1,ispec),'vp:',vpstore(1,1,1,ispec),'vs:',vsstore(1,1,1,ispec)
- enddo
-
-! checks material domains
- do ispec=1,nspec
- if( (ispec_is_acoustic(ispec) .eqv. .false.) &
- .and. (ispec_is_elastic(ispec) .eqv. .false.) &
- .and. (ispec_is_poroelastic(ispec) .eqv. .false.) ) then
- print*,'error material domain not assigned to element:',ispec
- print*,'acoustic: ',ispec_is_acoustic(ispec)
- print*,'elastic: ',ispec_is_elastic(ispec)
- print*,'poroelastic: ',ispec_is_poroelastic(ispec)
- stop 'error material domain index element'
- endif
- enddo
-
-
-! !! DK DK store the position of the six stations to be able to
-! !! DK DK exclude circles around each station to make sure they are on the bedrock
-! !! DK DK and not in the ice
-! utm_x_station(1) = 783500.6250000d0
-! utm_y_station(1) = -11828.7519531d0
-
-! utm_x_station(2) = 853644.5000000d0
-! utm_y_station(2) = -114.0138092d0
-
-! utm_x_station(3) = 863406.0000000d0
-! utm_y_station(3) = -53736.1640625d0
-
-! utm_x_station(4) = 823398.8125000d0
-! utm_y_station(4) = 29847.4511719d0
-
-! utm_x_station(5) = 863545.3750000d0
-! utm_y_station(5) = 19669.6621094d0
-
-! utm_x_station(6) = 817099.3750000d0
-! utm_y_station(6) = -24430.2871094d0
-
-! print*,myrank,'après store the position of the six stations'
-! call flush(6)
-
-! print*, myrank,minval(nodes_coords_ext_mesh(1,:))
-! call flush(6)
-
-
-! print*, myrank,maxval(nodes_coords_ext_mesh(1,:))
-! call flush(6)
-
-
-! do ispec = 1, nspec
-
-! zmesh = zstore(2,2,2,ispec)
-
-! ! if(doubling_index == IFLAG_ONE_LAYER_TOPOGRAPHY) then
-! if(any(ibelm_top == ispec)) then
-! doubling_value_found_for_Piero = IFLAG_ONE_LAYER_TOPOGRAPHY
-
-! else if(zmesh < Z_23p4km) then
-! doubling_value_found_for_Piero = IFLAG_MANTLE_BELOW_23p4km
-
-! else if(zmesh < Z_14km) then
-! doubling_value_found_for_Piero = IFLAG_14km_to_23p4km
-
-! else
-! doubling_value_found_for_Piero = IFLAG_BEDROCK_down_to_14km
-! endif
-! idoubling(ispec) = doubling_value_found_for_Piero
-
-! do k = 1, NGLLZ
-! do j = 1, NGLLY
-! do i = 1, NGLLX
-
-
-! if(idoubling(ispec) == IFLAG_ONE_LAYER_TOPOGRAPHY .or. &
-! idoubling(ispec) == IFLAG_BEDROCK_down_to_14km) then
-
-! ! since we have suppressed UTM projection for Piero Basini, UTMx is the same as long
-! ! and UTMy is the same as lat
-! long = xstore(i,j,k,ispec)
-! lat = ystore(i,j,k,ispec)
-
-! ! get coordinate of corner in model
-! icornerlong = int((long - ORIG_LONG_TOPO) / DEGREES_PER_CELL_TOPO) + 1
-! icornerlat = int((lat - ORIG_LAT_TOPO) / DEGREES_PER_CELL_TOPO) + 1
-
-! ! avoid edge effects and extend with identical point if outside model
-! if(icornerlong < 1) icornerlong = 1
-! if(icornerlong > NX_TOPO-1) icornerlong = NX_TOPO-1
-! if(icornerlat < 1) icornerlat = 1
-! if(icornerlat > NY_TOPO-1) icornerlat = NY_TOPO-1
-
-! ! compute coordinates of corner
-! long_corner = ORIG_LONG_TOPO + (icornerlong-1)*DEGREES_PER_CELL_TOPO
-! lat_corner = ORIG_LAT_TOPO + (icornerlat-1)*DEGREES_PER_CELL_TOPO
-
-! ! compute ratio for interpolation
-! ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO
-! ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO
-
-! ! avoid edge effects
-! if(ratio_xi < 0.) ratio_xi = 0.
-! if(ratio_xi > 1.) ratio_xi = 1.
-! if(ratio_eta < 0.) ratio_eta = 0.
-! if(ratio_eta > 1.) ratio_eta = 1.
-
-! ! interpolate elevation at current point
-! elevation_bedrock = &
-! ibedrock(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
-! ibedrock(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
-! ibedrock(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
-! ibedrock(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
-
-! !! DK DK exclude circles around each station to make sure they are on the bedrock
-! !! DK DK and not in the ice
-! is_around_a_station = .false.
-! do istation = 1,NUMBER_OF_STATIONS
-! if(sqrt((long - utm_x_station(istation))**2 + (lat - utm_y_station(istation))**2) < RADIUS_TO_EXCLUDE) then
-! is_around_a_station = .true.
-! exit
-! endif
-! enddo
-
-! ! define elastic parameters in the model
-
-! ! we are above the bedrock interface i.e. in the ice, and not too close to a station
-! if(zmesh >= elevation_bedrock .and. .not. is_around_a_station) then
-! vp = 3800.d0
-! vs = 1900.d0
-! rho = 900.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_ICE
-
-! ! we are below the bedrock interface i.e. in the bedrock, or close to a station
-! else
-! vp = 5800.d0
-! vs = 3200.d0
-! rho = 2600.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
-! endif
-
-! else if(idoubling(ispec) == IFLAG_14km_to_23p4km) then
-! vp = 6800.d0
-! vs = 3900.d0
-! rho = 2900.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
-
-! else if(idoubling(ispec) == IFLAG_MANTLE_BELOW_23p4km) then
-! vp = 8100.d0
-! vs = 4480.d0
-! rho = 3380.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
-
-! endif
-
-! !pll 8/06
-! if(CUSTOM_REAL == SIZE_REAL) then
-! rhostore(i,j,k,ispec) = sngl(rho)
-! vpstore(i,j,k,ispec) = sngl(vp)
-! vsstore(i,j,k,ispec) = sngl(vs)
-! else
-! rhostore(i,j,k,ispec) = rho
-! vpstore(i,j,k,ispec) = vp
-! vsstore(i,j,k,ispec) = vs
-! end if
-
-! kappastore(i,j,k,ispec) = rhostore(i,j,k,ispec)*(vpstore(i,j,k,ispec)*vpstore(i,j,k,ispec) - &
-! 4.d0*vsstore(i,j,k,ispec)*vsstore(i,j,k,ispec)/3.d0)
-! mustore(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)*&
-! vsstore(i,j,k,ispec)
-
-! ! Stacey, a completer par la suite
-! rho_vp(i,j,k,ispec) = rhostore(i,j,k,ispec)*vpstore(i,j,k,ispec)
-! rho_vs(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)
-! !end pll
-
-! ! kappastore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
-! ! (materials_ext_mesh(2,mat_ext_mesh(ispec))*materials_ext_mesh(2,mat_ext_mesh(ispec)) - &
-! ! 4.d0*materials_ext_mesh(3,mat_ext_mesh(ispec))*materials_ext_mesh(3,mat_ext_mesh(ispec))/3.d0)
-! ! mustore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
-! materials_ext_mesh(3,mat_ext_mesh(ispec))*&
-! ! x materials_ext_mesh(3,mat_ext_mesh(ispec))
-! enddo
-! enddo
-! enddo
-! enddo
-
-end subroutine crm_ext_determine_velocity
-
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
subroutine crm_ext_setup_indexing(ibool, &
xstore,ystore,zstore,nspec,nglob,npointot, &
nnodes_ext_mesh,nodes_coords_ext_mesh,myrank)
@@ -1140,1568 +775,40 @@
!-------------------------------------------------------------------------------------------------
!
-subroutine crm_ext_prepare_MPI(myrank,nglob,nspec,ibool, &
- nelmnts_ext_mesh,elmnts_ext_mesh, &
- my_nelmnts_neighbours_ext_mesh, my_interfaces_ext_mesh, &
- ibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh, &
- num_interfaces_ext_mesh,max_interface_size_ext_mesh, &
- my_neighbours_ext_mesh,NPROC)
-
-! sets up the MPI interface for communication between partitions
-
- use create_regions_mesh_ext_par
- implicit none
-
- integer :: myrank,nglob,nspec,NPROC
-
-! global indexing
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-
-! external mesh, element indexing
- integer :: nelmnts_ext_mesh
- integer, dimension(ESIZE,nelmnts_ext_mesh) :: elmnts_ext_mesh
-
- integer :: num_interfaces_ext_mesh,max_interface_size_ext_mesh
-
- integer, dimension(num_interfaces_ext_mesh) :: my_nelmnts_neighbours_ext_mesh
- integer, dimension(6,max_interface_size_ext_mesh,num_interfaces_ext_mesh) :: my_interfaces_ext_mesh
-
- integer, dimension(num_interfaces_ext_mesh) :: my_neighbours_ext_mesh
-
- integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh
- integer, dimension(NGLLX*NGLLX*max_interface_size_ext_mesh,num_interfaces_ext_mesh) :: ibool_interfaces_ext_mesh
-
-
- !integer :: nnodes_ext_mesh
- !double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
-
-!local parameters
- double precision, dimension(:), allocatable :: xp,yp,zp
- double precision, dimension(:), allocatable :: work_ext_mesh
-
- integer, dimension(:), allocatable :: locval
- integer, dimension(:), allocatable :: nibool_interfaces_ext_mesh_true
-
- ! for MPI buffers
- integer, dimension(:), allocatable :: reorder_interface_ext_mesh,ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh
- integer, dimension(:), allocatable :: ibool_interface_ext_mesh_dummy
- logical, dimension(:), allocatable :: ifseg
- integer :: iinterface,ilocnum
- integer :: num_points1, num_points2
-
- ! assembly test
- integer :: i,j,k,ispec,iglob,count,inum
- integer :: max_nibool_interfaces_ext_mesh
- integer,dimension(:),allocatable :: test_flag
- real(kind=CUSTOM_REAL), dimension(:),allocatable :: test_flag_cr
- integer, dimension(:,:), allocatable :: ibool_interfaces_dummy
-
-! gets global indices for points on MPI interfaces (defined by my_interfaces_ext_mesh) between different partitions
-! and stores them in ibool_interfaces_ext_mesh & nibool_interfaces_ext_mesh (number of total points)
- call prepare_assemble_MPI( nelmnts_ext_mesh,elmnts_ext_mesh, &
- ibool,nglob,ESIZE, &
- num_interfaces_ext_mesh, max_interface_size_ext_mesh, &
- my_nelmnts_neighbours_ext_mesh, my_interfaces_ext_mesh, &
- ibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh )
-
- allocate(nibool_interfaces_ext_mesh_true(num_interfaces_ext_mesh))
-
-! sorts ibool comm buffers lexicographically for all MPI interfaces
- num_points1 = 0
- num_points2 = 0
- do iinterface = 1, num_interfaces_ext_mesh
-
- allocate(xp(nibool_interfaces_ext_mesh(iinterface)))
- allocate(yp(nibool_interfaces_ext_mesh(iinterface)))
- allocate(zp(nibool_interfaces_ext_mesh(iinterface)))
- allocate(locval(nibool_interfaces_ext_mesh(iinterface)))
- allocate(ifseg(nibool_interfaces_ext_mesh(iinterface)))
- allocate(reorder_interface_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
- allocate(ibool_interface_ext_mesh_dummy(nibool_interfaces_ext_mesh(iinterface)))
- allocate(ind_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
- allocate(ninseg_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
- allocate(iwork_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
- allocate(work_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
-
- ! gets x,y,z coordinates of global points on MPI interface
- do ilocnum = 1, nibool_interfaces_ext_mesh(iinterface)
- xp(ilocnum) = xstore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
- yp(ilocnum) = ystore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
- zp(ilocnum) = zstore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
- enddo
-
- ! sorts (lexicographically?) ibool_interfaces_ext_mesh and updates value
- ! of total number of points nibool_interfaces_ext_mesh_true(iinterface)
- call sort_array_coordinates(nibool_interfaces_ext_mesh(iinterface),xp,yp,zp, &
- ibool_interfaces_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
- reorder_interface_ext_mesh,locval,ifseg,nibool_interfaces_ext_mesh_true(iinterface), &
- ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh,work_ext_mesh)
-
- ! checks that number of MPI points are still the same
- num_points1 = num_points1 + nibool_interfaces_ext_mesh(iinterface)
- num_points2 = num_points2 + nibool_interfaces_ext_mesh_true(iinterface)
- if( num_points1 /= num_points2 ) then
- write(*,*) 'error sorting MPI interface points:',myrank
- write(*,*) ' interface:',iinterface,num_points1,num_points2
- call exit_mpi(myrank,'error sorting MPI interface')
- endif
- !write(*,*) myrank,'intfc',iinterface,num_points2,nibool_interfaces_ext_mesh_true(iinterface)
-
- ! cleanup temporary arrays
- deallocate(xp)
- deallocate(yp)
- deallocate(zp)
- deallocate(locval)
- deallocate(ifseg)
- deallocate(reorder_interface_ext_mesh)
- deallocate(ibool_interface_ext_mesh_dummy)
- deallocate(ind_ext_mesh)
- deallocate(ninseg_ext_mesh)
- deallocate(iwork_ext_mesh)
- deallocate(work_ext_mesh)
-
- enddo
-
- ! cleanup
- deallocate(nibool_interfaces_ext_mesh_true)
-
- ! outputs total number of MPI interface points
- call sum_all_i(num_points2,ilocnum)
- if( myrank == 0 ) then
- write(IMAIN,*) ' total MPI interface points: ',ilocnum
- endif
-
-! checks with assembly of test fields
- allocate(test_flag(nglob),test_flag_cr(nglob))
- test_flag(:) = 0
- test_flag_cr(:) = 0._CUSTOM_REAL
- count = 0
- do ispec = 1, nspec
- ! sets flags on global points
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- ! global index
- iglob = ibool(i,j,k,ispec)
-
- ! counts number of unique global points to set
- if( test_flag(iglob) == 0 ) count = count+1
-
- ! sets identifier
- test_flag(iglob) = myrank + 1
- test_flag_cr(iglob) = myrank + 1.0
- enddo
- enddo
- enddo
- enddo
- call sync_all()
-
- ! collects contributions from different MPI partitions
- ! sets up MPI communications
- max_nibool_interfaces_ext_mesh = maxval( nibool_interfaces_ext_mesh(:) )
- allocate(ibool_interfaces_dummy(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
-
- count = 0
- do iinterface = 1, num_interfaces_ext_mesh
- ibool_interfaces_dummy(:,iinterface) = ibool_interfaces_ext_mesh(1:max_nibool_interfaces_ext_mesh,iinterface)
- count = count + nibool_interfaces_ext_mesh(iinterface)
- !write(*,*) myrank,'interfaces ',iinterface,nibool_interfaces_ext_mesh(iinterface),max_nibool_interfaces_ext_mesh
- enddo
- call sync_all()
-
- call sum_all_i(count,iglob)
- if( myrank == 0 ) then
- if( iglob /= ilocnum ) call exit_mpi(myrank,'error total global MPI interface points')
- endif
-
- ! adds contributions from different partitions to flag arrays
- ! integer arrays
- call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,test_flag, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_dummy,&
- my_neighbours_ext_mesh)
- ! custom_real arrays
- call assemble_MPI_scalar_ext_mesh(NPROC,nglob,test_flag_cr, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_dummy, &
- my_neighbours_ext_mesh)
-
- ! checks number of interface points
- i = 0
- j = 0
- do iglob=1,nglob
- ! only counts flags with MPI contributions
- if( test_flag(iglob) > myrank+1 ) i = i + 1
- if( test_flag_cr(iglob) > myrank+1.0) j = j + 1
- enddo
- call sum_all_i(i,inum)
- call sum_all_i(j,iglob)
- if( myrank == 0 ) then
- write(IMAIN,*) ' total assembled MPI interface points:',inum
- if( inum /= iglob .or. inum > ilocnum ) call exit_mpi(myrank,'error MPI assembly')
- endif
-
-end subroutine crm_ext_prepare_MPI
-
-
!
-!-------------------------------------------------------------------------------------------------
!
-
-subroutine crm_ext_create_mass_matrix(nglob,nspec,ibool)
-
-! returns precomputed mass matrix in rmass array
-
- use create_regions_mesh_ext_par
- implicit none
-
-! number of spectral elements in each block
- integer :: nspec
- integer :: nglob
-
-! arrays with the mesh global indices
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-
-! local parameters
- double precision :: weight
- real(kind=CUSTOM_REAL) :: jacobianl
- integer :: ispec,i,j,k,iglob,ier
-
-! allocates memory
- allocate(rmass(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
- allocate(rmass_acoustic(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
- allocate(rmass_solid_poroelastic(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
- allocate(rmass_fluid_poroelastic(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
-
-! creates mass matrix
- rmass(:) = 0._CUSTOM_REAL
- rmass_acoustic(:) = 0._CUSTOM_REAL
- rmass_solid_poroelastic(:) = 0._CUSTOM_REAL
- rmass_fluid_poroelastic(:) = 0._CUSTOM_REAL
-
- do ispec=1,nspec
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
-
- weight = wxgll(i)*wygll(j)*wzgll(k)
- jacobianl = jacobianstore(i,j,k,ispec)
-
-! acoustic mass matrix
- if( ispec_is_acoustic(ispec) ) then
- ! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
- sngl( dble(jacobianl) * weight / dble(kappastore(i,j,k,ispec)) )
- else
- rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
- jacobianl * weight / kappastore(i,j,k,ispec)
- endif
- endif
-
-! elastic mass matrix
- if( ispec_is_elastic(ispec) ) then
- if(CUSTOM_REAL == SIZE_REAL) then
- rmass(iglob) = rmass(iglob) + &
- sngl( dble(jacobianl) * weight * dble(rhostore(i,j,k,ispec)) )
- else
- rmass(iglob) = rmass(iglob) + &
- jacobianl * weight * rhostore(i,j,k,ispec)
- endif
- endif
-
-! poroelastic mass matrices
- if( ispec_is_poroelastic(ispec) ) then
-
- stop 'poroelastic mass matrices not implemented yet'
-
- !rho_solid = density(1,kmato(ispec))
- !rho_fluid = density(2,kmato(ispec))
- !phi = porosity(kmato(ispec))
- !tort = tortuosity(kmato(ispec))
- !rho_bar = (1._CUSTOM_REAL-phil)*rhol_s + phil*rhol_f
- !
- !if(CUSTOM_REAL == SIZE_REAL) then
- ! ! for the solid mass matrix
- ! rmass_solid_poroelastic(iglob) = rmass_solid_poroelastic(iglob) + &
- ! sngl( dble(jacobianl) * weight * dble(rho_bar - phi*rho_fluid/tort) )
- !
- ! ! for the fluid mass matrix
- ! rmass_fluid_poroelastic(iglob) = rmass_fluid_poroelastic(iglob) + &
- ! sngl( dble(jacobianl) * weight * dble(rho_bar*rho_fluid*tort - &
- ! phi*rho_fluid*rho_fluid)/dble(rho_bar*phi) )
- !else
- ! rmass_solid_poroelastic(iglob) = rmass_solid_poroelastic(iglob) + &
- ! jacobianl * weight * (rho_bar - phi*rho_fluid/tort)
- !
- ! rmass_fluid_poroelastic(iglob) = rmass_fluid_poroelastic(iglob) + &
- ! jacobianl * weight * (rho_bar*rho_fluid*tort - &
- ! phi*rho_fluid*rho_fluid) / (rho_bar*phi)
- !endif
- endif
-
- enddo
- enddo
- enddo
- enddo ! nspec
-
-end subroutine crm_ext_create_mass_matrix
-
-
+!subroutine bubble_sort( arr, ndim )
!
-!-------------------------------------------------------------------------------------------------
+!! sorts values in array arr[ndim] in increasing order
!
+! implicit none
+!
+! integer :: ndim
+! integer :: arr(ndim)
+!
+! logical :: swapped
+! integer :: j,tmp
+!
+! swapped = .true.
+! do while( swapped )
+! swapped = .false.
+! do j = 1, ndim-1
+! if( arr(j+1) < arr(j) ) then
+! tmp = arr(j)
+! arr(j) = arr(j+1)
+! arr(j+1) = tmp
+! swapped = .true.
+! endif
+! enddo
+! enddo
+!
+!end subroutine bubble_sort
-subroutine crm_ext_setup_abs_boundary(myrank,nspec,nglob,ibool, &
- nodes_coords_ext_mesh,nnodes_ext_mesh, &
- ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
- nodes_ibelm_xmin,nodes_ibelm_xmax,nodes_ibelm_ymin,nodes_ibelm_ymax, &
- nodes_ibelm_bottom,nodes_ibelm_top, &
- nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
- nspec2D_bottom,nspec2D_top)
-
-! determines absorbing boundaries/free-surface, 2D jacobians, face normals for Stacey conditions
-
- use create_regions_mesh_ext_par
- implicit none
-
-! number of spectral elements in each block
- integer :: myrank,nspec,nglob
-
-! arrays with the mesh
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-
-! data from the external mesh
- integer :: nnodes_ext_mesh
- double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
-
-! absorbing boundaries (as defined in CUBIT)
- integer :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, NSPEC2D_BOTTOM, NSPEC2D_TOP
- ! element indices containing a boundary
- integer, dimension(nspec2D_xmin) :: ibelm_xmin
- integer, dimension(nspec2D_xmax) :: ibelm_xmax
- integer, dimension(nspec2D_ymin) :: ibelm_ymin
- integer, dimension(nspec2D_ymax) :: ibelm_ymax
- integer, dimension(NSPEC2D_BOTTOM) :: ibelm_bottom
- integer, dimension(NSPEC2D_TOP) :: ibelm_top
-
- ! corner node indices of boundary faces coming from CUBIT
- integer, dimension(4,nspec2D_xmin) :: nodes_ibelm_xmin
- integer, dimension(4,nspec2D_xmax) :: nodes_ibelm_xmax
- integer, dimension(4,nspec2D_ymin) :: nodes_ibelm_ymin
- integer, dimension(4,nspec2D_ymax) :: nodes_ibelm_ymax
- integer, dimension(4,NSPEC2D_BOTTOM) :: nodes_ibelm_bottom
- integer, dimension(4,NSPEC2D_TOP) :: nodes_ibelm_top
-
-! local parameters
- logical, dimension(:,:),allocatable :: iboun ! pll
-
- ! (assumes NGLLX=NGLLY=NGLLZ)
- real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
- real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
- integer:: ijk_face(3,NGLLX,NGLLY)
-
- ! corner locations for faces
- real(kind=CUSTOM_REAL), dimension(:,:,:),allocatable :: xcoord_iboun,ycoord_iboun,zcoord_iboun
-
- ! face corner locations
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
- integer :: ispec,ispec2D,icorner,ier,iabs,iface,igll,i,j,igllfree,ifree
-
-! allocate temporary flag array
- allocate(iboun(6,nspec), &
- xcoord_iboun(NGNOD2D,6,nspec), &
- ycoord_iboun(NGNOD2D,6,nspec), &
- zcoord_iboun(NGNOD2D,6,nspec),stat=ier)
- if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
-
-! sets flag in array iboun for elements with an absorbing boundary faces
- iboun(:,:) = .false.
-
-! abs face counter
- iabs = 0
-
- ! xmin
- do ispec2D = 1, nspec2D_xmin
- ! sets element
- ispec = ibelm_xmin(ispec2D)
-
- !if(myrank == 0 ) print*,'xmin:',ispec2D,ispec
-
- ! looks for i,j,k indices of GLL points on boundary face
- ! determines element face by given CUBIT corners
- do icorner=1,NGNOD2D
- xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_xmin(icorner,ispec2D))
- ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_xmin(icorner,ispec2D))
- zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_xmin(icorner,ispec2D))
- !print*,'corner look:',icorner,xcoord(icorner),ycoord(icorner),zcoord(icorner)
- enddo
-
- ! sets face id of reference element associated with this face
- call get_element_face_id(ispec,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- iface)
-
- iboun(iface,ispec) = .true.
-
- ! ijk indices of GLL points for face id
- call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ)
-
- ! weighted jacobian and normal
- call get_jacobian_boundary_face(myrank,nspec, &
- xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
-
- ! normal convention: points away from element
- ! switch normal direction if necessary
- do j=1,NGLLZ
- do i=1,NGLLX
- call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal_face(:,i,j) )
- enddo
- enddo
-
- ! sets face infos
- iabs = iabs + 1
- abs_boundary_ispec(iabs) = ispec
-
- ! gll points -- assuming NGLLX = NGLLY = NGLLZ
- igll = 0
- do j=1,NGLLZ
- do i=1,NGLLX
- igll = igll+1
- abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
- abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
- abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
- enddo
- enddo
-
- enddo ! nspec2D_xmin
-
- ! xmax
- do ispec2D = 1, nspec2D_xmax
- ! sets element
- ispec = ibelm_xmax(ispec2D)
-
- ! looks for i,j,k indices of GLL points on boundary face
- ! determines element face by given CUBIT corners
- do icorner=1,NGNOD2D
- xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_xmax(icorner,ispec2D))
- ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_xmax(icorner,ispec2D))
- zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_xmax(icorner,ispec2D))
- enddo
-
- ! sets face id of reference element associated with this face
- call get_element_face_id(ispec,xcoord,ycoord,zcoord,&
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- iface )
- iboun(iface,ispec) = .true.
-
- ! ijk indices of GLL points on face
- call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ)
-
- ! weighted jacobian and normal
- call get_jacobian_boundary_face(myrank,nspec, &
- xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
-
- ! normal convention: points away from element
- ! switch normal direction if necessary
- do j=1,NGLLZ
- do i=1,NGLLX
- call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal_face(:,i,j) )
- enddo
- enddo
-
- ! sets face infos
- iabs = iabs + 1
- abs_boundary_ispec(iabs) = ispec
-
- ! gll points -- assuming NGLLX = NGLLY = NGLLZ
- igll = 0
- do j=1,NGLLZ
- do i=1,NGLLX
- igll = igll+1
- abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
- abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
- abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
- enddo
- enddo
-
- enddo
-
- ! ymin
- do ispec2D = 1, nspec2D_ymin
- ! sets element
- ispec = ibelm_ymin(ispec2D)
-
- ! looks for i,j,k indices of GLL points on boundary face
- ! determines element face by given CUBIT corners
- do icorner=1,NGNOD2D
- xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_ymin(icorner,ispec2D))
- ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_ymin(icorner,ispec2D))
- zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_ymin(icorner,ispec2D))
- enddo
-
- ! sets face id of reference element associated with this face
- call get_element_face_id(ispec,xcoord,ycoord,zcoord,&
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- iface )
- iboun(iface,ispec) = .true.
-
- ! ijk indices of GLL points on face
- call get_element_face_gll_indices(iface,ijk_face,NGLLY,NGLLZ)
-
- ! weighted jacobian and normal
- call get_jacobian_boundary_face(myrank,nspec, &
- xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
- ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ)
-
- ! normal convention: points away from element
- ! switch normal direction if necessary
- do j=1,NGLLZ
- do i=1,NGLLY
- call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal_face(:,i,j) )
- enddo
- enddo
-
- ! sets face infos
- iabs = iabs + 1
- abs_boundary_ispec(iabs) = ispec
-
- ! gll points -- assuming NGLLX = NGLLY = NGLLZ
- igll = 0
- do j=1,NGLLZ
- do i=1,NGLLY
- igll = igll+1
- abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
- abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
- abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
- enddo
- enddo
-
- enddo
-
- ! ymax
- do ispec2D = 1, nspec2D_ymax
- ! sets element
- ispec = ibelm_ymax(ispec2D)
-
- ! looks for i,j,k indices of GLL points on boundary face
- ! determines element face by given CUBIT corners
- do icorner=1,NGNOD2D
- xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_ymax(icorner,ispec2D))
- ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_ymax(icorner,ispec2D))
- zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_ymax(icorner,ispec2D))
- enddo
-
- ! sets face id of reference element associated with this face
- call get_element_face_id(ispec,xcoord,ycoord,zcoord,&
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- iface )
- iboun(iface,ispec) = .true.
-
- ! ijk indices of GLL points on face
- call get_element_face_gll_indices(iface,ijk_face,NGLLY,NGLLZ)
-
- ! weighted jacobian and normal
- call get_jacobian_boundary_face(myrank,nspec, &
- xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ)
-
- ! normal convention: points away from element
- ! switch normal direction if necessary
- do j=1,NGLLZ
- do i=1,NGLLY
- call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal_face(:,i,j) )
- enddo
- enddo
-
- ! sets face infos
- iabs = iabs + 1
- abs_boundary_ispec(iabs) = ispec
-
- ! gll points -- assuming NGLLX = NGLLY = NGLLZ
- igll = 0
- do j=1,NGLLY
- do i=1,NGLLX
- igll = igll+1
- abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
- abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
- abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
- enddo
- enddo
-
- enddo
-
- ! bottom
- do ispec2D = 1, NSPEC2D_BOTTOM
- ! sets element
- ispec = ibelm_bottom(ispec2D)
-
- ! looks for i,j,k indices of GLL points on boundary face
- ! determines element face by given CUBIT corners
- do icorner=1,NGNOD2D
- xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_bottom(icorner,ispec2D))
- ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_bottom(icorner,ispec2D))
- zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_bottom(icorner,ispec2D))
- enddo
-
- ! sets face id of reference element associated with this face
- call get_element_face_id(ispec,xcoord,ycoord,zcoord,&
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- iface )
- iboun(iface,ispec) = .true.
-
- ! ijk indices of GLL points on face
- call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY)
-
- ! weighted jacobian and normal
- call get_jacobian_boundary_face(myrank,nspec, &
- xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
-
- ! normal convention: points away from element
- ! switch normal direction if necessary
- do j=1,NGLLY
- do i=1,NGLLX
- call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal_face(:,i,j) )
- enddo
- enddo
-
- ! sets face infos
- iabs = iabs + 1
- abs_boundary_ispec(iabs) = ispec
-
- ! gll points -- assuming NGLLX = NGLLY = NGLLZ
- igll = 0
- do j=1,NGLLY
- do i=1,NGLLX
- igll = igll+1
- abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
- abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
- abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
- enddo
- enddo
-
- enddo
-
- ! top
- ! free surface face counter
- ifree = 0
- do ispec2D = 1, NSPEC2D_TOP
- ! sets element
- ispec = ibelm_top(ispec2D)
-
- ! looks for i,j,k indices of GLL points on boundary face
- ! determines element face by given CUBIT corners
- do icorner=1,NGNOD2D
- xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_top(icorner,ispec2D))
- ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_top(icorner,ispec2D))
- zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_top(icorner,ispec2D))
- enddo
-
- ! sets face id of reference element associated with this face
- call get_element_face_id(ispec,xcoord,ycoord,zcoord,&
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- iface )
- iboun(iface,ispec) = .true.
-
- ! ijk indices of GLL points on face
- call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY)
-
- ! weighted jacobian and normal
- call get_jacobian_boundary_face(myrank,nspec, &
- xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
-
- ! normal convention: points away from element
- ! switch normal direction if necessary
- do j=1,NGLLY
- do i=1,NGLLX
- call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal_face(:,i,j) )
- enddo
- enddo
-
- ! stores surface infos
- if( .not. ABSORB_FREE_SURFACE ) then
- ! store for free surface
- !jacobian2D_top(:,:,ispec2D) = jacobian2Dw_face(:,:)
- !normal_top(:,:,:,ispec2D) = normal_face(:,:,:)
-
- ! sets face infos
- ifree = ifree + 1
- free_surface_ispec(ifree) = ispec
-
- ! gll points -- assuming NGLLX = NGLLY = NGLLZ
- igllfree = 0
- do j=1,NGLLY
- do i=1,NGLLX
- igllfree = igllfree+1
- free_surface_ijk(:,igllfree,ifree) = ijk_face(:,i,j)
- free_surface_jacobian2Dw(igllfree,ifree) = jacobian2Dw_face(i,j)
- free_surface_normal(:,igllfree,ifree) = normal_face(:,i,j)
- enddo
- enddo
- else
- ! adds face infos to absorbing boundary surface
- iabs = iabs + 1
- abs_boundary_ispec(iabs) = ispec
-
- ! gll points -- assuming NGLLX = NGLLY = NGLLZ
- igll = 0
- do j=1,NGLLY
- do i=1,NGLLX
- igll = igll+1
- abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
- abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
- abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
- enddo
- enddo
-
- ! resets free surface
- ifree = 1
- free_surface_ispec(:) = 0
- free_surface_ijk(:,:,:) = 0
- free_surface_jacobian2Dw(:,:) = 0.0
- free_surface_normal(:,:,:) = 0.0
- endif
- enddo
-
-! checks counters
- if( ifree /= num_free_surface_faces ) then
- print*,'error number of free surface faces:',ifree,num_free_surface_faces
- stop 'error number of free surface faces'
- endif
-
- if( iabs /= num_abs_boundary_faces ) then
- print*,'error number of absorbing faces:',iabs,num_abs_boundary_faces
- stop 'error number of absorbing faces'
- endif
-
- call sum_all_i(num_abs_boundary_faces,iabs)
- if( myrank == 0 ) then
- write(IMAIN,*) ' absorbing boundary:'
- write(IMAIN,*) ' total number of faces = ',iabs
- if( ABSORB_FREE_SURFACE ) then
- write(IMAIN,*) ' absorbing boundary includes free surface'
- endif
- endif
-
-end subroutine crm_ext_setup_abs_boundary
-
-
!
!-------------------------------------------------------------------------------------------------
!
-subroutine crm_ext_detect_ac_el_surface(myrank, &
- nspec,nglob,ibool,NPROC, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
- num_interfaces_ext_mesh,max_interface_size_ext_mesh, &
- my_neighbours_ext_mesh)
-
-! determines coupling surface for acoustic-elastic domains
- use create_regions_mesh_ext_par
- implicit none
-
-! number of spectral elements in each block
- integer :: myrank,nspec,nglob,NPROC
-
-! arrays with the mesh
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-
-! MPI communication
- integer :: num_interfaces_ext_mesh,max_interface_size_ext_mesh
- integer, dimension(num_interfaces_ext_mesh) :: my_neighbours_ext_mesh
- integer, dimension(NGLLX*NGLLX*max_interface_size_ext_mesh,num_interfaces_ext_mesh) :: &
- ibool_interfaces_ext_mesh
- integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh
-
-! local parameters
- ! (assumes NGLLX=NGLLY=NGLLZ)
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
- real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
- real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
- real(kind=CUSTOM_REAL),dimension(:,:,:),allocatable :: tmp_normal
- real(kind=CUSTOM_REAL),dimension(:,:),allocatable :: tmp_jacobian2Dw
- integer :: ijk_face(3,NGLLX,NGLLY)
- integer,dimension(:,:,:),allocatable :: tmp_ijk
- integer,dimension(:),allocatable :: tmp_ispec
-
- integer,dimension(NGNOD2D) :: iglob_corners_ref !,iglob_corners
- integer :: ispec,i,j,k,igll,ier,iglob
- integer :: inum,iface_ref,icorner,iglob_midpoint ! iface,ispec_neighbor
- integer :: count_elastic,count_acoustic
-
- ! mpi interface communication
- integer, dimension(:), allocatable :: elastic_flag,acoustic_flag,test_flag
- integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh_dummy
- integer :: max_nibool_interfaces_ext_mesh
- logical, dimension(:), allocatable :: mask_ibool
-
- ! corners indices of reference cube faces
- integer,dimension(3,4),parameter :: iface1_corner_ijk = &
- reshape( (/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ /),(/3,4/)) ! xmin
- integer,dimension(3,4),parameter :: iface2_corner_ijk = &
- reshape( (/ NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),(/3,4/)) ! xmax
- integer,dimension(3,4),parameter :: iface3_corner_ijk = &
- reshape( (/ 1,1,1, 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,1,1 /),(/3,4/)) ! ymin
- integer,dimension(3,4),parameter :: iface4_corner_ijk = &
- reshape( (/ 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! ymax
- integer,dimension(3,4),parameter :: iface5_corner_ijk = &
- reshape( (/ 1,1,1, 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,1,1 /),(/3,4/)) ! bottom
- integer,dimension(3,4),parameter :: iface6_corner_ijk = &
- reshape( (/ 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! top
- integer,dimension(3,4,6),parameter :: iface_all_corner_ijk = &
- reshape( (/ iface1_corner_ijk,iface2_corner_ijk, &
- iface3_corner_ijk,iface4_corner_ijk, &
- iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
- ! midpoint indices for each face (xmin,xmax,ymin,ymax,zmin,zmax)
- integer,dimension(3,6),parameter :: iface_all_midpointijk = &
- reshape( (/ 1,2,2, NGLLX,2,2, 2,1,2, 2,NGLLY,2, 2,2,1, 2,2,NGLLZ /),(/3,6/)) ! top
-
-
- ! test vtk output
- !integer,dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: gll_data
- !character(len=256):: prname_file
-
-! allocates temporary arrays
- allocate(tmp_normal(NDIM,NGLLSQUARE,nspec*6))
- allocate(tmp_jacobian2Dw(NGLLSQUARE,nspec*6))
- allocate(tmp_ijk(3,NGLLSQUARE,nspec*6))
- allocate(tmp_ispec(nspec*6))
- tmp_ispec(:) = 0
- tmp_ijk(:,:,:) = 0
- tmp_normal(:,:,:) = 0.0
- tmp_jacobian2Dw(:,:) = 0.0
-
- ! sets flags for acoustic / elastic on global points
- allocate(elastic_flag(nglob),stat=ier)
- allocate(acoustic_flag(nglob),stat=ier)
- allocate(test_flag(nglob),stat=ier)
- allocate(mask_ibool(nglob),stat=ier)
- if( ier /= 0 ) stop 'error allocate flag array'
- elastic_flag(:) = 0
- acoustic_flag(:) = 0
- test_flag(:) = 0
- count_elastic = 0
- count_acoustic = 0
- do ispec = 1, nspec
- ! counts elements
- if( ispec_is_elastic(ispec) ) count_elastic = count_elastic + 1
- if( ispec_is_acoustic(ispec) ) count_acoustic = count_acoustic + 1
-
- ! sets flags on global points
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- ! global index
- iglob = ibool(i,j,k,ispec)
- ! sets elastic flag
- if( ispec_is_elastic(ispec) ) elastic_flag(iglob) = myrank+1
- ! sets acoustic flag
- if( ispec_is_acoustic(ispec) ) acoustic_flag(iglob) = myrank+1
- ! sets test flag
- test_flag(iglob) = myrank+1
- enddo
- enddo
- enddo
- enddo
- call sum_all_i(count_acoustic,inum)
- if( myrank == 0 ) then
- write(IMAIN,*) ' total acoustic elements:',inum
- endif
- call sum_all_i(count_elastic,inum)
- if( myrank == 0 ) then
- write(IMAIN,*) ' total elastic elements :',inum
- endif
-
- ! collects contributions from different MPI partitions
- ! sets up MPI communications
- max_nibool_interfaces_ext_mesh = maxval( nibool_interfaces_ext_mesh(:) )
- allocate(ibool_interfaces_ext_mesh_dummy(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
- if( ier /= 0 ) stop 'error allocating array'
- do i = 1, num_interfaces_ext_mesh
- ibool_interfaces_ext_mesh_dummy(:,i) = ibool_interfaces_ext_mesh(1:max_nibool_interfaces_ext_mesh,i)
- enddo
- ! sums elastic flags
- call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,elastic_flag, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
- my_neighbours_ext_mesh)
- ! sums acoustic flags
- call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,acoustic_flag, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
- my_neighbours_ext_mesh)
-
- ! sums test flags
- call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,test_flag, &
- num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
- nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
- my_neighbours_ext_mesh)
-
- ! loops over all element faces and
- ! counts number of coupling faces between acoustic and elastic elements
- mask_ibool(:) = .false.
- inum = 0
- do ispec=1,nspec
-
- ! loops over each face
- do iface_ref= 1, 6
-
- ! takes indices of corners of reference face
- do icorner = 1,NGNOD2D
- i = iface_all_corner_ijk(1,icorner,iface_ref)
- j = iface_all_corner_ijk(2,icorner,iface_ref)
- k = iface_all_corner_ijk(3,icorner,iface_ref)
- ! global reference indices
- iglob_corners_ref(icorner) = ibool(i,j,k,ispec)
-
- ! reference corner coordinates
- xcoord(icorner) = xstore_dummy(iglob_corners_ref(icorner))
- ycoord(icorner) = ystore_dummy(iglob_corners_ref(icorner))
- zcoord(icorner) = zstore_dummy(iglob_corners_ref(icorner))
- enddo
-
- ! checks if face has acoustic side
- if( acoustic_flag( iglob_corners_ref(1) ) >= 1 .and. &
- acoustic_flag( iglob_corners_ref(2) ) >= 1 .and. &
- acoustic_flag( iglob_corners_ref(3) ) >= 1 .and. &
- acoustic_flag( iglob_corners_ref(4) ) >= 1) then
- ! checks if face is has an elastic side
- if( elastic_flag( iglob_corners_ref(1) ) >= 1 .and. &
- elastic_flag( iglob_corners_ref(2) ) >= 1 .and. &
- elastic_flag( iglob_corners_ref(3) ) >= 1 .and. &
- elastic_flag( iglob_corners_ref(4) ) >= 1) then
-
- ! reference midpoint on face (used to avoid redundant face counting)
- i = iface_all_midpointijk(1,iface_ref)
- j = iface_all_midpointijk(2,iface_ref)
- k = iface_all_midpointijk(3,iface_ref)
- iglob_midpoint = ibool(i,j,k,ispec)
-
- ! checks if points on this face are masked already
- if( .not. mask_ibool(iglob_midpoint) ) then
-
- ! gets face GLL points i,j,k indices from element face
- call get_element_face_gll_indices(iface_ref,ijk_face,NGLLX,NGLLY)
-
- ! takes each element face only once, if it lies on an MPI interface
- ! note: this is not exactly load balanced
- ! lowest rank process collects as many faces as possible, second lowest as so forth
- if( (test_flag(iglob_midpoint) == myrank+1) .or. &
- (test_flag(iglob_midpoint) > 2*(myrank+1)) ) then
-
- ! gets face GLL 2Djacobian, weighted from element face
- call get_jacobian_boundary_face(myrank,nspec, &
- xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
-
- ! normal convention: points away from acoustic, reference element
- ! switch normal direction if necessary
- do j=1,NGLLY
- do i=1,NGLLX
- ! directs normals such that they point outwards of element
- call get_element_face_normal(ispec,iface_ref,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal_face(:,i,j) )
- ! makes sure that it always points away from acoustic element,
- ! otherwise switch direction
- if( ispec_is_elastic(ispec) ) normal_face(:,i,j) = - normal_face(:,i,j)
- enddo
- enddo
-
- ! stores informations about this face
- inum = inum + 1
- tmp_ispec(inum) = ispec
- igll = 0
- do j=1,NGLLY
- do i=1,NGLLX
- ! adds all gll points on this face
- igll = igll + 1
-
- ! do we need to store local i,j,k,ispec info? or only global indices iglob?
- tmp_ijk(:,igll,inum) = ijk_face(:,i,j)
-
- ! stores weighted jacobian and normals
- tmp_jacobian2Dw(igll,inum) = jacobian2Dw_face(i,j)
- tmp_normal(:,igll,inum) = normal_face(:,i,j)
-
- ! masks global points ( to avoid redundant counting of faces)
- iglob = ibool(ijk_face(1,i,j),ijk_face(2,i,j),ijk_face(3,i,j),ispec)
- mask_ibool(iglob) = .true.
- enddo
- enddo
- else
- ! assumes to be already collected by lower rank process, masks face points
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(ijk_face(1,i,j),ijk_face(2,i,j),ijk_face(3,i,j),ispec)
- mask_ibool(iglob) = .true.
- enddo
- enddo
- endif ! test_flag
- endif ! mask_ibool
- endif ! elastic_flag
- endif ! acoustic_flag
- enddo ! iface_ref
- enddo ! ispec
-
-! stores completed coupling face informations
-!
-! note: no need to store material parameters on these coupling points
-! for acoustic-elastic interface
- num_coupling_ac_el_faces = inum
- allocate(coupling_ac_el_normal(NDIM,NGLLSQUARE,num_coupling_ac_el_faces))
- allocate(coupling_ac_el_jacobian2Dw(NGLLSQUARE,num_coupling_ac_el_faces))
- allocate(coupling_ac_el_ijk(3,NGLLSQUARE,num_coupling_ac_el_faces))
- allocate(coupling_ac_el_ispec(num_coupling_ac_el_faces))
- do inum = 1,num_coupling_ac_el_faces
- coupling_ac_el_normal(:,:,inum) = tmp_normal(:,:,inum)
- coupling_ac_el_jacobian2Dw(:,inum) = tmp_jacobian2Dw(:,inum)
- coupling_ac_el_ijk(:,:,inum) = tmp_ijk(:,:,inum)
- coupling_ac_el_ispec(inum) = tmp_ispec(inum)
- enddo
-
-! user output
- call sum_all_i(num_coupling_ac_el_faces,inum)
- if( myrank == 0 ) then
- write(IMAIN,*) ' acoustic-elastic coupling:'
- write(IMAIN,*) ' total number of faces = ',inum
- endif
-
-end subroutine crm_ext_detect_ac_el_surface
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-subroutine crm_ext_create_ocean_load_mass(nglob,nspec,ibool,OCEANS,&
- UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION,NX_TOPO,NY_TOPO, &
- ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO, &
- itopo_bathy)
-
-! returns precomputed mass matrix in rmass array
-
- use create_regions_mesh_ext_par
- implicit none
-
-! number of spectral elements in each block
- integer :: nspec
- integer :: nglob
-
-! arrays with the mesh global indices
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
- logical :: OCEANS
-
-! use integer array to store topography values
- integer :: UTM_PROJECTION_ZONE
- logical :: SUPPRESS_UTM_PROJECTION
- integer :: NX_TOPO,NY_TOPO
- double precision :: ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO
- integer, dimension(NX_TOPO,NY_TOPO) :: itopo_bathy
-
-
-! local parameters
- double precision :: weight
- double precision :: xval,yval,long,lat,elevation
- double precision :: height_oceans
- double precision :: long_corner,lat_corner,ratio_xi,ratio_eta
- integer :: ix_oceans,iy_oceans,iz_oceans,ispec_oceans,ispec2D,igll,iglobnum
- integer :: icornerlong,icornerlat
-
-! creates ocean load mass matrix
- if(OCEANS) then
-
- ! adding ocean load mass matrix at ocean bottom
- NGLOB_OCEAN = nglob
- allocate(rmass_ocean_load(NGLOB_OCEAN))
-
- ! create ocean load mass matrix for degrees of freedom at ocean bottom
- rmass_ocean_load(:) = 0._CUSTOM_REAL
-
- ! add contribution of the oceans for surface elements exactly at ocean bottom
- do ispec2D = 1,num_free_surface_faces
-
- ispec_oceans = free_surface_ispec(ispec2D)
-
- ! only adds contribution if top boundary is elastic, no need to add this approximate calculation
- ! if top is already acoustic/poroelastic
- if( ispec_is_elastic(ispec_oceans) ) then
-
- do igll=1,NGLLSQUARE
- ix_oceans = free_surface_ijk(1,igll,ispec2D)
- iy_oceans = free_surface_ijk(1,igll,ispec2D)
- iz_oceans = free_surface_ijk(1,igll,ispec2D)
-
- iglobnum=ibool(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
-
- ! compute local height of oceans
-
- ! get coordinates of current point
- xval = xstore_dummy(iglobnum)
- yval = ystore_dummy(iglobnum)
-
- ! project x and y in UTM back to long/lat since topo file is in long/lat
- call utm_geo(long,lat,xval,yval,UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
-
- ! get coordinate of corner in bathy/topo model
- icornerlong = int((long - ORIG_LONG_TOPO) / DEGREES_PER_CELL_TOPO) + 1
- icornerlat = int((lat - ORIG_LAT_TOPO) / DEGREES_PER_CELL_TOPO) + 1
-
- ! avoid edge effects and extend with identical point if outside model
- if(icornerlong < 1) icornerlong = 1
- if(icornerlong > NX_TOPO-1) icornerlong = NX_TOPO-1
- if(icornerlat < 1) icornerlat = 1
- if(icornerlat > NY_TOPO-1) icornerlat = NY_TOPO-1
-
- ! compute coordinates of corner
- long_corner = ORIG_LONG_TOPO + (icornerlong-1)*DEGREES_PER_CELL_TOPO
- lat_corner = ORIG_LAT_TOPO + (icornerlat-1)*DEGREES_PER_CELL_TOPO
-
- ! compute ratio for interpolation
- ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO
- ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO
-
- ! avoid edge effects
- if(ratio_xi < 0.) ratio_xi = 0.
- if(ratio_xi > 1.) ratio_xi = 1.
- if(ratio_eta < 0.) ratio_eta = 0.
- if(ratio_eta > 1.) ratio_eta = 1.
-
- ! interpolate elevation at current point
- elevation = &
- itopo_bathy(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
- itopo_bathy(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
- itopo_bathy(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
- itopo_bathy(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
-
- ! suppress positive elevation, which means no oceans
- if(elevation >= - MINIMUM_THICKNESS_3D_OCEANS) then
- height_oceans = 0.d0
- else
- height_oceans = dabs(elevation)
- endif
-
- ! take into account inertia of water column
- weight = dble( free_surface_jacobian2Dw(igll,ispec2D)) &
- * dble(RHO_OCEANS) * height_oceans
-
- ! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- rmass_ocean_load(iglobnum) = rmass_ocean_load(iglobnum) + sngl(weight)
- else
- rmass_ocean_load(iglobnum) = rmass_ocean_load(iglobnum) + weight
- endif
-
- enddo ! igll
- endif ! ispec_is_elastic
- enddo ! num_free_surface_faces
-
- ! add regular mass matrix to ocean load contribution
- rmass_ocean_load(:) = rmass_ocean_load(:) + rmass(:)
-
- else
-
- ! allocate dummy array if no oceans
- NGLOB_OCEAN = 1
- allocate(rmass_ocean_load(NGLOB_OCEAN))
-
- endif
-
-end subroutine crm_ext_create_ocean_load_mass
-
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-subroutine create_regions_mesh_save_moho( myrank,nglob,nspec, &
- nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, &
- nodes_coords_ext_mesh,nnodes_ext_mesh,ibool )
-
- use create_regions_mesh_ext_par
- implicit none
-
- integer :: nspec2D_moho_ext
- integer, dimension(nspec2D_moho_ext) :: ibelm_moho
- integer, dimension(4,nspec2D_moho_ext) :: nodes_ibelm_moho
-
- integer :: myrank,nglob,nspec
-
- ! data from the external mesh
- integer :: nnodes_ext_mesh
- double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
-
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-
-! local parameters
- ! Moho mesh
- real(CUSTOM_REAL), dimension(:,:,:),allocatable :: normal_moho_top
- real(CUSTOM_REAL), dimension(:,:,:),allocatable :: normal_moho_bot
- integer,dimension(:,:,:),allocatable :: ijk_moho_top, ijk_moho_bot
- integer,dimension(:),allocatable :: ibelm_moho_top, ibelm_moho_bot
- integer :: NSPEC2D_MOHO
- logical, dimension(:),allocatable :: is_moho_top, is_moho_bot
-
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
- real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
- real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
- real(kind=CUSTOM_REAL),dimension(NDIM):: normal
- integer :: ijk_face(3,NGLLX,NGLLY)
-
- real(kind=CUSTOM_REAL),dimension(:,:),allocatable :: iglob_normals
- integer,dimension(:),allocatable:: iglob_is_surface
-
- integer :: imoho_bot,imoho_top
- integer :: ispec2D,ispec,icorner,iface,i,j,k,igll,iglob
- integer :: iglob_midpoint,idirect,counter
- integer :: imoho_top_all,imoho_bot_all,imoho_all
-
- ! corners indices of reference cube faces
- integer,dimension(3,4),parameter :: iface1_corner_ijk = &
- reshape( (/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ /),(/3,4/)) ! xmin
- integer,dimension(3,4),parameter :: iface2_corner_ijk = &
- reshape( (/ NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),(/3,4/)) ! xmax
- integer,dimension(3,4),parameter :: iface3_corner_ijk = &
- reshape( (/ 1,1,1, 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,1,1 /),(/3,4/)) ! ymin
- integer,dimension(3,4),parameter :: iface4_corner_ijk = &
- reshape( (/ 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! ymax
- integer,dimension(3,4),parameter :: iface5_corner_ijk = &
- reshape( (/ 1,1,1, 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,1,1 /),(/3,4/)) ! bottom
- integer,dimension(3,4),parameter :: iface6_corner_ijk = &
- reshape( (/ 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! top
- integer,dimension(3,4,6),parameter :: iface_all_corner_ijk = &
- reshape( (/ iface1_corner_ijk,iface2_corner_ijk, &
- iface3_corner_ijk,iface4_corner_ijk, &
- iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
- ! midpoint indices for each face (xmin,xmax,ymin,ymax,zmin,zmax)
- integer,dimension(3,6),parameter :: iface_all_midpointijk = &
- reshape( (/ 1,2,2, NGLLX,2,2, 2,1,2, 2,NGLLY,2, 2,2,1, 2,2,NGLLZ /),(/3,6/)) ! top
-
- ! temporary arrays for passing information
- allocate(iglob_is_surface(nglob))
- allocate(iglob_normals(NDIM,nglob))
- iglob_is_surface = 0
- iglob_normals = 0._CUSTOM_REAL
-
- ! loops over given moho surface elements
- do ispec2D=1, nspec2D_moho_ext
-
- ! gets element id
- ispec = ibelm_moho(ispec2D)
-
- ! looks for i,j,k indices of GLL points on boundary face
- ! determines element face by given CUBIT corners
- ! (note: uses point locations rather than point indices to find the element face,
- ! because the indices refer no more to the newly indexed ibool array )
- do icorner=1,NGNOD2D
- xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_moho(icorner,ispec2D))
- ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_moho(icorner,ispec2D))
- zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_moho(icorner,ispec2D))
- enddo
-
- ! sets face id of reference element associated with this face
- call get_element_face_id(ispec,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- iface)
-
- ! ijk indices of GLL points for face id
- call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ)
-
- ! weighted jacobian and normal
- call get_jacobian_boundary_face(myrank,nspec, &
- xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
-
- ! normal convention: points away from element
- ! switch normal direction if necessary
- do j=1,NGLLY
- do i=1,NGLLX
- call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal_face(:,i,j) )
- enddo
- enddo
-
- ! stores information on global points on moho surface
- igll = 0
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(ijk_face(1,i,j),ijk_face(2,i,j),ijk_face(3,i,j),ispec)
- ! sets flag
- iglob_is_surface(iglob) = ispec2D
- ! sets normals
- iglob_normals(:,iglob) = normal_face(:,i,j)
- enddo
- enddo
- enddo
-
- ! stores moho elements
- NSPEC2D_MOHO = nspec2D_moho_ext
-
- allocate(ibelm_moho_bot(NSPEC2D_MOHO))
- allocate(ibelm_moho_top(NSPEC2D_MOHO))
- allocate(normal_moho_top(NDIM,NGLLSQUARE,NSPEC2D_MOHO))
- allocate(normal_moho_bot(NDIM,NGLLSQUARE,NSPEC2D_MOHO))
- allocate(ijk_moho_bot(3,NGLLSQUARE,NSPEC2D_MOHO))
- allocate(ijk_moho_top(3,NGLLSQUARE,NSPEC2D_MOHO))
- ibelm_moho_bot = 0
- ibelm_moho_top = 0
-
- ! element flags
- allocate(is_moho_top(nspec))
- allocate(is_moho_bot(nspec))
- is_moho_top = .false.
- is_moho_bot = .false.
-
- ! finds spectral elements with moho surface
- imoho_top = 0
- imoho_bot = 0
- do ispec=1,nspec
-
- ! loops over each face
- do iface = 1,6
- ! checks if corners of face on surface
- counter = 0
- do icorner = 1,NGNOD2D
- i = iface_all_corner_ijk(1,icorner,iface)
- j = iface_all_corner_ijk(2,icorner,iface)
- k = iface_all_corner_ijk(3,icorner,iface)
- iglob = ibool(i,j,k,ispec)
-
- ! checks if point on surface
- if( iglob_is_surface(iglob) > 0 ) then
- counter = counter+1
-
- ! reference corner coordinates
- xcoord(icorner) = xstore_dummy(iglob)
- ycoord(icorner) = ystore_dummy(iglob)
- zcoord(icorner) = zstore_dummy(iglob)
- endif
- enddo
-
- ! stores moho informations
- if( counter == NGNOD2D ) then
-
- ! gets face GLL points i,j,k indices from element face
- call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY)
-
- ! re-computes face infos
- ! weighted jacobian and normal
- call get_jacobian_boundary_face(myrank,nspec, &
- xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
- ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
-
- ! normal convention: points away from element
- ! switch normal direction if necessary
- do j=1,NGLLZ
- do i=1,NGLLX
- call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal_face(:,i,j) )
- enddo
- enddo
-
- ! takes normal stored temporary on a face midpoint
- i = iface_all_midpointijk(1,iface)
- j = iface_all_midpointijk(2,iface)
- k = iface_all_midpointijk(3,iface)
- iglob_midpoint = ibool(i,j,k,ispec)
- normal(:) = iglob_normals(:,iglob_midpoint)
-
- ! determines whether normal points into element or not (top/bottom distinction)
- call get_element_face_normal_idirect(ispec,iface,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal,idirect )
-
- ! takes moho surface element id given by id on midpoint
- ispec2D = iglob_is_surface(iglob_midpoint)
-
- ! sets face infos for bottom (normal points away from element)
- if( idirect == 1 ) then
-
- ! checks validity
- if( is_moho_bot( ispec) .eqv. .true. ) then
- print*,'error: moho surface geometry bottom'
- print*,' does not allow for mulitple element faces in kernel computation'
- call exit_mpi(myrank,'error moho bottom elements')
- endif
-
- imoho_bot = imoho_bot + 1
- is_moho_bot(ispec) = .true.
- ibelm_moho_bot(ispec2D) = ispec
-
- ! stores on surface gll points (assuming NGLLX = NGLLY = NGLLZ)
- igll = 0
- do j=1,NGLLZ
- do i=1,NGLLX
- igll = igll+1
- ijk_moho_bot(:,igll,ispec2D) = ijk_face(:,i,j)
- normal_moho_bot(:,igll,ispec2D) = normal_face(:,i,j)
- enddo
- enddo
-
- ! sets face infos for top element
- else if( idirect == 2 ) then
-
- ! checks validity
- if( is_moho_top( ispec) .eqv. .true. ) then
- print*,'error: moho surface geometry top'
- print*,' does not allow for mulitple element faces kernel computation'
- call exit_mpi(myrank,'error moho top elements')
- endif
-
- imoho_top = imoho_top + 1
- is_moho_top(ispec) = .true.
- ibelm_moho_top(ispec2D) = ispec
-
- ! gll points
- igll = 0
- do j=1,NGLLZ
- do i=1,NGLLX
- igll = igll+1
- ijk_moho_top(:,igll,ispec) = ijk_face(:,i,j)
- ! note: top elements have normal pointing into element
- normal_moho_top(:,igll,ispec) = - normal_face(:,i,j)
- enddo
- enddo
- endif
-
- endif ! counter
-
- enddo ! iface
-
- ! checks validity of top/bottom distinction
- if( is_moho_top(ispec) .and. is_moho_bot(ispec) ) then
- print*,'error: moho surface elements confusing'
- print*,' element:',ispec,'has top and bottom surface'
- call exit_mpi(myrank,'error moho surface element')
- endif
-
- enddo ! ispec2D
-
- ! note: surface e.g. could be at the free-surface and have no top elements etc...
- ! user output
- call sum_all_i( imoho_top, imoho_top_all )
- call sum_all_i( imoho_bot, imoho_bot_all )
- call sum_all_i( NSPEC2D_MOHO, imoho_all )
- if( myrank == 0 ) then
- write(IMAIN,*) '********'
- write(IMAIN,*) 'Moho surface:'
- write(IMAIN,*) ' total surface elements: ',imoho_all
- write(IMAIN,*) ' top elements :',imoho_top_all
- write(IMAIN,*) ' bottom elements:',imoho_bot_all
- write(IMAIN,*) '********'
- endif
-
- ! saves moho files: total number of elements, corner points, all points
- open(unit=27,file=prname(1:len_trim(prname))//'ibelm_moho.bin',status='unknown',form='unformatted')
- write(27) NSPEC2D_MOHO
- write(27) ibelm_moho_top
- write(27) ibelm_moho_bot
- write(27) ijk_moho_top
- write(27) ijk_moho_bot
- close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'normal_moho.bin',status='unknown',form='unformatted')
- write(27) normal_moho_top
- write(27) normal_moho_bot
- close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'is_moho.bin',status='unknown',form='unformatted')
- write(27) is_moho_top
- write(27) is_moho_bot
- close(27)
-
-end subroutine create_regions_mesh_save_moho
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-subroutine bubble_sort( arr, ndim )
-
-! sorts values in array arr[ndim] in increasing order
-
- implicit none
-
- integer :: ndim
- integer :: arr(ndim)
-
- logical :: swapped
- integer :: j,tmp
-
- swapped = .true.
- do while( swapped )
- swapped = .false.
- do j = 1, ndim-1
- if( arr(j+1) < arr(j) ) then
- tmp = arr(j)
- arr(j) = arr(j+1)
- arr(j+1) = tmp
- swapped = .true.
- endif
- enddo
- enddo
-
-end subroutine bubble_sort
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-
!pll
! subroutine interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2010-02-25 21:58:34 UTC (rev 16342)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2010-02-25 23:07:29 UTC (rev 16343)
@@ -888,6 +888,9 @@
end do
end do
+
+ ! format:
+ ! # ispec_local # material_index_1 # material_index_2 # corner_id1 # corner_id2 # ... # corner_id8
write(IIN_database,*) glob2loc_elmnts(i)+1, num_modele(1,i+1), num_modele(2,i+1),(loc_nodes(k)+1, k=0,ngnod-1)
end if
end do
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90 2010-02-25 21:58:34 UTC (rev 16342)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90 2010-02-25 23:07:29 UTC (rev 16343)
@@ -593,7 +593,11 @@
read(IIN,*) nelmnts_ext_mesh
allocate(elmnts_ext_mesh(esize,nelmnts_ext_mesh))
allocate(mat_ext_mesh(2,nelmnts_ext_mesh))
+
+ ! reads in material association for each spectral element and corner node indices
do ispec = 1, nelmnts_ext_mesh
+ ! format:
+ ! # ispec_local # material_index_1 # material_index_2 # corner_id1 # corner_id2 # ... # corner_id8
read(IIN,*) dummy_elmnt, mat_ext_mesh(1,ispec),mat_ext_mesh(2,ispec), &
elmnts_ext_mesh(1,ispec), elmnts_ext_mesh(2,ispec), elmnts_ext_mesh(3,ispec), elmnts_ext_mesh(4,ispec), &
elmnts_ext_mesh(5,ispec), elmnts_ext_mesh(6,ispec), elmnts_ext_mesh(7,ispec), elmnts_ext_mesh(8,ispec)
Added: seismo/3D/SPECFEM3D_SESAME/trunk/get_MPI.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/get_MPI.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/get_MPI.f90 2010-02-25 23:07:29 UTC (rev 16343)
@@ -0,0 +1,229 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+ subroutine get_MPI(myrank,nglob,nspec,ibool, &
+ nelmnts_ext_mesh,elmnts_ext_mesh, &
+ my_nelmnts_neighbours_ext_mesh, my_interfaces_ext_mesh, &
+ ibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh, &
+ num_interfaces_ext_mesh,max_interface_size_ext_mesh, &
+ my_neighbours_ext_mesh,NPROC)
+
+! sets up the MPI interface for communication between partitions
+
+ use create_regions_mesh_ext_par
+ implicit none
+
+ integer :: myrank,nglob,nspec,NPROC
+
+! global indexing
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+! external mesh, element indexing
+ integer :: nelmnts_ext_mesh
+ integer, dimension(ESIZE,nelmnts_ext_mesh) :: elmnts_ext_mesh
+
+ integer :: num_interfaces_ext_mesh,max_interface_size_ext_mesh
+
+ integer, dimension(num_interfaces_ext_mesh) :: my_nelmnts_neighbours_ext_mesh
+ integer, dimension(6,max_interface_size_ext_mesh,num_interfaces_ext_mesh) :: my_interfaces_ext_mesh
+
+ integer, dimension(num_interfaces_ext_mesh) :: my_neighbours_ext_mesh
+
+ integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh
+ integer, dimension(NGLLX*NGLLX*max_interface_size_ext_mesh,num_interfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+
+
+ !integer :: nnodes_ext_mesh
+ !double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
+
+!local parameters
+ double precision, dimension(:), allocatable :: xp,yp,zp
+ double precision, dimension(:), allocatable :: work_ext_mesh
+
+ integer, dimension(:), allocatable :: locval
+ integer, dimension(:), allocatable :: nibool_interfaces_ext_mesh_true
+
+ ! for MPI buffers
+ integer, dimension(:), allocatable :: reorder_interface_ext_mesh,ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh
+ integer, dimension(:), allocatable :: ibool_interface_ext_mesh_dummy
+ logical, dimension(:), allocatable :: ifseg
+ integer :: iinterface,ilocnum
+ integer :: num_points1, num_points2
+
+ ! assembly test
+ integer :: i,j,k,ispec,iglob,count,inum
+ integer :: max_nibool_interfaces_ext_mesh
+ integer,dimension(:),allocatable :: test_flag
+ real(kind=CUSTOM_REAL), dimension(:),allocatable :: test_flag_cr
+ integer, dimension(:,:), allocatable :: ibool_interfaces_dummy
+
+! gets global indices for points on MPI interfaces (defined by my_interfaces_ext_mesh) between different partitions
+! and stores them in ibool_interfaces_ext_mesh & nibool_interfaces_ext_mesh (number of total points)
+ call prepare_assemble_MPI( nelmnts_ext_mesh,elmnts_ext_mesh, &
+ ibool,nglob,ESIZE, &
+ num_interfaces_ext_mesh, max_interface_size_ext_mesh, &
+ my_nelmnts_neighbours_ext_mesh, my_interfaces_ext_mesh, &
+ ibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh )
+
+ allocate(nibool_interfaces_ext_mesh_true(num_interfaces_ext_mesh))
+
+! sorts ibool comm buffers lexicographically for all MPI interfaces
+ num_points1 = 0
+ num_points2 = 0
+ do iinterface = 1, num_interfaces_ext_mesh
+
+ allocate(xp(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(yp(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(zp(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(locval(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ifseg(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(reorder_interface_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ibool_interface_ext_mesh_dummy(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ind_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ninseg_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(iwork_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(work_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+
+ ! gets x,y,z coordinates of global points on MPI interface
+ do ilocnum = 1, nibool_interfaces_ext_mesh(iinterface)
+ xp(ilocnum) = xstore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
+ yp(ilocnum) = ystore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
+ zp(ilocnum) = zstore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
+ enddo
+
+ ! sorts (lexicographically?) ibool_interfaces_ext_mesh and updates value
+ ! of total number of points nibool_interfaces_ext_mesh_true(iinterface)
+ call sort_array_coordinates(nibool_interfaces_ext_mesh(iinterface),xp,yp,zp, &
+ ibool_interfaces_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
+ reorder_interface_ext_mesh,locval,ifseg,nibool_interfaces_ext_mesh_true(iinterface), &
+ ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh,work_ext_mesh)
+
+ ! checks that number of MPI points are still the same
+ num_points1 = num_points1 + nibool_interfaces_ext_mesh(iinterface)
+ num_points2 = num_points2 + nibool_interfaces_ext_mesh_true(iinterface)
+ if( num_points1 /= num_points2 ) then
+ write(*,*) 'error sorting MPI interface points:',myrank
+ write(*,*) ' interface:',iinterface,num_points1,num_points2
+ call exit_mpi(myrank,'error sorting MPI interface')
+ endif
+ !write(*,*) myrank,'intfc',iinterface,num_points2,nibool_interfaces_ext_mesh_true(iinterface)
+
+ ! cleanup temporary arrays
+ deallocate(xp)
+ deallocate(yp)
+ deallocate(zp)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(reorder_interface_ext_mesh)
+ deallocate(ibool_interface_ext_mesh_dummy)
+ deallocate(ind_ext_mesh)
+ deallocate(ninseg_ext_mesh)
+ deallocate(iwork_ext_mesh)
+ deallocate(work_ext_mesh)
+
+ enddo
+
+ ! cleanup
+ deallocate(nibool_interfaces_ext_mesh_true)
+
+ ! outputs total number of MPI interface points
+ call sum_all_i(num_points2,ilocnum)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' total MPI interface points: ',ilocnum
+ endif
+
+! checks with assembly of test fields
+ allocate(test_flag(nglob),test_flag_cr(nglob))
+ test_flag(:) = 0
+ test_flag_cr(:) = 0._CUSTOM_REAL
+ count = 0
+ do ispec = 1, nspec
+ ! sets flags on global points
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ ! global index
+ iglob = ibool(i,j,k,ispec)
+
+ ! counts number of unique global points to set
+ if( test_flag(iglob) == 0 ) count = count+1
+
+ ! sets identifier
+ test_flag(iglob) = myrank + 1
+ test_flag_cr(iglob) = myrank + 1.0
+ enddo
+ enddo
+ enddo
+ enddo
+ call sync_all()
+
+ ! collects contributions from different MPI partitions
+ ! sets up MPI communications
+ max_nibool_interfaces_ext_mesh = maxval( nibool_interfaces_ext_mesh(:) )
+ allocate(ibool_interfaces_dummy(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
+
+ count = 0
+ do iinterface = 1, num_interfaces_ext_mesh
+ ibool_interfaces_dummy(:,iinterface) = ibool_interfaces_ext_mesh(1:max_nibool_interfaces_ext_mesh,iinterface)
+ count = count + nibool_interfaces_ext_mesh(iinterface)
+ !write(*,*) myrank,'interfaces ',iinterface,nibool_interfaces_ext_mesh(iinterface),max_nibool_interfaces_ext_mesh
+ enddo
+ call sync_all()
+
+ call sum_all_i(count,iglob)
+ if( myrank == 0 ) then
+ if( iglob /= ilocnum ) call exit_mpi(myrank,'error total global MPI interface points')
+ endif
+
+ ! adds contributions from different partitions to flag arrays
+ ! integer arrays
+ call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,test_flag, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_dummy,&
+ my_neighbours_ext_mesh)
+ ! custom_real arrays
+ call assemble_MPI_scalar_ext_mesh(NPROC,nglob,test_flag_cr, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_dummy, &
+ my_neighbours_ext_mesh)
+
+ ! checks number of interface points
+ i = 0
+ j = 0
+ do iglob=1,nglob
+ ! only counts flags with MPI contributions
+ if( test_flag(iglob) > myrank+1 ) i = i + 1
+ if( test_flag_cr(iglob) > myrank+1.0) j = j + 1
+ enddo
+ call sum_all_i(i,inum)
+ call sum_all_i(j,iglob)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' total assembled MPI interface points:',inum
+ if( inum /= iglob .or. inum > ilocnum ) call exit_mpi(myrank,'error MPI assembly')
+ endif
+
+ end subroutine get_MPI
Added: seismo/3D/SPECFEM3D_SESAME/trunk/get_absorbing_boundary.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/get_absorbing_boundary.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/get_absorbing_boundary.f90 2010-02-25 23:07:29 UTC (rev 16343)
@@ -0,0 +1,498 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+ subroutine get_absorbing_boundary(myrank,nspec,nglob,ibool, &
+ nodes_coords_ext_mesh,nnodes_ext_mesh, &
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+ nodes_ibelm_xmin,nodes_ibelm_xmax,nodes_ibelm_ymin,nodes_ibelm_ymax, &
+ nodes_ibelm_bottom,nodes_ibelm_top, &
+ nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
+ nspec2D_bottom,nspec2D_top)
+
+! determines absorbing boundaries/free-surface, 2D jacobians, face normals for Stacey conditions
+
+ use create_regions_mesh_ext_par
+ implicit none
+
+! number of spectral elements in each block
+ integer :: myrank,nspec,nglob
+
+! arrays with the mesh
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+! data from the external mesh
+ integer :: nnodes_ext_mesh
+ double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
+
+! absorbing boundaries (as defined in CUBIT)
+ integer :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, NSPEC2D_BOTTOM, NSPEC2D_TOP
+ ! element indices containing a boundary
+ integer, dimension(nspec2D_xmin) :: ibelm_xmin
+ integer, dimension(nspec2D_xmax) :: ibelm_xmax
+ integer, dimension(nspec2D_ymin) :: ibelm_ymin
+ integer, dimension(nspec2D_ymax) :: ibelm_ymax
+ integer, dimension(NSPEC2D_BOTTOM) :: ibelm_bottom
+ integer, dimension(NSPEC2D_TOP) :: ibelm_top
+
+ ! corner node indices of boundary faces coming from CUBIT
+ integer, dimension(4,nspec2D_xmin) :: nodes_ibelm_xmin
+ integer, dimension(4,nspec2D_xmax) :: nodes_ibelm_xmax
+ integer, dimension(4,nspec2D_ymin) :: nodes_ibelm_ymin
+ integer, dimension(4,nspec2D_ymax) :: nodes_ibelm_ymax
+ integer, dimension(4,NSPEC2D_BOTTOM) :: nodes_ibelm_bottom
+ integer, dimension(4,NSPEC2D_TOP) :: nodes_ibelm_top
+
+! local parameters
+ logical, dimension(:,:),allocatable :: iboun ! pll
+
+ ! (assumes NGLLX=NGLLY=NGLLZ)
+ real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
+ real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
+ integer:: ijk_face(3,NGLLX,NGLLY)
+
+ ! corner locations for faces
+ real(kind=CUSTOM_REAL), dimension(:,:,:),allocatable :: xcoord_iboun,ycoord_iboun,zcoord_iboun
+
+ ! face corner locations
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
+ integer :: ispec,ispec2D,icorner,ier,iabs,iface,igll,i,j,igllfree,ifree
+
+! allocate temporary flag array
+ allocate(iboun(6,nspec), &
+ xcoord_iboun(NGNOD2D,6,nspec), &
+ ycoord_iboun(NGNOD2D,6,nspec), &
+ zcoord_iboun(NGNOD2D,6,nspec),stat=ier)
+ if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
+
+! sets flag in array iboun for elements with an absorbing boundary faces
+ iboun(:,:) = .false.
+
+! abs face counter
+ iabs = 0
+
+ ! xmin
+ do ispec2D = 1, nspec2D_xmin
+ ! sets element
+ ispec = ibelm_xmin(ispec2D)
+
+ !if(myrank == 0 ) print*,'xmin:',ispec2D,ispec
+
+ ! looks for i,j,k indices of GLL points on boundary face
+ ! determines element face by given CUBIT corners
+ do icorner=1,NGNOD2D
+ xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_xmin(icorner,ispec2D))
+ ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_xmin(icorner,ispec2D))
+ zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_xmin(icorner,ispec2D))
+ !print*,'corner look:',icorner,xcoord(icorner),ycoord(icorner),zcoord(icorner)
+ enddo
+
+ ! sets face id of reference element associated with this face
+ call get_element_face_id(ispec,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ iface)
+
+ iboun(iface,ispec) = .true.
+
+ ! ijk indices of GLL points for face id
+ call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ)
+
+ ! weighted jacobian and normal
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
+
+ ! normal convention: points away from element
+ ! switch normal direction if necessary
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ enddo
+ enddo
+
+ ! sets face infos
+ iabs = iabs + 1
+ abs_boundary_ispec(iabs) = ispec
+
+ ! gll points -- assuming NGLLX = NGLLY = NGLLZ
+ igll = 0
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ igll = igll+1
+ abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
+ abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
+ abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
+ enddo
+ enddo
+
+ enddo ! nspec2D_xmin
+
+ ! xmax
+ do ispec2D = 1, nspec2D_xmax
+ ! sets element
+ ispec = ibelm_xmax(ispec2D)
+
+ ! looks for i,j,k indices of GLL points on boundary face
+ ! determines element face by given CUBIT corners
+ do icorner=1,NGNOD2D
+ xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_xmax(icorner,ispec2D))
+ ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_xmax(icorner,ispec2D))
+ zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_xmax(icorner,ispec2D))
+ enddo
+
+ ! sets face id of reference element associated with this face
+ call get_element_face_id(ispec,xcoord,ycoord,zcoord,&
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ iface )
+ iboun(iface,ispec) = .true.
+
+ ! ijk indices of GLL points on face
+ call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ)
+
+ ! weighted jacobian and normal
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
+
+ ! normal convention: points away from element
+ ! switch normal direction if necessary
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ enddo
+ enddo
+
+ ! sets face infos
+ iabs = iabs + 1
+ abs_boundary_ispec(iabs) = ispec
+
+ ! gll points -- assuming NGLLX = NGLLY = NGLLZ
+ igll = 0
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ igll = igll+1
+ abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
+ abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
+ abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
+ enddo
+ enddo
+
+ enddo
+
+ ! ymin
+ do ispec2D = 1, nspec2D_ymin
+ ! sets element
+ ispec = ibelm_ymin(ispec2D)
+
+ ! looks for i,j,k indices of GLL points on boundary face
+ ! determines element face by given CUBIT corners
+ do icorner=1,NGNOD2D
+ xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_ymin(icorner,ispec2D))
+ ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_ymin(icorner,ispec2D))
+ zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_ymin(icorner,ispec2D))
+ enddo
+
+ ! sets face id of reference element associated with this face
+ call get_element_face_id(ispec,xcoord,ycoord,zcoord,&
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ iface )
+ iboun(iface,ispec) = .true.
+
+ ! ijk indices of GLL points on face
+ call get_element_face_gll_indices(iface,ijk_face,NGLLY,NGLLZ)
+
+ ! weighted jacobian and normal
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ)
+
+ ! normal convention: points away from element
+ ! switch normal direction if necessary
+ do j=1,NGLLZ
+ do i=1,NGLLY
+ call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ enddo
+ enddo
+
+ ! sets face infos
+ iabs = iabs + 1
+ abs_boundary_ispec(iabs) = ispec
+
+ ! gll points -- assuming NGLLX = NGLLY = NGLLZ
+ igll = 0
+ do j=1,NGLLZ
+ do i=1,NGLLY
+ igll = igll+1
+ abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
+ abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
+ abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
+ enddo
+ enddo
+
+ enddo
+
+ ! ymax
+ do ispec2D = 1, nspec2D_ymax
+ ! sets element
+ ispec = ibelm_ymax(ispec2D)
+
+ ! looks for i,j,k indices of GLL points on boundary face
+ ! determines element face by given CUBIT corners
+ do icorner=1,NGNOD2D
+ xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_ymax(icorner,ispec2D))
+ ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_ymax(icorner,ispec2D))
+ zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_ymax(icorner,ispec2D))
+ enddo
+
+ ! sets face id of reference element associated with this face
+ call get_element_face_id(ispec,xcoord,ycoord,zcoord,&
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ iface )
+ iboun(iface,ispec) = .true.
+
+ ! ijk indices of GLL points on face
+ call get_element_face_gll_indices(iface,ijk_face,NGLLY,NGLLZ)
+
+ ! weighted jacobian and normal
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ)
+
+ ! normal convention: points away from element
+ ! switch normal direction if necessary
+ do j=1,NGLLZ
+ do i=1,NGLLY
+ call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ enddo
+ enddo
+
+ ! sets face infos
+ iabs = iabs + 1
+ abs_boundary_ispec(iabs) = ispec
+
+ ! gll points -- assuming NGLLX = NGLLY = NGLLZ
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ igll = igll+1
+ abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
+ abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
+ abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
+ enddo
+ enddo
+
+ enddo
+
+ ! bottom
+ do ispec2D = 1, NSPEC2D_BOTTOM
+ ! sets element
+ ispec = ibelm_bottom(ispec2D)
+
+ ! looks for i,j,k indices of GLL points on boundary face
+ ! determines element face by given CUBIT corners
+ do icorner=1,NGNOD2D
+ xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_bottom(icorner,ispec2D))
+ ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_bottom(icorner,ispec2D))
+ zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_bottom(icorner,ispec2D))
+ enddo
+
+ ! sets face id of reference element associated with this face
+ call get_element_face_id(ispec,xcoord,ycoord,zcoord,&
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ iface )
+ iboun(iface,ispec) = .true.
+
+ ! ijk indices of GLL points on face
+ call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY)
+
+ ! weighted jacobian and normal
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
+
+ ! normal convention: points away from element
+ ! switch normal direction if necessary
+ do j=1,NGLLY
+ do i=1,NGLLX
+ call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ enddo
+ enddo
+
+ ! sets face infos
+ iabs = iabs + 1
+ abs_boundary_ispec(iabs) = ispec
+
+ ! gll points -- assuming NGLLX = NGLLY = NGLLZ
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ igll = igll+1
+ abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
+ abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
+ abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
+ enddo
+ enddo
+
+ enddo
+
+ ! top
+ ! free surface face counter
+ ifree = 0
+ do ispec2D = 1, NSPEC2D_TOP
+ ! sets element
+ ispec = ibelm_top(ispec2D)
+
+ ! looks for i,j,k indices of GLL points on boundary face
+ ! determines element face by given CUBIT corners
+ do icorner=1,NGNOD2D
+ xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_top(icorner,ispec2D))
+ ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_top(icorner,ispec2D))
+ zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_top(icorner,ispec2D))
+ enddo
+
+ ! sets face id of reference element associated with this face
+ call get_element_face_id(ispec,xcoord,ycoord,zcoord,&
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ iface )
+ iboun(iface,ispec) = .true.
+
+ ! ijk indices of GLL points on face
+ call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY)
+
+ ! weighted jacobian and normal
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
+
+ ! normal convention: points away from element
+ ! switch normal direction if necessary
+ do j=1,NGLLY
+ do i=1,NGLLX
+ call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ enddo
+ enddo
+
+ ! stores surface infos
+ if( .not. ABSORB_FREE_SURFACE ) then
+ ! store for free surface
+ !jacobian2D_top(:,:,ispec2D) = jacobian2Dw_face(:,:)
+ !normal_top(:,:,:,ispec2D) = normal_face(:,:,:)
+
+ ! sets face infos
+ ifree = ifree + 1
+ free_surface_ispec(ifree) = ispec
+
+ ! gll points -- assuming NGLLX = NGLLY = NGLLZ
+ igllfree = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ igllfree = igllfree+1
+ free_surface_ijk(:,igllfree,ifree) = ijk_face(:,i,j)
+ free_surface_jacobian2Dw(igllfree,ifree) = jacobian2Dw_face(i,j)
+ free_surface_normal(:,igllfree,ifree) = normal_face(:,i,j)
+ enddo
+ enddo
+ else
+ ! adds face infos to absorbing boundary surface
+ iabs = iabs + 1
+ abs_boundary_ispec(iabs) = ispec
+
+ ! gll points -- assuming NGLLX = NGLLY = NGLLZ
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ igll = igll+1
+ abs_boundary_ijk(:,igll,iabs) = ijk_face(:,i,j)
+ abs_boundary_jacobian2Dw(igll,iabs) = jacobian2Dw_face(i,j)
+ abs_boundary_normal(:,igll,iabs) = normal_face(:,i,j)
+ enddo
+ enddo
+
+ ! resets free surface
+ ifree = 1
+ free_surface_ispec(:) = 0
+ free_surface_ijk(:,:,:) = 0
+ free_surface_jacobian2Dw(:,:) = 0.0
+ free_surface_normal(:,:,:) = 0.0
+ endif
+ enddo
+
+! checks counters
+ if( ifree /= num_free_surface_faces ) then
+ print*,'error number of free surface faces:',ifree,num_free_surface_faces
+ stop 'error number of free surface faces'
+ endif
+
+ if( iabs /= num_abs_boundary_faces ) then
+ print*,'error number of absorbing faces:',iabs,num_abs_boundary_faces
+ stop 'error number of absorbing faces'
+ endif
+
+ call sum_all_i(num_abs_boundary_faces,iabs)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' absorbing boundary:'
+ write(IMAIN,*) ' total number of faces = ',iabs
+ if( ABSORB_FREE_SURFACE ) then
+ write(IMAIN,*) ' absorbing boundary includes free surface'
+ endif
+ endif
+
+ end subroutine get_absorbing_boundary
+
Added: seismo/3D/SPECFEM3D_SESAME/trunk/get_coupling_surfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/get_coupling_surfaces.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/get_coupling_surfaces.f90 2010-02-25 23:07:29 UTC (rev 16343)
@@ -0,0 +1,308 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+ subroutine get_coupling_surfaces(myrank, &
+ nspec,nglob,ibool,NPROC, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ num_interfaces_ext_mesh,max_interface_size_ext_mesh, &
+ my_neighbours_ext_mesh)
+
+! determines coupling surface for acoustic-elastic domains
+
+ use create_regions_mesh_ext_par
+ implicit none
+
+! number of spectral elements in each block
+ integer :: myrank,nspec,nglob,NPROC
+
+! arrays with the mesh
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+! MPI communication
+ integer :: num_interfaces_ext_mesh,max_interface_size_ext_mesh
+ integer, dimension(num_interfaces_ext_mesh) :: my_neighbours_ext_mesh
+ integer, dimension(NGLLX*NGLLX*max_interface_size_ext_mesh,num_interfaces_ext_mesh) :: &
+ ibool_interfaces_ext_mesh
+ integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh
+
+! local parameters
+ ! (assumes NGLLX=NGLLY=NGLLZ)
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
+ real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
+ real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
+ real(kind=CUSTOM_REAL),dimension(:,:,:),allocatable :: tmp_normal
+ real(kind=CUSTOM_REAL),dimension(:,:),allocatable :: tmp_jacobian2Dw
+ integer :: ijk_face(3,NGLLX,NGLLY)
+ integer,dimension(:,:,:),allocatable :: tmp_ijk
+ integer,dimension(:),allocatable :: tmp_ispec
+
+ integer,dimension(NGNOD2D) :: iglob_corners_ref !,iglob_corners
+ integer :: ispec,i,j,k,igll,ier,iglob
+ integer :: inum,iface_ref,icorner,iglob_midpoint ! iface,ispec_neighbor
+ integer :: count_elastic,count_acoustic
+
+ ! mpi interface communication
+ integer, dimension(:), allocatable :: elastic_flag,acoustic_flag,test_flag
+ integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh_dummy
+ integer :: max_nibool_interfaces_ext_mesh
+ logical, dimension(:), allocatable :: mask_ibool
+
+ ! corners indices of reference cube faces
+ integer,dimension(3,4),parameter :: iface1_corner_ijk = &
+ reshape( (/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ /),(/3,4/)) ! xmin
+ integer,dimension(3,4),parameter :: iface2_corner_ijk = &
+ reshape( (/ NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),(/3,4/)) ! xmax
+ integer,dimension(3,4),parameter :: iface3_corner_ijk = &
+ reshape( (/ 1,1,1, 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,1,1 /),(/3,4/)) ! ymin
+ integer,dimension(3,4),parameter :: iface4_corner_ijk = &
+ reshape( (/ 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! ymax
+ integer,dimension(3,4),parameter :: iface5_corner_ijk = &
+ reshape( (/ 1,1,1, 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,1,1 /),(/3,4/)) ! bottom
+ integer,dimension(3,4),parameter :: iface6_corner_ijk = &
+ reshape( (/ 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! top
+ integer,dimension(3,4,6),parameter :: iface_all_corner_ijk = &
+ reshape( (/ iface1_corner_ijk,iface2_corner_ijk, &
+ iface3_corner_ijk,iface4_corner_ijk, &
+ iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
+ ! midpoint indices for each face (xmin,xmax,ymin,ymax,zmin,zmax)
+ integer,dimension(3,6),parameter :: iface_all_midpointijk = &
+ reshape( (/ 1,2,2, NGLLX,2,2, 2,1,2, 2,NGLLY,2, 2,2,1, 2,2,NGLLZ /),(/3,6/)) ! top
+
+
+ ! test vtk output
+ !integer,dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: gll_data
+ !character(len=256):: prname_file
+
+! allocates temporary arrays
+ allocate(tmp_normal(NDIM,NGLLSQUARE,nspec*6))
+ allocate(tmp_jacobian2Dw(NGLLSQUARE,nspec*6))
+ allocate(tmp_ijk(3,NGLLSQUARE,nspec*6))
+ allocate(tmp_ispec(nspec*6))
+ tmp_ispec(:) = 0
+ tmp_ijk(:,:,:) = 0
+ tmp_normal(:,:,:) = 0.0
+ tmp_jacobian2Dw(:,:) = 0.0
+
+ ! sets flags for acoustic / elastic on global points
+ allocate(elastic_flag(nglob),stat=ier)
+ allocate(acoustic_flag(nglob),stat=ier)
+ allocate(test_flag(nglob),stat=ier)
+ allocate(mask_ibool(nglob),stat=ier)
+ if( ier /= 0 ) stop 'error allocate flag array'
+ elastic_flag(:) = 0
+ acoustic_flag(:) = 0
+ test_flag(:) = 0
+ count_elastic = 0
+ count_acoustic = 0
+ do ispec = 1, nspec
+ ! counts elements
+ if( ispec_is_elastic(ispec) ) count_elastic = count_elastic + 1
+ if( ispec_is_acoustic(ispec) ) count_acoustic = count_acoustic + 1
+
+ ! sets flags on global points
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ ! global index
+ iglob = ibool(i,j,k,ispec)
+ ! sets elastic flag
+ if( ispec_is_elastic(ispec) ) elastic_flag(iglob) = myrank+1
+ ! sets acoustic flag
+ if( ispec_is_acoustic(ispec) ) acoustic_flag(iglob) = myrank+1
+ ! sets test flag
+ test_flag(iglob) = myrank+1
+ enddo
+ enddo
+ enddo
+ enddo
+ call sum_all_i(count_acoustic,inum)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' total acoustic elements:',inum
+ endif
+ call sum_all_i(count_elastic,inum)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' total elastic elements :',inum
+ endif
+
+ ! collects contributions from different MPI partitions
+ ! sets up MPI communications
+ max_nibool_interfaces_ext_mesh = maxval( nibool_interfaces_ext_mesh(:) )
+ allocate(ibool_interfaces_ext_mesh_dummy(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array'
+ do i = 1, num_interfaces_ext_mesh
+ ibool_interfaces_ext_mesh_dummy(:,i) = ibool_interfaces_ext_mesh(1:max_nibool_interfaces_ext_mesh,i)
+ enddo
+ ! sums elastic flags
+ call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,elastic_flag, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
+ my_neighbours_ext_mesh)
+ ! sums acoustic flags
+ call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,acoustic_flag, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
+ my_neighbours_ext_mesh)
+
+ ! sums test flags
+ call assemble_MPI_scalar_i_ext_mesh(NPROC,nglob,test_flag, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy,&
+ my_neighbours_ext_mesh)
+
+ ! loops over all element faces and
+ ! counts number of coupling faces between acoustic and elastic elements
+ mask_ibool(:) = .false.
+ inum = 0
+ do ispec=1,nspec
+
+ ! loops over each face
+ do iface_ref= 1, 6
+
+ ! takes indices of corners of reference face
+ do icorner = 1,NGNOD2D
+ i = iface_all_corner_ijk(1,icorner,iface_ref)
+ j = iface_all_corner_ijk(2,icorner,iface_ref)
+ k = iface_all_corner_ijk(3,icorner,iface_ref)
+ ! global reference indices
+ iglob_corners_ref(icorner) = ibool(i,j,k,ispec)
+
+ ! reference corner coordinates
+ xcoord(icorner) = xstore_dummy(iglob_corners_ref(icorner))
+ ycoord(icorner) = ystore_dummy(iglob_corners_ref(icorner))
+ zcoord(icorner) = zstore_dummy(iglob_corners_ref(icorner))
+ enddo
+
+ ! checks if face has acoustic side
+ if( acoustic_flag( iglob_corners_ref(1) ) >= 1 .and. &
+ acoustic_flag( iglob_corners_ref(2) ) >= 1 .and. &
+ acoustic_flag( iglob_corners_ref(3) ) >= 1 .and. &
+ acoustic_flag( iglob_corners_ref(4) ) >= 1) then
+ ! checks if face is has an elastic side
+ if( elastic_flag( iglob_corners_ref(1) ) >= 1 .and. &
+ elastic_flag( iglob_corners_ref(2) ) >= 1 .and. &
+ elastic_flag( iglob_corners_ref(3) ) >= 1 .and. &
+ elastic_flag( iglob_corners_ref(4) ) >= 1) then
+
+ ! reference midpoint on face (used to avoid redundant face counting)
+ i = iface_all_midpointijk(1,iface_ref)
+ j = iface_all_midpointijk(2,iface_ref)
+ k = iface_all_midpointijk(3,iface_ref)
+ iglob_midpoint = ibool(i,j,k,ispec)
+
+ ! checks if points on this face are masked already
+ if( .not. mask_ibool(iglob_midpoint) ) then
+
+ ! gets face GLL points i,j,k indices from element face
+ call get_element_face_gll_indices(iface_ref,ijk_face,NGLLX,NGLLY)
+
+ ! takes each element face only once, if it lies on an MPI interface
+ ! note: this is not exactly load balanced
+ ! lowest rank process collects as many faces as possible, second lowest as so forth
+ if( (test_flag(iglob_midpoint) == myrank+1) .or. &
+ (test_flag(iglob_midpoint) > 2*(myrank+1)) ) then
+
+ ! gets face GLL 2Djacobian, weighted from element face
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
+
+ ! normal convention: points away from acoustic, reference element
+ ! switch normal direction if necessary
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! directs normals such that they point outwards of element
+ call get_element_face_normal(ispec,iface_ref,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ ! makes sure that it always points away from acoustic element,
+ ! otherwise switch direction
+ if( ispec_is_elastic(ispec) ) normal_face(:,i,j) = - normal_face(:,i,j)
+ enddo
+ enddo
+
+ ! stores informations about this face
+ inum = inum + 1
+ tmp_ispec(inum) = ispec
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! adds all gll points on this face
+ igll = igll + 1
+
+ ! do we need to store local i,j,k,ispec info? or only global indices iglob?
+ tmp_ijk(:,igll,inum) = ijk_face(:,i,j)
+
+ ! stores weighted jacobian and normals
+ tmp_jacobian2Dw(igll,inum) = jacobian2Dw_face(i,j)
+ tmp_normal(:,igll,inum) = normal_face(:,i,j)
+
+ ! masks global points ( to avoid redundant counting of faces)
+ iglob = ibool(ijk_face(1,i,j),ijk_face(2,i,j),ijk_face(3,i,j),ispec)
+ mask_ibool(iglob) = .true.
+ enddo
+ enddo
+ else
+ ! assumes to be already collected by lower rank process, masks face points
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(ijk_face(1,i,j),ijk_face(2,i,j),ijk_face(3,i,j),ispec)
+ mask_ibool(iglob) = .true.
+ enddo
+ enddo
+ endif ! test_flag
+ endif ! mask_ibool
+ endif ! elastic_flag
+ endif ! acoustic_flag
+ enddo ! iface_ref
+ enddo ! ispec
+
+! stores completed coupling face informations
+!
+! note: no need to store material parameters on these coupling points
+! for acoustic-elastic interface
+ num_coupling_ac_el_faces = inum
+ allocate(coupling_ac_el_normal(NDIM,NGLLSQUARE,num_coupling_ac_el_faces))
+ allocate(coupling_ac_el_jacobian2Dw(NGLLSQUARE,num_coupling_ac_el_faces))
+ allocate(coupling_ac_el_ijk(3,NGLLSQUARE,num_coupling_ac_el_faces))
+ allocate(coupling_ac_el_ispec(num_coupling_ac_el_faces))
+ do inum = 1,num_coupling_ac_el_faces
+ coupling_ac_el_normal(:,:,inum) = tmp_normal(:,:,inum)
+ coupling_ac_el_jacobian2Dw(:,inum) = tmp_jacobian2Dw(:,inum)
+ coupling_ac_el_ijk(:,:,inum) = tmp_ijk(:,:,inum)
+ coupling_ac_el_ispec(inum) = tmp_ispec(inum)
+ enddo
+
+! user output
+ call sum_all_i(num_coupling_ac_el_faces,inum)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' acoustic-elastic coupling:'
+ write(IMAIN,*) ' total number of faces = ',inum
+ endif
+
+ end subroutine get_coupling_surfaces
+
Added: seismo/3D/SPECFEM3D_SESAME/trunk/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/get_model.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/get_model.f90 2010-02-25 23:07:29 UTC (rev 16343)
@@ -0,0 +1,411 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+ subroutine get_model(nspec,ibool,mat_ext_mesh,nelmnts_ext_mesh, &
+ materials_ext_mesh,nmat_ext_mesh, &
+ undef_mat_prop,nundefMat_ext_mesh, &
+ ANISOTROPY)
+
+ use create_regions_mesh_ext_par
+ implicit none
+
+ ! number of spectral elements in each block
+ integer :: nspec
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+ ! external mesh
+ integer :: nelmnts_ext_mesh
+ integer :: nmat_ext_mesh,nundefMat_ext_mesh
+
+ integer, dimension(2,nelmnts_ext_mesh) :: mat_ext_mesh
+ double precision, dimension(6,nmat_ext_mesh) :: materials_ext_mesh
+ character (len=30), dimension(6,nundefMat_ext_mesh):: undef_mat_prop
+
+ ! anisotropy
+ logical :: ANISOTROPY
+
+ ! local parameters
+ real(kind=CUSTOM_REAL) :: vp,vs,rho
+ real(kind=CUSTOM_REAL) :: c11,c12,c13,c14,c15,c16,c22,c23,c24,c25, &
+ c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+ integer :: ispec,i,j,k,iundef,iflag_atten
+ integer :: iflag,flag_below,flag_above
+ integer :: iflag_aniso,idomain_id,imaterial_id
+
+! ! Piero, read bedrock file
+! allocate(ibedrock(NX_TOPO_ANT,NY_TOPO_ANT))
+! if(myrank == 0) then
+! call read_bedrock_file(ibedrock)
+! ! write(IMAIN,*)
+! ! write(IMAIN,*) 'regional bedrock file read ranges in m from ',minval(ibedrock),' to ',maxval(ibedrock)
+! ! write(IMAIN,*)
+! endif
+! ! broadcast the information read on the master to the nodes
+! ! call MPI_BCAST(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT,MPI_REAL,0,MPI_COMM_WORLD,ier)
+! call bcast_all_cr(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT)
+
+ ispec_is_acoustic(:) = .false.
+ ispec_is_elastic(:) = .false.
+ ispec_is_poroelastic(:) = .false.
+
+ ! material properties on all GLL points: taken from material values defined for
+ ! each spectral element in input mesh
+ do ispec = 1, nspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ ! material index 1: associated material number
+ imaterial_id = mat_ext_mesh(1,ispec)
+
+ ! check if the material is known or unknown
+ if( imaterial_id > 0) then
+ ! gets velocity model as specified by (cubit) mesh files
+
+ ! density
+ ! materials_ext_mesh format:
+ ! #index1 = rho #index2 = vp #index3 = vs #index4 = Q_flag #index5 = 0
+ rho = materials_ext_mesh(1,imaterial_id)
+
+ ! isotropic values: vp, vs
+ vp = materials_ext_mesh(2,imaterial_id)
+ vs = materials_ext_mesh(3,imaterial_id)
+
+ ! attenuation
+ iflag_atten = materials_ext_mesh(4,imaterial_id)
+ !change for piero :
+ !if(mat_ext_mesh(1,ispec) == 1) then
+ ! iflag_attenuation_store(i,j,k,ispec) = 1
+ !else
+ ! iflag_attenuation_store(i,j,k,ispec) = 2
+ !endif
+
+ ! anisotropy
+ iflag_aniso = materials_ext_mesh(5,imaterial_id)
+
+ ! material domain_id
+ idomain_id = materials_ext_mesh(6,imaterial_id)
+
+ else if (mat_ext_mesh(2,ispec) == 1) then
+
+ stop 'material: interface not implemented yet'
+
+ do iundef = 1,nundefMat_ext_mesh
+ if(trim(undef_mat_prop(2,iundef)) == 'interface') then
+ read(undef_mat_prop(4,iundef),'(1i3)') flag_below
+ read(undef_mat_prop(5,iundef),'(1i3)') flag_above
+ endif
+ enddo
+
+ !call interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
+
+ iflag = 1
+ rho = materials_ext_mesh(1,iflag)
+ vp = materials_ext_mesh(2,iflag)
+ vs = materials_ext_mesh(3,iflag)
+ iflag_atten = materials_ext_mesh(4,iflag)
+ !change for piero :
+ ! if(iflag == 1) then
+ ! iflag_attenuation_store(i,j,k,ispec) = 1
+ ! else
+ ! iflag_attenuation_store(i,j,k,ispec) = 2
+ ! endif
+ iflag_aniso = materials_ext_mesh(5,iflag)
+ idomain_id = materials_ext_mesh(6,iflag)
+ else
+
+ stop 'material: tomography not implemented yet'
+ ! call tomography()
+
+ end if
+
+ ! adds/gets velocity model as specified in model_external_values.f90
+ if( USE_MODEL_EXTERNAL_VALUES ) then
+ call model_external_values(i,j,k,ispec,idomain_id,imaterial_id, &
+ nspec,ibool, &
+ iflag_aniso,iflag_atten, &
+ rho,vp,vs, &
+ c11,c12,c13,c14,c15,c16, &
+ c22,c23,c24,c25,c26,c33, &
+ c34,c35,c36,c44,c45,c46, &
+ c55,c56,c66,ANISOTROPY)
+ endif
+
+ ! adds anisotropic default model
+ if( ANISOTROPY .and. .not. USE_MODEL_EXTERNAL_VALUES ) then
+ call model_aniso(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
+ c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45, &
+ c46,c55,c56,c66)
+
+ endif
+
+ ! stores velocity model
+
+ ! density
+ rhostore(i,j,k,ispec) = rho
+
+ ! kappa, mu
+ kappastore(i,j,k,ispec) = rho*( vp*vp - FOUR_THIRDS*vs*vs )
+ mustore(i,j,k,ispec) = rho*vs*vs
+
+ ! attenuation
+ iflag_attenuation_store(i,j,k,ispec) = iflag_atten
+
+ ! Stacey, a completer par la suite
+ rho_vp(i,j,k,ispec) = rho*vp
+ rho_vs(i,j,k,ispec) = rho*vs
+ !end pll
+
+ ! adds anisotropic perturbation to vp, vs
+ if( ANISOTROPY ) then
+ !call model_aniso(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
+ ! c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
+ c11store(i,j,k,ispec) = c11
+ c12store(i,j,k,ispec) = c12
+ c13store(i,j,k,ispec) = c13
+ c14store(i,j,k,ispec) = c14
+ c15store(i,j,k,ispec) = c15
+ c16store(i,j,k,ispec) = c16
+ c22store(i,j,k,ispec) = c22
+ c23store(i,j,k,ispec) = c23
+ c24store(i,j,k,ispec) = c24
+ c25store(i,j,k,ispec) = c25
+ c26store(i,j,k,ispec) = c26
+ c33store(i,j,k,ispec) = c33
+ c34store(i,j,k,ispec) = c34
+ c35store(i,j,k,ispec) = c35
+ c36store(i,j,k,ispec) = c36
+ c44store(i,j,k,ispec) = c44
+ c45store(i,j,k,ispec) = c45
+ c46store(i,j,k,ispec) = c46
+ c55store(i,j,k,ispec) = c55
+ c56store(i,j,k,ispec) = c56
+ c66store(i,j,k,ispec) = c66
+ endif
+
+ ! material domain
+ !print*,'velocity model:',ispec,idomain_id
+ if( idomain_id == IDOMAIN_ACOUSTIC ) then
+ ispec_is_acoustic(ispec) = .true.
+ else if( idomain_id == IDOMAIN_ELASTIC ) then
+ ispec_is_elastic(ispec) = .true.
+ else if( idomain_id == IDOMAIN_POROELASTIC ) then
+ stop 'poroelastic material domain not implemented yet'
+ ispec_is_poroelastic(ispec) = .true.
+ else
+ stop 'error material domain index'
+ endif
+
+ enddo
+ enddo
+ enddo
+ !print*,myrank,'ispec:',ispec,'rho:',rhostore(1,1,1,ispec),'vp:',vpstore(1,1,1,ispec),'vs:',vsstore(1,1,1,ispec)
+ enddo
+
+ ! checks material domains
+ do ispec=1,nspec
+ if( (ispec_is_acoustic(ispec) .eqv. .false.) &
+ .and. (ispec_is_elastic(ispec) .eqv. .false.) &
+ .and. (ispec_is_poroelastic(ispec) .eqv. .false.) ) then
+ print*,'error material domain not assigned to element:',ispec
+ print*,'acoustic: ',ispec_is_acoustic(ispec)
+ print*,'elastic: ',ispec_is_elastic(ispec)
+ print*,'poroelastic: ',ispec_is_poroelastic(ispec)
+ stop 'error material domain index element'
+ endif
+ enddo
+
+
+! !! DK DK store the position of the six stations to be able to
+! !! DK DK exclude circles around each station to make sure they are on the bedrock
+! !! DK DK and not in the ice
+! utm_x_station(1) = 783500.6250000d0
+! utm_y_station(1) = -11828.7519531d0
+
+! utm_x_station(2) = 853644.5000000d0
+! utm_y_station(2) = -114.0138092d0
+
+! utm_x_station(3) = 863406.0000000d0
+! utm_y_station(3) = -53736.1640625d0
+
+! utm_x_station(4) = 823398.8125000d0
+! utm_y_station(4) = 29847.4511719d0
+
+! utm_x_station(5) = 863545.3750000d0
+! utm_y_station(5) = 19669.6621094d0
+
+! utm_x_station(6) = 817099.3750000d0
+! utm_y_station(6) = -24430.2871094d0
+
+! print*,myrank,'après store the position of the six stations'
+! call flush(6)
+
+! print*, myrank,minval(nodes_coords_ext_mesh(1,:))
+! call flush(6)
+
+
+! print*, myrank,maxval(nodes_coords_ext_mesh(1,:))
+! call flush(6)
+
+
+! do ispec = 1, nspec
+
+! zmesh = zstore(2,2,2,ispec)
+
+! ! if(doubling_index == IFLAG_ONE_LAYER_TOPOGRAPHY) then
+! if(any(ibelm_top == ispec)) then
+! doubling_value_found_for_Piero = IFLAG_ONE_LAYER_TOPOGRAPHY
+
+! else if(zmesh < Z_23p4km) then
+! doubling_value_found_for_Piero = IFLAG_MANTLE_BELOW_23p4km
+
+! else if(zmesh < Z_14km) then
+! doubling_value_found_for_Piero = IFLAG_14km_to_23p4km
+
+! else
+! doubling_value_found_for_Piero = IFLAG_BEDROCK_down_to_14km
+! endif
+! idoubling(ispec) = doubling_value_found_for_Piero
+
+! do k = 1, NGLLZ
+! do j = 1, NGLLY
+! do i = 1, NGLLX
+
+
+! if(idoubling(ispec) == IFLAG_ONE_LAYER_TOPOGRAPHY .or. &
+! idoubling(ispec) == IFLAG_BEDROCK_down_to_14km) then
+
+! ! since we have suppressed UTM projection for Piero Basini, UTMx is the same as long
+! ! and UTMy is the same as lat
+! long = xstore(i,j,k,ispec)
+! lat = ystore(i,j,k,ispec)
+
+! ! get coordinate of corner in model
+! icornerlong = int((long - ORIG_LONG_TOPO) / DEGREES_PER_CELL_TOPO) + 1
+! icornerlat = int((lat - ORIG_LAT_TOPO) / DEGREES_PER_CELL_TOPO) + 1
+
+! ! avoid edge effects and extend with identical point if outside model
+! if(icornerlong < 1) icornerlong = 1
+! if(icornerlong > NX_TOPO-1) icornerlong = NX_TOPO-1
+! if(icornerlat < 1) icornerlat = 1
+! if(icornerlat > NY_TOPO-1) icornerlat = NY_TOPO-1
+
+! ! compute coordinates of corner
+! long_corner = ORIG_LONG_TOPO + (icornerlong-1)*DEGREES_PER_CELL_TOPO
+! lat_corner = ORIG_LAT_TOPO + (icornerlat-1)*DEGREES_PER_CELL_TOPO
+
+! ! compute ratio for interpolation
+! ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO
+! ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO
+
+! ! avoid edge effects
+! if(ratio_xi < 0.) ratio_xi = 0.
+! if(ratio_xi > 1.) ratio_xi = 1.
+! if(ratio_eta < 0.) ratio_eta = 0.
+! if(ratio_eta > 1.) ratio_eta = 1.
+
+! ! interpolate elevation at current point
+! elevation_bedrock = &
+! ibedrock(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
+! ibedrock(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
+! ibedrock(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
+! ibedrock(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
+
+! !! DK DK exclude circles around each station to make sure they are on the bedrock
+! !! DK DK and not in the ice
+! is_around_a_station = .false.
+! do istation = 1,NUMBER_OF_STATIONS
+! if(sqrt((long - utm_x_station(istation))**2 + (lat - utm_y_station(istation))**2) < RADIUS_TO_EXCLUDE) then
+! is_around_a_station = .true.
+! exit
+! endif
+! enddo
+
+! ! define elastic parameters in the model
+
+! ! we are above the bedrock interface i.e. in the ice, and not too close to a station
+! if(zmesh >= elevation_bedrock .and. .not. is_around_a_station) then
+! vp = 3800.d0
+! vs = 1900.d0
+! rho = 900.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_ICE
+
+! ! we are below the bedrock interface i.e. in the bedrock, or close to a station
+! else
+! vp = 5800.d0
+! vs = 3200.d0
+! rho = 2600.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
+! endif
+
+! else if(idoubling(ispec) == IFLAG_14km_to_23p4km) then
+! vp = 6800.d0
+! vs = 3900.d0
+! rho = 2900.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
+
+! else if(idoubling(ispec) == IFLAG_MANTLE_BELOW_23p4km) then
+! vp = 8100.d0
+! vs = 4480.d0
+! rho = 3380.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
+
+! endif
+
+! !pll 8/06
+! if(CUSTOM_REAL == SIZE_REAL) then
+! rhostore(i,j,k,ispec) = sngl(rho)
+! vpstore(i,j,k,ispec) = sngl(vp)
+! vsstore(i,j,k,ispec) = sngl(vs)
+! else
+! rhostore(i,j,k,ispec) = rho
+! vpstore(i,j,k,ispec) = vp
+! vsstore(i,j,k,ispec) = vs
+! end if
+
+! kappastore(i,j,k,ispec) = rhostore(i,j,k,ispec)*(vpstore(i,j,k,ispec)*vpstore(i,j,k,ispec) - &
+! 4.d0*vsstore(i,j,k,ispec)*vsstore(i,j,k,ispec)/3.d0)
+! mustore(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)*&
+! vsstore(i,j,k,ispec)
+
+! ! Stacey, a completer par la suite
+! rho_vp(i,j,k,ispec) = rhostore(i,j,k,ispec)*vpstore(i,j,k,ispec)
+! rho_vs(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)
+! !end pll
+
+! ! kappastore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
+! ! (materials_ext_mesh(2,mat_ext_mesh(ispec))*materials_ext_mesh(2,mat_ext_mesh(ispec)) - &
+! ! 4.d0*materials_ext_mesh(3,mat_ext_mesh(ispec))*materials_ext_mesh(3,mat_ext_mesh(ispec))/3.d0)
+! ! mustore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
+! materials_ext_mesh(3,mat_ext_mesh(ispec))*&
+! ! x materials_ext_mesh(3,mat_ext_mesh(ispec))
+! enddo
+! enddo
+! enddo
+! enddo
+
+ end subroutine get_model
+
Added: seismo/3D/SPECFEM3D_SESAME/trunk/model_aniso.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/model_aniso.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/model_aniso.f90 2010-02-25 23:07:29 UTC (rev 16343)
@@ -0,0 +1,300 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+!=====================================================================
+! 07/09/04 Last changed by Min Chen
+! Users need to modify this subroutine to implement their own
+! anisotropic models.
+!=====================================================================
+
+ subroutine model_aniso(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
+ c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
+
+ implicit none
+
+ include "constants.h"
+
+! see for example:
+!
+! M. Chen & J. Tromp, 2006. Theoretical & numerical investigations
+! of global and regional seismic wave propagation in weakly anisotropic earth models,
+! GJI, 168, 1130-1152.
+
+!------------------------------------------------------------------------------
+! for anisotropy simulations in a halfspace model
+
+! only related to body waves
+! one-zeta term
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_CS1p_A = 0.2_CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_CS1sv_A = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_CS1sh_N = 0._CUSTOM_REAL
+! three-zeta term
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_CS3_L = 0._CUSTOM_REAL
+
+! Relative to Love wave
+! four-zeta term
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_N = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_E_N = 0._CUSTOM_REAL
+
+! Relative to Rayleigh wave
+! two-zeta term
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_A = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_C = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_F = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_H_F = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_B_A = 0._CUSTOM_REAL
+
+! Relative to both Love wave and Rayleigh wave
+! two-zeta term
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_L = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_G_L = 0._CUSTOM_REAL
+
+!------------------------------------------------------------------------------
+
+ !integer idoubling
+ integer iflag_aniso
+
+ !real(kind=CUSTOM_REAL) zmesh
+ real(kind=CUSTOM_REAL) rho,vp,vs
+ real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36, &
+ c44,c45,c46,c55,c56,c66
+
+! local parameters
+ real(kind=CUSTOM_REAL) vpv,vph,vsv,vsh,eta_aniso
+ real(kind=CUSTOM_REAL) aa,cc,nn,ll,ff
+ real(kind=CUSTOM_REAL) A,C,F,AL,AN,Bc,Bs,Gc,Gs,Hc,Hs,Ec,Es,C1p,C1sv,C1sh,C3,S1p,S1sv,S1sh,S3
+ real(kind=CUSTOM_REAL) d11,d12,d13,d14,d15,d16,d22,d23,d24,d25,d26,d33,d34,d35,d36, &
+ d44,d45,d46,d55,d56,d66
+
+! assumes vp,vs given in m/s, rho in kg/m**3
+ vph = vp
+ vpv = vp
+ vsh = vs
+ vsv = vs
+ eta_aniso = 1.0_CUSTOM_REAL
+
+
+! for definition, see for example:
+!
+! Dziewonski & Anderson, 1981. Preliminary reference earth model, PEPI, 25, 297-356.
+! page 305:
+ aa = rho*vph*vph
+ cc = rho*vpv*vpv
+ nn = rho*vsh*vsh
+ ll = rho*vsv*vsv
+ ff = eta_aniso*(aa - 2.*ll)
+
+! Add anisotropic perturbation
+
+! notation: see Chen & Tromp, 2006, appendix A, page 1151
+!
+! zeta-independant terms:
+! A = \delta A
+! C = \delta C
+! AN = \delta N
+! AL = \delta L
+! F = \delta F
+!
+! zeta-dependant terms:
+! C1p = J_c
+! C1sv = K_c
+! C1sh = M_c
+! S1p = J_s
+! S1sv = K_s
+! S1sh = M_s
+!
+! two-zeta dependant terms:
+! Gc = G_c
+! Gs = G_s
+! Bc = B_c
+! Bs = B_s
+! Hc = H_c
+! Hs = H_s
+!
+! three-zeta dependant terms:
+! C3 = D_c
+! S3 = D_s
+!
+! four-zeta dependant terms:
+! Ec = E_c
+! Es = E_s
+
+! no anisotropic perturbation
+ if( iflag_aniso <= 0 ) then
+ ! zeta-independant
+ A = aa
+ C = cc
+ AN = nn
+ AL = ll
+ F = ff
+
+ ! zeta-dependant terms
+ C1p = 0._CUSTOM_REAL
+ C1sv = 0._CUSTOM_REAL
+ C1sh = 0._CUSTOM_REAL
+ S1p = 0._CUSTOM_REAL
+ S1sv = 0._CUSTOM_REAL
+ S1sh = 0._CUSTOM_REAL
+
+ ! two-zeta dependant terms
+ Gc = 0._CUSTOM_REAL
+ Gs = 0._CUSTOM_REAL
+
+ Bc = 0._CUSTOM_REAL
+ Bs = 0._CUSTOM_REAL
+
+ Hc = 0._CUSTOM_REAL
+ Hs = 0._CUSTOM_REAL
+
+ ! three-zeta dependant terms
+ C3 = 0._CUSTOM_REAL
+ S3 = 0._CUSTOM_REAL
+
+ ! four-zeta dependant terms
+ Ec = 0._CUSTOM_REAL
+ Es = 0._CUSTOM_REAL
+ endif
+
+! perturbation model 1
+ if( iflag_aniso == IANISOTROPY_MODEL1 ) then
+ ! zeta-independant
+ A = aa*(1.0_CUSTOM_REAL + FACTOR_A)
+ C = cc*(1.0_CUSTOM_REAL + FACTOR_C)
+ AN = nn*(1.0_CUSTOM_REAL + FACTOR_N)
+ AL = ll*(1.0_CUSTOM_REAL + FACTOR_L)
+ F = ff*(1.0_CUSTOM_REAL + FACTOR_F)
+
+ ! zeta-dependant terms
+ C1p = FACTOR_CS1p_A*aa
+ C1sv = FACTOR_CS1sv_A*aa
+ C1sh = FACTOR_CS1sh_N*nn
+ S1p = 0._CUSTOM_REAL
+ S1sv = 0._CUSTOM_REAL
+ S1sh = 0._CUSTOM_REAL
+
+ ! two-zeta dependant terms
+ Gc = FACTOR_G_L*ll
+ Bc = FACTOR_B_A*aa
+ Hc = FACTOR_H_F*ff
+ Gs = 0._CUSTOM_REAL
+ Bs = 0._CUSTOM_REAL
+ Hs = 0._CUSTOM_REAL
+
+ ! three-zeta dependant terms
+ C3 = FACTOR_CS3_L*ll
+ S3 = 0._CUSTOM_REAL
+
+ ! four-zeta dependant terms
+ Ec = FACTOR_E_N*nn
+ Es = 0._CUSTOM_REAL
+ endif
+
+! perturbation model 2
+ if( iflag_aniso == IANISOTROPY_MODEL2 ) then
+ ! zeta-independant
+ A = aa*(1.0_CUSTOM_REAL + FACTOR_A + 0.1)
+ C = cc*(1.0_CUSTOM_REAL + FACTOR_C + 0.1)
+ AN = nn*(1.0_CUSTOM_REAL + FACTOR_N + 0.1)
+ AL = ll*(1.0_CUSTOM_REAL + FACTOR_L + 0.1)
+ F = ff*(1.0_CUSTOM_REAL + FACTOR_F + 0.1)
+
+ ! zeta-dependant terms
+ C1p = FACTOR_CS1p_A*aa
+ C1sv = FACTOR_CS1sv_A*aa
+ C1sh = FACTOR_CS1sh_N*nn
+ S1p = 0._CUSTOM_REAL
+ S1sv = 0._CUSTOM_REAL
+ S1sh = 0._CUSTOM_REAL
+
+ ! two-zeta dependant terms
+ Gc = FACTOR_G_L*ll
+ Bc = FACTOR_B_A*aa
+ Hc = FACTOR_H_F*ff
+ Gs = 0._CUSTOM_REAL
+ Bs = 0._CUSTOM_REAL
+ Hs = 0._CUSTOM_REAL
+
+ ! three-zeta dependant terms
+ C3 = FACTOR_CS3_L*ll
+ S3 = 0._CUSTOM_REAL
+
+ ! four-zeta dependant terms
+ Ec = FACTOR_E_N*nn
+ Es = 0._CUSTOM_REAL
+ endif
+
+
+! The mapping from the elastic coefficients to the elastic tensor elements
+! in the local Cartesian coordinate system (classical geographic) used in the
+! global code (1---South, 2---East, 3---up)
+! Always keep the following part when you modify this subroutine
+ d11 = A + Ec + Bc
+ d12 = A - 2.*AN - Ec
+ d13 = F + Hc
+ d14 = S3 + 2.*S1sh + 2.*S1p
+ d15 = 2.*C1p + C3
+ d16 = -Bs/2. - Es
+ d22 = A + Ec - Bc
+ d23 = F - Hc
+ d24 = 2.*S1p - S3
+ d25 = 2.*C1p - 2.*C1sh - C3
+ d26 = -Bs/2. + Es
+ d33 = C
+ d34 = 2.*(S1p - S1sv)
+ d35 = 2.*(C1p - C1sv)
+ d36 = -Hs
+ d44 = AL - Gc
+ d45 = -Gs
+ d46 = C1sh - C3
+ d55 = AL + Gc
+ d56 = S3 - S1sh
+ d66 = AN - Ec
+
+! The mapping to the global Cartesian coordinate system used in the code
+! (1---East, 2---North, 3---up)
+ c11 = d22
+ c12 = d12
+ c13 = d23
+ c14 = - d25
+ c15 = d24
+ c16 = - d26
+ c22 = d11
+ c23 = d13
+ c24 = - d15
+ c25 = d14
+ c26 = - d16
+ c33 = d33
+ c34 = - d35
+ c35 = d34
+ c36 = - d36
+ c44 = d55
+ c45 = - d45
+ c46 = d56
+ c55 = d44
+ c56 = - d46
+ c66 = d66
+
+ end subroutine model_aniso
+
Added: seismo/3D/SPECFEM3D_SESAME/trunk/model_external_values.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/model_external_values.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/model_external_values.f90 2010-02-25 23:07:29 UTC (rev 16343)
@@ -0,0 +1,211 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+! generic model file
+!
+! note: the idea is to super-impose velocity model values on the GLL points,
+! additional to the ones assigned on the CUBIT mesh
+
+ module external_model
+
+!---
+!
+! ADD YOUR MODEL HERE
+!
+!---
+
+ ! only a dummy variable used as example
+ integer :: idummy_size = 1
+
+ ! type model_external_variables
+ ! sequence
+ ! double precision dvs(0:dummy_size)
+ ! end type model_external_variables
+ !
+ ! type (model_external_variables) MEXT_V
+
+ end module external_model
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine model_external_values(i,j,k,ispec,idomain_id,imaterial_id,&
+ nspec,ibool, &
+ iflag_aniso,iflag_atten, &
+ rho,vp,vs, &
+ c11,c12,c13,c14,c15,c16, &
+ c22,c23,c24,c25,c26,c33, &
+ c34,c35,c36,c44,c45,c46, &
+ c55,c56,c66,ANISOTROPY)
+
+! given a GLL point, returns super-imposed velocity model values
+
+ use external_model
+ use create_regions_mesh_ext_par
+
+ implicit none
+
+ ! GLL point indices
+ integer :: i,j,k,ispec
+
+ ! acoustic/elastic/.. domain flag ( 1 = acoustic / 2 = elastic / ... )
+ integer :: idomain_id
+
+ ! associated material flag (in cubit, this would be the volume id number)
+ integer :: imaterial_id
+
+ ! local-to-global index array
+ integer :: nspec
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+ ! anisotropy flag
+ integer :: iflag_aniso
+
+ ! attenuation flag
+ integer :: iflag_atten
+
+ ! density, Vp and Vs
+ real(kind=CUSTOM_REAL) :: vp,vs,rho
+
+ ! all anisotropy coefficients
+ real(kind=CUSTOM_REAL) :: c11,c12,c13,c14,c15,c16,c22,c23,c24,c25, &
+ c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+ logical :: ANISOTROPY
+
+ ! local parameters
+ real(kind=CUSTOM_REAL) :: x,y,z
+ real(kind=CUSTOM_REAL) :: xmin,xmax,ymin,ymax,zmin,zmax
+ real(kind=CUSTOM_REAL) :: depth
+ integer :: iglob
+
+!---
+!
+! ADD YOUR MODEL HERE
+!
+!---
+
+ ! GLL point location
+ iglob = ibool(i,j,k,ispec)
+ x = xstore_dummy(iglob)
+ y = ystore_dummy(iglob)
+ z = zstore_dummy(iglob)
+
+ ! model dimensions
+ xmin = minval(xstore_dummy)
+ xmax = maxval(xstore_dummy)
+ ymin = minval(ystore_dummy)
+ ymax = maxval(ystore_dummy)
+ zmin = minval(zstore_dummy)
+ zmax = maxval(zstore_dummy)
+
+ ! depth in Z-direction
+ depth = zmax - z
+
+ ! normalizes depth between 0 and 1
+ if( abs( zmax - zmin ) > TINYVAL ) depth = depth / (zmax - zmin)
+
+
+ ! super-imposes values
+ !rho = 2.6910d0+0.6924d0*depth
+ !vp = 4.1875d0+3.9382d0*depth
+ !vs = 2.1519d0+2.3481d0*depth
+
+ ! adds a velocity depth gradient
+ ! (e.g. from PREM mantle gradients:
+ ! vp : 3.9382*6371/5.5
+ ! vs : 2.3481*6371/5.5
+ ! rho : 0.6924*6371/5.5 )
+ rho = rho + 802.d0 * depth
+ vp = vp + 4562.d0 * depth
+ vs = vs + 2720.d0 * depth
+
+ ! adds anisotropic velocity values
+ if( ANISOTROPY ) &
+ call model_aniso(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
+ c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45, &
+ c46,c55,c56,c66)
+
+ ! to avoid compiler warnings
+ idummy_size = imaterial_id
+ idummy_size = idomain_id
+ idummy_size = iflag_atten
+
+ end subroutine model_external_values
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine model_external_broadcast(myrank)
+
+! standard routine to setup model
+
+ use external_model
+
+ implicit none
+
+ include "constants.h"
+ ! standard include of the MPI library
+ include 'mpif.h'
+
+ integer :: myrank
+
+ ! local parameters
+ integer :: ier
+
+ ! dummy to ignore compiler warnings
+ ier = myrank
+
+!---
+!
+! ADD YOUR MODEL HERE
+!
+!---
+
+ ! the variables read are declared and stored in structure D3MM_V
+ !if(myrank == 0) call read_external_model()
+
+ ! broadcast the information read on the master to the nodes
+ !call MPI_BCAST(MEXT_V%dvs,size(MEXT_V%dvs),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+
+ end subroutine model_external_broadcast
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+!
+! subroutine read_external_model()
+!
+! use external_model
+!
+! implicit none
+!
+! include "constants.h"
+!
+! end subroutine read_external_model
+
\ No newline at end of file
Added: seismo/3D/SPECFEM3D_SESAME/trunk/save_moho_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/save_moho_arrays.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/save_moho_arrays.f90 2010-02-25 23:07:29 UTC (rev 16343)
@@ -0,0 +1,330 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+ subroutine save_moho_arrays( myrank,nglob,nspec, &
+ nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, &
+ nodes_coords_ext_mesh,nnodes_ext_mesh,ibool )
+
+ use create_regions_mesh_ext_par
+ implicit none
+
+ integer :: nspec2D_moho_ext
+ integer, dimension(nspec2D_moho_ext) :: ibelm_moho
+ integer, dimension(4,nspec2D_moho_ext) :: nodes_ibelm_moho
+
+ integer :: myrank,nglob,nspec
+
+ ! data from the external mesh
+ integer :: nnodes_ext_mesh
+ double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+! local parameters
+ ! Moho mesh
+ real(CUSTOM_REAL), dimension(:,:,:),allocatable :: normal_moho_top
+ real(CUSTOM_REAL), dimension(:,:,:),allocatable :: normal_moho_bot
+ integer,dimension(:,:,:),allocatable :: ijk_moho_top, ijk_moho_bot
+ integer,dimension(:),allocatable :: ibelm_moho_top, ibelm_moho_bot
+ integer :: NSPEC2D_MOHO
+ logical, dimension(:),allocatable :: is_moho_top, is_moho_bot
+
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
+ real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
+ real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
+ real(kind=CUSTOM_REAL),dimension(NDIM):: normal
+ integer :: ijk_face(3,NGLLX,NGLLY)
+
+ real(kind=CUSTOM_REAL),dimension(:,:),allocatable :: iglob_normals
+ integer,dimension(:),allocatable:: iglob_is_surface
+
+ integer :: imoho_bot,imoho_top
+ integer :: ispec2D,ispec,icorner,iface,i,j,k,igll,iglob
+ integer :: iglob_midpoint,idirect,counter
+ integer :: imoho_top_all,imoho_bot_all,imoho_all
+
+ ! corners indices of reference cube faces
+ integer,dimension(3,4),parameter :: iface1_corner_ijk = &
+ reshape( (/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ /),(/3,4/)) ! xmin
+ integer,dimension(3,4),parameter :: iface2_corner_ijk = &
+ reshape( (/ NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),(/3,4/)) ! xmax
+ integer,dimension(3,4),parameter :: iface3_corner_ijk = &
+ reshape( (/ 1,1,1, 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,1,1 /),(/3,4/)) ! ymin
+ integer,dimension(3,4),parameter :: iface4_corner_ijk = &
+ reshape( (/ 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! ymax
+ integer,dimension(3,4),parameter :: iface5_corner_ijk = &
+ reshape( (/ 1,1,1, 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,1,1 /),(/3,4/)) ! bottom
+ integer,dimension(3,4),parameter :: iface6_corner_ijk = &
+ reshape( (/ 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! top
+ integer,dimension(3,4,6),parameter :: iface_all_corner_ijk = &
+ reshape( (/ iface1_corner_ijk,iface2_corner_ijk, &
+ iface3_corner_ijk,iface4_corner_ijk, &
+ iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
+ ! midpoint indices for each face (xmin,xmax,ymin,ymax,zmin,zmax)
+ integer,dimension(3,6),parameter :: iface_all_midpointijk = &
+ reshape( (/ 1,2,2, NGLLX,2,2, 2,1,2, 2,NGLLY,2, 2,2,1, 2,2,NGLLZ /),(/3,6/)) ! top
+
+ ! temporary arrays for passing information
+ allocate(iglob_is_surface(nglob))
+ allocate(iglob_normals(NDIM,nglob))
+ iglob_is_surface = 0
+ iglob_normals = 0._CUSTOM_REAL
+
+ ! loops over given moho surface elements
+ do ispec2D=1, nspec2D_moho_ext
+
+ ! gets element id
+ ispec = ibelm_moho(ispec2D)
+
+ ! looks for i,j,k indices of GLL points on boundary face
+ ! determines element face by given CUBIT corners
+ ! (note: uses point locations rather than point indices to find the element face,
+ ! because the indices refer no more to the newly indexed ibool array )
+ do icorner=1,NGNOD2D
+ xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_moho(icorner,ispec2D))
+ ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_moho(icorner,ispec2D))
+ zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_moho(icorner,ispec2D))
+ enddo
+
+ ! sets face id of reference element associated with this face
+ call get_element_face_id(ispec,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ iface)
+
+ ! ijk indices of GLL points for face id
+ call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ)
+
+ ! weighted jacobian and normal
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
+
+ ! normal convention: points away from element
+ ! switch normal direction if necessary
+ do j=1,NGLLY
+ do i=1,NGLLX
+ call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ enddo
+ enddo
+
+ ! stores information on global points on moho surface
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(ijk_face(1,i,j),ijk_face(2,i,j),ijk_face(3,i,j),ispec)
+ ! sets flag
+ iglob_is_surface(iglob) = ispec2D
+ ! sets normals
+ iglob_normals(:,iglob) = normal_face(:,i,j)
+ enddo
+ enddo
+ enddo
+
+ ! stores moho elements
+ NSPEC2D_MOHO = nspec2D_moho_ext
+
+ allocate(ibelm_moho_bot(NSPEC2D_MOHO))
+ allocate(ibelm_moho_top(NSPEC2D_MOHO))
+ allocate(normal_moho_top(NDIM,NGLLSQUARE,NSPEC2D_MOHO))
+ allocate(normal_moho_bot(NDIM,NGLLSQUARE,NSPEC2D_MOHO))
+ allocate(ijk_moho_bot(3,NGLLSQUARE,NSPEC2D_MOHO))
+ allocate(ijk_moho_top(3,NGLLSQUARE,NSPEC2D_MOHO))
+ ibelm_moho_bot = 0
+ ibelm_moho_top = 0
+
+ ! element flags
+ allocate(is_moho_top(nspec))
+ allocate(is_moho_bot(nspec))
+ is_moho_top = .false.
+ is_moho_bot = .false.
+
+ ! finds spectral elements with moho surface
+ imoho_top = 0
+ imoho_bot = 0
+ do ispec=1,nspec
+
+ ! loops over each face
+ do iface = 1,6
+ ! checks if corners of face on surface
+ counter = 0
+ do icorner = 1,NGNOD2D
+ i = iface_all_corner_ijk(1,icorner,iface)
+ j = iface_all_corner_ijk(2,icorner,iface)
+ k = iface_all_corner_ijk(3,icorner,iface)
+ iglob = ibool(i,j,k,ispec)
+
+ ! checks if point on surface
+ if( iglob_is_surface(iglob) > 0 ) then
+ counter = counter+1
+
+ ! reference corner coordinates
+ xcoord(icorner) = xstore_dummy(iglob)
+ ycoord(icorner) = ystore_dummy(iglob)
+ zcoord(icorner) = zstore_dummy(iglob)
+ endif
+ enddo
+
+ ! stores moho informations
+ if( counter == NGNOD2D ) then
+
+ ! gets face GLL points i,j,k indices from element face
+ call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY)
+
+ ! re-computes face infos
+ ! weighted jacobian and normal
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
+
+ ! normal convention: points away from element
+ ! switch normal direction if necessary
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ enddo
+ enddo
+
+ ! takes normal stored temporary on a face midpoint
+ i = iface_all_midpointijk(1,iface)
+ j = iface_all_midpointijk(2,iface)
+ k = iface_all_midpointijk(3,iface)
+ iglob_midpoint = ibool(i,j,k,ispec)
+ normal(:) = iglob_normals(:,iglob_midpoint)
+
+ ! determines whether normal points into element or not (top/bottom distinction)
+ call get_element_face_normal_idirect(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal,idirect )
+
+ ! takes moho surface element id given by id on midpoint
+ ispec2D = iglob_is_surface(iglob_midpoint)
+
+ ! sets face infos for bottom (normal points away from element)
+ if( idirect == 1 ) then
+
+ ! checks validity
+ if( is_moho_bot( ispec) .eqv. .true. ) then
+ print*,'error: moho surface geometry bottom'
+ print*,' does not allow for mulitple element faces in kernel computation'
+ call exit_mpi(myrank,'error moho bottom elements')
+ endif
+
+ imoho_bot = imoho_bot + 1
+ is_moho_bot(ispec) = .true.
+ ibelm_moho_bot(ispec2D) = ispec
+
+ ! stores on surface gll points (assuming NGLLX = NGLLY = NGLLZ)
+ igll = 0
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ igll = igll+1
+ ijk_moho_bot(:,igll,ispec2D) = ijk_face(:,i,j)
+ normal_moho_bot(:,igll,ispec2D) = normal_face(:,i,j)
+ enddo
+ enddo
+
+ ! sets face infos for top element
+ else if( idirect == 2 ) then
+
+ ! checks validity
+ if( is_moho_top( ispec) .eqv. .true. ) then
+ print*,'error: moho surface geometry top'
+ print*,' does not allow for mulitple element faces kernel computation'
+ call exit_mpi(myrank,'error moho top elements')
+ endif
+
+ imoho_top = imoho_top + 1
+ is_moho_top(ispec) = .true.
+ ibelm_moho_top(ispec2D) = ispec
+
+ ! gll points
+ igll = 0
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ igll = igll+1
+ ijk_moho_top(:,igll,ispec) = ijk_face(:,i,j)
+ ! note: top elements have normal pointing into element
+ normal_moho_top(:,igll,ispec) = - normal_face(:,i,j)
+ enddo
+ enddo
+ endif
+
+ endif ! counter
+
+ enddo ! iface
+
+ ! checks validity of top/bottom distinction
+ if( is_moho_top(ispec) .and. is_moho_bot(ispec) ) then
+ print*,'error: moho surface elements confusing'
+ print*,' element:',ispec,'has top and bottom surface'
+ call exit_mpi(myrank,'error moho surface element')
+ endif
+
+ enddo ! ispec2D
+
+ ! note: surface e.g. could be at the free-surface and have no top elements etc...
+ ! user output
+ call sum_all_i( imoho_top, imoho_top_all )
+ call sum_all_i( imoho_bot, imoho_bot_all )
+ call sum_all_i( NSPEC2D_MOHO, imoho_all )
+ if( myrank == 0 ) then
+ write(IMAIN,*) '********'
+ write(IMAIN,*) 'Moho surface:'
+ write(IMAIN,*) ' total surface elements: ',imoho_all
+ write(IMAIN,*) ' top elements :',imoho_top_all
+ write(IMAIN,*) ' bottom elements:',imoho_bot_all
+ write(IMAIN,*) '********'
+ endif
+
+ ! saves moho files: total number of elements, corner points, all points
+ open(unit=27,file=prname(1:len_trim(prname))//'ibelm_moho.bin',status='unknown',form='unformatted')
+ write(27) NSPEC2D_MOHO
+ write(27) ibelm_moho_top
+ write(27) ibelm_moho_bot
+ write(27) ijk_moho_top
+ write(27) ijk_moho_bot
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'normal_moho.bin',status='unknown',form='unformatted')
+ write(27) normal_moho_top
+ write(27) normal_moho_bot
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'is_moho.bin',status='unknown',form='unformatted')
+ write(27) is_moho_top
+ write(27) is_moho_bot
+ close(27)
+
+end subroutine save_moho_arrays
More information about the CIG-COMMITS
mailing list