[cig-commits] r16434 - in seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D: . DATA
pieyre at geodynamics.org
pieyre at geodynamics.org
Fri Mar 19 13:11:02 PDT 2010
Author: pieyre
Date: 2010-03-19 13:11:02 -0700 (Fri, 19 Mar 2010)
New Revision: 16434
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/DATA/Par_file
seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/Makefile.in
seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/compute_parameters.f90
seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/constants.h.in
seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/create_regions_mesh.f90
seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/define_subregions.f90
seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/meshfem3D.f90
seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/read_parameter_file.f90
seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/save_databases.f90
Log:
added doublings layers option in internal mesher. Not validated yet.
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/DATA/Par_file
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/DATA/Par_file 2010-03-19 05:07:35 UTC (rev 16433)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/DATA/Par_file 2010-03-19 20:11:02 UTC (rev 16434)
@@ -21,6 +21,13 @@
NPROC_XI = 2
NPROC_ETA = 1
+# Regular/irregular mesh
+USE_REGULAR_MESH = .false.
+# Only for irregular meshes, number of doubling layers (1 or 2) and their position
+NDOUBLINGS = 2
+NZ_DOUGLING_1 = 8
+NZ_DOUGLING_2 = 4
+
# path to store the databases files
LOCAL_PATH = OUTPUT_FILES
@@ -31,7 +38,7 @@
# Q_flag : 0=no attenuation/1=IATTENUATION_SEDIMENTS_40, 2=..., 13=IATTENUATION_BEDROCK
# anisotropy_flag : 0=no anisotropy/ 1,2,.. check with implementation in aniso_model.f90
# domain_id : 1=acoustic / 2=elastic / 3=poroelastic
-1 1200 1500 750 1 0 1
+1 1200 1500 750 1 0 2
2 1100 1600 800 3 0 1
3 1000 1500 0 0 0 2
4 1300 1400 700 2 0 3
@@ -39,8 +46,8 @@
NREGIONS = 4
# define the different regions of the model as :
#NEX_XI_BEGIN #NEX_XI_END #NEX_ETA_BEGIN #NEX_ETA_END #NZ_BEGIN #NZ_END #material_id
-1 16 1 16 1 3 1
-1 16 1 16 4 7 2
-1 16 1 16 8 20 3
-14 16 9 16 6 20 4
+1 16 1 16 1 2 1
+1 16 1 16 3 3 2
+1 16 1 16 4 20 3
+14 16 7 16 7 10 4
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/Makefile.in 2010-03-19 05:07:35 UTC (rev 16433)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/Makefile.in 2010-03-19 20:11:02 UTC (rev 16434)
@@ -83,6 +83,7 @@
$O/create_regions_mesh.o \
$O/define_subregions.o \
$O/define_subregions_heuristic.o \
+ $O/define_superbrick.o \
$O/exit_mpi.o \
$O/get_MPI_cutplanes_eta.o \
$O/get_MPI_cutplanes_xi.o \
@@ -262,3 +263,6 @@
$O/define_subregions_heuristic.o: constants.h define_subregions_heuristic.f90
${FCCOMPILE_CHECK} -c -o $O/define_subregions_heuristic.o define_subregions_heuristic.f90
+
+$O/define_superbrick.o: constants.h define_superbrick.f90
+ ${FCCOMPILE_CHECK} -c -o $O/define_superbrick.o define_superbrick.f90
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/compute_parameters.f90 2010-03-19 05:07:35 UTC (rev 16433)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/compute_parameters.f90 2010-03-19 20:11:02 UTC (rev 16434)
@@ -23,29 +23,36 @@
!
!=====================================================================
- subroutine compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
- NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
- NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM, &
- NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,USE_REGULAR_MESH)
+subroutine compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
+ NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM, &
+ NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,&
+ USE_REGULAR_MESH,NDOUBLINGS,ner_doublings)
implicit none
include "constants.h"
-! parameters read from parameter file
+ ! parameters read from parameter file
integer NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA
integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM
+ integer NER_REGULAR1,NER_REGULAR2,NER_REGULAR3
+ integer NGLOB_DOUBLING,NGLOB_DOUBLING_1,NGLOB_DOUBLING_2,NGLOB_NO_DOUBLING
+ integer NDOUBLINGS
+ integer ner_doublings(2)
-! parameters to be computed based upon parameters above read from file
+ ! parameters to be computed based upon parameters above read from file
integer NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NER
+ integer NEX_DOUBLING_ABOVE_PER_PROC_XI,NEX_DOUBLING_ABOVE_PER_PROC_ETA
+ integer NEX_DOUBLING_ABOVE_XI,NEX_DOUBLING_ABOVE_ETA
integer NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB
integer NEX_DOUBLING_SEDIM_XI,NEX_DOUBLING_SEDIM_ETA
integer NEX_DOUBLING_SEDIM_PER_PROC_XI,NEX_DOUBLING_SEDIM_PER_PROC_ETA
@@ -63,200 +70,291 @@
logical USE_REGULAR_MESH
-!
-!--- case of a regular mesh
-!
- if(USE_REGULAR_MESH) then
-! number of elements horizontally in each slice (i.e. per processor)
-! these two values MUST be equal in all cases
+ ! number of elements horizontally in each slice (i.e. per processor)
+ ! these two values MUST be equal in all cases
NEX_PER_PROC_XI = NEX_XI / NPROC_XI
NEX_PER_PROC_ETA = NEX_ETA / NPROC_ETA
-! total number of processors in each of the six chunks
+ ! total number of processors in each of the six chunks
NPROC = NPROC_XI * NPROC_ETA
-! exact number of spectral elements without doubling layers
- NSPEC_NO_DOUBLING = NEX_XI*NEX_ETA*NER
+ !
+ !--- case of a regular mesh
+ !
+ if(USE_REGULAR_MESH) then
-! %%%%%%%%%%%%%% surface elements %%%%%%%%%%%%%%%%%%%
+ ! exact number of spectral elements without doubling layers
+ NSPEC_NO_DOUBLING = NEX_XI*NEX_ETA*NER
-! exact number of surface elements for a chunk without doubling layers
+ ! %%%%%%%%%%%%%% surface elements %%%%%%%%%%%%%%%%%%%
- NSPEC2D_NO_DOUBLING_XI = NEX_PER_PROC_XI*NER
+ ! exact number of surface elements for a chunk without doubling layers
- NSPEC2D_NO_DOUBLING_ETA = NEX_PER_PROC_ETA*NER
+ NSPEC2D_NO_DOUBLING_XI = NEX_PER_PROC_XI*NER
-! exact number of spectral elements
-! NSPEC_AB = NSPEC_NO_DOUBLING / NPROC
- NSPEC_AB = NER * NEX_PER_PROC_XI * NEX_PER_PROC_ETA
+ NSPEC2D_NO_DOUBLING_ETA = NEX_PER_PROC_ETA*NER
-! exact number of surface elements for faces A and B along XI and ETA
- NSPEC2D_A_XI = NSPEC2D_NO_DOUBLING_XI
- NSPEC2D_B_XI = NSPEC2D_NO_DOUBLING_XI
- NSPEC2D_A_ETA = NSPEC2D_NO_DOUBLING_ETA
- NSPEC2D_B_ETA = NSPEC2D_NO_DOUBLING_ETA
+ ! exact number of spectral elements
+ ! NSPEC_AB = NSPEC_NO_DOUBLING / NPROC
+ NSPEC_AB = NER * NEX_PER_PROC_XI * NEX_PER_PROC_ETA
-! exact number of surface elements on the bottom and top boundaries
-! and theoretical number of spectral elements in radial direction
+ ! exact number of surface elements for faces A and B along XI and ETA
+ NSPEC2D_A_XI = NSPEC2D_NO_DOUBLING_XI
+ NSPEC2D_B_XI = NSPEC2D_NO_DOUBLING_XI
+ NSPEC2D_A_ETA = NSPEC2D_NO_DOUBLING_ETA
+ NSPEC2D_B_ETA = NSPEC2D_NO_DOUBLING_ETA
- NSPEC2D_TOP = NEX_XI*NEX_ETA / NPROC
- NSPEC2D_BOTTOM = NSPEC2D_TOP
+ ! exact number of surface elements on the bottom and top boundaries
+ ! and theoretical number of spectral elements in radial direction
- NSPEC1D_RADIAL_BEDROCK = NER
+ NSPEC2D_TOP = NEX_XI*NEX_ETA / NPROC
+ NSPEC2D_BOTTOM = NSPEC2D_TOP
-! face with max number of elements is type B here
-! maximum number of surface elements on vertical boundaries of the slices
- NSPEC2DMAX_XMIN_XMAX = NSPEC2D_B_ETA
- NSPEC2DMAX_YMIN_YMAX = NSPEC2D_B_XI
+ NSPEC1D_RADIAL_BEDROCK = NER
-! theoretical number of Gauss-Lobatto points in radial direction
- NPOIN1D_RADIAL_BEDROCK = NSPEC1D_RADIAL_BEDROCK*(NGLLZ-1)+1
+ ! face with max number of elements is type B here
+ ! maximum number of surface elements on vertical boundaries of the slices
+ NSPEC2DMAX_XMIN_XMAX = NSPEC2D_B_ETA
+ NSPEC2DMAX_YMIN_YMAX = NSPEC2D_B_XI
-! 2-D addressing and buffers for summation between slices
-! we add one to number of points because of the flag after the last point
- NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX*NGLLY*NGLLZ + 1
- NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX*NGLLX*NGLLZ + 1
+ ! theoretical number of Gauss-Lobatto points in radial direction
+ NPOIN1D_RADIAL_BEDROCK = NSPEC1D_RADIAL_BEDROCK*(NGLLZ-1)+1
-! exact number of global points
- NGLOB_AB = (NEX_PER_PROC_XI*(NGLLX-1)+1) * (NEX_PER_PROC_ETA*(NGLLY-1)+1) * (NER*(NGLLZ-1)+1)
+ ! 2-D addressing and buffers for summation between slices
+ ! we add one to number of points because of the flag after the last point
+ NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX*NGLLY*NGLLZ + 1
+ NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX*NGLLX*NGLLZ + 1
-!
-!--- case of a non-regular mesh with mesh doublings
-!
- else
+ ! exact number of global points
+ NGLOB_AB = (NEX_PER_PROC_XI*(NGLLX-1)+1) * (NEX_PER_PROC_ETA*(NGLLY-1)+1) * (NER*(NGLLZ-1)+1)
-! total number of spectral elements along radius
- NER = NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT + NER_BASEMENT_SEDIM + NER_SEDIM
+ !
+ !--- case of a non-regular mesh with mesh doublings
+ !
+ else if(NDOUBLINGS == 1) then
-! number of elements horizontally in each slice (i.e. per processor)
-! these two values MUST be equal in all cases
- NEX_PER_PROC_XI = NEX_XI / NPROC_XI
- NEX_PER_PROC_ETA = NEX_ETA / NPROC_ETA
+ ! number of spectral elements at the bottom of the doubling below the moho
+ NEX_DOUBLING_ABOVE_XI=NEX_XI/2
+ NEX_DOUBLING_ABOVE_ETA=NEX_ETA/2
+ NEX_DOUBLING_ABOVE_PER_PROC_XI=NEX_PER_PROC_XI/2
+ NEX_DOUBLING_ABOVE_PER_PROC_ETA=NEX_PER_PROC_ETA/2
-! total number of processors in each of the six chunks
- NPROC = NPROC_XI * NPROC_ETA
+ ! exact number of spectral elements without doubling layers
+ NER_REGULAR2 = NER - ner_doublings(1)
+ NER_REGULAR1 = ner_doublings(1) - 2
-! number of spectral elements at the bottom of the doubling below the moho
- NEX_DOUBLING_SEDIM_XI=NEX_XI/2
- NEX_DOUBLING_SEDIM_ETA=NEX_ETA/2
- NEX_DOUBLING_SEDIM_PER_PROC_XI=NEX_PER_PROC_XI/2
- NEX_DOUBLING_SEDIM_PER_PROC_ETA=NEX_PER_PROC_ETA/2
+ NSPEC_NO_DOUBLING = &
+ (NEX_XI/2)*(NEX_ETA/2)*NER_REGULAR1 &
+ + NEX_XI*NEX_ETA*NER_REGULAR2
-! exact number of spectral elements without doubling layers
- NSPEC_NO_DOUBLING = &
- (NEX_DOUBLING_SEDIM_XI*NEX_DOUBLING_SEDIM_ETA*(NER_BASEMENT_SEDIM/2-3) &
- +(NEX_XI/4)*(NEX_ETA/4)*(NER_16_BASEMENT/2-3) &
- +(NEX_XI/4)*(NEX_ETA/4)*(NER_MOHO_16/2) &
- +(NEX_XI/4)*(NEX_ETA/4)*(NER_BOTTOM_MOHO/4)) + NEX_XI*NEX_ETA*NER_SEDIM
+ ! exact number of spectral elements in the doubling regions
-! exact number of spectral elements in the doubling regions
+ ! number of elementary bricks in the two regions with doubling
+ NUM_DOUBLING_BRICKS = (NEX_DOUBLING_ABOVE_XI*NEX_DOUBLING_ABOVE_ETA)/4
-! number of elementary bricks in the two regions with doubling
- NUM_DOUBLING_BRICKS = ((NEX_XI/4)*(NEX_ETA/4) &
- +NEX_DOUBLING_SEDIM_XI*NEX_DOUBLING_SEDIM_ETA)/4
+ ! for type AB, each doubling brick contains 32 elements on 2 levels
+ NSPEC_DOUBLING_AB=32*NUM_DOUBLING_BRICKS
-! for type AB, each doubling brick contains 40 elements on 3 levels
- NSPEC_DOUBLING_AB=40*NUM_DOUBLING_BRICKS
+ ! %%%%%%%%%%%%%% surface elements %%%%%%%%%%%%%%%%%%%
-! %%%%%%%%%%%%%% surface elements %%%%%%%%%%%%%%%%%%%
+ ! exact number of surface elements for a chunk without doubling layers
-! exact number of surface elements for a chunk without doubling layers
+ NSPEC2D_NO_DOUBLING_XI = &
+ (NEX_PER_PROC_XI/2)*NER_REGULAR1 &
+ + NEX_PER_PROC_XI*NER_REGULAR2
- NSPEC2D_NO_DOUBLING_XI = &
- NEX_DOUBLING_SEDIM_PER_PROC_XI*(NER_BASEMENT_SEDIM/2-3) &
- +(NEX_PER_PROC_XI/4)*(NER_16_BASEMENT/2-3) &
- +(NEX_PER_PROC_XI/4)*(NER_MOHO_16/2) &
- +(NEX_PER_PROC_XI/4)*(NER_BOTTOM_MOHO/4) + NEX_PER_PROC_XI*NER_SEDIM
+ NSPEC2D_NO_DOUBLING_ETA = &
+ (NEX_PER_PROC_ETA/2)*NER_REGULAR1 &
+ +NEX_PER_PROC_ETA*NER_REGULAR2
- NSPEC2D_NO_DOUBLING_ETA = &
- NEX_DOUBLING_SEDIM_PER_PROC_ETA*(NER_BASEMENT_SEDIM/2-3) &
- +(NEX_PER_PROC_ETA/4)*(NER_16_BASEMENT/2-3) &
- +(NEX_PER_PROC_ETA/4)*(NER_MOHO_16/2) &
- +(NEX_PER_PROC_ETA/4)*(NER_BOTTOM_MOHO/4) + NEX_PER_PROC_ETA*NER_SEDIM
+ ! exact number of surface elements in the doubling regions
-! exact number of surface elements in the doubling regions
+ ! number of elementary bricks in the two regions with doubling
+ NUM2D_DOUBLING_BRICKS_XI = (NEX_PER_PROC_XI/4 &
+ + NEX_DOUBLING_ABOVE_PER_PROC_XI)/2
-! number of elementary bricks in the two regions with doubling
- NUM2D_DOUBLING_BRICKS_XI = ((NEX_PER_PROC_XI/4) &
- +NEX_DOUBLING_SEDIM_PER_PROC_XI)/2
+ NUM2D_DOUBLING_BRICKS_ETA = (NEX_PER_PROC_ETA/2 &
+ + NEX_DOUBLING_ABOVE_PER_PROC_ETA)/2
- NUM2D_DOUBLING_BRICKS_ETA = ((NEX_PER_PROC_ETA/4) &
- +NEX_DOUBLING_SEDIM_PER_PROC_ETA)/2
+ ! for type A, each doubling brick contains 10 elements on 3 levels
+ NSPEC2D_DOUBLING_A_XI=10*NUM2D_DOUBLING_BRICKS_XI
+ NSPEC2D_DOUBLING_A_ETA=10*NUM2D_DOUBLING_BRICKS_ETA
-! for type A, each doubling brick contains 10 elements on 3 levels
- NSPEC2D_DOUBLING_A_XI=10*NUM2D_DOUBLING_BRICKS_XI
- NSPEC2D_DOUBLING_A_ETA=10*NUM2D_DOUBLING_BRICKS_ETA
+ ! for type B, each doubling brick contains 12 elements on 3 levels
+ NSPEC2D_DOUBLING_B_XI=10*NUM2D_DOUBLING_BRICKS_XI
+ NSPEC2D_DOUBLING_B_ETA=10*NUM2D_DOUBLING_BRICKS_ETA
-! for type B, each doubling brick contains 12 elements on 3 levels
- NSPEC2D_DOUBLING_B_XI=12*NUM2D_DOUBLING_BRICKS_XI
- NSPEC2D_DOUBLING_B_ETA=12*NUM2D_DOUBLING_BRICKS_ETA
+ ! exact number of spectral elements
+ NSPEC_AB = (NSPEC_NO_DOUBLING + NSPEC_DOUBLING_AB) / NPROC
-! exact number of spectral elements
- NSPEC_AB = (NSPEC_NO_DOUBLING + NSPEC_DOUBLING_AB) / NPROC
+ ! exact number of surface elements for faces A and B
+ ! along XI and ETA for doubling region
+ NSPEC2D_A_XI = NSPEC2D_NO_DOUBLING_XI + NSPEC2D_DOUBLING_A_XI
+ NSPEC2D_B_XI = NSPEC2D_NO_DOUBLING_XI + NSPEC2D_DOUBLING_B_XI
+ NSPEC2D_A_ETA = NSPEC2D_NO_DOUBLING_ETA + NSPEC2D_DOUBLING_A_ETA
+ NSPEC2D_B_ETA = NSPEC2D_NO_DOUBLING_ETA + NSPEC2D_DOUBLING_B_ETA
-! exact number of surface elements for faces A and B
-! along XI and ETA for doubling region
- NSPEC2D_A_XI = NSPEC2D_NO_DOUBLING_XI + NSPEC2D_DOUBLING_A_XI
- NSPEC2D_B_XI = NSPEC2D_NO_DOUBLING_XI + NSPEC2D_DOUBLING_B_XI
- NSPEC2D_A_ETA = NSPEC2D_NO_DOUBLING_ETA + NSPEC2D_DOUBLING_A_ETA
- NSPEC2D_B_ETA = NSPEC2D_NO_DOUBLING_ETA + NSPEC2D_DOUBLING_B_ETA
+ ! exact number of surface elements on the bottom and top boundaries
+ ! and theoretical number of spectral elements in radial direction
-! exact number of surface elements on the bottom and top boundaries
-! and theoretical number of spectral elements in radial direction
+ NSPEC2D_TOP = NEX_XI*NEX_ETA / NPROC
+ NSPEC2D_BOTTOM = (NEX_XI/2)*(NEX_ETA/2) / NPROC
- NSPEC2D_TOP = NEX_XI*NEX_ETA / NPROC
- NSPEC2D_BOTTOM = (NEX_XI/4)*(NEX_ETA/4) / NPROC
+ ! face with max number of elements is type B here
+ ! maximum number of surface elements on vertical boundaries of the slices
+ NSPEC2DMAX_XMIN_XMAX = NSPEC2D_B_ETA
+ NSPEC2DMAX_YMIN_YMAX = NSPEC2D_B_XI
- NSPEC1D_RADIAL_BEDROCK = (NER_BASEMENT_SEDIM+NER_16_BASEMENT+NER_MOHO_16)/2 + NER_BOTTOM_MOHO/4
+ ! theoretical number of Gauss-Lobatto points in radial direction
+ ! NPOIN1D_RADIAL_BEDROCK = NSPEC1D_RADIAL_BEDROCK*(NGLLZ-1)+1
-! face with max number of elements is type B here
-! maximum number of surface elements on vertical boundaries of the slices
- NSPEC2DMAX_XMIN_XMAX = NSPEC2D_B_ETA
- NSPEC2DMAX_YMIN_YMAX = NSPEC2D_B_XI
+ ! 2-D addressing and buffers for summation between slices
+ ! we add one to number of points because of the flag after the last point
+ NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX*2*2 + 1
+ NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX*2*2 + 1
-! theoretical number of Gauss-Lobatto points in radial direction
- NPOIN1D_RADIAL_BEDROCK = NSPEC1D_RADIAL_BEDROCK*(NGLLZ-1)+1
+ ! exact number of global points
-! 2-D addressing and buffers for summation between slices
-! we add one to number of points because of the flag after the last point
- NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX*NGLLY*NGLLZ + 1
- NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX*NGLLX*NGLLZ + 1
+ ! case of the doubling regions
-! exact number of global points
+ NGLOB_DOUBLING_1 = (NEX_PER_PROC_XI/4)*(NEX_PER_PROC_ETA/4)*(32*NGLLX**3 - 70*NGLLX**2 + 52*NGLLX - 13) &
+ - ((NEX_PER_PROC_XI/4-1)*(NEX_PER_PROC_ETA/4) + (NEX_PER_PROC_ETA/4-1)*(NEX_PER_PROC_XI/4)) * (8*NGLLX**2-11*NGLLX+4)&
+ + (NEX_PER_PROC_XI/4-1)*(NEX_PER_PROC_ETA/4-1)*(NGLLX+1)
-! case of the doubling regions
-! formulas computed using Mathematica for the basic 3D doubling brick
+ ! NGLOB_DOUBLING_1 = (NEX_PER_PROC_XI/2)*(NEX_PER_PROC_ETA/2)*(32*NGLLX**3 - 70*NGLLX**2 + 52*NGLLX - 13) &
+ ! - ((NEX_PER_PROC_XI/2-1)*(NEX_PER_PROC_ETA/2) + (NEX_PER_PROC_ETA/2-1)*(NEX_PER_PROC_XI/2)) * (8*NGLLX**2-11*NGLLX+4)&
+ ! + (NEX_PER_PROC_XI/2-1)*(NEX_PER_PROC_ETA/2-1)*NGLLX
-! compute number of points in blocks with no doubling
-! exclude the three surfaces in contact with the doubling regions
- nglob_no_doubling_volume = (4*(NGLLX-1)+1)*(4*(NGLLX-1)+1)*((NER_BASEMENT_SEDIM/2-3 )*(NGLLX-1)-1) &
- +(2*(NGLLX-1)+1)*(2*(NGLLX-1)+1)*(((NER_16_BASEMENT/2+NER_MOHO_16/2+NER_BOTTOM_MOHO/4)-3)*(NGLLX-1)+0)
+ ! NGLOB_DOUBLING_2 =(NEX_PER_PROC_XI/8)*(NEX_PER_PROC_ETA/8)*(32*NGLLX**3 - 70*NGLLX**2 + 52*NGLLX - 13) &
+ ! - ((NEX_PER_PROC_XI/8-1)*(NEX_PER_PROC_ETA/8) + (NEX_PER_PROC_ETA/8-1)*(NEX_PER_PROC_XI/8)) * (8*NGLLX**2-11*NGLLX+4)&
+ ! + (NEX_PER_PROC_XI/8-1)*(NEX_PER_PROC_ETA/8-1)*NGLLX
-! number of basic blocks in each slice
- nblocks_xi = NEX_PER_PROC_XI / 8
- nblocks_eta = NEX_PER_PROC_ETA / 8
+ NGLOB_DOUBLING = NGLOB_DOUBLING_1 !+ NGLOB_DOUBLING_2
- NGLOB_AB = nblocks_xi*nblocks_eta*(200*NGLLX**3 - 484*NGLLX**2 + 392*NGLLX - 106 + nglob_no_doubling_volume)
+ ! compute number of points in blocks with no doubling
+ ! exclude the four surfaces in contact with the doubling regions
-! same thing for 2D surfaces for the three types of faces
- nglob_no_doubling_surface = (4*(NGLLX-1)+1)*((NER_BASEMENT_SEDIM/2-3)*(NGLLX-1)-1) &
- +(2*(NGLLX-1)+1)*(((NER_16_BASEMENT/2+NER_MOHO_16/2+NER_BOTTOM_MOHO/4)-3)*(NGLLX-1)+0)
+ NGLOB_NO_DOUBLING = (NEX_PER_PROC_XI/2 + 1)*(NEX_PER_PROC_ETA/2 + 1)*NER_REGULAR1 &
+ + (NEX_PER_PROC_XI + 1)*(NEX_PER_PROC_ETA + 1)*NER_REGULAR2
- nglob_surface_typeA = 30*NGLLX**2 - 45 * NGLLX + 17
- nglob_surface_typeB = 36*NGLLX**2 - 57 * NGLLX + 23
+ NGLOB_AB = NGLOB_DOUBLING + NGLOB_NO_DOUBLING
-! final number of points in volume obtained by removing planes counted twice
- NGLOB_AB = NGLOB_AB &
- - (nblocks_xi-1)*nblocks_eta*(nglob_surface_typeA + nglob_no_doubling_surface) &
- - (nblocks_eta-1)*nblocks_xi*(nglob_surface_typeB + nglob_no_doubling_surface) &
- + (nblocks_eta-1)*(nblocks_xi-1)*NPOIN1D_RADIAL_BEDROCK
+ else if (NDOUBLINGS == 2) then
-! add number of points in the sediments
- NGLOB_AB = NGLOB_AB + (NEX_PER_PROC_XI*(NGLLX-1)+1) &
- *(NEX_PER_PROC_ETA*(NGLLY-1)+1)*(NER_SEDIM*(NGLLZ-1)+0)
+ ! number of spectral elements at the bottom of the doubling below the moho
+ NEX_DOUBLING_ABOVE_XI=NEX_XI/2
+ NEX_DOUBLING_ABOVE_ETA=NEX_ETA/2
+ NEX_DOUBLING_ABOVE_PER_PROC_XI=NEX_PER_PROC_XI/2
+ NEX_DOUBLING_ABOVE_PER_PROC_ETA=NEX_PER_PROC_ETA/2
+ ! exact number of spectral elements without doubling layers
+ NER_REGULAR3 = NER - ner_doublings(1)
+ NER_REGULAR2 = (ner_doublings(1) - 2) - ner_doublings(2)
+ NER_REGULAR1 = ner_doublings(2) - 2
+
+ NSPEC_NO_DOUBLING = &
+ (NEX_XI/4)*(NEX_ETA/4)*NER_REGULAR1 &
+ + (NEX_XI/2)*(NEX_ETA/2)*NER_REGULAR2 &
+ + NEX_XI*NEX_ETA*NER_REGULAR3
+
+ ! exact number of spectral elements in the doubling regions
+
+ ! number of elementary bricks in the two regions with doubling
+ NUM_DOUBLING_BRICKS = (NEX_DOUBLING_ABOVE_XI*NEX_DOUBLING_ABOVE_ETA)/4&
+ + (NEX_DOUBLING_ABOVE_XI*NEX_DOUBLING_ABOVE_ETA)/16
+
+ ! for type AB, each doubling brick contains 32 elements on 2 levels
+ NSPEC_DOUBLING_AB=32*NUM_DOUBLING_BRICKS
+
+
+ ! %%%%%%%%%%%%%% surface elements %%%%%%%%%%%%%%%%%%%
+
+ ! exact number of surface elements for a chunk without doubling layers
+
+ NSPEC2D_NO_DOUBLING_XI = &
+ (NEX_PER_PROC_XI/4)*NER_REGULAR1 &
+ + (NEX_PER_PROC_XI/2)*NER_REGULAR2 &
+ + NEX_PER_PROC_XI*NER_REGULAR3
+
+ NSPEC2D_NO_DOUBLING_ETA = &
+ (NEX_PER_PROC_ETA/4)*NER_REGULAR1 &
+ + (NEX_PER_PROC_ETA/4)*NER_REGULAR2 &
+ + NEX_PER_PROC_ETA*NER_REGULAR3
+
+ ! exact number of surface elements in the doubling regions
+
+ ! number of elementary bricks in the two regions with doubling
+ NUM2D_DOUBLING_BRICKS_XI = (NEX_PER_PROC_XI/2 &
+ + NEX_DOUBLING_ABOVE_PER_PROC_XI)/2
+
+ NUM2D_DOUBLING_BRICKS_ETA = (NEX_PER_PROC_ETA/2 &
+ + NEX_DOUBLING_ABOVE_PER_PROC_ETA)/2
+
+ ! for type A, each doubling brick contains 10 elements on 3 levels
+ NSPEC2D_DOUBLING_A_XI=7*NUM2D_DOUBLING_BRICKS_XI
+ NSPEC2D_DOUBLING_A_ETA=7*NUM2D_DOUBLING_BRICKS_ETA
+
+ ! for type B, each doubling brick contains 12 elements on 3 levels
+ NSPEC2D_DOUBLING_B_XI=10*NUM2D_DOUBLING_BRICKS_XI
+ NSPEC2D_DOUBLING_B_ETA=12*NUM2D_DOUBLING_BRICKS_ETA
+
+ ! exact number of spectral elements
+ NSPEC_AB = (NSPEC_NO_DOUBLING + NSPEC_DOUBLING_AB) / NPROC
+
+ ! exact number of surface elements for faces A and B
+ ! along XI and ETA for doubling region
+ NSPEC2D_A_XI = NSPEC2D_NO_DOUBLING_XI + NSPEC2D_DOUBLING_A_XI
+ NSPEC2D_B_XI = NSPEC2D_NO_DOUBLING_XI + NSPEC2D_DOUBLING_B_XI
+ NSPEC2D_A_ETA = NSPEC2D_NO_DOUBLING_ETA + NSPEC2D_DOUBLING_A_ETA
+ NSPEC2D_B_ETA = NSPEC2D_NO_DOUBLING_ETA + NSPEC2D_DOUBLING_B_ETA
+
+ ! exact number of surface elements on the bottom and top boundaries
+ ! and theoretical number of spectral elements in radial direction
+
+ NSPEC2D_TOP = NEX_XI*NEX_ETA / NPROC
+ NSPEC2D_BOTTOM = (NEX_XI/4)*(NEX_ETA/4) / NPROC
+
+ ! NSPEC1D_RADIAL_BEDROCK = (NER_BASEMENT_SEDIM+NER_16_BASEMENT+NER_MOHO_16)/2 + NER_BOTTOM_MOHO/4
+
+ ! face with max number of elements is type B here
+ ! maximum number of surface elements on vertical boundaries of the slices
+ NSPEC2DMAX_XMIN_XMAX = NSPEC2D_B_ETA
+ NSPEC2DMAX_YMIN_YMAX = NSPEC2D_B_XI
+
+ ! theoretical number of Gauss-Lobatto points in radial direction
+ ! NPOIN1D_RADIAL_BEDROCK = NSPEC1D_RADIAL_BEDROCK*(NGLLZ-1)+1
+
+ ! 2-D addressing and buffers for summation between slices
+ ! we add one to number of points because of the flag after the last point
+ NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX*2*2 + 1
+ NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX*2*2 + 1
+
+ ! exact number of global points
+
+ ! case of the doubling regions
+
+ NGLOB_DOUBLING_1 = (NEX_PER_PROC_XI/8)*(NEX_PER_PROC_ETA/8)*(32*NGLLX**3 - 70*NGLLX**2 + 52*NGLLX - 13) &
+ - ((NEX_PER_PROC_XI/8-1)*(NEX_PER_PROC_ETA/8) + (NEX_PER_PROC_ETA/8-1)*(NEX_PER_PROC_XI/8)) * (8*NGLLX**2-11*NGLLX+4)&
+ + (NEX_PER_PROC_XI/8-1)*(NEX_PER_PROC_ETA/8-1)*(NGLLX+1)
+
+ NGLOB_DOUBLING_2 = (NEX_PER_PROC_XI/4)*(NEX_PER_PROC_ETA/4)*(32*NGLLX**3 - 70*NGLLX**2 + 52*NGLLX - 13) &
+ - ((NEX_PER_PROC_XI/4-1)*(NEX_PER_PROC_ETA/4) + (NEX_PER_PROC_ETA/4-1)*(NEX_PER_PROC_XI/4)) * (8*NGLLX**2-11*NGLLX+4)&
+ + (NEX_PER_PROC_XI/4-1)*(NEX_PER_PROC_ETA/4-1)*(NGLLX+1)
+
+ NGLOB_DOUBLING = NGLOB_DOUBLING_1 + NGLOB_DOUBLING_2
+
+ ! compute number of points in blocks with no doubling
+ ! exclude the four surfaces in contact with the doubling regions
+
+ NGLOB_NO_DOUBLING = (NEX_PER_PROC_XI/4 + 1)*(NEX_PER_PROC_ETA/4 + 1)*NER_REGULAR1 &
+ + (NEX_PER_PROC_XI/2 + 1)*(NEX_PER_PROC_ETA/2 + 1)*(NER_REGULAR2 - 1) &
+ + (NEX_PER_PROC_XI + 1)*(NEX_PER_PROC_ETA + 1)*NER_REGULAR3
+
+ NGLOB_AB = NGLOB_DOUBLING + NGLOB_NO_DOUBLING
+
+
endif ! end of section for non-regular mesh with doublings
- end subroutine compute_parameters
+end subroutine compute_parameters
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/constants.h.in 2010-03-19 05:07:35 UTC (rev 16433)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/constants.h.in 2010-03-19 20:11:02 UTC (rev 16434)
@@ -320,3 +320,9 @@
integer, parameter :: ITYPE_UNUSUAL_4 = 3
integer, parameter :: ITYPE_UNUSUAL_4p = 4
+! define number of spectral elements and points in basic symmetric mesh doubling superbrick
+ integer, parameter :: NSPEC_DOUBLING_SUPERBRICK = 32
+ integer, parameter :: NGLOB_DOUBLING_SUPERBRICK = 67
+ integer, parameter :: NSPEC_SUPERBRICK_1L = 28
+ integer, parameter :: NGLOB_SUPERBRICK_1L = 58
+ integer, parameter :: NGNOD_EIGHT_CORNERS = 8
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/create_regions_mesh.f90 2010-03-19 05:07:35 UTC (rev 16433)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/create_regions_mesh.f90 2010-03-19 20:11:02 UTC (rev 16434)
@@ -23,32 +23,7 @@
!
!=====================================================================
-! subroutine create_regions_mesh(xgrid,ygrid,zgrid,ibool,idoubling, &
-! xstore,ystore,zstore,npx,npy,iproc_xi,iproc_eta,addressing,nspec, &
-! volume_local,area_local_bottom,area_local_top, &
-! NGLOB_AB,npointot, &
-! NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER, &
-! NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
-! NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
-! HARVARD_3D_GOCAD_MODEL,NPROC_XI,NPROC_ETA,NSPEC2D_A_XI,NSPEC2D_B_XI, &
-! NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
-! myrank,LOCAL_PATH,UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK,UTM_PROJECTION_ZONE, &
-! HAUKSSON_REGIONAL_MODEL,OCEANS, &
-! VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
-! IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,MOHO_MAP_LUPEI, &
-! ANISOTROPY,SAVE_MESH_FILES,SUPPRESS_UTM_PROJECTION, &
-! ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO,NX_TOPO,NY_TOPO,USE_REGULAR_MESH)
-! subroutine create_regions_mesh(xgrid,ygrid,zgrid,ibool,idoubling, &
-! xstore,ystore,zstore,npx,npy,iproc_xi,iproc_eta,addressing,nspec, &
-! NGLOB_AB,npointot, &
-! NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER, &
-! NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
-! NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
-! HARVARD_3D_GOCAD_MODEL,NPROC_XI,NPROC_ETA,NSPEC2D_A_XI,NSPEC2D_B_XI, &
-! NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
-! myrank,LOCAL_PATH,UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
-! HAUKSSON_REGIONAL_MODEL,USE_REGULAR_MESH)
module createRegMesh
contains
@@ -63,7 +38,7 @@
NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
nsubregions,subregions,nblayers,ner_layer,NMATERIALS,material_properties, &
myrank,LOCAL_PATH,UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
- USE_REGULAR_MESH)
+ USE_REGULAR_MESH,NDOUBLINGS,ner_doublings)
! create the different regions of the mesh
@@ -74,26 +49,23 @@
! number of spectral elements in each block
integer nspec
- integer NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NER!,UTM_PROJECTION_ZONE
-! integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER
+ integer NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NER
integer NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP
integer NPROC_XI,NPROC_ETA,NSPEC2D_A_XI,NSPEC2D_B_XI
integer NSPEC2D_A_ETA,NSPEC2D_B_ETA
-! integer NX_TOPO,NY_TOPO
integer npx,npy
integer npointot
-! logical HARVARD_3D_GOCAD_MODEL,HAUKSSON_REGIONAL_MODEL
- logical USE_REGULAR_MESH!,OCEANS,IMPOSE_MINIMUM_VP_GOCAD
-! logical MOHO_MAP_LUPEI,SUPPRESS_UTM_PROJECTION
+ logical USE_REGULAR_MESH
+ integer NDOUBLINGS
+ integer, dimension(2) :: ner_doublings
+
double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
-! double precision VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM
- double precision horiz_size,vert_size!,THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR
-! double precision ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO
+ double precision horiz_size,vert_size
character(len=150) LOCAL_PATH
@@ -104,22 +76,12 @@
! real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: nodes_coords
double precision, dimension(:,:), allocatable :: nodes_coords
-! double precision xstore_local(NGLLX,NGLLY,NGLLZ)
-! double precision ystore_local(NGLLX,NGLLY,NGLLZ)
-! double precision zstore_local(NGLLX,NGLLY,NGLLZ)
-
double precision xgrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA)
double precision ygrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA)
double precision zgrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA)
integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
-! use integer array to store topography values
-! integer icornerlat,icornerlong
-! double precision lat,long,elevation
-! double precision long_corner,lat_corner,ratio_xi,ratio_eta
-! integer itopo_bathy(NX_TOPO,NY_TOPO)
-
! auxiliary variables to generate the mesh
integer ix,iy,ir,ir1,ir2,dir
integer ix1,ix2,dix,iy1,iy2,diy
@@ -140,17 +102,6 @@
integer nblayers
integer ner_layer(nblayers)
-! Gauss-Lobatto-Legendre points and weights of integration
-! double precision, dimension(:), allocatable :: xigll,yigll,zigll,wxgll,wygll,wzgll
-
-! 3D shape functions and their derivatives
-! double precision, dimension(:,:,:,:), allocatable :: shape3D
-! double precision, dimension(:,:,:,:,:), allocatable :: dershape3D
-
-! 2D shape functions and their derivatives
-! double precision, dimension(:,:,:), allocatable :: shape2D_x,shape2D_y,shape2D_bottom,shape2D_top
-! double precision, dimension(:,:,:,:), allocatable :: dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top
-
! topology of the elements
integer iaddx(NGNOD)
integer iaddy(NGNOD)
@@ -160,39 +111,12 @@
double precision yelm(NGNOD)
double precision zelm(NGNOD)
-
-! ! for model density
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhostore,kappastore,mustore
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: c11store,c12store,c13store,c14store,c15store,c16store,&
-! c22store,c23store,c24store,c25store,c26store,c33store,c34store,c35store,c36store,c44store,c45store,c46store,&
-! c55store,c56store,c66store
-
-! ! the jacobian
-! real(kind=CUSTOM_REAL) jacobianl
-
! boundary locator
logical, dimension(:,:), allocatable :: iboun
-! arrays with mesh parameters
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: xixstore,xiystore,xizstore, &
-! etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,jacobianstore
-
-! mass matrix and bathymetry for ocean load
-! integer ix_oceans,iy_oceans,iz_oceans,ispec_oceans
-! integer ispec2D_ocean_bottom
-! integer nglob_oceans
-! double precision xval,yval
-! double precision height_oceans
-! real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_ocean_load
-
! proc numbers for MPI
integer myrank
-! check area and volume of the final mesh
-! double precision weight
-! double precision area_local_bottom,area_local_top
-! double precision volume_local
-
! variables for creating array ibool (some arrays also used for AVS or DX files)
integer, dimension(:), allocatable :: iglob,locval
logical, dimension(:), allocatable :: ifseg
@@ -201,33 +125,9 @@
integer nglob,NGLOB_AB
integer ieoff,ilocnum
-! mass matrix
-! real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass
-
! boundary parameters locator
integer, dimension(:), allocatable :: ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top
-! ---- Moho Vars here ------
-! ! Moho boundary locator
-! integer, dimension(:), allocatable :: ibelm_moho_top, ibelm_moho_bot
-! logical, dimension(:), allocatable :: is_moho_top, is_moho_bot
-
-! ! 2-D jacobian and normals
-! real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: jacobian2D_moho
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_moho
-
-! ! number of elements on the boundaries
-! integer nspec_moho_top, nspec_moho_bottom
-! ---------------------------
-
-! 2-D jacobians and normals
-! real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
-! jacobian2D_xmin,jacobian2D_xmax, &
-! jacobian2D_ymin,jacobian2D_ymax,jacobian2D_bottom,jacobian2D_top
-
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
-! normal_xmin,normal_xmax,normal_ymin,normal_ymax,normal_bottom,normal_top
-
! MPI cut-planes parameters along xi and along eta
logical, dimension(:,:), allocatable :: iMPIcut_xi,iMPIcut_eta
@@ -237,32 +137,9 @@
! number of elements on the boundaries
integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax
- integer i,j,k,ia,ispec,itype_element,ipoin
+ integer i,j,k,ia,ispec,ispec_superbrick,itype_element,ipoin
integer iproc_xi,iproc_eta
-! double precision rho,vp,vs
-! double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
-
-! for the Harvard 3-D basin model
-! double precision vp_block_gocad_MR(0:NX_GOCAD_MR-1,0:NY_GOCAD_MR-1,0:NZ_GOCAD_MR-1)
-! double precision vp_block_gocad_HR(0:NX_GOCAD_HR-1,0:NY_GOCAD_HR-1,0:NZ_GOCAD_HR-1)
-! integer irecord,nrecord,i_vp
-! character(len=150) BASIN_MODEL_3D_MEDIUM_RES_FILE,BASIN_MODEL_3D_HIGH_RES_FILE
-
-! for the harvard 3D salton sea model
-! real :: vp_st_gocad(GOCAD_ST_NU,GOCAD_ST_NV,GOCAD_ST_NW)
-! double precision :: umesh, vmesh, wmesh, vp_st, vs_st, rho_st
-
-! for Hauksson's model
-! double precision, dimension(NLAYERS_HAUKSSON,NGRID_NEW_HAUKSSON,NGRID_NEW_HAUKSSON) :: vp_hauksson,vs_hauksson
-! integer ilayer
-! character(len=150 ) HAUKSSON_REGIONAL_MODEL_FILE
-
-! Stacey put back
-! indices for Clayton-Engquist absorbing conditions
-! integer, dimension(:,:), allocatable :: nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_vp,rho_vs
-
! flag indicating whether point is in the sediments
! logical point_is_in_sediments
logical, dimension(:,:,:,:), allocatable :: flag_sediments
@@ -274,68 +151,19 @@
! second dimension : #rho #vp #vs #Q_flag #anisotropy_flag #domain_id
double precision , dimension(NMATERIALS,6) :: material_properties
+! doublings zone
+ integer :: nspec_sb
+ integer :: offset_x,offset_y,offset_z
+ logical, dimension(NSPEC_DOUBLING_SUPERBRICK,6) :: iboun_sb
+ integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC_DOUBLING_SUPERBRICK) :: ibool_superbrick
+ double precision, dimension(NGLOB_DOUBLING_SUPERBRICK) :: x_superbrick,y_superbrick,z_superbrick
+
! **************
! create the name for the database of the current slide and region
call create_name_database(prname,myrank,LOCAL_PATH)
-! Gauss-Lobatto-Legendre points of integration
-! allocate(xigll(NGLLX))
-! allocate(yigll(NGLLY))
-! allocate(zigll(NGLLZ))
-
-! Gauss-Lobatto-Legendre weights of integration
-! allocate(wxgll(NGLLX))
-! allocate(wygll(NGLLY))
-! allocate(wzgll(NGLLZ))
-
-! 3D shape functions and their derivatives
- ! allocate(shape3D(NGNOD,NGLLX,NGLLY,NGLLZ))
-! allocate(dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ))
-
-! 2D shape functions and their derivatives
-! allocate(shape2D_x(NGNOD2D,NGLLY,NGLLZ))
-! allocate(shape2D_y(NGNOD2D,NGLLX,NGLLZ))
-! allocate(shape2D_bottom(NGNOD2D,NGLLX,NGLLY))
-! allocate(shape2D_top(NGNOD2D,NGLLX,NGLLY))
-! allocate(dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ))
-! allocate(dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ))
-! allocate(dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY))
-! allocate(dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY))
-
-! array with model density
-! allocate(rhostore(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(kappastore(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(mustore(NGLLX,NGLLY,NGLLZ,nspec))
-
-! array with anisotropy
-! allocate(c11store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c12store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c13store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c14store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c15store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c16store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c22store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c23store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c24store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c25store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c26store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c33store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c34store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c35store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c36store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c44store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c45store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c46store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c55store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c56store(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(c66store(NGLLX,NGLLY,NGLLZ,nspec))
-
-! Stacey
-! allocate(rho_vp(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(rho_vs(NGLLX,NGLLY,NGLLZ,nspec))
-
! flag indicating whether point is in the sediments
allocate(flag_sediments(NGLLX,NGLLY,NGLLZ,nspec))
allocate(not_fully_in_bedrock(nspec))
@@ -343,18 +171,6 @@
! boundary locator
allocate(iboun(6,nspec))
-! arrays with mesh parameters
-! allocate(xixstore(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(xiystore(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(xizstore(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(etaxstore(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(etaystore(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(etazstore(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(gammaxstore(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(gammaystore(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(gammazstore(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(jacobianstore(NGLLX,NGLLY,NGLLZ,nspec))
-
! boundary parameters locator
allocate(ibelm_xmin(NSPEC2DMAX_XMIN_XMAX))
allocate(ibelm_xmax(NSPEC2DMAX_XMIN_XMAX))
@@ -363,67 +179,10 @@
allocate(ibelm_bottom(NSPEC2D_BOTTOM))
allocate(ibelm_top(NSPEC2D_TOP))
-! 2-D jacobians and normals
-! allocate(jacobian2D_xmin(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
-! allocate(jacobian2D_xmax(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
-! allocate(jacobian2D_ymin(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
-! allocate(jacobian2D_ymax(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
-! allocate(jacobian2D_bottom(NGLLX,NGLLY,NSPEC2D_BOTTOM))
-! allocate(jacobian2D_top(NGLLX,NGLLY,NSPEC2D_TOP))
-
-! allocate(normal_xmin(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
-! allocate(normal_xmax(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
-! allocate(normal_ymin(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
-! allocate(normal_ymax(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
-! allocate(normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM))
-! allocate(normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP))
-
-! Moho boundary parameters, 2-D jacobians and normals
-! if (SAVE_MOHO_MESH) then
-! allocate(ibelm_moho_top(NSPEC2D_BOTTOM))
-! allocate(ibelm_moho_bot(NSPEC2D_BOTTOM))
-! allocate(is_moho_top(nspec))
-! allocate(is_moho_bot(nspec))
-! is_moho_top = .false.
-! is_moho_bot = .false.
-! nspec_moho_top = 0
-! nspec_moho_bottom = 0
-! allocate(jacobian2D_moho(NGLLX,NGLLY,NSPEC2D_BOTTOM))
-! allocate(normal_moho(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM))
-! endif
-
-
-! Stacey put back
-! allocate(nimin(2,NSPEC2DMAX_YMIN_YMAX))
-! allocate(nimax(2,NSPEC2DMAX_YMIN_YMAX))
-! allocate(njmin(2,NSPEC2DMAX_XMIN_XMAX))
-! allocate(njmax(2,NSPEC2DMAX_XMIN_XMAX))
-! allocate(nkmin_xi(2,NSPEC2DMAX_XMIN_XMAX))
-! allocate(nkmin_eta(2,NSPEC2DMAX_YMIN_YMAX))
-
! MPI cut-planes parameters along xi and along eta
allocate(iMPIcut_xi(2,nspec))
allocate(iMPIcut_eta(2,nspec))
-! set up coordinates of the Gauss-Lobatto-Legendre points
-! call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
-! call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
-! call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
-
-! if number of points is odd, the middle abscissa is exactly zero
-! if(mod(NGLLX,2) /= 0) xigll((NGLLX-1)/2+1) = ZERO
-! if(mod(NGLLY,2) /= 0) yigll((NGLLY-1)/2+1) = ZERO
-! if(mod(NGLLZ,2) /= 0) zigll((NGLLZ-1)/2+1) = ZERO
-
-! get the 3-D shape functions
-! call get_shape3D(myrank,shape3D,dershape3D,xigll,yigll,zigll)
-
-! get the 2-D shape functions
-! call get_shape2D(myrank,shape2D_x,dershape2D_x,yigll,zigll,NGLLY,NGLLZ)
-! call get_shape2D(myrank,shape2D_y,dershape2D_y,xigll,zigll,NGLLX,NGLLZ)
-! call get_shape2D(myrank,shape2D_bottom,dershape2D_bottom,xigll,yigll,NGLLX,NGLLY)
-! call get_shape2D(myrank,shape2D_top,dershape2D_top,xigll,yigll,NGLLX,NGLLY)
-
! allocate memory for arrays
allocate(iglob(npointot))
allocate(locval(npointot))
@@ -432,208 +191,144 @@
allocate(yp(npointot))
allocate(zp(npointot))
-! !--- read Hauksson's model
-! if(HAUKSSON_REGIONAL_MODEL) then
-! call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
-! 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
-! 'DATA/hauksson_model/hauksson_final_grid_smooth.dat')
-! ! call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
-! ! 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
-! ! 'DATA/lin_model/lin_final_grid_smooth.dat')
-! open(unit=14,file=HAUKSSON_REGIONAL_MODEL_FILE,status='old',action='read')
-! do iy = 1,NGRID_NEW_HAUKSSON
-! do ix = 1,NGRID_NEW_HAUKSSON
-! read(14,*) (vp_hauksson(ilayer,ix,iy),ilayer=1,NLAYERS_HAUKSSON), &
-! (vs_hauksson(ilayer,ix,iy),ilayer=1,NLAYERS_HAUKSSON)
-! enddo
-! enddo
-! close(14)
-! vp_hauksson(:,:,:) = vp_hauksson(:,:,:) * 1000.d0
-! vs_hauksson(:,:,:) = vs_hauksson(:,:,:) * 1000.d0
-! endif
-! !--- read the Harvard 3-D basin model
-! if(HARVARD_3D_GOCAD_MODEL) then
-! ! read medium-resolution model
-
-! ! initialize array to undefined values everywhere
-! vp_block_gocad_MR(:,:,:) = 20000.
-
-! ! read Vp from extracted text file
-! call get_value_string(BASIN_MODEL_3D_MEDIUM_RES_FILE, &
-! 'model.BASIN_MODEL_3D_MEDIUM_RES_FILE', &
-! 'DATA/la_3D_block_harvard/la_3D_medium_res/LA_MR_voxet_extracted.txt')
-! open(unit=27,file=BASIN_MODEL_3D_MEDIUM_RES_FILE,status='old',action='read')
-! read(27,*) nrecord
-! do irecord = 1,nrecord
-! read(27,*) ix,iy,iz,i_vp
-! if(ix<0 .or. ix>NX_GOCAD_MR-1 .or. iy<0 .or. iy>NY_GOCAD_MR-1 .or. iz<0 .or. iz>NZ_GOCAD_MR-1) &
-! stop 'wrong array index read in Gocad medium-resolution file'
-! vp_block_gocad_MR(ix,iy,iz) = dble(i_vp)
-! enddo
-! close(27)
-
-! ! read high-resolution model
-
-! ! initialize array to undefined values everywhere
-! vp_block_gocad_HR(:,:,:) = 20000.
-
-! ! read Vp from extracted text file
-! call get_value_string(BASIN_MODEL_3D_HIGH_RES_FILE, &
-! 'model.BASIN_MODEL_3D_HIGH_RES_FILE', &
-! 'DATA/la_3D_block_harvard/la_3D_high_res/LA_HR_voxet_extracted.txt')
-! open(unit=27,file=BASIN_MODEL_3D_HIGH_RES_FILE,status='old',action='read')
-! read(27,*) nrecord
-! do irecord = 1,nrecord
-! read(27,*) ix,iy,iz,i_vp
-! if(ix<0 .or. ix>NX_GOCAD_HR-1 .or. iy<0 .or. iy>NY_GOCAD_HR-1 .or. iz<0 .or. iz>NZ_GOCAD_HR-1) &
-! stop 'wrong array index read in Gocad high-resolution file'
-! vp_block_gocad_HR(ix,iy,iz) = dble(i_vp)
-! enddo
-! close(27)
-
-! ! read Salton Trough model
-! call read_salton_sea_model(vp_st_gocad)
-
-! endif
-
!--- apply heuristic rule to modify doubling regions to balance angles
- if(APPLY_HEURISTIC_RULE .and. .not. USE_REGULAR_MESH) then
+! if(APPLY_HEURISTIC_RULE .and. .not. USE_REGULAR_MESH) then
- stop 'pas encore implemente'
+! stop 'pas encore implemente'
-! define number of subregions affected by heuristic rule in doubling regions
- nsubregions = 8
+! ! define number of subregions affected by heuristic rule in doubling regions
+! nsubregions = 8
- do isubregion = 1,nsubregions
+! do isubregion = 1,nsubregions
-! define shape of elements for heuristic
- call define_subregions_heuristic(myrank,isubregion,iaddx,iaddy,iaddz, &
- nblayers,ner_layer, &
- ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar, &
- itype_element,npx,npy)
- ! NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM)
+! ! define shape of elements for heuristic
+! call define_subregions_heuristic(myrank,isubregion,iaddx,iaddy,iaddz, &
+! nblayers,ner_layer, &
+! ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar, &
+! itype_element,npx,npy)
+! ! NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM)
-! loop on all the mesh points in current subregion
- do ir = ir1,ir2,dir
- do iy = iy1,iy2,diy
- do ix = ix1,ix2,dix
+! ! loop on all the mesh points in current subregion
+! do ir = ir1,ir2,dir
+! do iy = iy1,iy2,diy
+! do ix = ix1,ix2,dix
-! this heuristic rule is only valid for 8-node elements
-! it would not work in the case of 27 nodes
+! ! this heuristic rule is only valid for 8-node elements
+! ! it would not work in the case of 27 nodes
-!----
- if(itype_element == ITYPE_UNUSUAL_1) then
+! !----
+! if(itype_element == ITYPE_UNUSUAL_1) then
-! side 1
- horiz_size = xgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) &
- - xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- xgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
- xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * MAGIC_RATIO
+! ! side 1
+! horiz_size = xgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) &
+! - xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
+! xgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
+! xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * MAGIC_RATIO
- vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
- - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
- zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
+! vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
+! - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
+! zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
+! zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
-! side 2
- horiz_size = xgrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
- - xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
- xgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
- xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + horiz_size * MAGIC_RATIO
+! ! side 2
+! horiz_size = xgrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
+! - xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
+! xgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
+! xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + horiz_size * MAGIC_RATIO
- vert_size = zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) &
- - zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
- zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
- zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + vert_size * MAGIC_RATIO / 0.50
+! vert_size = zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) &
+! - zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
+! zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
+! zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + vert_size * MAGIC_RATIO / 0.50
-!----
- else if(itype_element == ITYPE_UNUSUAL_1p) then
+! !----
+! else if(itype_element == ITYPE_UNUSUAL_1p) then
-! side 1
- horiz_size = xgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) &
- - xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- xgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
- xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * (1. - MAGIC_RATIO)
+! ! side 1
+! horiz_size = xgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) &
+! - xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
+! xgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
+! xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * (1. - MAGIC_RATIO)
- vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
- - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
- zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
+! vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
+! - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
+! zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
+! zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
-! side 2
- horiz_size = xgrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
- - xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
- xgrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
- xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + horiz_size * (1. - MAGIC_RATIO)
+! ! side 2
+! horiz_size = xgrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
+! - xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
+! xgrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
+! xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + horiz_size * (1. - MAGIC_RATIO)
- vert_size = zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) &
- - zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
- zgrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
- zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + vert_size * MAGIC_RATIO / 0.50
+! vert_size = zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) &
+! - zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
+! zgrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
+! zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + vert_size * MAGIC_RATIO / 0.50
-!----
- else if(itype_element == ITYPE_UNUSUAL_4) then
+! !----
+! else if(itype_element == ITYPE_UNUSUAL_4) then
-! side 1
- horiz_size = ygrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
- - ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
- ygrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
- ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + horiz_size * (1. - MAGIC_RATIO)
+! ! side 1
+! horiz_size = ygrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
+! - ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
+! ygrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
+! ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + horiz_size * (1. - MAGIC_RATIO)
- vert_size = zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) &
- - zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
- zgrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
- zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + vert_size * MAGIC_RATIO / 0.50
+! vert_size = zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) &
+! - zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
+! zgrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
+! zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + vert_size * MAGIC_RATIO / 0.50
-! side 2
- horiz_size = ygrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) &
- - ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- ygrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
- ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * (1. - MAGIC_RATIO)
+! ! side 2
+! horiz_size = ygrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) &
+! - ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
+! ygrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
+! ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * (1. - MAGIC_RATIO)
- vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
- - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
- zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
+! vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
+! - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
+! zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
+! zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
-!----
- else if(itype_element == ITYPE_UNUSUAL_4p) then
+! !----
+! else if(itype_element == ITYPE_UNUSUAL_4p) then
-! side 1
- horiz_size = ygrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
- - ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
- ygrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
- ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + horiz_size * MAGIC_RATIO
+! ! side 1
+! horiz_size = ygrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
+! - ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
+! ygrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
+! ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + horiz_size * MAGIC_RATIO
- vert_size = zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) &
- - zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
- zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
- zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + vert_size * MAGIC_RATIO / 0.50
+! vert_size = zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) &
+! - zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
+! zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
+! zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + vert_size * MAGIC_RATIO / 0.50
-! side 2
- horiz_size = ygrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) &
- - ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- ygrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
- ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * MAGIC_RATIO
+! ! side 2
+! horiz_size = ygrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) &
+! - ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
+! ygrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
+! ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * MAGIC_RATIO
- vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
- - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
- zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
+! vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
+! - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
+! zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
+! zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
- endif
+! endif
-enddo
-enddo
-enddo
+! enddo
+! enddo
+! enddo
-enddo
+! enddo
-endif
+! endif
!---
@@ -651,7 +346,7 @@
do iy = iy1,iy2,diy
do ix = ix1,ix2,dix
- material_num(ir,ix,ir) = imaterial_number
+ material_num(ir,ix,iy) = imaterial_number
! end of loop on all the mesh points in current subregion
enddo
@@ -665,13 +360,23 @@
if(USE_REGULAR_MESH) then
nmeshregions = 1
else
- ! TO DO
+ call define_superbrick(x_superbrick,y_superbrick,z_superbrick,ibool_superbrick,iboun_sb)
+ !call define_superbrick_one_layer(x_superbrick,y_superbrick,z_superbrick,ibool_superbrick,iboun_sb)
+ nspec_sb = NSPEC_DOUBLING_SUPERBRICK
+ !nspec_sb = NSPEC_SUPERBRICK_1L
+ if(NDOUBLINGS == 1) then
+ nmeshregions = 3
+ else if(NDOUBLINGS == 2) then
+ nmeshregions = 5
+ else
+ stop 'Wrong number of mesh regions'
+ end if
endif
do isubregion = 1,nmeshregions
! define shape of elements
call define_mesh_regions(USE_REGULAR_MESH,isubregion,NER,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,iproc_xi,iproc_eta,&
- nblayers,ner_layer, &
+ nblayers,ner_layer,NDOUBLINGS,ner_doublings,&
iaddx,iaddy,iaddz,ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar)
! loop on all the mesh points in current subregion
@@ -679,18 +384,20 @@
do iy = iy1,iy2,diy
do ix = ix1,ix2,dix
-! loop over the NGNOD nodes
- do ia=1,NGNOD
- xelm(ia) = xgrid(ir+iar*iaddz(ia),ix+iax*iaddx(ia),iy+iay*iaddy(ia))
- yelm(ia) = ygrid(ir+iar*iaddz(ia),ix+iax*iaddx(ia),iy+iay*iaddy(ia))
- zelm(ia) = zgrid(ir+iar*iaddz(ia),ix+iax*iaddx(ia),iy+iay*iaddy(ia))
- enddo
+ if(modulo(isubregion,2) == 1) then ! Regular subregion case
+
+ ! loop over the NGNOD nodes
+ do ia=1,NGNOD
+ xelm(ia) = xgrid(ir+iar*iaddz(ia),ix+iax*iaddx(ia),iy+iay*iaddy(ia))
+ yelm(ia) = ygrid(ir+iar*iaddz(ia),ix+iax*iaddx(ia),iy+iay*iaddy(ia))
+ zelm(ia) = zgrid(ir+iar*iaddz(ia),ix+iax*iaddx(ia),iy+iay*iaddy(ia))
+ enddo
! add one spectral element to the list and store its material number
ispec = ispec + 1
- if(ispec > nspec) then
- call exit_MPI(myrank,'ispec greater than nspec in mesh creation')
- end if
+ if(ispec > nspec) then
+ call exit_MPI(myrank,'ispec greater than nspec in mesh creation')
+ end if
! A revoir
if((ir == ir2) .and. (isubregion == nmeshregions)) then
@@ -699,214 +406,8 @@
doubling_index = IFLAG_BASEMENT_TOPO
endif
- true_material_num(ispec) = material_num(ir,ix,ir)
+ true_material_num(ispec) = material_num(ir,ix,iy)
-! ! assign Moho surface element
-! if (SAVE_MOHO_MESH) then
-! if (isubregion == 15 .and. ir == ir1) then
-! nspec_moho_top = nspec_moho_top + 1
-! if (nspec_moho_top > NSPEC2D_BOTTOM) call exit_mpi(myrank,"Error counting moho top elements")
-! ibelm_moho_top(nspec_moho_top) = ispec
-! call compute_jacobian_2D(myrank,nspec_moho_top,xelm(1:NGNOD2D),yelm(1:NGNOD2D),zelm(1:NGNOD2D), &
-! dershape2D_bottom,jacobian2D_moho,normal_moho,NGLLX,NGLLY,NSPEC2D_BOTTOM)
-! is_moho_top(ispec) = .true.
-! else if (isubregion == 28 .and. ir+dir > ir2) then
-! nspec_moho_bottom = nspec_moho_bottom + 1
-! if (nspec_moho_bottom > NSPEC2D_BOTTOM) call exit_mpi(myrank,"Error counting moho bottom elements")
-! ibelm_moho_bot(nspec_moho_bottom) = ispec
-! is_moho_bot(ispec) = .true.
-! endif
-! endif
-
-! ! initialize flag indicating whether element is in sediments
-! not_fully_in_bedrock(ispec) = .false.
-
-! ! create mesh element
-! do k=1,NGLLZ
-! do j=1,NGLLY
-! do i=1,NGLLX
-
-! ! compute mesh coordinates
-! xmesh = ZERO
-! ymesh = ZERO
-! zmesh = ZERO
-! do ia=1,NGNOD
-! xmesh = xmesh + shape3D(ia,i,j,k)*xelm(ia)
-! ymesh = ymesh + shape3D(ia,i,j,k)*yelm(ia)
-! zmesh = zmesh + shape3D(ia,i,j,k)*zelm(ia)
-! enddo
-
-! xstore_local(i,j,k) = xmesh
-! ystore_local(i,j,k) = ymesh
-! zstore_local(i,j,k) = zmesh
-
-! ! initialize flag indicating whether point is in the sediments
-! point_is_in_sediments = .false.
-
-! if(ANISOTROPY) then
-! call aniso_model(doubling_index,zmesh,rho,vp,vs,c11,c12,c13,c14,c15,c16,&
-! c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
-! else
-! ! get the regional model parameters
-! if(HAUKSSON_REGIONAL_MODEL) then
-! ! get density from socal model
-! call socal_model(doubling_index,rho,vp,vs)
-! ! get vp and vs from Hauksson
-! call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp,vs,MOHO_MAP_LUPEI)
-! ! if Moho map is used, then assume homogeneous medium below the Moho
-! ! and use bottom layer of Hauksson's model in the halfspace
-! if(MOHO_MAP_LUPEI .and. doubling_index == IFLAG_HALFSPACE_MOHO) &
-! call socal_model(IFLAG_HALFSPACE_MOHO,rho,vp,vs)
-! else
-! call socal_model(doubling_index,rho,vp,vs)
-! ! include attenuation in first SoCal layer if needed
-! ! uncomment line below to include attenuation in the 1D case
-! ! if(zmesh >= DEPTH_5p5km_SOCAL) point_is_in_sediments = .true.
-! endif
-
-! ! get the Harvard 3-D basin model
-! if(HARVARD_3D_GOCAD_MODEL .and. &
-! (doubling_index == IFLAG_ONE_LAYER_TOPOGRAPHY &
-! .or. doubling_index == IFLAG_BASEMENT_TOPO) &
-! .and. xmesh >= ORIG_X_GOCAD_MR &
-! .and. xmesh <= END_X_GOCAD_MR &
-! .and. ymesh >= ORIG_Y_GOCAD_MR &
-! .and. ymesh <= END_Y_GOCAD_MR) then
-
-! ! use medium-resolution model first
-! call interpolate_gocad_block_MR(vp_block_gocad_MR, &
-! xmesh,ymesh,zmesh,rho,vp,vs,point_is_in_sediments, &
-! VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
-! IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_MR, &
-! vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL,&
-! MOHO_MAP_LUPEI)
-
-! ! then superimpose high-resolution model
-! if(xmesh >= ORIG_X_GOCAD_HR &
-! .and. xmesh <= END_X_GOCAD_HR &
-! .and. ymesh >= ORIG_Y_GOCAD_HR &
-! .and. ymesh <= END_Y_GOCAD_HR) &
-! call interpolate_gocad_block_HR(vp_block_gocad_HR,vp_block_gocad_MR,&
-! xmesh,ymesh,zmesh,rho,vp,vs,point_is_in_sediments, &
-! VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
-! IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR, &
-! vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL, &
-! MOHO_MAP_LUPEI)
-
-! endif
-! ! get the Harvard Salton Trough model
-! if (HARVARD_3D_GOCAD_MODEL) then
-! call vx_xyz2uvw(xmesh, ymesh, zmesh, umesh, vmesh, wmesh)
-! if (umesh >= 0 .and. umesh <= GOCAD_ST_NU-1 .and. &
-! vmesh >= 0 .and. vmesh <= GOCAD_ST_NV-1 .and. &
-! wmesh >= 0 .and. wmesh <= GOCAD_ST_NW-1) then
-! call vx_xyz_interp(umesh,vmesh,wmesh, vp_st, vs_st, rho_st, vp_st_gocad)
-! if (abs(vp_st - GOCAD_ST_NO_DATA_VALUE) > 1.0d-3) then
-! vp = vp_st
-! vs = vs_st
-! rho = rho_st
-! endif
-! endif
-! endif
-! endif
-! ! store flag indicating whether point is in the sediments
-! flag_sediments(i,j,k,ispec) = point_is_in_sediments
-! if(point_is_in_sediments) not_fully_in_bedrock(ispec) = .true.
-
-! ! define elastic parameters in the model
-! ! distinguish between single and double precision for reals
-! if(ANISOTROPY) then
-
-! if(CUSTOM_REAL == SIZE_REAL) then
-! rhostore(i,j,k,ispec) = sngl(rho)
-! kappastore(i,j,k,ispec) = sngl(rho*(vp*vp - 4.d0*vs*vs/3.d0))
-! mustore(i,j,k,ispec) = sngl(rho*vs*vs)
-! c11store(i,j,k,ispec) = sngl(c11)
-! c12store(i,j,k,ispec) = sngl(c12)
-! c13store(i,j,k,ispec) = sngl(c13)
-! c14store(i,j,k,ispec) = sngl(c14)
-! c15store(i,j,k,ispec) = sngl(c15)
-! c16store(i,j,k,ispec) = sngl(c16)
-! c22store(i,j,k,ispec) = sngl(c22)
-! c23store(i,j,k,ispec) = sngl(c23)
-! c24store(i,j,k,ispec) = sngl(c24)
-! c25store(i,j,k,ispec) = sngl(c25)
-! c26store(i,j,k,ispec) = sngl(c26)
-! c33store(i,j,k,ispec) = sngl(c33)
-! c34store(i,j,k,ispec) = sngl(c34)
-! c35store(i,j,k,ispec) = sngl(c35)
-! c36store(i,j,k,ispec) = sngl(c36)
-! c44store(i,j,k,ispec) = sngl(c44)
-! c45store(i,j,k,ispec) = sngl(c45)
-! c46store(i,j,k,ispec) = sngl(c46)
-! c55store(i,j,k,ispec) = sngl(c55)
-! c56store(i,j,k,ispec) = sngl(c56)
-! c66store(i,j,k,ispec) = sngl(c66)
-! ! Stacey
-! rho_vp(i,j,k,ispec) = sngl(rho*vp)
-! rho_vs(i,j,k,ispec) = sngl(rho*vs)
-! else
-! rhostore(i,j,k,ispec) = rho
-! kappastore(i,j,k,ispec) = rho*(vp*vp - 4.d0*vs*vs/3.d0)
-! mustore(i,j,k,ispec) = rho*vs*vs
-! 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
-! ! Stacey
-! rho_vp(i,j,k,ispec) = rho*vp
-! rho_vs(i,j,k,ispec) = rho*vs
-! endif
-
-
-! else
-! if(CUSTOM_REAL == SIZE_REAL) then
-! rhostore(i,j,k,ispec) = sngl(rho)
-! kappastore(i,j,k,ispec) = sngl(rho*(vp*vp - 4.d0*vs*vs/3.d0))
-! mustore(i,j,k,ispec) = sngl(rho*vs*vs)
-
-! ! Stacey
-! rho_vp(i,j,k,ispec) = sngl(rho*vp)
-! rho_vs(i,j,k,ispec) = sngl(rho*vs)
-! else
-! rhostore(i,j,k,ispec) = rho
-! kappastore(i,j,k,ispec) = rho*(vp*vp - 4.d0*vs*vs/3.d0)
-! mustore(i,j,k,ispec) = rho*vs*vs
-
-! ! Stacey
-! rho_vp(i,j,k,ispec) = rho*vp
-! rho_vs(i,j,k,ispec) = rho*vs
-! endif
-! endif
-
-! enddo
-! enddo
-! enddo
-
-! ! compute coordinates and jacobian
-! call calc_jacobian(myrank,xixstore,xiystore,xizstore, &
-! etaxstore,etaystore,etazstore, &
-! gammaxstore,gammaystore,gammazstore,jacobianstore, &
-! xstore,ystore,zstore, &
-! xelm,yelm,zelm,shape3D,dershape3D,ispec,nspec)
-
! store coordinates
call store_coords(xstore,ystore,zstore,xelm,yelm,zelm,ispec,nspec)
@@ -916,20 +417,59 @@
iboun,iMPIcut_xi,iMPIcut_eta,NPROC_XI,NPROC_ETA, &
UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK)
-! end of loop on all the mesh points in current subregion
- enddo
- enddo
- enddo
+ else ! Irregular subregion case
+
+ ! loop on all the elements in the mesh doubling superbrick
+ do ispec_superbrick = 1,nspec_sb
+ ! loop on all the corner nodes of this element
+ do ia = 1,NGNOD_EIGHT_CORNERS
+ ! define topological coordinates of this mesh point
+ offset_x = ix + iax*x_superbrick(ibool_superbrick(ia,ispec_superbrick))
+ offset_y = iy + iay*y_superbrick(ibool_superbrick(ia,ispec_superbrick))
+ offset_z = ir + iar*z_superbrick(ibool_superbrick(ia,ispec_superbrick))
-! end of loop on all the subregions of the current region the mesh
- enddo
+ xelm(ia) = xgrid(offset_z,offset_x,offset_y)
+ yelm(ia) = ygrid(offset_z,offset_x,offset_y)
+ zelm(ia) = zgrid(offset_z,offset_x,offset_y)
+ enddo
+
+ ! add one spectral element to the list and store its material number
+ ispec = ispec + 1
+ if(ispec > nspec) then
+ call exit_MPI(myrank,'ispec greater than nspec in mesh creation')
+ end if
+ ! A revoir
+ if((ir == ir2) .and. (isubregion == nmeshregions)) then
+ doubling_index = IFLAG_ONE_LAYER_TOPOGRAPHY
+ else
+ doubling_index = IFLAG_BASEMENT_TOPO
+ endif
+
+ true_material_num(ispec) = material_num(ir,ix,iy)
+
+ ! store coordinates
+ call store_coords(xstore,ystore,zstore,xelm,yelm,zelm,ispec,nspec)
+
+ ! detect mesh boundaries
+ call get_flags_boundaries(nspec,iproc_xi,iproc_eta,ispec,doubling_index, &
+ xstore(:,:,:,ispec),ystore(:,:,:,ispec),zstore(:,:,:,ispec), &
+ iboun,iMPIcut_xi,iMPIcut_eta,NPROC_XI,NPROC_ETA, &
+ UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK)
+ end do
+ end if
+ ! end of loop on all the mesh points in current subregion
+ enddo
+ enddo
+ enddo
+
+ ! end of loop on all the subregions of the current region the mesh
+enddo
-
-
! check total number of spectral elements created
+
if(ispec /= nspec) call exit_MPI(myrank,'ispec should equal nspec')
do ispec=1,nspec
@@ -951,6 +491,7 @@
! put in classical format
allocate(nodes_coords(nglob,3))
+
do ispec=1,nspec
ieoff = NGLLCUBE*(ispec-1)
ilocnum = 0
@@ -968,6 +509,7 @@
enddo
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1071,92 +613,18 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(minval(ibool(:,:,:,:)) /= 1 .or. maxval(ibool(:,:,:,:)) /= NGLOB_AB) then
+ print*,maxval(ibool(:,:,:,:)) ,NGLOB_AB
call exit_MPI(myrank,'incorrect global numbering')
end if
-! ! create a new indirect addressing array instead, to reduce cache misses
-! ! in memory access in the solver
-! allocate(copy_ibool_ori(NGLLX,NGLLY,NGLLZ,nspec))
-! allocate(mask_ibool(nglob))
-! mask_ibool(:) = -1
-! copy_ibool_ori(:,:,:,:) = ibool(:,:,:,:)
-! inumber = 0
-! do ispec=1,nspec
-! do k=1,NGLLZ
-! do j=1,NGLLY
-! do i=1,NGLLX
-! if(mask_ibool(copy_ibool_ori(i,j,k,ispec)) == -1) then
-! ! create a new point
-! inumber = inumber + 1
-! ibool(i,j,k,ispec) = inumber
-! mask_ibool(copy_ibool_ori(i,j,k,ispec)) = inumber
-! else
-! ! use an existing point created previously
-! ibool(i,j,k,ispec) = mask_ibool(copy_ibool_ori(i,j,k,ispec))
-! endif
-! enddo
-! enddo
-! enddo
-! enddo
-! deallocate(copy_ibool_ori)
-! deallocate(mask_ibool)
-! ! creating mass matrix (will be fully assembled with MPI in the solver)
-! allocate(rmass(nglob))
-! rmass(:) = 0._CUSTOM_REAL
-
-! do ispec=1,nspec
-! do k=1,NGLLZ
-! do j=1,NGLLY
-! do i=1,NGLLX
-! weight=wxgll(i)*wygll(j)*wzgll(k)
-! iglobnum=ibool(i,j,k,ispec)
-
-! jacobianl=jacobianstore(i,j,k,ispec)
-
-! ! distinguish between single and double precision for reals
-! if(CUSTOM_REAL == SIZE_REAL) then
-! rmass(iglobnum) = rmass(iglobnum) + &
-! sngl(dble(rhostore(i,j,k,ispec)) * dble(jacobianl) * weight)
-! else
-! rmass(iglobnum) = rmass(iglobnum) + rhostore(i,j,k,ispec) * jacobianl * weight
-! endif
-
-! enddo
-! enddo
-! enddo
-! enddo
-
-! call get_jacobian_boundaries(myrank,iboun,nspec,xstore,ystore,zstore, &
-! dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
-! ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
-! nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
-! jacobian2D_xmin,jacobian2D_xmax, &
-! jacobian2D_ymin,jacobian2D_ymax, &
-! jacobian2D_bottom,jacobian2D_top, &
-! normal_xmin,normal_xmax, &
-! normal_ymin,normal_ymax, &
-! normal_bottom,normal_top, &
-! NSPEC2D_BOTTOM,NSPEC2D_TOP, &
-! NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX)
-
call store_boundaries(myrank,iboun,nspec, &
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
NSPEC2D_BOTTOM,NSPEC2D_TOP, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX)
-! create MPI buffers
-! arrays locval(npointot) and ifseg(npointot) used to save memory
-! call get_MPI_cutplanes_xi(myrank,prname,nspec,iMPIcut_xi,ibool, &
-! xstore,ystore,zstore,ifseg,npointot, &
-! NSPEC2D_A_ETA,NSPEC2D_B_ETA)
-! call get_MPI_cutplanes_eta(myrank,prname,nspec,iMPIcut_eta,ibool, &
-! xstore,ystore,zstore,ifseg,npointot, &
-! NSPEC2D_A_XI,NSPEC2D_B_XI)
-
-
call save_databases(prname,nspec,nglob,iproc_xi,iproc_eta,NPROC_XI,NPROC_ETA,addressing,iMPIcut_xi,iMPIcut_eta,&
ibool,nodes_coords,true_material_num,nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,&
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,&
@@ -1165,3 +633,5 @@
end subroutine create_regions_mesh
end module createRegMesh
+
+
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/define_subregions.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/define_subregions.f90 2010-03-19 05:07:35 UTC (rev 16433)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/define_subregions.f90 2010-03-19 20:11:02 UTC (rev 16434)
@@ -87,7 +87,7 @@
subroutine define_mesh_regions(USE_REGULAR_MESH,isubregion,NER,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,iproc_xi,iproc_eta,&
- nblayers,ner_layer,&
+ nblayers,ner_layer,ndoublings,ner_doublings,&
iaddx,iaddy,iaddz,ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar)
implicit none
@@ -104,13 +104,15 @@
integer iax,iay,iar
integer nblayers
integer num_material
+ integer ndoublings
! topology of the elements
integer iaddx(NGNOD)
integer iaddy(NGNOD)
integer iaddz(NGNOD)
integer ner_layer(nblayers)
-
+ integer ner_doublings(2)
+
! **************
!
@@ -138,9 +140,173 @@
iar=1
else
+ if(ndoublings == 1) then
- ! not iplemented yet
+ select case (isubregion)
+ case (1)
+ ix1=0
+ ix2=2*(NEX_PER_PROC_XI - 1)
+ dix=4
+
+ iy1=0
+ iy2=2*(NEX_PER_PROC_ETA - 1)
+ diy=4
+
+ ir1=0
+ ir2=2*(ner_doublings(1) -2 -1)
+ dir=2
+
+ iax=2
+ iay=2
+ iar=1
+
+ case (2)
+
+ ix1=0
+ ix2=2*(NEX_PER_PROC_XI - 1)
+ dix=8
+
+ iy1=0
+ iy2=2*(NEX_PER_PROC_ETA - 1)
+ diy=8
+
+ ir1=2*(ner_doublings(1) - 2)
+ ir2=2*(ner_doublings(1) - 2)
+ dir=2
+
+ iax=4
+ iay=4
+ iar=2
+
+ case (3)
+
+ ix1=0
+ ix2=2*(NEX_PER_PROC_XI - 1)
+ dix=2
+
+ iy1=0
+ iy2=2*(NEX_PER_PROC_ETA - 1)
+ diy=2
+
+ ir1=2*(ner_doublings(1))
+ ir2=2*(NER - 1)
+ dir=2
+
+ iax=1
+ iay=1
+ iar=1
+
+ case default
+ stop 'Wrong number of subregions'
+
+ end select
+
+ else if(ndoublings == 2) then
+
+ select case (isubregion)
+
+ case (1)
+ ix1=0
+ ix2=2*(NEX_PER_PROC_XI - 1)
+ dix=8
+
+ iy1=0
+ iy2=2*(NEX_PER_PROC_ETA - 1)
+ diy=8
+
+ ir1=0
+ ir2=2*(ner_doublings(2) -2 -1)
+ dir=2
+
+ iax=4
+ iay=4
+ iar=1
+
+ case (2)
+
+ ix1=0
+ ix2=2*(NEX_PER_PROC_XI - 1)
+ dix=16
+
+ iy1=0
+ iy2=2*(NEX_PER_PROC_ETA - 1)
+ diy=16
+
+ ir1=2*(ner_doublings(2) - 2)
+ ir2=2*(ner_doublings(2) - 2)
+ dir=2
+
+ iax=8
+ iay=8
+ iar=2
+
+
+ case (3)
+
+ ix1=0
+ ix2=2*(NEX_PER_PROC_XI - 1)
+ dix=4
+
+ iy1=0
+ iy2=2*(NEX_PER_PROC_ETA - 1)
+ diy=4
+
+ ir1=2*ner_doublings(2)
+ ir2=2*(ner_doublings(1) -2 -1)
+ dir=2
+
+ iax=2
+ iay=2
+ iar=1
+
+ case (4)
+
+ ix1=0
+ ix2=2*(NEX_PER_PROC_XI - 1)
+ dix=8
+
+ iy1=0
+ iy2=2*(NEX_PER_PROC_ETA - 1)
+ diy=8
+
+ ir1=2*(ner_doublings(1) - 2)
+ ir2=2*(ner_doublings(1) - 2)
+ dir=2
+
+ iax=4
+ iay=4
+ iar=2
+
+ case (5)
+
+ ix1=0
+ ix2=2*(NEX_PER_PROC_XI - 1)
+ dix=2
+
+ iy1=0
+ iy2=2*(NEX_PER_PROC_ETA - 1)
+ diy=2
+
+ ir1=2*(ner_doublings(1))
+ ir2=2*(NER - 1)
+ dir=2
+
+ iax=1
+ iay=1
+ iar=1
+
+ case default
+ stop 'Wrong number of subregions'
+
+ end select
+
+
+ else
+ stop 'Wrong number of doublings'
+
+ end if
+
end if
@@ -151,841 +317,840 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine define_subregions_old(myrank,isubregion,iaddx,iaddy,iaddz, &
- ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar, &
- doubling_index,npx,npy, &
- NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER,USE_REGULAR_MESH)
+ ! subroutine define_subregions_old(myrank,isubregion,iaddx,iaddy,iaddz, &
+! ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar, &
+! doubling_index,npx,npy, &
+! NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER,USE_REGULAR_MESH)
-! define shape of elements in current subregion of the mesh
+! ! define shape of elements in current subregion of the mesh
- implicit none
+! implicit none
- include "constants.h"
+! include "constants.h"
- integer myrank
- integer ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir
- integer iax,iay,iar
- integer isubregion,doubling_index
- integer npx,npy
+! integer myrank
+! integer ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir
+! integer iax,iay,iar
+! integer isubregion,doubling_index
+! integer npx,npy
- integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER
+! integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER
- logical USE_REGULAR_MESH
+! logical USE_REGULAR_MESH
-! topology of the elements
- integer iaddx(NGNOD)
- integer iaddy(NGNOD)
- integer iaddz(NGNOD)
+! ! topology of the elements
+! integer iaddx(NGNOD)
+! integer iaddy(NGNOD)
+! integer iaddz(NGNOD)
-! **************
+! ! **************
-!
-!--- case of a regular mesh
-!
- if(USE_REGULAR_MESH) then
+! !
+! !--- case of a regular mesh
+! !
+! if(USE_REGULAR_MESH) then
-! use two layers even for a regular mesh, because the algorithm detects the top of the mesh
-! (the "topography") based on one layer of elements with flag IFLAG_ONE_LAYER_TOPOGRAPHY
- if(isubregion == 2) then
+! ! use two layers even for a regular mesh, because the algorithm detects the top of the mesh
+! ! (the "topography") based on one layer of elements with flag IFLAG_ONE_LAYER_TOPOGRAPHY
+! if(isubregion == 2) then
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+! call usual_hex_nodes(iaddx,iaddy,iaddz)
- iy1=0
- iy2=npy-2
- diy=2
+! iy1=0
+! iy2=npy-2
+! diy=2
- ix1=0
- ix2=npx-2
- dix=2
+! ix1=0
+! ix2=npx-2
+! dix=2
- ir1=0
- ir2=2*(NER - 2)
- dir=2
+! ir1=0
+! ir2=2*(NER - 2)
+! dir=2
- iax=1
- iay=1
- iar=1
+! iax=1
+! iay=1
+! iar=1
- doubling_index = IFLAG_BASEMENT_TOPO
+! doubling_index = IFLAG_BASEMENT_TOPO
- else if(isubregion == 1) then
+! else
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+! iy1=0
+! iy2=npy-2
+! diy=2
- iy1=0
- iy2=npy-2
- diy=2
+! ix1=0
+! ix2=npx-2
+! dix=2
- ix1=0
- ix2=npx-2
- dix=2
+! ir1=2*(NER - 1)
+! ir2=ir1
+! dir=2
- ir1=2*(NER - 1)
- ir2=ir1
- dir=2
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! doubling_index = IFLAG_ONE_LAYER_TOPOGRAPHY
- doubling_index = IFLAG_ONE_LAYER_TOPOGRAPHY
+! else
- else
+! call exit_MPI(myrank,'incorrect subregion code')
- call exit_MPI(myrank,'incorrect subregion code')
+! endif
- endif
+! !
+! !--- case of a non-regular mesh with mesh doublings
+! !
+! else
-!
-!--- case of a non-regular mesh with mesh doublings
-!
- else
+! ! this last region only defined when NER_SEDIM > 1
+! if(isubregion == 30) then
-! this last region only defined when NER_SEDIM > 1
- if(isubregion == 30) then
+! call usual_hex_nodes(iaddx,iaddy,iaddz)
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+! iy1=0
+! iy2=npy-2
+! diy=2
- iy1=0
- iy2=npy-2
- diy=2
+! ix1=0
+! ix2=npx-2
+! dix=2
- ix1=0
- ix2=npx-2
- dix=2
+! ir1=2*(NER - NER_SEDIM)
+! ir2=2*(NER - 2)
+! dir=2
- ir1=2*(NER - NER_SEDIM)
- ir2=2*(NER - 2)
- dir=2
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! doubling_index = IFLAG_BASEMENT_TOPO
- doubling_index = IFLAG_BASEMENT_TOPO
+! else if(isubregion == 29) then
- else if(isubregion == 29) then
+! call usual_hex_nodes(iaddx,iaddy,iaddz)
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+! iy1=0
+! iy2=npy-2
+! diy=2
- iy1=0
- iy2=npy-2
- diy=2
+! ix1=0
+! ix2=npx-2
+! dix=2
- ix1=0
- ix2=npx-2
- dix=2
+! ir1=2*(NER - 1)
+! ir2=ir1
+! dir=2
- ir1=2*(NER - 1)
- ir2=ir1
- dir=2
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! doubling_index = IFLAG_ONE_LAYER_TOPOGRAPHY
- doubling_index = IFLAG_ONE_LAYER_TOPOGRAPHY
+! else if(isubregion == 28) then
- else if(isubregion == 28) then
+! call usual_hex_nodes(iaddx,iaddy,iaddz)
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+! iy1=0
+! iy2=npy-8
+! diy=8
- iy1=0
- iy2=npy-8
- diy=8
+! ix1=0
+! ix2=npx-8
+! dix=8
- ix1=0
- ix2=npx-8
- dix=8
+! ir1= 0
+! ir2= 2*NER_BOTTOM_MOHO-8
+! dir=8
- ir1= 0
- ir2= 2*NER_BOTTOM_MOHO-8
- dir=8
+! iax=4
+! iay=4
+! iar=4
- iax=4
- iay=4
- iar=4
+! doubling_index= IFLAG_HALFSPACE_MOHO
- doubling_index= IFLAG_HALFSPACE_MOHO
+! else if(isubregion == 27) then
- else if(isubregion == 27) then
+! call unusual_hex_nodes1(iaddx,iaddy,iaddz)
- call unusual_hex_nodes1(iaddx,iaddy,iaddz)
+! ! generating stage 2 of the mesh doubling below 670
-! generating stage 2 of the mesh doubling below 670
+! iy1=0
+! iy2=npy-8
+! diy=8
- iy1=0
- iy2=npy-8
- diy=8
+! ix1=0
+! ix2=npx-16
+! dix=16
- ix1=0
- ix2=npx-16
- dix=16
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 26) then
- else if(isubregion == 26) then
+! call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
- call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
+! ! generating stage 3 of the mesh doubling below 670
-! generating stage 3 of the mesh doubling below 670
+! iy1=0
+! iy2=npy-8
+! diy=8
- iy1=0
- iy2=npy-8
- diy=8
+! ix1=8
+! ix2=npx-8
+! dix=16
- ix1=8
- ix2=npx-8
- dix=16
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 25) then
- else if(isubregion == 25) then
+! call unusual_hex_nodes2(iaddx,iaddy,iaddz)
- call unusual_hex_nodes2(iaddx,iaddy,iaddz)
+! ! generating stage 4 of the mesh doubling below 670
-! generating stage 4 of the mesh doubling below 670
+! iy1=0
+! iy2=npy-8
+! diy=8
- iy1=0
- iy2=npy-8
- diy=8
+! ix1=0
+! ix2=npx-16
+! dix=16
- ix1=0
- ix2=npx-16
- dix=16
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 24) then
- else if(isubregion == 24) then
+! call unusual_hex_nodes2p(iaddx,iaddy,iaddz)
- call unusual_hex_nodes2p(iaddx,iaddy,iaddz)
+! ! generating stage 5 of the mesh doubling below 670
-! generating stage 5 of the mesh doubling below 670
+! iy1=0
+! iy2=npy-8
+! diy=8
- iy1=0
- iy2=npy-8
- diy=8
+! ix1=12
+! ix2=npx-4
+! dix=16
- ix1=12
- ix2=npx-4
- dix=16
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 23) then
- else if(isubregion == 23) then
+! call unusual_hex_nodes3(iaddx,iaddy,iaddz)
- call unusual_hex_nodes3(iaddx,iaddy,iaddz)
+! ! generating stage 6 of the mesh doubling below 670
-! generating stage 6 of the mesh doubling below 670
+! iy1=0
+! iy2=npy-8
+! diy=8
- iy1=0
- iy2=npy-8
- diy=8
+! ix1=4
+! ix2=npx-12
+! dix=16
- ix1=4
- ix2=npx-12
- dix=16
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 22) then
- else if(isubregion == 22) then
+! call unusual_hex_nodes3(iaddx,iaddy,iaddz)
- call unusual_hex_nodes3(iaddx,iaddy,iaddz)
+! ! generating stage 7 of the mesh doubling below 670
-! generating stage 7 of the mesh doubling below 670
+! iy1=0
+! iy2=npy-8
+! diy=8
- iy1=0
- iy2=npy-8
- diy=8
+! ix1=8
+! ix2=npx-8
+! dix=16
- ix1=8
- ix2=npx-8
- dix=16
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 21) then
- else if(isubregion == 21) then
+! call unusual_hex_nodes4(iaddx,iaddy,iaddz)
- call unusual_hex_nodes4(iaddx,iaddy,iaddz)
+! ! generating stage 8 of the mesh doubling below 670
-! generating stage 8 of the mesh doubling below 670
+! iy1=8
+! iy2=npy-8
+! diy=16
- iy1=8
- iy2=npy-8
- diy=16
+! ix1=0
+! ix2=npx-4
+! dix=4
- ix1=0
- ix2=npx-4
- dix=4
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 20) then
- else if(isubregion == 20) then
+! call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
- call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
+! ! generating stage 9 of the mesh doubling below 670
-! generating stage 9 of the mesh doubling below 670
+! iy1=0
+! iy2=npy-16
+! diy=16
- iy1=0
- iy2=npy-16
- diy=16
+! ix1=0
+! ix2=npx-4
+! dix=4
- ix1=0
- ix2=npx-4
- dix=4
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 19) then
- else if(isubregion == 19) then
+! call usual_hex_nodes(iaddx,iaddy,iaddz)
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+! ! generating stage 10 of the mesh doubling below 670
-! generating stage 10 of the mesh doubling below 670
+! iy1=8
+! iy2=npy-8
+! diy=16
- iy1=8
- iy2=npy-8
- diy=16
+! ix1=0
+! ix2=npx-4
+! dix=4
- ix1=0
- ix2=npx-4
- dix=4
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 18) then
- else if(isubregion == 18) then
+! call usual_hex_nodes(iaddx,iaddy,iaddz)
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+! ! generating stage 11 of the mesh doubling below 670
-! generating stage 11 of the mesh doubling below 670
+! iy1=4
+! iy2=npy-12
+! diy=16
- iy1=4
- iy2=npy-12
- diy=16
+! ix1=0
+! ix2=npx-4
+! dix=4
- ix1=0
- ix2=npx-4
- dix=4
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 17) then
- else if(isubregion == 17) then
+! call unusual_hex_nodes6(iaddx,iaddy,iaddz)
- call unusual_hex_nodes6(iaddx,iaddy,iaddz)
+! ! generating stage 12 of the mesh doubling below 670
-! generating stage 12 of the mesh doubling below 670
+! iy1=12
+! iy2=npy-4
+! diy=16
- iy1=12
- iy2=npy-4
- diy=16
+! ix1=0
+! ix2=npx-4
+! dix=4
- ix1=0
- ix2=npx-4
- dix=4
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 16) then
- else if(isubregion == 16) then
+! call unusual_hex_nodes6p(iaddx,iaddy,iaddz)
- call unusual_hex_nodes6p(iaddx,iaddy,iaddz)
+! ! generating stage 13 of the mesh doubling below 670
-! generating stage 13 of the mesh doubling below 670
+! iy1=0
+! iy2=npy-16
+! diy=16
- iy1=0
- iy2=npy-16
- diy=16
+! ix1=0
+! ix2=npx-4
+! dix=4
- ix1=0
- ix2=npx-4
- dix=4
+! dir=4
- dir=4
+! iax=2
+! iay=2
+! iar=1
- iax=2
- iay=2
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
- ir2=ir1
+! doubling_index=IFLAG_16km_BASEMENT
- doubling_index=IFLAG_16km_BASEMENT
+! else if(isubregion == 15) then
- else if(isubregion == 15) then
+! call usual_hex_nodes(iaddx,iaddy,iaddz)
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+! iy1=0
+! iy2=npy-8
+! diy=8
- iy1=0
- iy2=npy-8
- diy=8
+! ix1=0
+! ix2=npx-8
+! dix=8
- ix1=0
- ix2=npx-8
- dix=8
+! ! honor So-Cal model discontinuity at 16 km for accuracy
+! ir1=2*NER_BOTTOM_MOHO
+! ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16) - 4
+! dir=4
-! honor So-Cal model discontinuity at 16 km for accuracy
- ir1=2*NER_BOTTOM_MOHO
- ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16) - 4
- dir=4
+! iax=4
+! iay=4
+! iar=2
- iax=4
- iay=4
- iar=2
+! doubling_index = IFLAG_MOHO_16km
- doubling_index = IFLAG_MOHO_16km
+! else if(isubregion == 14) then
- else if(isubregion == 14) then
+! call usual_hex_nodes(iaddx,iaddy,iaddz)
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+! iy1=0
+! iy2=npy-8
+! diy=8
- iy1=0
- iy2=npy-8
- diy=8
+! ix1=0
+! ix2=npx-8
+! dix=8
- ix1=0
- ix2=npx-8
- dix=8
+! ! honor So-Cal model discontinuity at 16 km for accuracy
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16)
+! ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT)-12
+! dir=4
-! honor So-Cal model discontinuity at 16 km for accuracy
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16)
- ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT)-12
- dir=4
+! iax=4
+! iay=4
+! iar=2
- iax=4
- iay=4
- iar=2
+! doubling_index = IFLAG_16km_BASEMENT
- doubling_index = IFLAG_16km_BASEMENT
+! else if(isubregion == 13) then
- else if(isubregion == 13) then
+! call usual_hex_nodes(iaddx,iaddy,iaddz)
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+! ! generating stage 1 of the mesh doubling below the Moho
-! generating stage 1 of the mesh doubling below the Moho
+! iy1=0
+! iy2=npy-4
+! diy=4
- iy1=0
- iy2=npy-4
- diy=4
+! ix1=0
+! ix2=npx-4
+! dix=4
- ix1=0
- ix2=npx-4
- dix=4
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT)
+! ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-12
+! dir=4
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT)
- ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-12
- dir=4
+! iax=2
+! iay=2
+! iar=2
- iax=2
- iay=2
- iar=2
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 12) then
- else if(isubregion == 12) then
+! call unusual_hex_nodes1(iaddx,iaddy,iaddz)
- call unusual_hex_nodes1(iaddx,iaddy,iaddz)
+! ! generating stage 2 of the mesh doubling below the Moho
-! generating stage 2 of the mesh doubling below the Moho
+! iy1=0
+! iy2=npy-4
+! diy=4
- iy1=0
- iy2=npy-4
- diy=4
+! ix1=0
+! ix2=npx-8
+! dix=8
- ix1=0
- ix2=npx-8
- dix=8
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 11) then
- else if(isubregion == 11) then
+! call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
- call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
+! ! generating stage 3 of the mesh doubling below the Moho
-! generating stage 3 of the mesh doubling below the Moho
+! iy1=0
+! iy2=npy-4
+! diy=4
- iy1=0
- iy2=npy-4
- diy=4
+! ix1=4
+! ix2=npx-4
+! dix=8
- ix1=4
- ix2=npx-4
- dix=8
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 10) then
- else if(isubregion == 10) then
+! call unusual_hex_nodes2(iaddx,iaddy,iaddz)
- call unusual_hex_nodes2(iaddx,iaddy,iaddz)
+! ! generating stage 4 of the mesh doubling below the Moho
-! generating stage 4 of the mesh doubling below the Moho
+! iy1=0
+! iy2=npy-4
+! diy=4
- iy1=0
- iy2=npy-4
- diy=4
+! ix1=0
+! ix2=npx-8
+! dix=8
- ix1=0
- ix2=npx-8
- dix=8
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 9) then
- else if(isubregion == 9) then
+! call unusual_hex_nodes2p(iaddx,iaddy,iaddz)
- call unusual_hex_nodes2p(iaddx,iaddy,iaddz)
+! ! generating stage 5 of the mesh doubling below the Moho
-! generating stage 5 of the mesh doubling below the Moho
+! iy1=0
+! iy2=npy-4
+! diy=4
- iy1=0
- iy2=npy-4
- diy=4
+! ix1=6
+! ix2=npx-2
+! dix=8
- ix1=6
- ix2=npx-2
- dix=8
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 8) then
- else if(isubregion == 8) then
+! call unusual_hex_nodes3(iaddx,iaddy,iaddz)
- call unusual_hex_nodes3(iaddx,iaddy,iaddz)
+! ! generating stage 6 of the mesh doubling below the Moho
-! generating stage 6 of the mesh doubling below the Moho
+! iy1=0
+! iy2=npy-4
+! diy=4
- iy1=0
- iy2=npy-4
- diy=4
+! ix1=2
+! ix2=npx-6
+! dix=8
- ix1=2
- ix2=npx-6
- dix=8
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 7) then
- else if(isubregion == 7) then
+! call unusual_hex_nodes3(iaddx,iaddy,iaddz)
- call unusual_hex_nodes3(iaddx,iaddy,iaddz)
+! ! generating stage 7 of the mesh doubling below the Moho
-! generating stage 7 of the mesh doubling below the Moho
+! iy1=0
+! iy2=npy-4
+! diy=4
- iy1=0
- iy2=npy-4
- diy=4
+! ix1=4
+! ix2=npx-4
+! dix=8
- ix1=4
- ix2=npx-4
- dix=8
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 6) then
- else if(isubregion == 6) then
+! call unusual_hex_nodes4(iaddx,iaddy,iaddz)
- call unusual_hex_nodes4(iaddx,iaddy,iaddz)
+! ! generating stage 8 of the mesh doubling below the Moho
-! generating stage 8 of the mesh doubling below the Moho
+! iy1=4
+! iy2=npy-4
+! diy=8
- iy1=4
- iy2=npy-4
- diy=8
+! ix1=0
+! ix2=npx-2
+! dix=2
- ix1=0
- ix2=npx-2
- dix=2
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 5) then
- else if(isubregion == 5) then
+! call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
- call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
+! ! generating stage 9 of the mesh doubling below the Moho
-! generating stage 9 of the mesh doubling below the Moho
+! iy1=0
+! iy2=npy-8
+! diy=8
- iy1=0
- iy2=npy-8
- diy=8
+! ix1=0
+! ix2=npx-2
+! dix=2
- ix1=0
- ix2=npx-2
- dix=2
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 4) then
- else if(isubregion == 4) then
+! call usual_hex_nodes(iaddx,iaddy,iaddz)
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+! ! generating stage 10 of the mesh doubling below the Moho
-! generating stage 10 of the mesh doubling below the Moho
+! iy1=4
+! iy2=npy-4
+! diy=8
- iy1=4
- iy2=npy-4
- diy=8
+! ix1=0
+! ix2=npx-2
+! dix=2
- ix1=0
- ix2=npx-2
- dix=2
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 3) then
- else if(isubregion == 3) then
+! call usual_hex_nodes(iaddx,iaddy,iaddz)
- call usual_hex_nodes(iaddx,iaddy,iaddz)
+! ! generating stage 11 of the mesh doubling below the Moho
-! generating stage 11 of the mesh doubling below the Moho
+! iy1=2
+! iy2=npy-6
+! diy=8
- iy1=2
- iy2=npy-6
- diy=8
+! ix1=0
+! ix2=npx-2
+! dix=2
- ix1=0
- ix2=npx-2
- dix=2
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 2) then
- else if(isubregion == 2) then
+! call unusual_hex_nodes6(iaddx,iaddy,iaddz)
- call unusual_hex_nodes6(iaddx,iaddy,iaddz)
+! ! generating stage 12 of the mesh doubling below the Moho
-! generating stage 12 of the mesh doubling below the Moho
+! iy1=6
+! iy2=npy-2
+! diy=8
- iy1=6
- iy2=npy-2
- diy=8
+! ix1=0
+! ix2=npx-2
+! dix=2
- ix1=0
- ix2=npx-2
- dix=2
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else if(isubregion == 1) then
- else if(isubregion == 1) then
+! call unusual_hex_nodes6p(iaddx,iaddy,iaddz)
- call unusual_hex_nodes6p(iaddx,iaddy,iaddz)
+! ! generating stage 13 of the mesh doubling below the Moho
-! generating stage 13 of the mesh doubling below the Moho
+! iy1=0
+! iy2=npy-8
+! diy=8
- iy1=0
- iy2=npy-8
- diy=8
+! ix1=0
+! ix2=npx-2
+! dix=2
- ix1=0
- ix2=npx-2
- dix=2
+! dir=4
- dir=4
+! iax=1
+! iay=1
+! iar=1
- iax=1
- iay=1
- iar=1
+! ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
+! ir2=ir1
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
- ir2=ir1
+! doubling_index=IFLAG_BASEMENT_TOPO
- doubling_index=IFLAG_BASEMENT_TOPO
+! else
- else
+! call exit_MPI(myrank,'incorrect subregion code')
- call exit_MPI(myrank,'incorrect subregion code')
+! endif
- endif
+! endif
- endif
+! end subroutine define_subregions_old
- end subroutine define_subregions_old
-
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/meshfem3D.f90 2010-03-19 05:07:35 UTC (rev 16433)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/meshfem3D.f90 2010-03-19 20:11:02 UTC (rev 16434)
@@ -247,6 +247,9 @@
logical SUPPRESS_UTM_PROJECTION,USE_REGULAR_MESH
! integer NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
+ integer NDOUBLINGS
+ integer, dimension(2) :: ner_doublings
+
character(len=150) OUTPUT_FILES,LOCAL_PATH,MODEL
! parameters deduced from parameters read from file
@@ -339,48 +342,24 @@
write(IMAIN,*)
endif
-! read the parameter file
-! call read_parameter_file (LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX, &
-! UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
-! NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT,NER_MOHO_16,NER_BOTTOM_MOHO, &
-! NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA,UTM_PROJECTION_ZONE, &
-! HARVARD_3D_GOCAD_MODEL,TOPOGRAPHY,LOCAL_PATH, &
-! THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
-! IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL, &
-! BASEMENT_MAP,MOHO_MAP_LUPEI, &
-! SUPPRESS_UTM_PROJECTION,MODEL,&
-! interfacesfile,NSUBREGIONS,subregions,NMATERIALS,material_properties,&
-! USE_REGULAR_MESH)
- USE_REGULAR_MESH = .true.
-
! nullify(subregions,material_properties)
call read_parameter_file(LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX, &
UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA,UTM_PROJECTION_ZONE, &
LOCAL_PATH,SUPPRESS_UTM_PROJECTION,&
- INTERFACES_FILE,NSUBREGIONS,subregions,NMATERIALS,material_properties)
+ INTERFACES_FILE,NSUBREGIONS,subregions,NMATERIALS,material_properties, &
+ USE_REGULAR_MESH,NDOUBLINGS,ner_doublings)
if (sizeprocs == 1 .and. (NPROC_XI /= 1 .or. NPROC_ETA /= 1)) then
stop 'must have NPROC_XI = NPROC_ETA = 1 for a serial run'
endif
-! ! compute other parameters based upon values read
-! call compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
-! NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
-! NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM, &
-! NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
-! NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
-! NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
-! NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,USE_REGULAR_MESH)
-
-! ! check that the code is running with the requested nb of processes
-! if(sizeprocs /= NPROC) call exit_MPI(myrank,'wrong number of MPI processes')
-
-
-!! PLL ici ouverture fichier interface
! get interface data from external file to count the spectral elements along Z
- write(IMAIN,*) 'Reading interface data from file DATA/',INTERFACES_FILE(1:len_trim(INTERFACES_FILE)), &
- ' to count the spectral elements'
+ if(myrank == 0) then
+ write(IMAIN,*) 'Reading interface data from file DATA/',INTERFACES_FILE(1:len_trim(INTERFACES_FILE)), &
+ ' to count the spectral elements'
+ end if
+
open(unit=IIN,file='DATA/'//INTERFACES_FILE,status='old')
max_npx_interface = -1
@@ -388,7 +367,7 @@
! read number of interfaces
call read_value_integer(DONT_IGNORE_JUNK,number_of_interfaces,'NINTERFACES')
- if(number_of_interfaces < 2) stop 'not enough interfaces (minimum is 2)'
+ if(number_of_interfaces < 1) stop 'not enough interfaces (minimum is 1, for topography)'
! loop on all the interfaces
do interface_current = 1,number_of_interfaces
@@ -432,7 +411,8 @@
NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,USE_REGULAR_MESH)
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,&
+ USE_REGULAR_MESH,NDOUBLINGS,ner_doublings)
! check that the code is running with the requested nb of processes
if(sizeprocs /= NPROC) call exit_MPI(myrank,'wrong number of MPI processes')
@@ -525,11 +505,11 @@
if(mod(NEX_XI,NPROC_XI) /= 0) call exit_MPI(myrank,'NEX_XI must be a multiple of NPROC_XI for a regular mesh')
if(mod(NEX_ETA,NPROC_ETA) /= 0) call exit_MPI(myrank,'NEX_ETA must be a multiple of NPROC_ETA for a regular mesh')
else
- if(mod(NEX_XI,8) /= 0) call exit_MPI(myrank,'NEX_XI must be a multiple of 8 for a non-regular mesh')
- if(mod(NEX_ETA,8) /= 0) call exit_MPI(myrank,'NEX_ETA must be a multiple of 8 for a non-regular mesh')
+ !if(mod(NEX_XI,8) /= 0) call exit_MPI(myrank,'NEX_XI must be a multiple of 8 for a non-regular mesh')
+ !if(mod(NEX_ETA,8) /= 0) call exit_MPI(myrank,'NEX_ETA must be a multiple of 8 for a non-regular mesh')
- if(mod(NEX_XI/8,NPROC_XI) /= 0) call exit_MPI(myrank,'NEX_XI must be a multiple of 8*NPROC_XI for a non-regular mesh')
- if(mod(NEX_ETA/8,NPROC_ETA) /= 0) call exit_MPI(myrank,'NEX_ETA must be a multiple of 8*NPROC_ETA for a non-regular mesh')
+ !if(mod(NEX_XI/8,NPROC_XI) /= 0) call exit_MPI(myrank,'NEX_XI must be a multiple of 8*NPROC_XI for a non-regular mesh')
+ !if(mod(NEX_ETA/8,NPROC_ETA) /= 0) call exit_MPI(myrank,'NEX_ETA must be a multiple of 8*NPROC_ETA for a non-regular mesh')
endif
if(myrank == 0) then
@@ -659,10 +639,12 @@
min_elevation = +HUGEVAL
max_elevation = -HUGEVAL
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'Reading interface data from file DATA/',INTERFACES_FILE(1:len_trim(INTERFACES_FILE))
+ write(IMAIN,*)
+ end if
- write(IMAIN,*)
- write(IMAIN,*) 'Reading interface data from file DATA/',INTERFACES_FILE(1:len_trim(INTERFACES_FILE))
- write(IMAIN,*)
open(unit=IIN,file='DATA/'//INTERFACES_FILE,status='old')
allocate(interface_bottom(max_npx_interface,max_npy_interface))
@@ -1056,7 +1038,8 @@
NPROC_XI,NPROC_ETA,NSPEC2D_A_XI,NSPEC2D_B_XI, &
NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
NSUBREGIONS,subregions,number_of_layers,ner_layer,NMATERIALS,material_properties, &
- myrank,LOCAL_PATH,UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK,USE_REGULAR_MESH)
+ myrank,LOCAL_PATH,UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK,&
+ USE_REGULAR_MESH,NDOUBLINGS,ner_doublings)
! ! print min and max of topography included
! if(TOPOGRAPHY) then
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/read_parameter_file.f90 2010-03-19 05:07:35 UTC (rev 16433)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/read_parameter_file.f90 2010-03-19 20:11:02 UTC (rev 16434)
@@ -30,8 +30,8 @@
UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA,UTM_PROJECTION_ZONE, &
LOCAL_PATH,SUPPRESS_UTM_PROJECTION,&
- INTERFACES_FILE,NSUBREGIONS,subregions,NMATERIALS,material_properties)!,&
- !USE_REGULAR_MESH)
+ INTERFACES_FILE,NSUBREGIONS,subregions,NMATERIALS,material_properties,&
+ USE_REGULAR_MESH,NDOUBLINGS,ner_doublings)
implicit none
@@ -42,8 +42,11 @@
double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, UTM_MAX
double precision LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX
- logical SUPPRESS_UTM_PROJECTION!,USE_REGULAR_MESH
+ logical SUPPRESS_UTM_PROJECTION,USE_REGULAR_MESH
+ integer NDOUBLINGS
+ integer, dimension(2) :: ner_doublings
+
character(len=150) LOCAL_PATH
character(len=50) INTERFACES_FILE
@@ -123,6 +126,16 @@
NEX_MAX = max(NEX_XI,NEX_ETA)
UTM_MAX = max(UTM_Y_MAX-UTM_Y_MIN, UTM_X_MAX-UTM_X_MIN)/1000.0 ! in KM
+ call read_value_logical(IGNORE_JUNK,USE_REGULAR_MESH, 'mesher.USE_REGULAR_MESH')
+ if(err_occurred() /= 0) return
+ call read_value_integer(IGNORE_JUNK,NDOUBLINGS, 'mesher.NDOUBLINGS')
+ if(err_occurred() /= 0) return
+ call read_value_integer(IGNORE_JUNK,ner_doublings(1), 'mesher.NZ_DOUGLING_1')
+ if(err_occurred() /= 0) return
+ call read_value_integer(IGNORE_JUNK,ner_doublings(2), 'mesher.NZ_DOUGLING_2')
+ if(err_occurred() /= 0) return
+
+
! file in which we store the databases
call read_value_string(IGNORE_JUNK,LOCAL_PATH, 'LOCAL_PATH')
if(err_occurred() /= 0) return
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/save_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/save_databases.f90 2010-03-19 05:07:35 UTC (rev 16433)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/meshfem3D/save_databases.f90 2010-03-19 20:11:02 UTC (rev 16434)
@@ -100,7 +100,7 @@
write(15,*) nspec
do ispec=1,nspec
- write(15,'(11i8)') ispec,true_material_num(ispec),1,ibool(1,1,1,ispec),ibool(2,1,1,ispec),&
+ write(15,'(11i14)') ispec,true_material_num(ispec),1,ibool(1,1,1,ispec),ibool(2,1,1,ispec),&
ibool(2,2,1,ispec),ibool(1,2,1,ispec),ibool(1,1,2,ispec),&
ibool(2,1,2,ispec),ibool(2,2,2,ispec),ibool(1,2,2,ispec)
end do
More information about the CIG-COMMITS
mailing list